* Date: Tue, 24 Nov 92 10:20:47 -0700
* From: seeger@gem.LANL.GOV
SUBROUTINE PGPLOT10
C
C Subroutines to implement high-level functions from the author's
C PLOT-10 library, and to interpret calls to PLOT-10 primatives
C as corresponding entries in the PGPLOT library. When logarithmic
C axes have been defined by calling GRAPH, the PLOT-10 move and draw
C calls will automatically take logarithms before plotting. Additional
C high-level calls (CONTOUR, CURVE, ERRBAR) which duplicate existing
C PGPLOT functions (but which support logarithmic scaling) are included
C at the end of this file.
C
C P. A. Seeger, Los Alamos National Laboratory, Oct. 6, 1992
C
C Entries:
C ANCHO DASHA DASHR DRAWA DRAWR DRWABS DRWREL
C DSHABS DSHREL DWINDO ERASE FINITT GRAPH HDCOPY
C HLABEL INITT MOVABS MOVEA MOVER MOVREL OSTRING
C PNTABS PNTREL POINTA POINTR RESET SCURSR SEEDW
C SEELOC SEETW SETCOLOR SWINDO SYMBOL TERM TWINDO
C VCURSR VLABEL VWINDO
C
C Auxiliary Externals:
C CHECK_TEXT LINAXIS LOGAXIS
C
C PGPLOT Externals:
C PGBEGIN PGBOX PGCURSE PGDRAW PGEND PGETXT PGMOVE
C PGMTEXT PGNUMB PGPOINT PGQINF PGQPOS PGQVP PGSCH
C PGSCI PGSLS PGSLW PGVPORT PGWINDOW
C
IMPLICIT NONE
C
C Variables which may occur in calling sequences of entry points
CHARACTER TITLE*(*),SUBTITLE*(*),XLABEL*(*),YLABEL*(*),CH*1
REAL XMN,XMX,X,YMN,YMX,Y
INTEGER LOGX,LOGY,L,IX,IY,LX,LY
C
C Local temporary variables:
CHARACTER XOPTION*7,YOPTION*8,XFORMAT*6,YFORMAT*6,STRING*12,
1 NEWSTRING*255
REAL XX,XMAJOR,XMINOR,YY,YMAJOR,YMINOR
INTEGER I,J,IXMAJOR,IXTICK,IYMAJOR,IYTICK,NC,LL,IX2,IY2
LOGICAL LPOINT
C
C Local variables to be saved between entries:
CHARACTER TERMTYPE*40,HDCPYTYPE*40
REAL ASPECT,X1,X2,XMIN,XMAX,XPERPIX,X0,Y1,Y2,YMIN,YMAX,YPERPIX,Y0
INTEGER LINE,ISYMB,ITERM,IHDCPY
LOGICAL XLOG,YLOG,PLOT10
SAVE PLOT10,TERMTYPE,HDCPYTYPE,ITERM,IHDCPY,ASPECT,LINE,ISYMB,
1 X1,X2,XMIN,XMAX,XLOG,XPERPIX,X0,
2 Y1,Y2,YMIN,YMAX,YLOG,YPERPIX,Y0
DATA PLOT10, TERMTYPE,HDCPYTYPE,ITERM,IHDCPY,ASPECT,LINE,ISYMB
1 /.FALSE.,'?', '?', 1, 1, 0.75, -1, -1/
DATA X1,X2, XPERPIX,X0,Y1,Y2, YPERPIX,Y0
1 /0.,1023.,1., 0.,0.,767.,1., 0./
C
ENTRY INITT(X)
C Initialize graphics terminal device
plot10 = .true.
call pgbegin(0, termtype(1:iterm), 1, 1)
if (termtype.eq.'?') then
call pgqinf('dev/type',termtype,iterm)
iterm = min0(iterm,index(termtype,'/')+3)
end if
return
C
ENTRY HDCOPY
C Initialize hardcopy output device
call pgbegin(0, hdcpytype(1:ihdcpy), 1, 1)
if (termtype.eq.'?') then
call pgqinf('dev/type',hdcpytype,ihdcpy)
ihdcpy = min0(ihdcpy,index(hdcpytype,'/')+3)
end if
return
C
ENTRY ERASE
ENTRY RESET
C Erase the screen
call pgpage
ENTRY TERM(IX,IY)
return
C
ENTRY FINITT(X,Y)
C Exit graphics
call pgend
return
C
C Define graph area in world co-ordinates, to fill the available
C viewport area. Place two lines of titles at top, and axis labels
C on bottom and left. Draw box with tick marks and labels; logarithmic
C if the corresponding LOGX or LOGY value is non-zero. XMN, XMX, YMN,
C and YMX will be changed to rounded values.
C
ENTRY GRAPH(TITLE, SUBTITLE, XLABEL, YLABEL,
1 XMN, XMX, LOGX, YMN, YMX, LOGY)
call pgetxt
call pgqvp(3,x1,x2,y1,y2)
aspect = (y2-y1+1.)/(x2-x1+1.)
call pgvport(0.20*aspect, 1.0-0.05*aspect, 0.10, 0.85 )
C
call pgsch(2.0)
call pgslw(3)
call check_text(title,newstring,j,plot10)
call pgmtext('T', 2.0, 0.5, 0.5, newstring(1:j))
C
call pgsch(1.5)
call pgslw(2)
call check_text(subtitle,newstring,j,plot10)
call pgmtext('T', 1.2, 0.5, 0.5, newstring(1:j))
C
call check_text(xlabel,newstring,j,plot10)
call pgmtext('B', 2.16, 0.5, 0.5, newstring(1:j))
C
call check_text(ylabel,newstring,j,plot10)
call pgmtext('L', 4.0, 0.5, 0.5, newstring(1:j))
call pgsch(1.0)
call pgslw(1)
C
xlog = logx.gt.0
ylog = logy.gt.0
if (xlog) then
call logaxis(xmn, xmx, ixmajor, ixtick, xminor, xformat)
xmin = alog10(xmn)
xmax = alog10(xmx)
xmajor = 1.
ixtick = 9
xoption = 'BCLNTS'
if (.not.ylog) xoption = 'ABCLNTS'
else
call linaxis(xmn, xmx, ixmajor, ixtick, xformat)
xmin = xmn
xmax = xmx
xmajor = (xmax-xmin)/float(ixmajor)
xoption = 'BCNTS'
if (.not.ylog) xoption = 'ABCNTS'
end if
C
if (ylog) then
call logaxis(ymn, ymx, iymajor, iytick, yminor, yformat)
ymin = alog10(ymn)
ymax = alog10(ymx)
ymajor = 1.
iytick = 9
yoption = 'BCLNTSV'
if (.not.xlog) yoption = 'ABCLNTSV'
else
call linaxis(ymn, ymx, iymajor, iytick, yformat)
ymin = ymn
ymax = ymx
ymajor = (ymax-ymin)/float(iymajor)
yoption = 'BCNTSV'
if (.not.xlog) yoption = 'ABCNTSV'
end if
C
call pgwindow(xmin, xmax, ymin, ymax)
call pgbox(xoption, xmajor, ixtick, yoption, ymajor, iytick)
xperpix = (xmax-xmin)/(x2-x1)/(1.-0.25*aspect)
x0 = x1 + 0.20*aspect*(x2-x1)
yperpix = (ymax-ymin)/(y2-y1)/0.75
y0 = y1 + 0.10*(y2-y1)
C
if (xlog) then
C Make sure ends of X-axis are labeled
if (amod(abs(xmin)+0.01,1.) .gt. 0.02) then
i = nint(xmin-0.5)
j = nint(xmn/10.**i)
call pgnumb(j, i, 0, string, nc)
call pgmtext('B', 1.25, 0., 0.5, string(1:nc))
end if
if (amod(abs(xmax)+0.01,1.) .gt. 0.02) then
i = nint(xmax-0.5)
j = nint(xmx/10.**i)
call pgnumb(j, i, 0, string, nc)
call pgmtext('B', 1.25, 1., 0.5, string(1:nc))
end if
end if
C
if (ylog) then
C Make sure ends of Y-axis are labeled
if (amod(abs(ymin)+0.01,1.) .gt. 0.02) then
i = nint(ymin-0.5)
j = nint(ymn/10.**i)
call pgnumb(j, i, 0, string, nc)
call pgmtext('LV', aspect, 0., 1., string(1:nc))
end if
if (amod(abs(ymax)+0.01,1.) .gt. 0.02) then
i = nint(ymax-0.5)
j = nint(ymx/10.**i)
call pgnumb(j, i, 0, string, nc)
call pgmtext('LV', aspect, 1., 1., string(1:nc))
end if
end if
return
C
C Following entries either move to a new point without drawing, or
C move to a new point and draw a single dot or a symbol, accounting
C for possibility of log scales.
C
ENTRY MOVEA(X,Y)
lpoint = .false.
xx = x
yy = y
go to 100
C
ENTRY MOVABS(IX,IY)
lpoint = .false.
xx = xmin + xperpix*(float(ix)-x0)
yy = ymin + yperpix*(float(iy)-y0)
go to 110
C
ENTRY POINTA(X,Y)
lpoint = .true.
ll = -1
xx = x
yy = y
go to 100
C
ENTRY SYMBOL(X,Y,L)
lpoint = .true.
ll = l
xx = x
yy = y
go to 100
C
ENTRY ANCHO(L)
lpoint = .true.
ll = l
call pgqpos(xx,yy)
go to 110
C
ENTRY PNTABS(IX,IY)
lpoint = .true.
ll = -1
xx = xmin + xperpix*(float(ix)-x0)
yy = ymin + yperpix*(float(iy)-y0)
go to 110
C
ENTRY MOVREL(IX,IY)
lpoint = .false.
go to 80
C
ENTRY PNTREL(IX,IY)
lpoint = .true.
ll = -1
80 call pgqpos(xx,yy)
xx = xx + xperpix*float(ix)
yy = yy + yperpix*float(iy)
go to 110
C
ENTRY MOVER(X,Y)
lpoint = .false.
go to 90
C
ENTRY POINTR(X,Y)
lpoint = .true.
ll = -1
90 continue
call pgqpos(xx,yy)
if (xlog) then
xx = (10.**xx) + x
else
xx = xx + x
end if
if (ylog) then
yy = (10.**yy) + y
else
yy = yy + y
end if
100 continue
if (xlog) then
if (xx.le.0.) then
xx = -38.
else
xx = alog10(xx)
end if
end if
if (ylog) then
if (yy.le.0.) then
yy = -38.
else
yy = alog10(yy)
end if
end if
110 if (lpoint) then
call pgpoint(1, xx, yy, ll)
else
call pgmove(xx,yy)
end if
return
C
C Following entries draw a solid or dashed line from the current
C location to a new point, accounting for possibility of log scales.
C
ENTRY DRAWA(X,Y)
ll = 0
xx = x
yy = y
go to 200
C
ENTRY DRWABS(IX,IY)
ll = 0
xx = xmin + xperpix*(float(ix)-x0)
yy = ymin + yperpix*(float(iy)-y0)
go to 210
C
ENTRY DASHA(X,Y,L)
ll = l
xx = x
yy = y
go to 200
C
ENTRY DSHABS(IX,IY,L)
ll = l
xx = xmin + xperpix*(float(ix)-x0)
yy = ymin + yperpix*(float(iy)-y0)
go to 210
C
ENTRY DRWREL(IX,IY)
ll = 0
go to 180
C
ENTRY DSHREL(IX,IY,L)
ll = l
180 call pgqpos(xx,yy)
xx = xx + xperpix*float(ix)
yy = yy + yperpix*float(iy)
go to 210
C
ENTRY DRAWR(X,Y)
ll = 0
go to 190
C
ENTRY DASHR(X,Y,L)
ll = l
190 continue
call pgqpos(xx,yy)
if (xlog) then
xx = (10.**xx) + x
else
xx = xx + x
end if
if (ylog) then
yy = (10.**yy) + y
else
yy = yy + y
end if
200 continue
if (xlog) then
if (xx.le.0.) then
xx = -38.
else
xx = alog10(xx)
end if
end if
if (ylog) then
if (yy.le.0.) then
yy = -38.
else
yy = alog10(yy)
end if
end if
210 if (ll.lt.0) then
call pgmove(xx,yy)
else
if (ll.ne.line) then
line = mod(ll,5)
call pgsls(line+1)
end if
call pgdraw(xx,yy)
end if
return
C
C Set boundaries of viewport or of world co-ordinate system
C
ENTRY VWINDO(XMN, X, YMN, Y)
xmax = xmn+x
ymax = ymn+y
go to 250
C
ENTRY DWINDO(XMN, XMX, YMN, YMX)
xmax = xmx
ymax = ymx
250 xmin = xmn
ymin = ymn
call pgwindow(xmin, xmax, ymin, ymax)
xperpix = (xmax-xmin)/(x2-x1+1.)
x0 = x1
yperpix = (ymax-ymin)/(y2-y1+1.)
y0 = y1
return
C
ENTRY SWINDOW(IX, LX, IY, LY)
ix2 = x+lx
iy2 = y+ly
go to 260
C
ENTRY TWINDO(IX, LX, IY, LY)
ix2 = lx
iy2 = ly
260 call pgqvp(3, x1, x2, y1, y2)
call pgvport(float(ix)/(x2-x1), float(ix2)/(x2-x1),
1 float(iy)/(y2-y1), float(iy2)/(y2-y1))
call pgqvp(3, x1, x2, y1, y2)
xperpix = (xmax-xmin)/(x2-x1+1.)
x0 = x1
yperpix = (ymax-ymin)/(y2-y1+1.)
y0 = y1
return
C
C Write horizontal or vertical text strings
C
ENTRY HLABEL(TITLE, IX, IY, LX, LY)
string = 'B'
xx = xperpix*(float(ix)-x0)/(xmax-xmin)
ly = nint((y2-y1)/40.)
lx = (3*ly)/4
yy = -(float(iy)-y0)/((y2-y1)/40.)
go to 300
C
ENTRY OSTRING(TITLE,L)
call pgqpos(xx,yy)
if (l.eq.0) then
string = 'B'
xx = (xx-xmin)/(xmax-xmin)
yy = (yy-ymin)/yperpix/(y2-y1)*40.
else
string = 'L'
xx = (yy-ymin)/(ymax-ymin)
yy = (xx-xmin)/yperpix/(y2-y1)*40.
end if
go to 300
C
ENTRY VLABEL(TITLE, IX, IY, LX, LY)
string = 'L'
xx = yperpix*(float(iy)-y0)/(ymax-ymin)
ly = nint((y2-y1)/40.)
lx = (3*ly)/4
yy = -(float(iy)-y0)/((y2-y1)/40.)
300 continue
call check_text(title,newstring,j,plot10)
call pgmtext(string, yy, xx, 0.5, newstring(1:j))
return
C
C Find ("see") various window parameters
C
ENTRY SEETW(IX, LX, IY, LY)
ix = x1
lx = x2
iy = y1
ly = y2
return
C
ENTRY SEEDW(XMN,XMX,YMN,YMX)
if (xlog) then
xmn = 10.**xmin
xmx = 10.**xmax
else
xmn = xmin
xmx = xmax
end if
if (ylog) then
ymn = 10.**ymin
ymx = 10.**ymax
else
ymn = ymin
ymx = ymax
end if
return
C
ENTRY SEELOC(IX,IY)
call pgqpos(xx,yy)
ix = x0+(xx-xmin)/xperpix
iy = y0+(yy-ymin)/yperpix
return
C
ENTRY SETCOLOR(L)
call pgsci(l)
return
C
ENTRY SCURSR(L, IX, IY)
xx = xmin+xperpix*(float(ix)-x0)
yy = ymin+yperpix*(float(iy)-y0)
call pgcurse(xx, yy, ch)
l = ichar(ch)
if (l.gt.0) then
ix = x0+(xx-xmin)/xperpix
iy = y0+(yy-ymin)/yperpix
end if
return
C
ENTRY VCURSR(L, X, Y)
if (xlog) then
if (x.le.0.) then
xx = xmin
else
xx = alog10(x)
end if
else
xx = x
end if
if (ylog) then
if (y.le.0.) then
yy = ymin
else
yy = alog10(y)
end if
else
yy = y
end if
call pgcurse(xx,yy,ch)
l = ichar(ch)
if (l.ne.0) then
if (xlog) then
x = 10.**xx
else
x = xx
end if
if (ylog) then
y = 10.**yy
else
y = yy
end if
end if
return
C
END
C
SUBROUTINE CHECK_TEXT(IN,OUT,KOUT,PLOT10_FLAG)
C
C Decode special characters in PLOT10 string to PGPLOT
C
IMPLICIT NONE
CHARACTER IN*(*),OUT*(*)
INTEGER KOUT
LOGICAL PLOT10_FLAG
C
CHARACTER SPECIAL(5)*1,CH2*2,SYMBOL(8)*6
INTEGER I,J,K,IIN,JIN,IOUT,MODE,KMAX
DATA SPECIAL/'<', '>', '?', '#', '&'/
DATA SYMBOL/'\(845)', '\(847)', '\(840)', '\(846)', '\(841)',
1 '\(842)', '\(843)', '\(852)'/
C
C Omit trailing blanks and nulls
do 2 jin=len(in),3,-1
if (in(jin:jin).ne.' ' .and. in(jin:jin).ne.char(0)) go to 3
2 continue
C Omit terminal '$'
3 if (in(jin:jin).eq.'$') jin = jin-1
C
kmax = len(out)
j = 0
k = 0
mode = 2
10 continue
j = j+1
if (plot10_flag) then
C Look for PLOT10 special characters
do 80 i=1,5
if (in(j:j).eq.special(i)) then
if ((i.eq.1.or.i.eq.2).and.(i.ne.mode)) then
k = k+3
out(k-2:k) = '\f1'
else if (i.eq.5) then
C Convert next 2 characters to lower case for testing
ch2(1:1) = char(ior(ichar(in(j+1:j+1)),32))
ch2(2:2) = char(ior(ichar(in(j+2:j+2)),32))
if (ch2.eq.'ex') then
C End of superscript
j = j+2
k = k+2
out(k-1:k) = '\d'
else if (ch2(1:1).eq.'e') then
C Beginning of superscript
j = j+1
k = k+2
out(k-1:k) = '\u'
else if (ch2.eq.'lx') then
C End of subscript
j = j+2
k = k+2
out(k-1:k) = '\u'
else if (ch2(1:1).eq.'l') then
C Beginning of subscript
j = j+1
k = k+2
out(k-1:k) = '\d'
end if
end if
mode = i
go to 100
end if
80 continue
C Not a special case, just keep the character
if (mode.eq.3 .or. mode.eq.4) then
C Character is Greek, must be preceded with flag
k = k+2
out(k-1:k) = '\g'
end if
k = k+1
if (mode.eq.1 .or. mode.eq.4) then
C Character must be lower case
out(k:k) = char(ior(ichar(in(j:j)),32))
else
C No case modification
out(k:k) = in(j:j)
end if
else
C
C Look for PGPLOT symbol numbers
if (in(j:j).eq.'\') then
i = ichar(in(j+1:j+1)) - ichar('0')
if (i.ge.1 .and. i.le.8) then
j = j+1
k = k+6
out(k-5:k) = symbol(i)
go to 100
end if
end if
C Not symbol number, copy unmodified character to output
k = k+1
out(k:k) = in(j:j)
end if
C
100 if (j.lt.jin .and. k.lt.kmax) go to 10
kout = min0(k,kmax)
return
END
C
SUBROUTINE LINAXIS(XONE,XTWO,MAJOR,MINOR,FORMAT)
C
C SCALE ENDS OF LINEAR AXIS TO ROUNDED NUMBERS
C
C An axis is defined to include XONE and XTWO, with each end rounded to
C be an integer number of units of the form (1, 2, or 5) x 10**n. The
C resulting number MAJOR of major divisions will be between 4 and 10, with
C MINOR minor divisions in each step. A character string (FORMAT = 'Fww.dd')
C is generated to use as the format for labeling the axis.
C
C P. A. Seeger, Los Alamos National Laboratory, May 24, 1986
C Modified to prevent log of X=0., Aug. 31, 1986
C Use E instead of F format if width > 10, Nov. 21, 1987
C Guard against axis ends > 10**38, Feb. 5, 1990
C
C No Externals
C
IMPLICIT NONE
REAL*4 XONE,XTWO,XMIN,XMAX,AXLEN,UNIT
INTEGER MAJOR,MINOR,IPOWR,NN,IW,ID
CHARACTER*6 FORMAT
C
XMIN = AMAX1(-0.7E38,AMIN1(XONE,XTWO))
XMAX = AMIN1( 0.7E38,AMAX1(XONE,XTWO))
IF (XMAX.EQ.XMIN) THEN
IF (XMAX.LT.0.) THEN
XMAX = 0.
ELSE IF (XMIN.GT.0.) THEN
XMIN = 0.
ELSE
XMAX = 0.001
END IF
END IF
AXLEN = XMAX-XMIN
IPOWR = NINT(ALOG10(AXLEN)-1.05)
UNIT = 10.**IPOWR
NN = AXLEN/UNIT
IF (NN.GE.15) THEN
UNIT=UNIT*5.
MINOR=1
IF (NN.LE.30) MINOR=5
ELSE IF (NN.GE.7) THEN
UNIT=UNIT*2.
MINOR=2
IF (NN.LE.10) MINOR=4
ELSE
MINOR=2
IF (NN.LE.4) MINOR=5
END IF
C
IF (XMIN.GE.0.) THEN
XMIN = UNIT*AINT((XMIN+0.01*AXLEN)/UNIT)
ELSE
XMIN = UNIT*AINT((XMIN+0.01*AXLEN)/UNIT-1.)
END IF
IF (XMAX.GE.0.) THEN
XMAX = UNIT*AINT((XMAX-0.01*AXLEN)/UNIT+1.)
ELSE
XMAX = UNIT*AINT((XMAX-0.01*AXLEN)/UNIT)
END IF
C
MAJOR = NINT((XMAX-XMIN)/UNIT)
IF (XTWO.GT.XONE) THEN
XONE = XMIN
XTWO = XMAX
ELSE
XONE = XMAX
XTWO = XMIN
END IF
C
ID = MAX0(0,-IPOWR)
IF (XMIN.LT.0.) XMIN=-10.*XMIN
IF (XMAX.LT.0.) XMAX=-10.*XMAX
IPOWR = ALOG10(AMAX1(XMIN,XMAX,1.))+0.001
IW = IPOWR+ID+2
IF (IW.LT.10) THEN
WRITE (FORMAT,100) IW,ID
100 FORMAT ('F',I2,'.',I2)
ELSE
FORMAT = 'E 9. 2'
END IF
C
RETURN
END
C
SUBROUTINE LOGAXIS(XONE,XTWO,MAJOR,MINOR,UNIT,FORMAT)
C
C SCALE ENDS OF LOGARITHMIC AXIS TO ROUNDED NUMBERS
C
C An axis is defined to include XONE and XTWO, with each end rounded to
C be (1, 2, 3, or 5) x some power of 10. The resulting axis spans MAJOR
C decades, including partial decades at one or both ends. The lowest decade
C on the axis has MINOR units of size UNIT; if it is a full decade, MINOR = 9
C and UNIT = min(X). A character string (FORMAT = 'Fww.dd') is generated to
C use as the format for labeling the axis.
C
C P. A. Seeger, Los Alamos National Laboratory, May 20, 1986
C Modified to avoid log(X) when X.LE.0., June 17, 1986
C Corrected FORMAT when XMAX<1., June 19, 1987
C Use E instead of F format when width > 10, Nov. 21, 1987
C Guard against axis ends > 10**38, Feb. 5, 1990
C
C No Externals
C
IMPLICIT NONE
REAL*4 XONE,XTWO,UNIT,XMIN,XMAX,XLOG,LOG2,LOG3,LOG5,DELTA
INTEGER MAJOR,MINOR,IMIN,IMAX,IW,ID
CHARACTER*6 FORMAT
PARAMETER (LOG2=0.3010300,LOG3=0.4771213,LOG5=0.6989700)
C
XMIN = AMIN1(XONE,XTWO)
XMAX = AMAX1(XONE,XTWO)
IF (XMAX.LE.0.) XMAX = 1.
IF (XMIN.LE.0.) THEN
XMIN = XMAX/1000.
ELSE IF (XMAX.GE.1.E38) THEN
XMAX = XMIN*1.E6
ELSE
XMIN = AMAX1(XMIN,XMAX/1.E24)
END IF
IF (XMIN.EQ.XMAX) THEN
XMAX = XMAX*2.
XMIN = XMIN/2.
END IF
DELTA = 0.01*ALOG10(XMAX/XMIN)
C
C Look at small end of axis first
XLOG = ALOG10(XMIN)+DELTA
IMIN = NINT(XLOG-0.5)
XMIN = 10.**IMIN
UNIT = XMIN
XLOG = XLOG-FLOAT(IMIN)
IF (XLOG.GT.LOG5) THEN
XMIN = 5.*XMIN
MINOR = 5
ELSE IF (XLOG.GT.LOG3) THEN
XMIN = 3.*XMIN
MINOR = 7
ELSE IF (XLOG.GT.LOG2) THEN
XMIN = 2.*XMIN
MINOR = 8
ELSE
MINOR = 9
END IF
ID = MAX0(0,-IMIN)
C
C Now do "same" thing at larger end
XLOG = ALOG10(XMAX)-DELTA
IMAX = NINT(XLOG-0.5)
C Compute total number of decades included
MAJOR = IMAX-IMIN+1
XMAX = 10.**IMAX
XLOG = XLOG-FLOAT(IMAX)
IF (XLOG.LT.LOG2) THEN
XMAX = 2.*XMAX
ELSE IF (XLOG.LT.LOG3) THEN
XMAX = 3.*XMAX
ELSE IF (XLOG.LT.LOG5) THEN
XMAX = 5.*XMAX
ELSE
XMAX = 10.*XMAX
IMAX = IMAX+1
END IF
IF (MAJOR.EQ.1) MINOR=NINT((XMAX-XMIN)/UNIT)
C
IW = MAX0(0,IMAX)+ID+2
IF (IW.LT.10) THEN
WRITE (FORMAT,100) IW,ID
100 FORMAT ('F',I2,'.',I2)
ELSE
FORMAT = 'E 9. 2'
END IF
IF (XONE.LT.XTWO) THEN
XONE = XMIN
XTWO = XMAX
ELSE
XONE = XMAX
XTWO = XMIN
END IF
C
RETURN
END
C
C***********************************************************************
C
SUBROUTINE CONTOUR(Z,NRZ,X,NX,Y,NY,CV,NCV,LINE,ZMAX,BITMAP)
C
C DRAW CONTOURS THROUGH EQUAL VALUES IN AN ARRAY
C
C From Collected Algorithms from ACM #531 "Contour Plotting [J6]"
C by: William V. Snyder, 1978
C Copied from GSAS; transferred to IBM-PC, Dec. 8, 1987 (P.A.Seeger)
C Added X, Y, and LINE to calling sequence, Dec. 9, 1987 (PAS)
C Restructured, closer to "standard", Dec. 11, 1987 (PAS)
C Changed I and J to II and JJ in innermost loop, Dec. 13, 1987 (PAS)
C Revised scan pattern from spiral to linear, Aug. 13, 1988 (PAS)
C
C Externals
C DASHA FILL0 LGETMARK MOVEA
C
C Arguments in calling sequence:
C Z R(NRZ,*) Input Array of values to be contoured; nodes must
C lie on a topologically rectangular grid.
C NRZ I Input Number of rows declared for array Z
C X R(*) Input Values of X at grid points of array Z
C NX I Input Limit for 1st subscript of Z and X
C Y R(*) Input Values of Y at grid points of array Z
C NY I Input Limit for 2nd subscript of Z and Y
C CV R(*) Input Values of contour levels
C NCV I Input No. of contour levels
C LINE I(*) Input Line style for each contour level
C ZMAX R Input Maximum Z to be considered; grid lines at
C a node with value above ZMAX are excluded
C BITMAP I(*) Scratch Work area of size (2*NX*NY*NCV+7)/8 bytes
C
IMPLICIT NONE
INTEGER NRZ,NX,NY,NCV,LINE(*)
REAL*4 Z(NRZ,*),X(*),Y(*),CV(*),ZMAX
INTEGER*4 BITMAP(*)
C
LOGICAL LGETMARK
C
C Local variables in CONTOUR:
C CVAL R Contour-line value being traced
C DELZ R Change in Z when moving 1 cell right or up
C IADD I Bit address in BITMAP, starting at zero
C IBORDER I Edge of current cell which is on a border
C ICV I Index of contour line being traced
C IEDGE I Edge of new cell where contour line enters
C IFLAG I 1=continue, 2=start at boundary, 3=start in interior,
C 4=end at boundary, 5=close contour, 6=none found yet
C I,J I Subscripts for search
C IJ I(2) Equivalent to I,J
C IJMAX I(2) Local copies of NX and NY
C I1(2),I2(2),I3(6) Used for subscript computations
C II,JJ I Cell with continuation of contour line being traced
C INDEX I L-1 + 2*(I-1 + NX*(J-1))
C K I Index of cell edges: 1=bottom, 2=left, 3=top, 4=right
C KS I Next cell boundary to cross
C L,LL I Orientation flags: 1 is horizontal line, 2 vertical
C LBORDER L Flag to show that only border lines are to be done
C NI I Number of edges of cell which contour crosses
C XINT R(4) Intersections of contour with edges of cell,
C in order bottom, left, top, right
C XX,YY R Plotting coordinates
C Z1,Z2 R Smaller and larger values at segment ends
C Z3,Z4 R Values at ends of segment for continuation of contour
C
REAL*4 XINT(4),Z1,Z2,Z3,Z4,DELZ,CVAL,XX,YY
INTEGER I1(2),I2(2),I3(6),I,J,IJ(2),IJMAX(2),K,L,LL,II,JJ,ICV,
, IBORDER,IEDGE,IFLAG,NI,KS
INTEGER*4 INDEX,MAXINDX,NBITS,IADD
LOGICAL LBORDER
EQUIVALENCE (IJ(1),I),(IJ(2),J)
DATA I1/1,0/, I2/1,-1/, I3/1,0,0,1,1,0/
C
IJMAX(1) = NX
IJMAX(2) = NY
C Clear bit map
NBITS = 2*NX*NY*NCV
CALL FILL0(BITMAP,NBITS)
IFLAG = 6
C
C Search every cell in rectangular array for a line segment such that:
C 1. the end points are not excluded because Z > ZMAX
C 2. current contour <= Z at one end and > Z at other end
C 3. no mark has been recorded for this contour on this segment
C
C Set start at lower left corner of array; do borders first, then interior
C
LBORDER = .TRUE.
MAXINDX = 2*NX*NY-2
100 CONTINUE
DO 500 INDEX=0,MAXINDX
L = MOD(INDEX,2)+1
I = MOD(INDEX/2,NX)+1
J = (INDEX/2)/NX+1
IF (Z(I,J).LE.ZMAX) THEN
C Node itself is within non-excluded range
II = I+I1(L)
JJ = J+I1(3-L)
IF (II.LE.NX .AND. JJ.LE.NY .AND. Z(II,JJ).LE.ZMAX) THEN
C Both ends of grid line within range; test if "border" segment
IBORDER = 0
IF (IJ(3-L).EQ.1 .OR.
1 Z(I-I1(3-L),J-I1(L)).GT.ZMAX .OR.
2 Z(I+I2(L),J+I2(3-L)).GT.ZMAX) IBORDER = 1
IF (IJ(3-L).GE.IJMAX(3-L) .OR.
1 Z(I+I1(3-L),J+I1(L)).GT.ZMAX .OR.
2 Z(I+1,J+1).GT.ZMAX) IBORDER = IBORDER+2
C 1st time, do ONLY borders (including edges of exclusions)
IF (IBORDER.NE.3 .AND. (LBORDER.EQV.(IBORDER.NE.0))) THEN
C Examine this line segment this pass
Z1 = AMIN1(Z(I,J),Z(II,JJ))
Z2 = AMAX1(Z(I,J),Z(II,JJ))
DELZ = Z(II,JJ)-Z(I,J)
DO 400 ICV=1,NCV
C Test for all possible contours crossing this grid line;
C first check if already done, and set bit to show done
IADD = ICV-1 + NCV*INDEX
IF (.NOT.LGETMARK(BITMAP,IADD) .AND.
1 CV(ICV).GT.Z1 .AND. CV(ICV).LE.Z2) THEN
C Found one we haven't done yet!!! Interpolate.
CVAL = CV(ICV)
C Decide which edge we are entering cell from
IEDGE = L
IF (IBORDER.EQ.2) IEDGE = IEDGE+2
XINT(IEDGE) = (CVAL-Z(I,J))/DELZ
C Move "pen" to starting point of contour
XX = X(I)+XINT(IEDGE)*(X(II)-X(I))
YY = Y(J)+XINT(IEDGE)*(Y(JJ)-Y(J))
IFLAG = 3
IF (LBORDER) IFLAG = 2
CALL MOVEA(XX,YY)
C
C Follow this contour until it hits boundary or closes on itself
II = I
JJ = J
IFLAG = 1
C "DO WHILE (IFLAG.LT.4)"
200 CONTINUE
C If haven't moved to next cell yet, do so
IF (IEDGE.EQ.3) JJ = JJ-1
IF (IEDGE.EQ.4) II = II-1
NI = 1
DO 300 K = 1,4
C Test interpolation on other 3 edges
IF (K.NE.IEDGE) THEN
Z3 = Z(II+I3(K),JJ+I3(K+1))
Z4 = Z(II+I3(K+1),JJ+I3(K+2))
IF (CVAL.GT.AMIN1(Z3,Z4) .AND.
. CVAL.LE.AMAX1(Z3,Z4)) THEN
C The contour also crosses this edge
IF (K.EQ.1 .OR. K.EQ.4) THEN
XINT(K) = (CVAL-Z4)/(Z3-Z4)
ELSE
XINT(K) = (CVAL-Z3)/(Z4-Z3)
END IF
C Count how many crossings
NI = NI+1
KS = K
END IF
END IF
300 CONTINUE
C
IF (NI.NE.2) THEN
C The contour crosses all four edges of the cell being examined. Choose the
C lines top-to-left and bottom-to-right if the interpolation point on the top
C edge is less than the interpolation point on the bottom edge. Otherwise,
C choose the other pair. This method produces the same result if the axes are
C reversed. The contour may close at any edge, but must not cross itself
C inside any cell.
IF (XINT(3).GE.XINT(1)) THEN
KS = 3-IEDGE
IF (KS.LE.0) KS = KS+4
ELSE
KS = 5-IEDGE
END IF
END IF
C
C Determine if the contour will close or run into a boundary at edge KS of the
C current cell.
IF (KS.LE.2) THEN
LL = KS
IEDGE = KS+2
ELSE
C Must move to adjacent cell before test for closure
II = II+I3(KS)
JJ = JJ+I3(KS+2)
LL = KS-2
IEDGE = KS-2
END IF
IADD = ICV-1+NCV*(LL-1+2*(II-1+NX*(JJ-1)))
IF (LGETMARK(BITMAP,IADD)) THEN
C We've already been here; contour has closed
IFLAG = 5
ELSE IF (LL.EQ.2.AND.(II.EQ.1.OR.II.GE.NX)
1 .OR. LL.EQ.1.AND.(JJ.EQ.1.OR.JJ.GE.NY)) THEN
C Segment is actual boundary of plot
IFLAG = 4
ELSE IF (Z(II-I1(3-LL),JJ-I1(LL)).GT.ZMAX
1 .OR. Z(II+I2(LL),JJ+I2(3-LL)).GT.ZMAX
2 .OR. Z(II+I1(3-LL),JJ+I1(LL)).GT.ZMAX
3 .OR. Z(II+1,JJ+1).GT.ZMAX) THEN
C Segment is boundary of an excluded cell
IFLAG = 4
END IF
C
C Draw piece of contour
XINT(IEDGE) = XINT(KS)
IF (LL.EQ.1) THEN
XX = X(II)+XINT(IEDGE)*(X(II+1)-X(II))
YY = Y(JJ)
ELSE
XX = X(II)
YY = Y(JJ)+XINT(IEDGE)*(Y(JJ+1)-Y(JJ))
END IF
CALL DASHA(XX,YY,LINE(ICV))
C
IF (IFLAG.LT.4) GO TO 200
C Reset II and JJ before looking for next contour
II = I+I1(L)
JJ = J+I1(3-L)
END IF
400 CONTINUE
END IF
END IF
END IF
500 CONTINUE
C
IF (.NOT.LBORDER) RETURN
LBORDER = .FALSE.
GO TO 100
END
C
C CONTOUR SUBROUTINES FOR BIT MANIPULATION
C
SUBROUTINE FILL0(BITMAP,N)
C Fills entire BITMAP with zeros
IMPLICIT NONE
INTEGER N,I,LOOP
INTEGER*4 BITMAP(*)
C
LOOP = (N-1)/32+1
DO 10 I=1,LOOP
BITMAP(I) = 0
10 CONTINUE
RETURN
END
C
LOGICAL FUNCTION LGETMARK(BITMAP,N)
C Tests bit in BITMAP, and then sets it to one
IMPLICIT NONE
INTEGER N,NWORD,NBIT
INTEGER*4 BITMAP(*)
LOGICAL BTEST
C
NWORD = N/32+1
NBIT = MOD(N,32)
LGETMARK = BTEST(BITMAP(NWORD),NBIT)
IF (.NOT.LGETMARK) BITMAP(NWORD) = IBSET(BITMAP(NWORD),NBIT)
RETURN
END
C
C******************************************************************************
C
SUBROUTINE CURVE(X,Y,N,LDASH,ISYM)
C
C PLOT N POINTS AT (X(I),Y(I)), USING SYMBOL CORRESPONDING TO ISYM, AND
C CONNECTING WITH LINE DESCRIBED BY LDASH.
C
C P. A. Seeger, Los Alamos National Laboratory, May 24, 1986
C
INTEGER N,LDASH,ISYM,I
REAL*4 X(*),Y(*)
C
C Externals
C DASHA MOVEA SYMBOL
C
CALL MOVEA(X(1),Y(1))
DO 100 I=1,N
CALL DASHA(X(I),Y(I),LDASH)
CALL SYMBOL(X(I),Y(I),ISYM)
100 CONTINUE
C
RETURN
END
C
C******************************************************************************
C
SUBROUTINE ERRBARS(X,Y,DY,N,INC,LOGY)
C
C PLOT ERROR BARS ON EVERY INCth POINT OF AN ARRAY OF N POINTS
C (X(I),Y(I)+-DY(I)). LOGY=1 IF ORDINATE IS LOGARITHMIC.
C
C P. A. Seeger, Los Alamos National Laboratory, Nov. 24, 1987
C
IMPLICIT NONE
INTEGER N,INC,LOGY,I
REAL*4 X(*),Y(*),DY(*),D,R,Y1,Y2
C
C Externals
C DRAWA MOVEA
C
DO 100 I=1,N,INC
D = ABS(DY(I))
IF (LOGY.EQ.1) THEN
IF (Y(I).LE.0.) THEN
Y1 = 2.E-38
Y2 = 2.E-38
ELSE IF (D.LT.0.05*Y(I)) THEN
Y1 = Y(I)-D
Y2 = Y(I)+D
ELSE IF (D.LT.88.*Y(I)) THEN
R = EXP(D/Y(I))
Y1 = Y(I)/R
Y2 = Y(I)*R
ELSE
Y1 = 2.E-38
Y2 = 1.E+38
END IF
ELSE
Y1 = Y(I)-D
Y2 = Y(I)+D
END IF
CALL MOVEA(X(I),Y1)
CALL DRAWA(X(I),Y2)
100 CONTINUE
C
RETURN
END