aboutsummaryrefslogtreecommitdiff
path: root/plotlib/plt_util.f
diff options
context:
space:
mode:
Diffstat (limited to 'plotlib/plt_util.f')
-rw-r--r--plotlib/plt_util.f930
1 files changed, 930 insertions, 0 deletions
diff --git a/plotlib/plt_util.f b/plotlib/plt_util.f
new file mode 100644
index 0000000..692f491
--- /dev/null
+++ b/plotlib/plt_util.f
@@ -0,0 +1,930 @@
+C***********************************************************************
+C Module: plt_util.f
+C
+C Copyright (C) 1996 Harold Youngren, Mark Drela
+C
+C This library is free software; you can redistribute it and/or
+C modify it under the terms of the GNU Library General Public
+C License as published by the Free Software Foundation; either
+C version 2 of the License, or (at your option) any later version.
+C
+C This library is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+C Library General Public License for more details.
+C
+C You should have received a copy of the GNU Library General Public
+C License along with this library; if not, write to the Free
+C Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+C
+C Report problems to: guppy@maine.com
+C or drela@mit.edu
+C***********************************************************************
+C
+***********************************************************************
+C --- Xplot11 utility routines
+C
+C Version 4.46 11/28/01
+C
+C Note: These are additional routines that supply additional plot
+C functionality. Included are routines for axis scaling,
+C axis plotting, line or curve plotting, contours, etc.
+C***********************************************************************
+
+
+ subroutine XAXIS(X1,Y1,XAXT,DXANN,FANN,DANN,CHT,NDIG)
+C.......................................................
+C
+C X1,Y1 starting point of x axis
+C XAXT length of x axis ( - = suppress zero annotation)
+C DXANN distance between annotations
+C FANN first annotation value
+C DANN delta annotation value
+C CHT character width ( - = annotation above axis)
+C NDIG number of digits to right of decimal point
+C = -1 no decimal point
+C = -2 number of digits determined internally
+C <= -3 no axis annotation (just hash marks)
+C.......................................................
+C
+ XAX = ABS(XAXT)
+ CH = ABS(CHT)
+C
+ IF(NDIG.LE.-2) THEN
+ ADANN = ABS(DANN)
+ ND = MAX( 1 , INT(-LOG10(ADANN)) )
+ IF(ADANN*10**ND - AINT(ADANN*10**ND+0.01) .GT. 0.01) ND = ND + 1
+ IF(ADANN*10**ND - AINT(ADANN*10**ND+0.01) .GT. 0.01) ND = ND + 1
+ ELSE
+ ND = NDIG
+ ENDIF
+C
+ CALL GETFACTORS(XFAC,YFAC)
+ CHX = CH
+ CHY = CH*XFAC/YFAC
+C
+C---- x-axis
+ CALL PLOT(X1,Y1,3)
+ CALL PLOT(X1+XAX,Y1,2)
+C
+C---- annotate x-axis
+ DO 10 NT=1, 12345
+ XT = X1 + DXANN*FLOAT(NT-1)
+ IF(XT-X1.GT.XAX+0.5*DXANN) GO TO 11
+C---- hash marks
+ CALL PLOT(XT,Y1-0.2*CHY,3)
+ CALL PLOT(XT,Y1+0.2*CHY,2)
+C---- Numeric annotations
+ IF(NDIG.LE.-3) GO TO 10
+ RN = FANN + DANN*FLOAT(NT-1)
+ IF(ABS(RN).LT.1.0E-5 .AND. XAXT.LT.0.0) GO TO 10
+ GRN = 0.
+ IF(RN.NE.0.0) GRN = LOG10(ABS(RN)+0.5/10.0**MAX(0,ND))
+ GRN = MAX(GRN,0.0)
+ NABC = INT(GRN) + 2 + ND
+ WIDTH = CHX*FLOAT(NABC)
+ IF(RN.LT.0.0) WIDTH = WIDTH + CHX
+ XNUM = XT - 0.5*WIDTH + 0.1*CHX
+ YNUM = Y1 - 2.1*CHY
+ IF(CHT.LT.0.0) YNUM = Y1 + 0.9*CHY
+ CALL PLNUMB(XNUM,YNUM,CH,RN,0.0,ND)
+ 10 CONTINUE
+ 11 CONTINUE
+C
+ RETURN
+ END
+
+
+ subroutine YAXIS(X1,Y1,YAXT,DYANN,FANN,DANN,CHT,NDIG)
+C.......................................................
+C
+C X1,Y1 starting point of y axis
+C YAXT length of y axis ( - = suppress zero annotation)
+C DYANN distance between annotations
+C FANN first annotation value
+C DANN delta annotation value
+C CHT character width ( - = annotation on right side )
+C NDIG number of digits to right of decimal point
+C = -1 no decimal point
+C = -2 number of digits determined internally
+C <= -3 no axis annotation (just hash marks)
+C.......................................................
+C
+ YAX = ABS(YAXT)
+ CH = ABS(CHT)
+C
+ IF(NDIG.LE.-2) THEN
+ ADANN = ABS(DANN)
+ ND = MAX( 1 , INT(-LOG10(ADANN)) )
+ IF(ADANN*10**ND - AINT(ADANN*10**ND+0.01) .GT. 0.01) ND = ND + 1
+ IF(ADANN*10**ND - AINT(ADANN*10**ND+0.01) .GT. 0.01) ND = ND + 1
+ ELSE
+ ND = NDIG
+ ENDIF
+C
+ CALL GETFACTORS(XFAC,YFAC)
+ CHX = CH
+ CHY = CH*XFAC/YFAC
+C
+C---- y-axis
+ CALL PLOT(X1,Y1,3)
+ CALL PLOT(X1,Y1+YAX,2)
+C
+C---- annotate y-axis
+ DO 10 NT=1, 12345
+ YT = Y1 + DYANN*FLOAT(NT-1)
+ IF(YT-Y1.GT.YAX+0.5*DYANN) GO TO 11
+C---- hash marks
+ CALL PLOT(X1-0.2*CHX,YT,3)
+ CALL PLOT(X1+0.2*CHX,YT,2)
+C---- Numeric annotations
+ IF(NDIG.LE.-3) GO TO 10
+ RN = FANN + DANN*FLOAT(NT-1)
+ IF(ABS(RN).LT.1.0E-5 .AND. YAXT.LT.0.0) GO TO 10
+ GRN = 0.
+ IF(RN.NE.0.0) GRN = LOG10(ABS(RN)+0.5/10.0**MAX(0,ND))
+ GRN = MAX(GRN,0.0)
+ NABC = INT(GRN) + 2 + ND
+ WIDTH = CHX*FLOAT(NABC)
+ IF(RN.LT.0.0) WIDTH = WIDTH + CHX
+ XT = X1 - (0.6*CHX + WIDTH)
+ IF(CHT.LT.0.0) XT = X1 + CHX
+ CALL PLNUMB(XT,YT-0.5*CHY,CH,RN,0.0,ND)
+ 10 CONTINUE
+ 11 CONTINUE
+C
+ RETURN
+ END
+
+
+
+ subroutine XYLINE(N,X,Y,XOFF,XWT,YOFF,YWT,ILIN)
+C....................................................................
+C
+C...General XY polyline plotting routine with offsets and scaling
+C
+C...INPUT X, Y Input arrays of length N
+C XOFF,XWT Offset and scale factor for X array...
+C YOFF,YWT Offset and scale factor for Y array...
+C Xplot = XWT*(X-XOFF)
+C Yplot = YWT*(Y-YOFF)
+C ILIN Selects line pattern
+C
+C...8 line patterns are available (repeat for ILIN>8)
+C 1 ***************************** SOLID
+C 2 **** **** **** **** **** **** LONG DASHED
+C 3 ** ** ** ** ** ** ** ** ** ** SHORT DASHED
+C 4 * * * * * * * * * * * * * * * DOTTED
+C 5 ***** * ***** * ***** * ***** DASH-DOT
+C 6 ***** * * ***** * * ***** * * DASH-DOT-DOT
+C 7 ***** * * * ***** * * * ***** DASH-DOT-DOT-DOT
+C 8 **** **** * * **** **** * * DASH-DASH-DOT-DOT
+C
+C....................................................................
+ DIMENSION X(N), Y(N)
+C
+ DIMENSION NMOV(7), SMOV(8,7)
+C
+ DATA NPAT / 8 /
+ DATA SCL1 / 0.125 /
+ DATA NMOV / 2, 2, 2, 4, 6, 8, 8 /
+ DATA SMOV /1.2, -.4, 0., 0., 0., 0., 0., 0.,
+ & .5, -.4, 0., 0., 0., 0., 0., 0.,
+ & .2, -.4, 0., 0., 0., 0., 0., 0.,
+ & 1.4, -.4, .2, -.4, 0., 0., 0., 0.,
+ & 1.4, -.4, .2, -.4, .2, -.4, 0., 0.,
+ & 1.4, -.4, .2, -.4, .2, -.4, .2, -.4,
+ & 1.2, -.4, 1.2, -.4, .2, -.4, .2, -.4 /
+C
+ IF(N.LE.1) RETURN
+C
+C---- set line pattern scale based on current user scaling factors
+ CALL GETFACTORS(XSCALE,YSCALE)
+ SCL = SCL1 / SQRT(XSCALE*YSCALE)
+C
+ NLIN = MAX(ILIN,1)
+ IPAT = MOD(NLIN-1,NPAT) + 1
+C
+ X2 = XWT*(X(1)-XOFF)
+ Y2 = YWT*(Y(1)-YOFF)
+ CALL PLOT(X2,Y2,3)
+C
+ IF (IPAT.EQ.1) THEN
+C...Plot using continuous line
+ DO 10 I=2, N
+ X1 = X2
+ Y1 = Y2
+ X2 = XWT*(X(I)-XOFF)
+ Y2 = YWT*(Y(I)-YOFF)
+ CALL PLOT(X2,Y2,2)
+ 10 CONTINUE
+C
+ ELSE
+C...Plot using stored patterns for lines
+ I = 1
+ S1 = 0.
+ S2 = 0.
+ S0 = 0.
+C
+ 20 DO 40 II=1, 99999
+C
+C...Pattern specifies pen up or down
+ IM = MOD(II+1,NMOV(IPAT-1)) + 1
+ IPEN = 3
+ IF(SMOV(IM,IPAT-1).GT.0.) IPEN = 2
+C
+ DS = SCL*ABS(SMOV(IM,IPAT-1))
+ SPAT = S0 + DS
+C
+C...Find data interval containing pattern point
+ 30 IF (SPAT.GE.S2 .AND. I+1.LE.N) THEN
+ I = I + 1
+ CALL PLOT(X2,Y2,IPEN)
+ X1 = X2
+ Y1 = Y2
+ S1 = S2
+ X2 = XWT*(X(I)-XOFF)
+ Y2 = YWT*(Y(I)-YOFF)
+ DS = SQRT((X2-X1)**2 + (Y2-Y1)**2)
+ S2 = S1 + DS
+ GO TO 30
+ ENDIF
+C
+C...Find point on interval using linear interpolation
+ IF (SPAT.GT.S2) SPAT = S2
+ IF(S2 .EQ. S1) THEN
+ FRAC = 0.0
+ ELSE
+ FRAC = (SPAT-S1)/(S2-S1)
+ ENDIF
+ XX = X1 + FRAC*(X2-X1)
+ YY = Y1 + FRAC*(Y2-Y1)
+C
+C...Move to new point using stored pattern to specify pen up or down
+ CALL PLOT(XX,YY,IPEN)
+ IF (I.GE.N .AND. SPAT.GE.S2) GO TO 50
+ S0 = SPAT
+ 40 CONTINUE
+C
+ 50 CONTINUE
+ ENDIF
+C
+ RETURN
+ END
+
+
+
+ subroutine XYSYMB(N,X,Y,XOFF,XWT,YOFF,YWT,SH,ISYM)
+C.............................................................
+C
+C...GENERAL XY MULTIPLE-SYMBOL PLOTTING ROUTINE
+C (useful for overplotting XYLINE plot with point symbols)
+C
+C...INPUT X, Y Input arrays of length N
+C XOFF,XWT Offset and scale factor for X array
+C YOFF,YWT Offset and scale factor for Y array
+C SH Symbol size
+C ISYM Selects symbol type
+C if ISYM < 0 ... no plotting
+C.............................................................
+C
+ DIMENSION X(N), Y(N)
+C
+ IF(ISYM.LT.0) RETURN
+C
+ DO 10 I=1, N
+ XPLT = XWT*(X(I)-XOFF)
+ YPLT = YWT*(Y(I)-YOFF)
+ CALL PLSYMB(XPLT,YPLT,SH,ISYM,0.0,0)
+ 10 CONTINUE
+C
+ RETURN
+ END
+
+
+ subroutine CONT_GRID(IX,JX,II,JJ,X,Y,F,FCON,XOFF,YOFF,XWT,YWT)
+ DIMENSION X(IX,JX), Y(IX,JX), F(IX,JX)
+ CALL CONTGRID(IX,JX,II,JJ,X,Y,F,FCON,XOFF,YOFF,XWT,YWT)
+ RETURN
+ END
+
+ subroutine CONTGRID(IX,JX,II,JJ,X,Y,F,FCON,XOFF,YOFF,XWT,YWT)
+ DIMENSION X(IX,JX), Y(IX,JX), F(IX,JX)
+C--------------------------------------------------------------------------
+C
+C Plots one contour of a function F on a logically rectangular grid.
+C (normally called repeatedly if a number of contours is to be drawn)
+C
+C IX JX dimensions of arrays X, Y, F
+C II JJ array limits of arrays X, Y, F
+C X(i,j) independent coordinates of point (i,j)
+C Y(i,j)
+C F(i,j) function value at point (i,j)
+C FCON value of F on the contour to be drawn
+C XOFF offset for X
+C YOFF offset for Y
+C XWT scaling factor for X
+C YWT scaling factor for Y
+C
+C XPLOT = (X - XOFF)*XWT
+C YPLOT = (Y - YOFF)*YWT
+C--------------------------------------------------------------------------
+C
+ LOGICAL FOUND
+C
+C---- go over all cells and draw contour in any cell which contains the contour
+ DO 10 IO=1, II-1
+ IP = IO+1
+C
+ DO 110 JO=1, JJ-1
+ JP = JO+1
+C
+ FOUND = .FALSE.
+C
+C op 3 pp
+C
+C 4 2
+C
+C oo 1 po
+C
+ XOO = X(IO,JO)
+ XOP = X(IO,JP)
+ XPO = X(IP,JO)
+ XPP = X(IP,JP)
+C
+ YOO = Y(IO,JO)
+ YOP = Y(IO,JP)
+ YPO = Y(IP,JO)
+ YPP = Y(IP,JP)
+C
+ FOO = F(IO,JO)
+ FOP = F(IO,JP)
+ FPO = F(IP,JO)
+ FPP = F(IP,JP)
+C
+C-------- bottom edge (side 1)
+ IF(FCON.GE.FOO .AND. FCON.LT.FPO .OR.
+ & FCON.LT.FOO .AND. FCON.GE.FPO ) THEN
+ XCON = XOO + (FCON-FOO)*(XPO-XOO)/(FPO-FOO)
+ YCON = YOO + (FCON-FOO)*(YPO-YOO)/(FPO-FOO)
+ IF(FOUND) THEN
+ CALL PLOT(XWT*(XCON-XOFF),YWT*(YCON-YOFF),2)
+ ELSE
+ CALL PLOT(XWT*(XCON-XOFF),YWT*(YCON-YOFF),3)
+ ENDIF
+ FOUND = .NOT.FOUND
+ ENDIF
+C
+C-------- left edge (side 4)
+ IF(FCON.GE.FOO .AND. FCON.LT.FOP .OR.
+ & FCON.LT.FOO .AND. FCON.GE.FOP ) THEN
+ XCON = XOO + (FCON-FOO)*(XOP-XOO)/(FOP-FOO)
+ YCON = YOO + (FCON-FOO)*(YOP-YOO)/(FOP-FOO)
+ IF(FOUND) THEN
+ CALL PLOT(XWT*(XCON-XOFF),YWT*(YCON-YOFF),2)
+ ELSE
+ CALL PLOT(XWT*(XCON-XOFF),YWT*(YCON-YOFF),3)
+ ENDIF
+ FOUND = .NOT.FOUND
+ ENDIF
+C
+C-------- right edge (side 2)
+ IF(FCON.GE.FPO .AND. FCON.LT.FPP .OR.
+ & FCON.LT.FPO .AND. FCON.GE.FPP ) THEN
+ XCON = XPO + (FCON-FPO)*(XPP-XPO)/(FPP-FPO)
+ YCON = YPO + (FCON-FPO)*(YPP-YPO)/(FPP-FPO)
+ IF(FOUND) THEN
+ CALL PLOT(XWT*(XCON-XOFF),YWT*(YCON-YOFF),2)
+ ELSE
+ CALL PLOT(XWT*(XCON-XOFF),YWT*(YCON-YOFF),3)
+ ENDIF
+ FOUND = .NOT.FOUND
+ ENDIF
+C
+C-------- top edge (side 3)
+ IF(FCON.GE.FOP .AND. FCON.LT.FPP .OR.
+ & FCON.LT.FOP .AND. FCON.GE.FPP ) THEN
+ XCON = XOP + (FCON-FOP)*(XPP-XOP)/(FPP-FOP)
+ YCON = YOP + (FCON-FOP)*(YPP-YOP)/(FPP-FOP)
+ IF(FOUND) THEN
+ CALL PLOT(XWT*(XCON-XOFF),YWT*(YCON-YOFF),2)
+ ELSE
+ CALL PLOT(XWT*(XCON-XOFF),YWT*(YCON-YOFF),3)
+ ENDIF
+ FOUND = .NOT.FOUND
+ ENDIF
+C
+ 110 CONTINUE
+ 10 CONTINUE
+C
+ RETURN
+ END
+
+
+
+
+ subroutine CONTQUAD(X,Y,F,
+ & FUPR,FLWR,
+ & NCU,XCU,YCU,
+ & NCL,XCL,YCL,
+ & NA,NE,NV,XP,YP)
+C--------------------------------------------------------------------------
+C Contour a single quadrilateral element between an upper and
+C a lower contour limit. The output from this routine is both the
+C line segments defining the upper and lower contours and a set of
+C polygons that define the area of the element lying between the
+C contour limits. These may be used to shade the area on a plot.
+C
+C X, Y Arrays containing quadrilateral coordinate data
+C F Array containing the quantity F(x,y) to be contoured,
+C ie. F(i) is defined at X(i),Y(i), i=1->4
+C FUPR, FLWR Upper and lower contour limits
+C NCU Number of upper limit contour line points
+C XCU, YCU Arrays of x,y points (in pairs) on upper contour
+C NCL Number of lower limit contour line points
+C XCL, YCL Arrays of x,y points (in pairs) on lower contour
+C NA Number of polygon areas in XP, YP, NE arrays
+C NE Array of numbers of points in each polygon area
+C NV Vertex count of points in XP,YP arrays
+C XP, YP Coordinate points of contour polygons
+C
+C Note: the output contour line points on the upper contour (XCU,YCU)
+C or lower contour (XCL,YCL) are only valid in point pairs, ie.
+C points 1 and 2 define a segment, points 3 and 4 define the next,
+C etc. In general, there is no guarantee that the cross-pair points
+C (like points 2 and 3) will be contiguous. DO NOT PLOT THESE AS A
+C CONTIGUOUS ARRAY OF POINTS.
+C
+C Note: the output areas (NA polygon areas) are stored with the points
+C for all polygons in one big XP,YP array. The number of points
+C in each polygon are stored in NE, i.e. NE(1) is the number of
+C points stored in XP,YP for the first polygon, NE(2) is the number
+C of vertices stored following those for the second polygon, etc.
+C The total number of vertices is NV=sum[NE(1)+NE(2)..+NE(NA)].
+C
+C Note: the NCU,NCL,NA,NV counters are cumulative in this routine!
+C If you want to contour each quadrilateral without accumulating
+C contour points or polygons reset NCU=0,NCL=0,NA=0,NV=0 before
+C each call to CONTQUAD.
+C
+C--------------------------------------------------------------------------
+ DIMENSION X(4), Y(4), F(4)
+ DIMENSION XP(1), YP(1), NE(1)
+ DIMENSION XCL(1), YCL(1), XCU(1), YCU(1)
+C
+ DIMENSION IANG(4)
+ DIMENSION XTMP(3), YTMP(3), FTMP(3)
+C
+C---- NCU tracks the number of upper contour points
+C NCL tracks the number of lower contour points
+C NA tracks the number of polygon areas
+C NV tracks the total number of stored vertices in the XP,YP arrays
+C---- Uncomment these if you want to reset the counters each time you
+C contour a quadrilateral
+C NA = 0
+C NV = 0
+C NCU = 0
+C NCL = 0
+C
+ FHI = FUPR
+ FLO = FLWR
+ IF (FHI.LT.FLO) THEN
+ FHI = FLWR
+ FLO = FUPR
+ ENDIF
+C
+C---- Extrema
+ FMAX = AMAX1(F(1),F(2),F(3),F(4))
+ FMIN = AMIN1(F(1),F(2),F(3),F(4))
+C
+C---- If cell is above contour band or below contour band skip it
+ IF (FMAX.LE.FLO .OR. FMIN.GE.FHI) RETURN
+C
+C---- If cell is totally within contour band there are no contour lines,
+C and cell can be shaded directly
+ IF (FMAX.LE.FHI .AND. FMIN.GE.FLO) THEN
+ DO 2 I = 1, 4
+ NV = NV + 1
+ XP(NV) = X(I)
+ YP(NV) = Y(I)
+ 2 CONTINUE
+ NA = NA + 1
+ NE(NA) = 4
+ GO TO 100
+ ENDIF
+C
+C----Check for convex or concave quadrilaterals
+ ISUM = 0
+ DO 3 J = 1, 4
+ JM = MOD(J+2,4) + 1
+ JP = MOD(J, 4) + 1
+ IANG(J) = 1
+ IF ( (X(J)-X(JM))*(Y(JP)-Y(J)) .LT.
+ & (X(JP)-X(J))*(Y(J)-Y(JM)) ) IANG(J) = -1
+ ISUM = ISUM + IANG(J)
+ 3 CONTINUE
+ DO 4 J = 1, 4
+ IF (IANG(J)*ISUM.LT.0) GO TO 10
+ 4 CONTINUE
+C
+C----All angles < 180 deg., split into 4 triangles with average center pt
+ XTMP(3) = 0.25*(X(1)+X(2)+X(3)+X(4))
+ YTMP(3) = 0.25*(Y(1)+Y(2)+Y(3)+Y(4))
+ FTMP(3) = 0.25*(F(1)+F(2)+F(3)+F(4))
+C
+ DO 5 I = 1, 4
+ IP = MOD(I,4) + 1
+ XTMP(1) = X(I)
+ YTMP(1) = Y(I)
+ FTMP(1) = F(I)
+ XTMP(2) = X(IP)
+ YTMP(2) = Y(IP)
+ FTMP(2) = F(IP)
+ CALL CONTTRI(XTMP,YTMP,FTMP,FHI,FLO,
+ & NCU,XCU,YCU,NCL,XCL,YCL,NA,NE,NV,XP,YP)
+ 5 CONTINUE
+ GO TO 100
+C
+C----Quadrilaterals with an angle > 180, two triangles
+ 10 XTMP(3) = X(J)
+ YTMP(3) = Y(J)
+ FTMP(3) = F(J)
+C
+ JP1 = MOD(J,4) + 1
+ XTMP(1) = X(JP1)
+ YTMP(1) = Y(JP1)
+ FTMP(1) = F(JP1)
+ JP2 = MOD(JP1,4) + 1
+ XTMP(2) = X(JP2)
+ YTMP(2) = Y(JP2)
+ FTMP(2) = F(JP2)
+ CALL CONTTRI(XTMP,YTMP,FTMP,FHI,FLO,
+ & NCU,XCU,YCU,NCL,XCL,YCL,NA,NE,NV,XP,YP)
+C
+ JM2 = MOD(JM+1,4) + 1
+ XTMP(1) = X(JM2)
+ YTMP(1) = Y(JM2)
+ FTMP(1) = F(JM2)
+ JM1 = MOD(JM2,4) + 1
+ XTMP(2) = X(JM1)
+ YTMP(2) = Y(JM1)
+ FTMP(2) = F(JM1)
+ CALL CONTTRI(XTMP,YTMP,FTMP,FHI,FLO,
+ & NCU,XCU,YCU,NCL,XCL,YCL,NA,NE,NV,XP,YP)
+C
+ 100 RETURN
+ END
+
+
+
+ subroutine CONTTRI(X,Y,F,
+ & FUPR,FLWR,
+ & NCU,XCU,YCU,
+ & NCL,XCL,YCL,
+ & NA,NE,NV,XP,YP)
+C
+C Contour a single triangular element between an upper and
+C a lower contour limit. The output from this routine is both the line
+C segments defining the upper and lower contours and a set of
+C polygons that define the area of the element lying between the
+C contour limits. These may be used to shade the area on a plot.
+C
+C X, Y, Arrays containing triangular element points
+C F Array containing the quantity F(x,y) to be contoured,
+C ie. F(i) is defined at X(i),Y(i), i=1->3
+C FUPR, FLWR Upper and lower contour limits
+C NCU Number of upper limit contour line points
+C XCU, YCU Arrays of x,y points on upper contour
+C NCL Number of lower limit contour line points
+C XCL, YCL Arrays of x,y points on lower contour
+C NA Number of polygon areas in XP, YP, NE arrays
+C NE Array of numbers of points in each polygon area
+C NV Vertex count of points in XP,YP arrays
+C XP, YP Coordinate points of contour polygons
+C
+C Note: the output contour line points on the upper contour (XCU,YCU)
+C or lower contour (XCL,YCL) are only valid in point pairs, ie.
+C points 1 and 2 define a segment, points 3 and 4 define the next,
+C etc. In general, there is no guarantee that the cross-pair points
+C (like points 2 and 3) will be contiguous. DO NOT PLOT THESE AS A
+C CONTIGUOUS ARRAY OF POINTS.
+C
+C Note: the output areas (NA polygons) are stored with the points
+C for all polys in one big XP,YP array. The number of points
+C in each polygon are stored in NE, i.e. NE(1) is the number of
+C points stored in XP,YP for the first polygon, NE(2) is the number
+C of vertices stored following those for the second polygon, etc.
+C The total number of vertices is NV=sum[NE(1)+NE(2)..+NE(NA)].
+C
+C Note: the NCU,NCL,NA,NV counters are cumulative in this routine!
+C If you want to contour each triangle without accumulating
+C contour points or polygons reset NCU=0,NCL=0,NA=0,NV=0 before
+C each call to CONTTRI.
+C
+ DIMENSION X(3), Y(3), F(3), FH(3), FL(3),
+ & XCU(1), YCU(1), XCL(1), YCL(1), XP(1), YP(1), NE(1)
+C
+C---- NCU tracks the number of upper contour points
+C NCL tracks the number of lower contour points
+C NA tracks the number of polygon areas
+C NV tracks the total number of stored vertices in the XP,YP arrays
+C---- Uncomment these if you want to reset the counters each time you
+C contour a triangle
+C NA = 0
+C NV = 0
+C NCU = 0
+C NCL = 0
+C
+ FHI = FUPR
+ FLO = FLWR
+ IF (FHI.LT.FLO) THEN
+ FHI = FLWR
+ FLO = FUPR
+ ENDIF
+
+ NVFRST = NV
+C
+ EPS = 0.0001*(FHI-FLO)
+C
+C---- Temporary values
+ DO 1 I = 1, 3
+ FH(I) = F(I) - FHI
+ FL(I) = F(I) - FLO
+ IF (FH(I).GE.0. .AND. FH(I).LT. EPS) FH(I) = EPS
+ IF (FH(I).LE.0. .AND. FH(I).GT.-EPS) FH(I) = -EPS
+ IF (FL(I).GE.0. .AND. FL(I).LT. EPS) FL(I) = EPS
+ IF (FL(I).LE.0. .AND. FL(I).GT.-EPS) FL(I) = -EPS
+ 1 CONTINUE
+C
+C----Check point by point for points in contour limits
+ DO 50 I = 1, 3
+C
+C----Inside contour limits
+ IF (FH(I).LT.0. .AND. FL(I).GT.0.) THEN
+ NV = NV + 1
+ XP(NV) = X(I)
+ YP(NV) = Y(I)
+C
+ ELSE
+C
+ IF (FH(I).GE.0.) THEN
+C----Check for intersections with previous and next point
+ IP = MOD(I,3) + 1
+ IM = MOD(I+1,3) + 1
+ IF (FH(IM).LT.0.) THEN
+ ETA = -(FH(IM)+FH(I))/(FH(IM)-FH(I))
+ NCU = NCU + 1
+ XCU(NCU) = 0.5*(X(I) + X(IM) + ETA*(X(IM)-X(I)) )
+ YCU(NCU) = 0.5*(Y(I) + Y(IM) + ETA*(Y(IM)-Y(I)) )
+ NV = NV + 1
+ XP(NV) = XCU(NCU)
+ YP(NV) = YCU(NCU)
+ ENDIF
+ IF (FH(IP).LT.0.) THEN
+ ETA = -(FH(IP)+FH(I))/(FH(IP)-FH(I))
+ NCU = NCU + 1
+ XCU(NCU) = 0.5*(X(I) + X(IP) + ETA*(X(IP)-X(I)) )
+ YCU(NCU) = 0.5*(Y(I) + Y(IP) + ETA*(Y(IP)-Y(I)) )
+ NV = NV + 1
+ XP(NV) = XCU(NCU)
+ YP(NV) = YCU(NCU)
+ ENDIF
+ ENDIF
+C
+ IF (FL(I).LE.0.) THEN
+C----Check for intersections with previous and next point
+ IP = MOD(I,3) + 1
+ IM = MOD(I+1,3) + 1
+ IF (FL(IM).GT.0.) THEN
+ ETA = -(FL(IM)+FL(I))/(FL(IM)-FL(I))
+ NCL = NCL + 1
+ XCL(NCL) = 0.5*(X(I) + X(IM) + ETA*(X(IM)-X(I)) )
+ YCL(NCL) = 0.5*(Y(I) + Y(IM) + ETA*(Y(IM)-Y(I)) )
+ NV = NV + 1
+ XP(NV) = XCL(NCL)
+ YP(NV) = YCL(NCL)
+ ENDIF
+ IF (FL(IP).GT.0.) THEN
+ ETA = -(FL(IP)+FL(I))/(FL(IP)-FL(I))
+ NCL = NCL + 1
+ XCL(NCL) = 0.5*(X(I) + X(IP) + ETA*(X(IP)-X(I)) )
+ YCL(NCL) = 0.5*(Y(I) + Y(IP) + ETA*(Y(IP)-Y(I)) )
+ NV = NV + 1
+ XP(NV) = XCL(NCL)
+ YP(NV) = YCL(NCL)
+ ENDIF
+ ENDIF
+C
+ ENDIF
+C
+ 50 CONTINUE
+C
+ IF (NV.GT.NVFRST+2) THEN
+ NA = NA + 1
+ NE(NA) = NV - NVFRST
+ ENDIF
+C
+ RETURN
+ END
+
+
+
+ subroutine AXISADJ(xmin,xmax,xspan,deltax,ntics)
+C...Make scaled axes with engineering increments between tics
+C
+C Input: xmin, xmax - input range for which scaled axis is desired
+C
+C Output: xmin, xmax - adjusted range for scaled axis
+C xspan - adjusted span of scaled axis
+C deltax - increment to be used for scaled axis
+C nincr - number of tics to be used on axis
+C note that ntics=1+(xspan/deltax)
+C
+ real xmin,xmax,xspan,deltax,xinc,xinctbl(5)
+ integer ntics,i
+ data xinctbl / 0.1, 0.2, 0.25, 0.5, 1. /
+c
+ xspan1 = xmax-xmin
+ if (xspan1.eq.0.) xspan1 = 1.
+c
+ xpon = ifix(log10(xspan1))
+ xspan = xspan1 / 10.**xpon
+c
+ do i = 1, 5
+ xinc = xinctbl(i)
+ ntics = 1 + ifix(xspan/xinc + 0.1)
+ if (ntics.LE.6) go to 1
+ end do
+c
+ 1 deltax = xinc*10.**xpon
+ xmin = deltax* ifloor(xmin/deltax)
+ xmax = deltax*iceiling(xmax/deltax)
+ xspan = xmax - xmin
+ ntics = 1 + ifix(xspan/deltax + 0.1)
+ return
+ end
+
+ function iceiling(x)
+c--- returns next highest integer value if fraction is non-zero
+ integer iceiling
+ real x
+ i = ifix(x)
+ if(x-i.GT.0.) i = i+1
+ iceiling = i
+ return
+ end
+
+ function ifloor(x)
+c--- returns next lowest integer value if fraction is negative, non-zero
+ integer ifloor
+ real x
+ i = ifix(x)
+ if(x-i.LT.0.) i = i-1
+ ifloor = i
+ return
+ end
+
+
+
+
+
+ subroutine ANNOT(CH)
+C------------------------------------------------------
+C Interactive annotation menu for adding custom
+C ornaments to an active plot (before PLEND call).
+C------------------------------------------------------
+ CHARACTER*80 AA
+ CHARACTER*1 OPT, KCHAR
+C
+ SAVE CHF, ISYMB
+ DATA CHF, ISYMB / 1.0, 0 /
+C
+ 900 CONTINUE
+C
+ 1000 FORMAT(A)
+ 1010 FORMAT(A,$)
+ 1020 FORMAT(A,F7.3,A,$)
+ 1030 FORMAT(A,I2 ,A,$)
+C
+ WRITE(*,1050)
+ 1050 FORMAT(/' C haracters | '
+ & /' S lant characters | '
+ & /' M ath characters | plot '
+ & /' P oint symbol | '
+ & /' L ine | '
+ & /' A rrow | '
+ & /' '
+ & /' W idth of characters | modify'
+ & /' T ype of point symbol | ')
+C
+ 905 WRITE(*,*)
+ WRITE(*,1010) ' Select option or <return>: '
+ READ(*,1000) OPT
+ IF(OPT.EQ.' ') RETURN
+C
+ CHI = CHF*CH
+C
+C-------------------------------------------------------------
+ IF(INDEX('CcSsMm',OPT).NE.0) THEN
+C
+ WRITE(*,*) 'Click on lower left point of character string...'
+ CALL GETCURSORXY(XX,YY,KCHAR)
+ WRITE(*,1010) ' Enter character string: '
+ READ (*,1000) AA
+C
+C---- find index of last non-blank character
+ DO 112 NA=80, 1, -1
+ IF(AA(NA:NA).NE.' ') GO TO 113
+ 112 CONTINUE
+ 113 CONTINUE
+C
+ CALL NEWPEN(3)
+ IF(INDEX('Cc',OPT).NE.0) CALL PLCHAR(XX,YY,CHI,AA,0.0,NA)
+ IF(INDEX('Ss',OPT).NE.0) CALL PLSLAN(XX,YY,CHI,AA,0.0,NA)
+ IF(INDEX('Mm',OPT).NE.0) CALL PLMATH(XX,YY,CHI,AA,0.0,NA)
+ CALL PLFLUSH
+C
+C-------------------------------------------------------------
+ ELSE IF(INDEX('Pp',OPT).NE.0) THEN
+C
+ WRITE(*,*) 'Click on symbol locations ...'
+ CALL GETCURSORXY(XX,YY,KCHAR)
+ CALL NEWPEN(2)
+ CALL PLSYMB(XX,YY,CHI,ISYMB,0.0,0)
+ CALL PLFLUSH
+C
+C-------------------------------------------------------------
+ ELSE IF(INDEX('LlAa',OPT).NE.0) THEN
+C
+ WRITE(*,*) 'Click on line points, twice on last point...'
+C
+ CALL NEWPEN(1)
+ CALL GETCURSORXY(XXM,YYM,KCHAR)
+ CALL PLOT(XXM,YYM,3)
+ CALL PLOT(XXM,YYM,2)
+ XXL = XXM
+ YYL = YYM
+ CALL PLFLUSH
+ DO 131 IP=1, 12345
+ CALL GETCURSORXY(XX,YY,KCHAR)
+ CALL PLOT(XX,YY,2)
+
+ CALL PLFLUSH
+
+ IF(XXM.EQ.XX .AND. YYM.EQ.YY) GO TO 132
+ XXL = XXM
+ YYL = YYM
+ XXM = XX
+ YYM = YY
+ 131 CONTINUE
+ 132 CONTINUE
+C
+ IF(INDEX('Aa',OPT).NE.0) THEN
+C------ add arrowhead
+ DX = XX - XXL
+ DY = YY - YYL
+ DS = SQRT(DX**2 + DY**2)
+ ARLEN = 1.5*CHI
+ HAR = 0.1
+cc IF(DS .GT. ARLEN) THEN
+ CALL PLOT(XX-ARLEN*(DX+HAR*DY)/DS,YY-ARLEN*(DY-HAR*DX)/DS,2)
+ CALL PLOT(XX-ARLEN*(DX-HAR*DY)/DS,YY-ARLEN*(DY+HAR*DX)/DS,2)
+ CALL PLOT(XX,YY,2)
+cc ENDIF
+ ENDIF
+C
+ XX = XX + 0.7*CHI
+ YY = YY - 0.5*CHI
+ CALL PLFLUSH
+C
+C-------------------------------------------------------------
+ ELSE IF(INDEX('Ww',OPT).NE.0) THEN
+C
+ 140 WRITE(*,1020)
+ & ' Enter new character width factor (currently =',CHF,'): '
+ READ (*,*,ERR=140) CHF
+C
+C-------------------------------------------------------------
+ ELSE IF(INDEX('Tt',OPT).NE.0) THEN
+C
+ WRITE(*,*)
+ WRITE(*,*) ' 0 square 7 Y '
+ WRITE(*,*) ' 1 circle 8 flipped Y '
+ WRITE(*,*) ' 2 triangle 9 * '
+ WRITE(*,*) ' 3 + 10 flipped * '
+ WRITE(*,*) ' 4 x 11 hourglass '
+ WRITE(*,*) ' 5 diamond 12 bowtie '
+ WRITE(*,*) ' 6 yield sign 13 star '
+ WRITE(*,*)
+ 160 WRITE(*,1030) ' Enter new symbol type (currently =',ISYMB,'): '
+ READ (*,*,ERR=160) ISYMB
+C
+C-------------------------------------------------------------
+ ELSE
+C
+ GO TO 900
+C
+ ENDIF
+C
+ GO TO 905
+ END ! ANNOT