
      FUNCTION     CARG (Z)

C

C      [COMPLEX ARGUMENT]
C     [05-MAR-74]

      COMPLEX      Z

      CARG=ATAN2(AIMAG(Z),REAL(Z))
      RETURN

      END

      SUBROUTINE  KONIT (I,J,K)
 
C     [INITIAL TRIANGLE]
C     [16-NOV-74]
 
      COMMON/KON/ M1(2),M2(2),M3(2),Z1,Z2,Z3
 
      GO TO (10,20),K
C

   10 M1(1)=I
      M1(2)=J
      M2(1)=I+1
      M2(2)=J
      M3(1)=I
      M3(2)=J+1
      RETURN
 
   20 M1(1)=I+1
      M1(2)=J+1
      M2(1)=I+1
      M2(2)=J
      M3(1)=I
      M3(2)=J+1
      RETURN
 
      END
      SUBROUTINE  KONNC
 
C     [NEXT, CONSTANT]
C     SELECT P3 TO DEFINE THE NEXT TRIANGLE ALONG A CONSTANT CONTOUR.
C     POINTS P1 AND P2 OF THE NEW TRIANGLE WILL BE POINTS OF THE OLD
C     TRIANGLE, WHILE P3 WILL BE A NEW POINT GOTTEN BY REFLECTION OF
C     THE MISSING POINT IN THE OPPOSITE EDGE.
C     [24-MAY-73]
 
      COMMON/KON/ M1(2),M2(2),M3(2),Z1,Z2,Z3
 
      IF (SIGN(1.0,Z1).EQ.SIGN(1.0,Z3)) GO TO 10
      I=M1(1)-M2(1)+M3(1)
      J=M1(2)-M2(2)+M3(2)
      CALL KONXV (2,3)
      M3(1)=I
      M3(2)=J
      RETURN
 
   10 IF (SIGN(1.0,Z2).EQ.SIGN(1.0,Z3)) RETURN
      I=-M1(1)+M2(1)+M3(1)
      J=-M1(2)+M2(2)+M3(2)
      CALL KONXV (1,3)
      M3(1)=I
      M3(2)=J
      RETURN
 
      END
      SUBROUTINE  KONRE
 
C     [RESTORE]
C     RESTORE THE INITIAL POINT OF THE CONTOUR
C     [14-FEB-74]
 
      COMMON/KON/ M(6),Z(3)
      COMMON/KQN/ N(6),W(3)
 
      DO 10 I=1,6
   10 M(I)=N(I)
      DO 20 I=1,3
   20 Z(I)=W(I)
      RETURN
 
      END
      SUBROUTINE  KONSA
 
C     [SAVE]
C     SAVE THE INITIAL POINT OF THE CONTOUR
C     [24-MAY-73]
 
      COMMON/KON/ M1(2),M2(2),M3(2),Z1,Z2,Z3
      COMMON/KQN/ N1(2),N2(2),N3(2),W1,W2,W3
 
      N1(1)=M1(1)
      N1(2)=M1(2)
      N2(1)=M3(1)
      N2(2)=M3(2)
      N3(1)=M2(1)
      N3(2)=M2(2)
      W1=Z1
      W2=Z3
      W3=Z2
      RETURN
 
      END
      SUBROUTINE  KONSC (Z0,XF,YF,IA,IB,JA,JB,ZE,NX,NY,PL)
 
C     [SINGLE CONTOUR]
C     A GENERAL PURPOSE SUBROUTINE WHICH MAY BE USED TO GENERATE SIMPLE
C     CONTOURS, OR CONTOURS OF APHIC RELIEF.
C     Z0	CONTOUR LEVEL SOUGHT
C     (XF,YF)	LIGHTING DIREATION FOR ORTHOGRAPHIC RELIEF
C     (IA,IB)	X-INTERVAL TO BE CONTOURED
C     (JA,JB)	Y-INTERVAL TO BED
C     ZE(NX,NY) ARRAY OF FUNCTION VALUES
C     PL	PEN MOVEMENT SUBROUTINE
C     [06-JAN-75]
 
      LOGICAL	  FE(35,35,2)
      DIMENSION   ZE(1)
      COMMON/KON/ I1,J1,I2,J2,I3,J3,Z1,Z2,Z3
 
      U(I,J)=ZE(I+NX*(J-1))-Z0+XF*FLOAT(I-1)+YF*FLOAT(J-1)
      ZP(I1,I2)=FLOAT(I1-1)-Z1*(FLOAT(I2-I1)/(Z2-Z1))
 
      IF ((IB-IA).GT.35) RETURN
      IF ((JB-JA).GT.35) RETURN
      XS=1.0/FLOAT(NX-1)
      YS=1.0/FLOAT(NY-1)
      II=MAX0(IA,IB-1)
      JJ=MAX0(JA,JB-1)
      DO 10 I=IA,II
      DO 10 J=JA,JJ
      Z11=U(I,J)
      Z12=U(I,J+1)
      Z21=U(I+1,J)
      Z22=U(I+1,J+1)
      ZP1=AMAX1(Z11,Z12,Z21)
      ZM1=AMIN1(Z11,Z12,Z21)
      ZP2=AMAX1(Z12,Z21,Z22)
      ZM2=AMIN1(Z12,Z21,Z22)
      FE(I-IA+1,J-JA+1,1)=(ZP1.LT.0.0).OR.(ZM1.GT.0.0)
   10 FE(I-IA+1,J-JA+1,2)=(ZP2.LT.0.0).OR.(ZM2.GT.0.0)
      DO 40 K=1,2
      DO 40 I=IA,II
      DO 40 J=JA,JJ
      IF (FE(I-IA+1,J-JA+1,K)) GO TO 40
      CALL KONIT (I,J,K)
      Z1=U(I1,J1)
      Z2=U(I2,J2)
      Z3=U(I3,J3)
      IF (SIGN(1.0,Z1).EQ.SIGN(1.0,Z2)) CALL KONXV (1,3)
      IF (SIGN(1.0,Z1).EQ.SIGN(1.0,Z3)) CALL KONXV (1,2)
      CALL KONSA
      DO 30 L=1,2
      CALL PL (XS*ZP(I1,I2),YS*ZP(J1,J2),.FALSE.)
   20 CALL KONNC
      CALL PL (XS*ZP(I1,I2),YS*ZP(J1,J2),.TRUE.)
      I0=MIN0(I1,I2,I3)-IA+1
      J0=MIN0(J1,J2,J3)-JA+1
      K0=MOD(I1+I2+I3,3)
      IF (FE(I0,J0,K0)) GO TO 30
      FE(I0,J0,K0)=.TRUE.
      IF ((I3.LT.IA).OR.(I3.GT.IB).OR.(J3.LT.JA).OR.(J3.GT.JB)) GO TO 30
      Z3=U(I3,J3)
      GO TO 20
   30 CALL KONRE
   40 CONTINUE
      RETURN
      END
      SUBROUTINE  KONSK (Z0,IA,IB,JA,JB,ZE,NX,NY,FU,PL)
 
C     [SINGLE COMPLEX CONTOUR]
C     Z0	CONTOUR LEVEL SOUGHT
C     (IA,IB)	X-INTERVAL TO BE CONTOURED
C     (JA,JB)	Y-INTERVAL TO BE CONTOURED
C     ZE(NX,NY) ARRAY OF FUNCTION VALUES
C     PL	PEN MOVEMENT SUBROUTINE
C     [16-NOV-74]
 
      LOGICAL	  FE(34,34,2)
      COMPLEX	  ZE(1)
      COMMON/KON/ I1,J1,I2,J2,I3,J3,Z1,Z2,Z3
 
      U(I,J)=FU(ZE(I+NX*(J-1)))-Z0
      ZP(I1,I2)=FLOAT(I1-1)-Z1*(FLOAT(I2-I1)/(Z2-Z1))
 
      XS=1.0/FLOAT(NX-1)
      YS=1.0/FLOAT(NY-1)
      II=MAX0(IA,IB-1)
      JJ=MAX0(JA,JB-1)
      DO 10 I=IA,II
      DO 10 J=JA,JJ
      Z11=U(I,J)
      Z12=U(I,J+1)
      Z21=U(I+1,J)
      Z22=U(I+1,J+1)
      ZP1=AMAX1(Z11,Z12,Z21)
      ZM1=AMIN1(Z11,Z12,Z21)
      ZP2=AMAX1(Z12,Z21,Z22)
      ZM2=AMIN1(Z12,Z21,Z22)
      FE(I-IA+1,J-JA+1,1)=(ZP1.LT.0.0).OR.(ZM1.GT.0.0)
   10 FE(I-IA+1,J-JA+1,2)=(ZP2.LT.0.0).OR.(ZM2.GT.0.0)
      DO 40 K=1,2
      DO 40 I=IA,II
      DO 40 J=JA,JJ
      IF (FE(I-IA+1,J-JA+1,K)) GO TO 40
      CALL KONIT (I,J,K)
      Z1=U(I1,J1)
      Z2=U(I2,J2)
      Z3=U(I3,J3)
      IF (SIGN(1.0,Z1).EQ.SIGN(1.0,Z2)) CALL KONXV (1,3)
      IF (SIGN(1.0,Z1).EQ.SIGN(1.0,Z3)) CALL KONXV (1,2)
      CALL KONSA
      DO 30 L=1,2
	CALL PL (XS*ZP(I1,I2),YS*ZP(J1,J2),.FALSE.)
   20 CALL KONNC
      CALL PL (XS*ZP(I1,I2),YS*ZP(J1,J2),.TRUE.)
      I0=MIN0(I1,I2,I3)-IA+1
      J0=MIN0(J1,J2,J3)-JA+1
      K0=MOD(I1+I2+I3,3)
      IF (FE(I0,J0,K0)) GO TO 30
      FE(I0,J0,K0)=.TRUE.
      IF ((I3.LT.IA).OR.(I3.GT.IB).OR.(J3.LT.JA).OR.(J3.GT.JB)) GO TO 30
      Z3=U(I3,J3)
      GO TO 20
   30 CALL KONRE
   40 CONTINUE
      RETURN
      END
      SUBROUTINE  KONXV (I,J)
 
C     [EXCHANGE VECTORS]
C     KONXV (I,J) EXCHANGES THE ITH AND JTH VECTORS IN COMMON.
C     [24-MAY-73]
 
      COMMON/KON/ MM(2,3),Z(3)
 
      DO 10 L=1,2
      N=MM(L,I)
      MM(L,I)=MM(L,J)
   10 MM(L,J)=N
      T=Z(I)
      Z(I)=Z(J)
      Z(J)=T
      RETURN
 
      END
      SUBROUTINE  PLT00
 
C     [INITIALIZATION]
C     INITIALIZING SUBROUTINE TO START OFF A SERIES OF GRAPHS.
C     CALLS <PLOTS>, THEN MOVES THE PEN AT MOST 11" TO THE RIGHT
C     TO INSURE ITS PROPER POSITIONING.
C     [18-FEB-73]
 
      CALL PLOTS (I)
      CALL PLOT  (0.0,-11.0,-3)
      RETURN
 
      END
      SUBROUTINE  PLTAX (X,Y,HE,NC,SZ,TH,V0,DV,L)
 
C     (X,Y)   POINT FROM WHICH AXIS ORIGINATES
C     HE      HEADING TO BE PLACED UNDER GRAPH
C     NC      NUMBER OF CHARACTERS IN HEADING
C     SZ      LENGTH OF AXIS, IN INCHES
C     TH      COUNTERCLOCKWISE ANGLE OF INCLINATION, DEGREES
C     V0      STARTING VALUE OF VARIABLE ALONG AXIS
C     DV      INCREMENT OF VARIABLE, PER INCH
C     L       =1, LETTERING ABOVE; =-1, LETTERING BELOW
C     [18-NOV-74]
 
      DIMENSION   HE(1)
 
      S=FLOAT(L)
      N=IFIX(SZ+0.5)
      CTH=COSD(TH)
      STH=SIND(TH)
      XB=X
      YB=Y
      XA=X-0.1*S*STH
      YA=Y+0.1*S*CTH
      CALL PLOT (YA,-XA,3)
      DO 20 I=1,N
      CALL PLOT (YB,-XB,2)
      XC=XB+CTH
      YC=YB+STH
      CALL PLOT (YC,-XC,2)
      XA=XA+CTH
      YA=YA+STH
      CALL PLOT (YA,-XA,2)
      XB=XC
   20 YB=YC
      IX=0
      NT=IFIX(ALOG10(DV)+0.001)
      IF (NT.LT.-1.OR.NT.GT.1) IX=NT
      ADV=DV*10.0**(-IX)
      ABV=V0*10.0**(-IX)+FLOAT(N)*ADV
      XA=XB-(0.20*S-0.05)*STH-0.0857*CTH
      YA=YB+(0.20*S-0.05)*CTH-0.0857*STH
      N=N+1
      DO 30 I=1,N
      CALL NUMBER (YA,-XA,0.1,ABV,TH-90.0,2)
      ABV=ABV-ADV
      XA=XA-CTH
   30 YA=YA-STH
      TA=FLOAT(NC+7)
      XA=X+(SZ/2.0-0.06*TA)*CTH-(-0.07+S*0.36)*STH
      YA=Y+(SZ/2.0-0.06*TA)*STH+(-0.07+S*0.36)*CTH
      IF (NC.NE.0) CALL SYMBOL (YA,-XA,0.12,HE,TH-90.0,NC)
      IF (IX.EQ.0) RETURN
      XA=XA+((TA-6.0)*0.12)*CTH
      YA=YA+((TA-6.0)*0.12)*STH
      CALL SYMBOL (YA,-XA,0.12,'(X10  )',TH-90.0,7)
      XA=XA+0.48*CTH-0.07*STH
      YA=YA+0.48*STH+0.07*CTH
      CALL NUMBER (YA,-XA,0.1,FLOAT(IX),TH-90.0,-1)
      RETURN
 
      END
      SUBROUTINE  PLTBH (X,Y,P)
 
C     [BOTTOM HALF]
C     SCALE THE CARTESIAN COORDINATES X,Y SO AS TO PLACE A GRAPH
C     IN THE LOWER HALF OF THE PLOTTER PAGE.
C     [20-APR-74]
 
      LOGICAL	  P
      DATA	  HX,HY/4.50,3.25/
 
      CALL PLTMS (HX*(Y-1.0),2.0*HY*(0.5-X),P)
      RETURN
 
      END
      SUBROUTINE  PLTBO
 
C     [BORDER]
C     SET UP AN 8-1/2" X 11" FRAME WITH AN INNER FRAME 1" INSIDE
C     OF IT AND DEFINE THE ORIGIN AT THE PAGE CENTER.
C     [15-NOV-73]
 
      DIMENSION   IH(10),IJOB(3),IDATE(2)
      EQUIVALENCE (IJOB(1),IH(3))
      EQUIVALENCE (IDATE(1),IH(7))
      EQUIVALENCE (ITIME,IH(10))
 
      IH(1)='ESFM:'
      IH(2)='	  '
      IH(6)='	  '
      IH(9)='	  '
      CALL PLOT   (0.0, 0.0, 3)
      CALL PLOT   (0.0,11.0, 2)
      CALL PLOT   (8.5,11.0, 1)
      CALL PLOT   (8.5, 0.0, 1)
      CALL PLOT   (0.0, 0.0, 1)
      CALL SYSJO  (IJOB)
      CALL DATE   (IDATE)
      CALL TIME   (ITIME)
      CALL SYMBOL (0.1,4.5,0.08,IH,-90.0,50)
      CALL PLOT   (1.0, 1.0, 3)
      CALL PLOT   (1.0,10.0, 2)
      CALL PLOT   (7.5,10.0, 1)
      CALL PLOT   (7.5, 1.0, 1)
      CALL PLOT   (1.0, 1.0, 1)
      CALL PLOT   (4.25,5.50,-3)
      RETURN
 
      END
      SUBROUTINE  PLTBS
 
C     [BACK SPACE]
C     PARTICULARLY FOR MAKING COLOR COMPOSITES, IT IS SOMETIMES REQUIRED
C     TO NEGATE THE EFFECT OF PLTEJ IN SUCH A WAY THAT AN INTERMEDIATE
C     PLT00 CAN BE EXECUTED.  AT THE SAME TIME, THE PEN CREEP OCCASIONED
C     BY THE PLOTTER SPOOLER CAN BE NULLIFIED. THEREFORE, PLTBS MUST NOT
C     BE USED WITHOUT THE PLOTTER SPOOLER.  SINCE IT DRAWS NO MARGINS,
C     IT AVOIDS SUPERIMPOSING COLORED VERSIONS OF THE IDENTIFICATION.
C     [11-JAN-75]
 
      DATA	  SX,SY/5.50,4.25/
 
      CALL PLOT (-SY-0.02,SX,-3)
      RETURN
 
      END
      SUBROUTINE  PLTBV (Z1,ZE,Z2,NX,MX,NY)
 
C     [BIRDSEYE VIEW]
C     [18-DEC-74]
 
      DIMENSION   ZE(1)
      DATA	  HX,HY/4.50,3.25/
 
      F(I,J)=(R*ZS)/(Z0-ZE(I+MX*(J-1)))
 
      K=1
      R=5.0
      Z0=1.3*Z2
      ZS=0.125*(Z2-Z1)
      DX=(1.75*HX)/FLOAT(MX-1)
      DY=(1.75*HY)/FLOAT(NY-1)
      X=-0.875*HX
      Y=-0.875*HY
      DO 20 J=1,NY
      I1=((NX+1)-K*(NX-1))/2
      I2=((NX+1)+K*(NX-1))/2
      EF=F(I1,J)
      CALL PLTMS (EF*X,EF*Y,.FALSE.)
      DO 10 I=I1,I2,K
      EF=F(I,J)
      CALL PLTMS (EF*X,EF*Y,.TRUE.)
   10 X=X+DX
      DX=-DX
      X=X+DX
      K=-K
   20 Y=Y+DY
      DX=-(1.75*HX)/FLOAT(MX-1)
      DY=-(1.75*HY)/FLOAT(NY-1)
      X= 0.875*HX
      Y= 0.875*HY
      K=-1
      DO 40 I=1,MX
      J1=((NY+1)-K*(NY-1))/2
      J2=((NY+1)+K*(NY-1))/2
      EF=F(MX-I+1,J1)
      CALL PLTMS (EF*X,EF*Y,.FALSE.)
      DO 30 J=J1,J2,K
      EF=F(MX-I+1,J)
      CALL PLTMS (EF*X,EF*Y,.TRUE.)
   30 Y=Y+DY
      DY=-DY
      Y=Y+DY
      K=-K
   40 X=X+DX
      RETURN
 
      END
      SUBROUTINE  PLTCA (X,Y,P)
 
C     [CARTESIAN]
C     SCALE (X,Y) TO THE PAGE WIDTH, ALLOWING THE COORDINATES TO BE
C     GENERATED IN THE INTERVAL (0.0.LE.X,Y.LE.1.0).
C     [12-FEB-74]
 
      LOGICAL	  P
      DATA	  HX,HY/4.50,3.25/
 
      EX=2.0*HX*X-HX
      WY=2.0*HY*Y-HY
      CALL PLTMS (EX,WY,P)
      RETURN
 
      END
      SUBROUTINE  PLTCI (X,Y,R)
 
C     [CIRCLE]
C     DRAW A CIRCLE ON THE PLOTTER WITH CENTER (X,Y), RADIUS R.
C     [28-SEP-73]
 
      N=60
      DT=6.28318/FLOAT(N)
      TH=DT
      CALL PLTMC (X+R,Y,.FALSE.)
      DO 10 I=1,N
      CALL PLTMC (X+R*COS(TH),Y+R*SIN(TH),.TRUE.)
   10 TH=TH+DT
 
      RETURN
      END
      SUBROUTINE  PLTEJ
 
C     [EJECT]
C     EJECT A PAGE ON THE PLOTTER, SUPPOSING THE ORIGIN AT PAGE CENTER.
C     [05-OCT-73]
 
      DATA	  SX,SY/5.50,4.25/
 
      CALL PLOT  (SY,-SX,-3)
 
      RETURN
      END
      SUBROUTINE  PLTEL (XI,ETA,P)
 
C     [ELLIPTICAL]
C     CHANGE (XI,ETA) INTO (X,Y) SO THAT PEN MOVEMENTS CAN BE
C     DEFINED DIRECTLY IN ELLIPTICAL COORDINATES.
C     0.0 .LE. XI  .LE. 1.0
C     0.0 .LE. ETA .LE. 1.0
C     S*COSH(XI=1)=HX
C     1.76=ARCCOSH(3.0)
C     [14-FEB-74]
 
      LOGICAL	  P
      DATA	  HX,HY/4.50,3.25/
 
      S=1.50
      E=6.28318*ETA
      X=XI*1.76
      CALL PLTMS (S*COSH(XI)*COS(E),S*SINH(XI)*SIN(E),P)
      RETURN
 
      END
      SUBROUTINE  PLTEV (Z1,ZE,Z2,NX,NE)
 
C     [ELLIPTICAL VIEW]
C     PROGRAM TO PRODUCE A PERSPECTIVE DRAWING OF A SINGLE VALUED
C     FUNCTION DEFINED OVER ELLIPTICAL COORDINATES, IN SUCH A WAY
C     AS TO EXHIBIT THE ARCS CORRESPONDING TO CONSTANT XI AND ETA.
C     X=COSH(XI)*COS(ETA)
C     Y=SINH(XI)*SIN(ETA)
C     (Z1,Z2)	RANGE OF Z VALUES
C     ZE(NX,NE) ARRAY OF FUNCTION VALUES
C     NX	NUMBER OF XI VALUES
C     NE	NUMBER (=4*N+1) OF ETA VALUES
C     [10-NOV-74]
 
      DIMENSION   ZE(1)
      DATA	  Q0,Q1,Q2,Q3,Q4/0.000,1.571,3.142,4.713,6.283/
      DATA	  S1,S3/1.570,4.712/
      DATA	  X1,X2/0.01,1.76/
 
      NQ=1+(NE-1)/4
      NN=NX*(NQ-1)
      CALL PLTFR
      CALL VISNH
      CALL VISES (Z1,ZE(3*NN+1),Z2,X1,X2,NX,Q3,Q4,NQ,-1, 1)
      CALL VISES (Z1,ZE        ,Z2,X1,X2,NX,Q0,S1,NQ, 1, 1)
      CALL VISNH
      CALL VISES (Z1,ZE(2*NN+1),Z2,X1,X2,NX,Q2,S3,NQ,-1,-1)
      CALL VISES (Z1,ZE(  NN+1),Z2,X1,X2,NX,Q1,Q2,NQ, 1,-1)
      CALL PLTEJ
      RETURN
 
      END
      SUBROUTINE  PLTFI (Y1,WY,Y2,N)
 
C     [FUNCTION OF INTEGERS]
C     PLOT A GRAPH IN RECTANGULAR COORDINATES, BY CONNECTING SUCCESSIVE
C     DATA POINTS BY STRAIGHT LINES. THE POINTS DEFINING THE GRAPH ARE
C     TAKEN FROM AN ARRAY OF Y-VALUES. THE X-VALUES ARE INTEGERS, LYING
C     BETWEEN 1 AND N. THE RESPECTIVE SCALES ARE INDICATED BY THE VALUES
C     TO BE ASSIGNED TO THE MARGINS OF THE GRAPH. ORDINARILY THE MARGINS
C     WOULD BE GIVEN ROUNDED VALUES SLIGHTLY LARGER THAN THE EXTREME
C     DATA VALUES.  HOWEVER, THE GRAPH MAY BE CENTERED IN VARIOUS WAYS
C     BY ASSIGNING THE Y-MARGINS CONSIDERABLY LARGER VALUES. LIKEWISE
C     EXCERPTS FROM THE GRAPH MAY BE CHOSEN BY GIVING THE Y-MARGINS
C     LESSER VALUES THAN THE EXTREMES.	THE X-RANGE CANNOT BE ALTERED,
C     THE MARGINS BEING FIXED AT 1 AND N.  HOWEVER, THE SUBROUTINE CAN
C     BE CALLED USING A SUBARRAY OF WY AS ITS ARGUMENT, OR WY COULD BE
C     EMBEDDED IN A LARGER ARRAY USING AN EQUIVALENCE. ON THE OTHER HAND
C     PLTRG OR PLTRI SHOULD PROBABLY BE USED WHEN IT IS NOT SATISFACTORY
C     TO GRAPH THE ARRAY AS IT STANDS.
C     Y1     Y LOWER LIMIT
C     WY(N)  ARRAY OF Y VALUES
C     Y2     Y UPPER LIMIT
C     N      NUMBER OF POINTS
C     [16-JAN-75]
 
      DIMENSION   WY(1)
 
      IF (N.LT.2) RETURN
      CALL PLTRI (1.0,Y1,1)
      CALL PLTRI (FLOAT(N),Y2,2)
      CALL PLTRI (1.0,WY(1),3)
      DO 10 I=2,N
   10 CALL PLTRI (FLOAT(I),WY(I),4)
      RETURN
 
      END
      SUBROUTINE  PLTFM (X,Y,R,PL)
 
C     [FIDUCIAL MARK]
C     (X,Y)	CENTER OF MARK
C     R 	RADIUS OF MARK
C     PL	PEN MOVEMENT SUBROUTINE
C     [25-APR-74]
 
      CALL PL (X  ,Y  ,.FALSE.)
      CALL PL (X-R,Y  ,.TRUE.)
      CALL PL (X+R,Y  ,.TRUE.)
      CALL PL (X  ,Y  ,.TRUE.)
      CALL PL (X  ,Y-R,.TRUE.)
      CALL PL (X  ,Y+R,.TRUE.)
      CALL PL (X  ,Y  ,.TRUE.)
      RETURN

      END
      SUBROUTINE  PLTFR
 
C     [FRAME]
C     SET UP AN 8-1/2" X 11" FRAME AND LOCATE THE ORIGIN AT THE
C     CENTER OF THE PAGE.
C     [15-NOV-73]
 
      DIMENSION   IH(10),IJOB(3),IDATE(2)
      EQUIVALENCE (IJOB(1),IH(3))
      EQUIVALENCE (IDATE(1),IH(7))
      EQUIVALENCE (ITIME,IH(10))
 
      IH(1)='INEN:'
      IH(2)='	  '
      IH(6)='	  '
      IH(9)='	  '
      CALL PLOT   (0.0, 0.0, 3)
      CALL PLOT   (0.0,11.0, 2)
      CALL PLOT   (8.5,11.0, 1)
      CALL PLOT   (8.5, 0.0, 1)
      CALL PLOT   (0.0, 0.0, 1)
      CALL SYSJO  (IJOB)
      CALL DATE   (IDATE)
      CALL TIME   (ITIME)
      CALL SYMBOL (0.1,4.5,0.08,IH,-90.0,50)
      CALL PLOT   (4.25,5.50,-3)
      RETURN
 
      END
      SUBROUTINE  PLTHP (NR,NT)
 
C     [HYPERBOLIC POLAR]
C     NR   NUMBER OF CIRCLES OF CONSTANT RADIUS
C     NT   NUMBER OF RADII OF CONSTANT ANGLE
C     POLAR COORDINATE GRID WITH RADIAL HYPERBOLIC TANGENT DISTORTION.
C     [14-MAR-73]
 
      DATA	  HX,HY/4.50,3.25/
      DATA	  RR,UU/3.25,3.00/
 
      CALL PLTBO
      CALL PLTRI (-HX,-HY,1)
      CALL PLTRI ( HX, HY,2)
      DT=6.28318/FLOAT(NT)
 
      IF (NR.LT.1) GO TO 11
      DO 10 I=1,NR
   10 CALL PLTCI (0.0,0.0,RR*TANH(FLOAT(I)/UU))
   11 CALL PLTCI (0.0,0.0,RR)
 
      T=0.0
      SS=(0.25*RR)/UU
      DO 20 I=1,NT,2
      CALL PLTRI (RR*COS(T),RR*SIN(T),3)
      CALL PLTRI (SS*COS(T),SS*SIN(T),4)
      T=T+DT
      CALL PLTRI (SS*COS(T),SS*SIN(T),3)
      CALL PLTRI (RR*COS(T),RR*SIN(T),4)
   20 T=T+DT
      RETURN
 
      END
      SUBROUTINE  PLTIL (X1,Y1,Z1,X2,Y2,Z2,PL)
 
C     [INTERRUPTED LINE]
C     DRAW A LINE FROM (X1,Y1) TO (X2,Y2) SHOWING ONLY THAT PORTION
C     WHERE Z IS POSITIVE, THIS REGION BEING DETERMINED BY LINEAR
C     INTERPOLATION FROM Z1 AND Z2.
C     PL IS THE PEN MOVEMENT SUBROUTINE, PERHAPS PLTCA.
C     [05-JAN-75]
 
      LOGICAL	  P1,P2
 
      ZP(W1,W2)=W1-Z1*((W2-W1)/(Z2-Z1))
 
      P1=Z1.GE.0.0
      P2=Z2.GE.0.0
      IF (P1.EQ.P2) GO TO 10
      CALL PL (ZP(X1,X2),ZP(Y1,Y2),P1)
   10 CALL PL (X2,Y2,P2)
      RETURN
 
      END
      SUBROUTINE  PLTIR
 
C     [INCH RETICULE]
C     IT IS ASSUMED THAT THE PEN HAS ALREADY BEEN PLACED AT PAGE CENTER.
C     [25-APR-74]
 
      EXTERNAL	  PLTMC
      X= 4.0
      Y= 3.0
      R= 0.125
      DX=-1.0
      DY=-1.0
      DO 20 I=1,7
      DO 10 J=1,9
      CALL PLTFM (X,Y,R,PLTMC)
   10 X=X+DX
      DX=-DX
      X=X+DX
   20 Y=Y+DY
 
      RETURN
      END
      SUBROUTINE  PLTKB (Z1,ZE,Z2,NZ,NX,NY,PL)
 
C     [CONTOUR BORDER]
C     ZE(NX,NY)  ARRAY FROM WHICH BORDER VALUES ARE TAKEN
C     (Z1,Z2)	 INTERVAL DEFINING SCALE
C     NZ	 NUMBER OF Z-INTERVALS TO BE MARKED
C     PL	 PEN MOVEMENT SUBROUTINE, PERHAPS PLTCA
C     [05-JAN-75]
 
      EXTERNAL	  PL
      DIMENSION   ZE(1)
      DATA	  EP/0.01/
 
      IX(I,J)=I+NX*(J-1)
 
      DX=1.0/FLOAT(NX-1)
      DY=1.0/FLOAT(NY-1)
      DZ=(Z2-Z1)/FLOAT(NZ-1)
      Z=Z1
      DO 50 K=1,NZ
      FK=FLOAT(K)
      X=-HX+DX
      Y=-HY+FK*EP
      CALL PL (X-DX,Y,.FALSE.)
      DO 10 I=2,NX
      I1=IX(I-1,1)
      I2=IX(I,1)
      CALL PLTIL (X-DX,Y,ZE(I1)-Z,X,Y,ZE(I2)-Z,PL)
   10 X=X+DX
      X=HX-FK*EP
      Y=-HY+DY
      CALL PL (X,Y-DY,.FALSE.)
      DO 20 I=2,NY
      I1=IX(NX,I-1)
      I2=IX(NX,I)
      CALL PLTIL (X,Y-DY,ZE(I1)-Z,X,Y,ZE(I2)-Z,PL)
   20 Y=Y+DY
      X=HX
      Y=HY-FK*EP
      CALL PL (X,Y,.FALSE.)
      DO 30 I=2,NX
      I1=IX(NX-I+2,NY)
      I2=IX(NX-I+1,NY)
      CALL PLTIL (X,Y,ZE(I1)-Z,X-DX,Y,ZE(I2)-Z,PL)
   30 X=X-DX
      X=-HX+FK*EP
      Y=HY
      CALL PL (X,Y,.FALSE.)
      DO 40 I=2,NY
      I1=IX(1,NY-I+2)
      I2=IX(1,NY-I+1)
      CALL PLTIL (X,Y,ZE(I1)-Z,X,Y-DY,ZE(I2)-Z,PL)
   40 Y=Y-DY
   50 Z=Z+DZ
      RETURN
 
      END
      SUBROUTINE  PLTKC (Z1,ZE,Z2,NZ,KX,NX,KY,NY,PL)
 
C     [CONTOUR A COMPLEX FUNCTION]
C     ZE(NX,NY)  ARRAY OF FUNCTION VALUES
C     Z1,Z2	 RANGE OF ABSOLUTE VALUES
C     NZ	 NUMBER OF CONTOURS OF ABSOLUTE VALUE
C     (THE ARGUMENT IS CONTOURED AT 30 DEGREE INTERVALS
C		 FROM -180 DEGREES TO 180 DEGREES)
C     KX,KY	 NUMBER OF BLOCKS IN X AND Y DIRECTIONS
C     PL	 PEN MOVEMENT SUBROUTINE
C     [05-MAR-74]
 
      EXTERNAL	  CABS,CARG,PL
      COMPLEX	  ZE(1)
      DATA	  PI,NA/3.14159,13/
 
      I0=MAX0(NX/(2*KX),1)
      J0=MAX0(NY/(2*KY),1)
      IX=2*I0
      IY=2*J0
      DZ=(Z2-Z1)/FLOAT(NZ-1)
      DA=PI/6.0
 
      II=I0+1
      JJ=J0+1
   10 JA=MAX0(JJ-J0,1)
      JB=MIN0(JJ+J0,NY)
   20 IA=MAX0(II-I0,1)
      IB=MIN0(II+I0,NX)
      Z=Z1
      DO 30 K=1,NZ
      CALL KONSK (Z,IA,IB,JA,JB,ZE,NX,NY,CABS,PL)
   30 Z=Z+DZ
      A=-PI
      DO 40 K=1,NA
      CALL KONSK (A,IA,IB,JA,JB,ZE,NX,NY,CARG,PL)
   40 A=A+DA
      II=II+IX
      IF ((II-I0.LT.NX).AND.(II+I0.GT.1)) GO TO 20
      IX=-IX
      II=II+IX
      JJ=JJ+IY
      IF (JJ-J0.LT.NY) GO TO 10
      RETURN
 
      END
      SUBROUTINE  PLTKO (Z1,ZE,Z2,NZ,KX,NX,KY,NY)
 
C     [CONTOUR PLOT FROM FUNCTION ARRAY]
C     ZE(NX,NY)  ARRAY TO BE CONTOURED
C     (Z1,Z2)	 CONTOURING INTERVAL
C     NZ	 NUMBER OF CONTOUR LEVELS
C     (KX,KY)	 NUMBER OF X- AND Y- SUBDIVISIONS
C     [03-JAN-75]
 
      EXTERNAL	  PLTCA
      DIMENSION   ZE(1)
 
      CALL PLTBO
      CALL PLTIR
      CALL PLTKB (Z1,ZE,Z2,NZ,NX,NY)
      CALL PLTKP (Z1,ZE,Z2,NZ,KX,NX,KY,NY,PLTCA)
      CALL PLTEJ
      RETURN
 
      END
      SUBROUTINE  PLTKP (Z1,ZE,Z2,NZ,KX,NX,KY,NY,PL)
 
C     [CONTOUR PLOT]
C     FUNCTION VALUES TAKEN FROM AN ARRAY
C     ESSENTIALLY PLTKO, BUT WITHOUT FRAMING AND EJECT
C     ZE(NX,NY)  ARRAY TO BE CONTOURED
C     (KX,KY)	 NUMBER OF X- AND Y-SUBDIVISIONS
C     (Z1,Z2)	 CONTOURING INTERVAL
C     NZ	 NUMBER OF CONTOUR LEVELS
C     PL	 COORDINATE CONVERSION SUBROUTINE
C     [16-FEB-74]
 
      EXTERNAL	  PL
      DIMENSION   ZE(1)
 
      I0=MAX0(NX/(2*KX),1)
      J0=MAX0(NY/(2*KY),1)
      IX=2*I0
      IY=2*J0
      DZ=(Z2-Z1)/FLOAT(NZ-1)
 
      II=I0+1
      JJ=J0+1
   10 JA=MAX0(JJ-J0,1)
      JB=MIN0(JJ+J0,NY)
   20 IA=MAX0(II-I0,1)
      IB=MIN0(II+I0,NX)
      Z=Z1
      DO 30 K=1,NZ
      CALL KONSC (Z,0.0,0.0,IA,IB,JA,JB,ZE,NX,NY,PL)
   30 Z=Z+DZ
      II=II+IX
      IF ((II-I0.LT.NX).AND.(II+I0.GT.1)) GO TO 20
      IX=-IX
      II=II+IX
      JJ=JJ+IY
      IF (JJ-J0.LT.NY) GO TO 10
      RETURN
 
      END
      SUBROUTINE  PLTKX (Z1,ZE,Z2,NX,NY,PL)
 
C     [X-CROSSHATCHING]
C     ZE(NX,NY)   DATA ARRAY
C     (Z1,Z2)	  INTERVAL TO BE MARKED
C     PL	  PEN MOVEMENT SUBROUTINE, PERHAPS PLTCA
C     [05-JAN-75]
 
      EXTERNAL	  PL
      DIMENSION   ZE(1)
 
      IX(I,J)=I+NX*(J-1)
      VI(Z)=(Z-Z1)*(Z2-Z)
 
      X=0.0
      Y=0.0
      DX=1.0/FLOAT(NX-1)
      DY=1.0/FLOAT(NY-1)
      DO 30 J=1,NY,2
      CALL PL (X,Y,.FALSE.)
      DO 10 I=2,NX
      I1=IX(I-1,J)
      I2=IX(I,J)
      CALL PLTIL (X,Y,VI(ZE(I1)),X+DX,Y,VI(ZE(I2)),PL)
   10 X=X+DX
      DX=-DX
      Y=Y+DY
      IF (J.GE.NY) RETURN
      CALL PL (X,Y,.FALSE.)
      DO 20 I=2,NX
      I1=IX(NX-I+2,J+1)
      I2=IX(NX-I+1,J+1)
      CALL PLTIL (X,Y,VI(ZE(I1)),X+DX,Y,VI(ZE(I2)),PL)
   20 X=X+DX
      DX=-DX
   30 Y=Y+DY
      RETURN
 
      END
      SUBROUTINE  PLTKY (Z1,ZE,Z2,NX,NY,PL)
 
C     [Y CROSSHATCHING]
C     ZE(NX,NY)   DATA ARRAY
C     (Z1,Z2)	  INTERVAL TO BE MARKED
C     PL	  PEN MOVEMENT SUBROUTINE, PERHAPS PLTCA
C     [05-JAN-75]
 
      EXTERNAL	  PL
      DIMENSION   ZE(1)
 
      IX(I,J)=I+NX*(J-1)
      VI(Z)=(Z-Z1)*(Z2-Z)
 
      X=1.0
      Y=1.0
      DX=-1.0/FLOAT(NX-1)
      DY=-1.0/FLOAT(NY-1)
      DO 30 I=1,NX,2
      CALL PL (X,Y,.FALSE.)
      DO 10 J=2,NY
      I1=IX(NX-I+1,NY-J+2)
      I2=IX(NX-I+1,NY-J+1)
      CALL PLTIL (X,Y,VI(ZE(I1)),X,Y+DY,VI(ZE(I2)),PL)
   10 Y=Y+DY
      DY=-DY
      X=X+DX
      IF (I.GE.NX) RETURN
      CALL PL (X,Y,.FALSE.)
      DO 20 J=2,NY
      I1=IX(NX-I,J-1)
      I2=IX(NX-I,J)
      CALL PLTIL (X,Y,VI(ZE(I1)),X,Y+DY,VI(ZE(I2)),PL)
   20 Y=Y+DY
      DY=-DY
   30 X=X+DX
      RETURN
 
      END
      SUBROUTINE  PLTLA (I)
 
C     [LABEL]
C     SOMETIMES IT IS CONVENIENT TO IDENTIFY PLOTTER SHEETS WITH
C     A SHORT LABEL IN THE LOWER MARGIN (MAXIMUM OF 5 LETTERS).
C     PLOTTER ORIGIN IS TAKEN AS PAGE CENTER, SO PLTLA MUST BE
C     CALLED AFTER CALLING PLTFR OR PLTBO, AND BEFORE CALLING
C     PLTEJ.
C     [21-MAR-74]
 
      CALL SYMBOL (-3.85,4.30,0.2,I,-90.0,5)
      RETURN
 
      END
      SUBROUTINE  PLTLH (X,Y,P)
 
C     [LEFT HALF]
C     SCALE THE CARTESIAN COORDINATES (X,Y) SO AS TO POSITION A GRAPH
C     IN THE LEFT HALF OF A PLOTTER PAGE.
C     [20-APR-74]
 
      LOGICAL	  P
      DATA	  HX,HY/4.50,3.25/
 
      CALL PLTMS (HX*(X-1.0),HY*(Y-0.5),P)
      RETURN
 
      END
      SUBROUTINE  PLTMA (X,Y,X0,Y0)
 
C     [MARGIN ADJUSTER]
C     PLTMA ADJUSTS THE COORDINATES (X0,Y0) TO PROVIDE A HYPERBOLIC
C     TANGENTIAL DISTORTION (X,Y) OF THE OUTER ONE INCH MARGIN OF A
C     GRAPH, SO THE PEN WILL NEVER RUN OVER THE EDGE OF THE PAGE.
C     [02-NOV-73]
 
      DATA	 HX,HY/4.50,3.25/
 
      X=X0
      IF (X.GT. HX) X= HX+TANH(X-HX)

      IF (X.LT.-HX) X=-HX+TANH(X+HX)
      Y=Y0
      IF (Y.GT. HY) Y= HY+TANH(Y-HY)
      IF (Y.LT.-HY) Y=-HY+TANH(Y+HY)
      RETURN
 
      END
      SUBROUTINE  PLTMC (X,Y,S)
 
C     [MARGIN COMPRESSED]
C     PLTMC (X,Y,S) MOVES THE PEN TO THE POINT (X,Y) ON AN 8-1/2" X 11"
C     SHEET.  HOWEVER, A BORDER IS MAINTAINED 1" FROM THE EDGE, BEYOND
C     WHICH NO PEN MOVEMENT IS PERMITTED.  LINEAR INTERPOLATION IS USED
C     TO OBTAIN AN ACCURATE INTERSECTION OF ANY LINE SEGMENT WHITH THIS
C     MARGIN. IF S=.TRUE. THE PEN IS LOWERED, OTHERWISE THE NEW PEN
C     POSITION IS ONLY NOTED AS (X0,Y0). THE PEN IS ALLOWED TO WRITE IN
C     THIS MARGIN, WHICH IS SUBJECT TO A HYPERBOLIC TANGENT DISTORTION.
C     [02-NOV-73]
 
      LOGICAL	 S,S0,P1,P2,Q
 
      IF (.NOT.S) GO TO 30
      X1=X0
      Y1=Y0
      X2=X
      Y2=Y
      CALL PLTMT (X1,Y1,P1,X2,Y2,P2,Q)
      IF (S0) GO TO 10
      CALL PLTMA (EX,WY,X0,Y0)
      CALL PLOT  (WY,-EX,3)
   10 IF (Q) GO TO 20
      CALL PLTME (X0,Y0,X,Y)
      GO TO 30
   20 IF (P1) CALL PLTME (X0,Y0,X1,Y1)
      CALL PLOT  (Y2,-X2,2)
      IF (P2) CALL PLTME (X2,Y2,X,Y)
   30 S0=S
      X0=X
      Y0=Y
      RETURN
 
      END
      SUBROUTINE  PLTME (X1,Y1,X2,Y2)
 
C     [MARGINAL EXCURSION]
C     N POINTS OF A DISTORTED LINE IN THE MARGIN, AS REQUIRED BY PLTMC.
C     [03-NOV-73]
 
      N=11
      EN=FLOAT(N-1)
      DX=(X2-X1)/EN
      DY=(Y2-Y1)/EN
      EX=X1
      WY=Y1
      DO 10 I=1,N
      CALL PLTMA (XX,YY,EX,WY)
      CALL PLOT  (YY,-XX,2)
      EX=EX+DX
   10 WY=WY+DY
      RETURN
 
      END
      SUBROUTINE  PLTMS (X,Y,S)
 
C     [MARGIN SUPPRESSED]
C     PLTMS (X,Y,S) MOVES THE PEN TO THE POINT (X,Y) ON AN 8-1/2" X 11"
C     SHEET.  HOWEVER, A BORDER IS MAINTAINED 1" FROM THE EDGE, BEYOND
C     WHICH NO PEN MOVEMENT IS PERMITTED. LINEAR INTERPOLATION IS USED
C     TO OBTAIN AN ACCURATE INTERSECTION OF ANY LINE SEGMENT WHITH THIS
C     MARGIN. IF S=.TRUE. THE PEN IS MOVED IN THE LOWERED POSITION,
C     OTHERWISE THE NEW PEN POSITION IS MERELY NOTED AS (X0,Y0). NULL
C     INTERVALS ARE INVISIBLE, SO THE PEN IS NOT RAISED OR LOWERED EVEN
C     THOUGH THIS MOVEMENT HAS BEEN CALLED FOR.
C     [15-MAY-74]
 
      LOGICAL	 S,S0,P0,P1,Q,EQ
 
      EQ(X,Y)=ABS(X-Y).LE.(1.0E-5)
 
      IF (S) GO TO 20
   10 S0=.TRUE.
   11 X0=X
      Y0=Y
      RETURN
   20 IF (EQ(X,X0).AND.EQ(Y,Y0)) GO TO 10
      X1=X
      Y1=Y
      CALL PLTMT (X0,Y0,P0,X1,Y1,P1,Q)
      IF (.NOT.Q) GO TO 10
      IF (S0.OR.P0) CALL PLOT (Y0,-X0,3)
      CALL PLOT  (Y1,-X1,2)
      S0=P1
      GO TO 11
 
      END
      SUBROUTINE  PLTMT (X1,Y1,P1,X2,Y2,P2,Q)
 
C     [MARGIN TRIMMER]
C     THE POINTS Z1=(X1,Y1) AND Z2=(X2,Y2) ARE ADJUSTED TO ENCOMPASS
C     ONLY THAT PART OF THE LINE JOINING THEM WHICH LIES WITHIN A
C     PLOTTER PAGE WHORE HALFWIDTHS ARE HX AND HY.
C     P1=.TRUE.  Z1 WAS MOVED
C     P2=.TRUE.  Z2 WAS MOVED
C     Q=.TRUE.	 SOME PORTION OF THE LINE IS VISIBLE
C     [02-NOV-73]
 
      LOGICAL	  P1,P2,Q
      DATA	  HX,HY/4.50,3.25/
 
      EX(Y)=X1+(Y-Y1)*((X2-X1)/(Y2-Y1))
      WY(X)=Y1+(X-X1)*((Y2-Y1)/(X2-X1))
 
      P1=.FALSE.
      IF (ABS(X1).LE.HX) GO TO 5
      IF (X1.EQ.X2) GO TO 20
      X0=SIGN(HX,X1)
      Y1=WY(X0)
      X1=X0
      P1=.TRUE.
    5 IF (ABS(Y1).LE.HY) GO TO 10
      IF (Y1.EQ.Y2) GO TO 20
      Y0=SIGN(HY,Y1)
      X1=EX(Y0)
      Y1=Y0
      P1=.TRUE.
   10 P2=.FALSE.
      IF (ABS(X2).LE.HX) GO TO 15
      IF (X1.EQ.X2) GO TO 20
      X0=SIGN(HX,X2)
      Y2=WY(X0)
      X2=X0
      P2=.TRUE.
   15 IF (ABS(Y2).LE.HY) GO TO 20
      IF (Y1.EQ.Y2) GO TO 20
      Y0=SIGN(HY,Y2)
      X2=EX(Y0)
      Y2=Y0
      P2=.TRUE.
   20 X0=0.5*(X1+X2)
      Y0=0.5*(Y1+Y2)
      Q=(ABS(X0).LT.HX).AND.(ABS(Y0).LT.HY)
      RETURN
 
      END
      SUBROUTINE  PLTOR (Z1,ZE,Z2,NZ,KX,NX,KY,NY,PL)
 
C     [ORTHOGONAL RELIEF]
C     ZE(NX,NY)  ARRAY TO BE CONTOURED
C     (KX,KY)	 NUMBER OF X- AND Y-SUBDIVISIONS
C     (Z1,Z2)	 CONTOURING INTERVAL
C     NZ	 NUMBER OF CONTOUR LEVELS
C     PL	 COORDINATE CONVERSION SUBROUTINE
C     [20-NOV-74]
 
      EXTERNAL	  PL
      DIMENSION   ZE(1)
 
      I0=MAX0(NX/(2*KX),1)
      J0=MAX0(NY/(2*KY),1)
      IX=2*I0
      IY=2*J0
      ZZ=Z2-Z1
      XF=ZZ/FLOAT(NX-1)
      YF=ZZ/FLOAT(NY-1)
      DZ=(3.0*ZZ)/FLOAT(NZ-1)
 
      II=I0+1
      JJ=J0+1
   10 JA=MAX0(JJ-J0,1)
      JB=MIN0(JJ+J0,NY)
   20 IA=MAX0(II-I0,1)
      IB=MIN0(II+I0,NX)
      Z=Z1-ZZ
      DO 30 K=1,NZ
      CALL KONSC (Z,-XF,YF,IA,IB,JA,JB,ZE,NX,NY,PL)
   30 Z=Z+DZ
      II=II+IX
      IF ((II-I0.LT.NX).AND.(II+I0.GT.1)) GO TO 20
      IX=-IX
      II=II+IX
      JJ=JJ+IY
      IF (JJ-J0.LT.NY) GO TO 10
      RETURN
 
      END
      SUBROUTINE  PLTPO (T,R,P)
 
C     [POLAR]
C     CHANGE (T,R) TO (X,Y) PERMITTING PEN MOVEMENTS TO BE DEFINED
C     DIRECTLY IN POLAR COORDINATES.
C     0.0.LE.T.LE.1.0 FILLS THE PAGE
C     0.0.LE.R.LE.1.0 FILLS THE PAGE
C     [19-MAY-74]
 
      LOGICAL	  P
      DATA	  HX,HY/4.50,3.25/
 
      S=HY*R
      U=6.28318*T
      CALL PLTMS (S*COS(U),S*SIN(U),P)
      RETURN
 
      END
      SUBROUTINE  PLTPV (Z1,ZE,Z2,NR,NP)
 
C     [POLAR VIEW]
C     PROGRAM TO PRODUCE A PERSPECTIVE DRAWING OF A SINGLE VALUED
C     FUNCTION DEFINED IN POLAR COORDINATES, IN SUCH A WAY AS TO
C     EXHIBIT THE RADIAL AND CIRCUMFERENTIAL ARCS.
C     (Z1,Z2)	RANGE OF FUNCTION VALUES
C     ZE(NR,NP) ARRAY OF FUNCTION VALUES
C     NR	NUMBER OF POINTS ON ONE RADIUS
C     NP	NUMBER (=4*N+1) OF ANGLES
C     [10-NOV-74]
 
      DIMENSION   ZE(1)
      DATA	  Q0,Q1,Q2,Q3,Q4/0.000,1.571,3.142,4.713,6.283/
      DATA	  S1,S3/1.570,4.712/
 
      NQ=1+(NP-1)/4
      NN=NR*(NQ-1)
      R1=0.02
      R2=1.00
      CALL PLTFR
      CALL VISNH
	CALL VISPS(Z1,ZE(3*NN+1),Z2,R1,R2,NR,Q3,Q4,NQ,-1, 1)
      CALL VISPS (Z1,ZE        ,Z2,R1,R2,NR,Q0,S1,NQ, 1, 1)
      CALL VISNH
      CALL VISPS (Z1,ZE(2*NN+1),Z2,R1,R2,NR,Q2,S3,NQ,-1,-1)
      CALL VISPS (Z1,ZE(  NN+1),Z2,R1,R2,NR,Q1,Q2,NQ, 1,-1)
      CALL PLTEJ
      RETURN
 
      END
      SUBROUTINE  PLTQ1 (X,Y,P)
 
C     [QUADRANT #1]
C     SCALE X,Y TO FIT INTO THE FIRST QUADRANT ASSUMING THAT THEY
C     ALREADY LIE IN THE RANGE (0.0 .LE. X,Y .LE. 1.0).
C     [15-FEB-74]
 
      LOGICAL	  P
      DATA	  HX,HY/4.50,3.25/
 
      CALL PLTMS (HX*X,HY*Y,P)
      RETURN
 
      END
      SUBROUTINE  PLTQ2 (X,Y,P)
 
C     [QUADRANT #2]
C     SCALE X,Y TO FIT INTO THE SECOND QUADRANT ASSUMING THAT THEY
C     ALREADY LIE IN THE RANGE (0.0 .LE. X,Y .LE. 1.0).
C     [15-FEB-74]
 
      LOGICAL	  P
      DATA	  HX,HY/4.50,3.25/
 
      CALL PLTMS (HX*(X-1.0),HY*Y,P)
      RETURN
 
      END
      SUBROUTINE  PLTQ3 (X,Y,P)
 
C     [QUADRANT #3]
C     SCALE X,Y TO FIT INTO THE THIRD QUADRANT ASSUMING THAT THEY
C     ALREADY LIE IN THE RANGE (0.0 .LE. X,Y .LE. 1.0).
C     [15-FEB-74]
 
      LOGICAL	  P
      DATA	  HX,HY/4.50,3.25/
 
      CALL PLTMS (HX*(X-1.0),HY*(Y-1.0),P)
      RETURN
 
      END
      SUBROUTINE  PLTQ4 (X,Y,P)
 
C     [QUADRANT #4]
C     SCALE X,Y TO FIT INTO THE FOURTH QUADRANT ASSUMING THAT THEY
C     ALREADY LIE IN THE RANGE (0.0 .LE. X,Y .LE. 1.0).
C     [15-FEB-74]

      LOGICAL	  P
      DATA	  HX,HY/4.50,3.25/
 
      CALL PLTMS (HX*X,HY*(Y-1.0),P)
      RETURN
 
      END
      SUBROUTINE  PLTRG (X1,X,X2,Y1,Y,Y2,N)
 
C     [RECTANGULAR GRAPH]
C     PLOT A GRAPH IN RECTANGULAR COORDINATES, BY CONNECTING SUCCESSIVE
C     DATA POINTS BY STRAIGHT LINES. THE POINTS DEFINING THE GRAPH ARE
C     TAKEN FROM TWO ARRAYS, ONE CONTAINING THE X-VALUES AND THE OTHER
C     CONTAINING THE Y-VALUES.	THE RESPECTIVE SCALES ARE INDICATED BY
C     THE VALUES TO BE ASSIGNED TO THE MARGINS OF THE GRAPH. ORDINARILY
C     THE MARGINS WOULD BE GIVEN ROUNDED VALUES SLIGHTLY LARGER THAN
C     THE EXTREME DATA VALUES.	HOWEVER, THE GRAPH MAY BE CENTERED IN
C     VARIOUS WAYS BY ASSIGNING ONE OR MORE MARGINS CONSIDERABLY LARGER
C     VALUES.  LIKEWISE EXCERPTS FROM THE GRAPH MAY BE CHOSEN BY GIVING
C     THE MARGINS LESSER VALUES THAN THE EXTREMES.
C     X1     X LOWER LIMIT
C     X(N)   ARRAY OF X VALUES
C     X2     X UPPER LIMIT
C     Y1     Y LOWER LIMIT
C     Y(N)   ARRAY OF Y VALUES
C     Y2     Y UPPER LIMIT
C     N      NUMBER OF POINTS
C     [12-OCT-74]
 
      DIMENSION   X(1),Y(1)
      DATA	  HX,HY/4.50,3.25/
 
      EX(X)=(X-X1)*SCX-HX
      WY(Y)=(Y-Y1)*SCY-HY
 
      IF (N.LT.2) RETURN
      SCX=(2.0*HX)/(X2-X1)
      SCY=(2.0*HY)/(Y2-Y1)
      CALL PLTMS (EX(X(1)),WY(Y(1)),.FALSE.)
      DO 10 I=2,N
   10 CALL PLTMS (EX(X(I)),WY(Y(I)),.TRUE.)
      RETURN
 
      END
      SUBROUTINE  PLTRH (X,Y,P)
 
C     [RIGHT HALF]
C     SCALE THE CARTESIAN COORDINATES (X,Y) SO AS TO POSITION A
C     GRAPH IN THE RIGHT HALF OF A PLOTTER PAGE.
C     [20-APR-74]
 
      LOGICAL	  P
      DATA	  HX,HY/4.50,3.25/
 
      CALL PLTMS (HX*X,HY*(Y-0.5),P)
      RETURN
 
      END
      SUBROUTINE  PLTRI (X,Y,L)
 
C     [RECTANGULAR INCREMENTAL GRAPH]
C     PLOT A RECTANGULAR GRAPH POINT BY POINT. PEN ORIGIN IS ASSUMED
C     TO BE THE CENTER OF A PAGE WHOSE DIMENSIONS ARE 2*HX X 2*HY.THE
C     OPTIONS AFFORDED BY L ARE:
C	L=1   (X1,Y1) RESPECTIVE LOWER LIMITS
C	L=2   (X2,Y2) RESPECTIVE UPPER LIMITS
C	L=3   FIRST POINT OF A SERIES
C	L=4   SUBSEQUENT POINTS
C	L=5   RECTANGULAR AXES THROUGH (X,Y)
C	L=6   TICK MARK AT (X,Y)
C	L=7   LARGER TICK MARK
C     AS LONG AS OPTIONS 6 AND 7 CALL PLTMC AND THE OTHERS CALL PLTMS,
C     TICK MARKS MAY BE PLACED IN ANY ORDER-BUT NOT BEFORE THE LIMITS
C     HAVE BEEN ESTABLISHED. THE LIMITS MUST BE DEFINED BEFORE STARTING
C     THE GRAPH.  THE INITIAL POINT SHOULD ONLY BE ENTERED BY OPTION
C     3, WHICH MAY ALSO BE USED TO CREATE GAPS IN THE GRAPH, OR TO
C     INIATE A NEW CURVE.  TO SUPPRESS ONE OF THE AXES DRAWN BY OPTION
C     5, CHOOSE A CROSSING POINT OUTSIDE OF THE RANGE OF THE GRAPH.
C     [10-NOV-74]
 
      EXTERNAL	   PLTMC
      DATA	   HX,HY/4.50,3.25/
 
      EX(X)=(X-X1)*SX-HX
      WY(Y)=(Y-Y1)*SY-HY
 
      GO TO (10,20,5,40,5,5,5),L
    5 SX=(2.0*HX)/(X2-X1)
      SY=(2.0*HY)/(Y2-Y1)
      GO TO (7,7,30,7,50,60,70),L
    7 RETURN
 
   10 X1=X
      Y1=Y
      RETURN
   20 X2=X
      Y2=Y
      RETURN
   30 CALL PLTMS (EX(X),WY(Y),.FALSE.)
      RETURN
   40 CALL PLTMS (EX(X),WY(Y),.TRUE.)
      RETURN
   50 CALL PLTMS (-HX,WY(Y),.FALSE.)
      CALL PLTMS ( HX,WY(Y),.TRUE.)
      CALL PLTMS (EX(X), HY,.FALSE.)
      CALL PLTMS (EX(X),-HY,.TRUE.)
      RETURN
   60 CALL PLTFM (EX(X),WY(Y),0.04,PLTMC)
      RETURN
   70 CALL PLTFM (EX(X),WY(Y),0.06,PLTMC)
      RETURN
 
      END
      SUBROUTINE  PLTRV (Z1,ZE,Z2,NX,MX,NY,TH)
 
C     [ROTATED VIEW]
C     PROGRAM TO PRODUCE A PERSPECTIVE DRAWING OF A SINGLE VALUED
C     FUNCTION DEFINED IN CARTESIAN COORDINATES, EXHIBITING ARCS
C     ON THE SURFACE PARALLEL TO THE COORDINATE AXES. FOR GREATER
C     VARIETY IN PRESENTATION, THE ENTIRE FIGURE MAY BE ROTATED
C     THROUGH AN ANGLE, WHICH SHOULD BE SPECIFIED IN DEGREES.
C     ZE     ARRAY OF FUNCTION VALUES
C     Z1,Z2  RANGE OF FUNCTION VALUES
C     NX     ACTUAL  LENGTH OF COLUMNS
C     MX     MAXIMUM LENGTH OF COLUMNS
C     NY     ACTUAL  LENGTH OF ROWS
C     TH     ANGLE OF ROTATION (-90.0.LT.TH.LT.90.0) IN DEGREES
C     [22-NOV-74]
 
      EXTERNAL	  PLTCA
      DIMENSION   ZE(1)
 
      IF ((TH.LE.-90.0).OR.(TH.GE.90.0)) RETURN
      CALL PLTFR
      CALL VISNH
      CALL VISRS (Z1,ZE,Z2,NX,MX,NY,NY,TH,PLTCA)
      CALL PLTEJ
      RETURN
 
      END
	SUBROUTINE  PLTSE (Z1,ZE,Z2,NX,MX,NY)
 
C     [SOUTHEAST VIEW]
C     PROGRAM TO PRODUCE A PERSPECTIVE DRAWING OF A SINGLE VALUED
C     FUNCTION DEFINED IN CARTESIAN COORDINATES, IN SUCH A WAY AS
C     TO EXHIBIT ARCS ON THE SURFACE PARALLEL TO THE COORDINATE
C     AXES. FOR GREATER CLARITY IN PRESENTATION, THE ENTIRE FIGURE
C     MAY BE SHEARED, HORIZONTALLY WHICH WILL GIVE THE ILLUSION OF
C     A SIDEWISE PERSPECTIVE, AND VERTICALLY TO GIVE THE ILLUSION OF
C     DEPTH AND TO EXPOSE THE REMOTER DETAILS WHICH WOULD OTHERWISE
C     BE HIDDEN.  SHEARING IS PREFERABLE TO ROTATION WHENEVER IT IS
C     DESIRED TO MAINTAIN HORIZONTAL LINES HORIZONTAL.
C     ZE     ARRAY OF FUNCTION VALUES
C     Z1,Z2  RANGE OF FUNCTION VALUES
C     NX     ACTUAL  LENGTH OF COLUMNS
C     MX     MAXIMUM LENGTH OF COLUMNS
C     NY     ACTUAL  LENGTH OF ROWS
C     [22-NOV-74]
 
      EXTERNAL	  PLTCA
      DIMENSION   ZE(1)
 
      CALL PLTFR
      CALL VISNH
      CALL VISDS (Z1,ZE,Z2,1,NX,MX,1,NY,NY,0.2,0.2,-1,1,PLTCA)
      CALL PLTEJ
      RETURN
 
      END
      SUBROUTINE  PLTSP (PH,TH,P)
 
C     [SPHERICAL POLAR]
C     CHANGE THE ANGULAR VARIABLES PH,TH TO THE CARTESIAN COORDINATES
C     X,Z SO AS TO DEFINE DIRECTLY IN SPHERICAL POLAR COORDINATES POINTS
C     WHICH LIE UPON THE SURFACE OF A CONSTANT SPHERE AND GRAPH THEIR
C     PROJECTION ON THE X-Z PLANE. PH,TH ARE BOTH SUPPOSED TO LIE IN
C     THE RANGE 0.0 .LE. PH,TH .LE. 1.0, SINCE THIS IS THE RANGE ASSUMED
C     BY SUCH SUBROUTINES AS THE CONTOURING PROGRAMS.
C     [16-NOV-74]
 
      LOGICAL	  P
      DATA	  HX,HY/4.50,3.25/
 
      R=HY
      THE=3.14159*TH
      PHI=6.2831*PH
      X=R*SIN(THE)*COS(PHI)
      Y=R*SIN(THE)*SIN(PHI)
      Z=R*COS(THE)
      CALL PLTMS (X,Z,(P.AND.(Y.GT.0.0)))
      RETURN
 
      END
      SUBROUTINE  PLTSS (Z1,ZE,Z2,NX,MX,NY)
 
C     [SOUTHERN STEREOPAIR]
C     PROGRAM TO PRODUCE TWO PERSPECTIVE DRAWINGS OF A FUNCTION
C     DEFINED IN A RECTANGULAR ARRAY.  THE DRAWINGS ARE OFFSET
C     IN A MANNER SUITABLE FOR THEIR USE AS A STEREOPAIR.
C     ZE     ARRAY OF FUNCTION VALUES
C     Z1,Z2  RANGE OF FUNCTION VALUES
C     NX     ACTUAL  LENGTH OF COLUMNS
C     MX     MAXIMUM LENGTH OF COLUMNS
C     NY     ACTUAL  LENGTH OF ROWS
C     [06-OCT-74]
 
      EXTERNAL	  PLTLH,PLTRH
      DIMENSION   ZE(1)
 
      CALL PLTFR
      CALL VISNH
      CALL VISDS (Z1,ZE,Z2,1,NX,MX,1,NY,NY,0.1,0.25,-1,1,PLTLH)
      CALL VISNH
      CALL VISDS (Z1,ZE,Z2,1,NX,MX,1,NY,NY,0.1,0.25, 1,1,PLTRH)
      CALL PLTEJ
      RETURN
      END
      SUBROUTINE  PLTSV (FU,NP,NT,O,PR,PL)
 
C     [SPHERICAL VIEW]
C     PROGRAM TO PRODUCE A PERSPECTIVE DRAWING OF A SINGLE VALUED
C     FUNCTION DEFINED OVER A SPHERICAL SURFACE, SO AS TO EXHIBIT
C     THE ARCS OF LATITUDE AND LONGITUDE.
C     FU(NP,NT) ARRAY OF FUNCTION VALUES
C     NP	NUMBER (=2*N) OF POINTS ON ONE LATITUDE
C     NT	NUMBER OF POINTS ON ONE LONGITUDE
C     O(3,3)	ORTHOGONAL ROTATION MATRIX
C     PR	PROJECTION SUBROUTINE
C     PL	PEN MOVEMENT SUBROUTINE
C     [22-NOV-74]
 
      EXTERNAL	  PR,PL
      LOGICAL	  B,C
      DIMENSION   FU(1),O(3,3)
 
      NH=NP/2
      CALL VISNP (PH,TH,JP,IT,NP,NT,O)
      IF (TH.GT.(1.57079)) GO TO 10
      I1=1
      I2=IT
      I3=IT
      I4=NT
      S1= 1.0
      S2=-1.0
      GO TO 12
   10 I1=IT
      I2=NT
      I3=1
      I4=IT
      S1=-1.0
      S2= 1.0
   12 J1=JP
      J2=JP+NH
      J3=JP-NH
      J4=JP
      CALL PR (R,P,0.1,TH,PH+0.05,O)
      B=((-0.25).LT.P).AND.(P.LE.(0.25))
      C=.NOT.B
      CALL VISNH
      CALL VISSS (FU,J1,J2,NP,I1,I2,NT,1,-1,S1,B,O,PR,PL)
      CALL VISSS (FU,J1,J2,NP,I1,I2,NT,1, 1,S2,B,O,PR,PL)
      CALL VISSS (FU,J1,J2,NP,I3,I4,NT,1, 1,S2,B,O,PR,PL)
      CALL VISSS (FU,J1,J2,NP,I3,I4,NT,1,-1,S1,B,O,PR,PL)
      CALL VISNH
      CALL VISSS (FU,J3,J4,NP,I1,I2,NT,-1,-1,S1,C,O,PR,PL)
      CALL VISSS (FU,J3,J4,NP,I1,I2,NT,-1, 1,S2,C,O,PR,PL)
      CALL VISSS (FU,J3,J4,NP,I3,I4,NT,-1, 1,S2,C,O,PR,PL)
      CALL VISSS (FU,J3,J4,NP,I3,I4,NT,-1,-1,S1,C,O,PR,PL)
      RETURN
 
      END
      SUBROUTINE  PLTSW (Z1,ZE,Z2,NX,MX,NY)
 
C     [SOUTHWEST VIEW]
C     ZE     ARRAY OF FUNCTION VALUES
C     Z1,Z2  RANGE OF FUNCTION VALUES
C     NX     ACTUAL  LENGTH OF COLUMNS
C     MX     MAXIMUM LENGTH OF COLUMNS
C     NY     ACTUAL  LENGTH OF ROWS
C     [22-NOV-74]
 
      EXTERNAL	  PLTCA
      DIMENSION   ZE(1)
 
      CALL PLTFR
      CALL VISNH
      CALL VISDS (Z1,ZE,Z2,1,NX,MX,1,NY,NY,0.2,0.2,1,1,PLTCA)
      CALL PLTEJ
      RETURN
 
      END
      SUBROUTINE  PLTTG (N)
 
C     [TRIANGULAR GRID]
C     PLTTG (N) SETS UP A TRIANGULAR GRID WITH N GRID INTERVALS.
C     [14-APR-73]
 
      DATA	  Z,U/0.0,1.0/
 
      IF (N.LE.0) RETURN
      CALL PLTTP (U,Z,Z,.FALSE.)
      I=0
   10 A=FLOAT(N-I)
      B=FLOAT(I)
      CALL PLTTP (A,B,Z,.TRUE.)
      CALL PLTTP (A,Z,B,.TRUE.)
      CALL PLTTP (Z,A,B,.TRUE.)
      CALL PLTTP (B,A,Z,.TRUE.)
      CALL PLTTP (B,Z,A,.TRUE.)
      CALL PLTTP (Z,B,A,.TRUE.)
      CALL PLTTP (A,B,Z,.TRUE.)
      I=I+1
      IF (N-I.GE.I) GO TO 10
      CALL PLTTP (U,U,U,.FALSE.)
      RETURN
 
      END
      SUBROUTINE  PLTTH (X,Y,P)
 
C     [TOP HALF]
C     SCALE THE CARTESIAN COORDINATES X,Y SO AS TO PLACE A GRAPH
C     IN THE TOP HALF OF A PLOTTER PAGE.
C     [20-APR-74]
 
      LOGICAL	  P
      DATA	  HX,HY/4.50,3.25/
 
      CALL PLTMS (HX*Y,2.0*HY*(0.5-X),P)
      RETURN
 
      END
      SUBROUTINE  PLTTP (X,Y,Z,P)
 
C     [TRIANGULAR POINT]
C     PLTTP (X,Y,Z,P) INSERTS A POINT ON A TRIANGULAR GRAPH
C     [14-APR-73]
 
      LOGICAL	  P
 
      S=X+Y+Z
      EX=(3.5*(Y-X))/S
      WY=(-2.02*(X+Y)+4.04*Z)/S
      CALL PLTMC (EX,WY,P)
 
      RETURN
      END
      SUBROUTINE  PLTTV (Z1,ZE,Z2,N,M)
 
C     [TRIANGULAR VIEW]
C     PROGRAM TO PRODUCE A PERSPECTIVE DRAWING OF A FUNCTION DEFINED
C     OTER THE UNIT TRIANGLE IN TERMS OF HOMOGENEOUS COORDINATES. THE
C     FUNCTION VALUES OCCUPY THE UPPER TRIANGULAR PORTION OF THE SQUARE
C     MATRIX ZE(M,M), OF WHICH ONLY THE FRAGMENT ZE(N,N) IS TO BE DRAWN.
C     IN GENERATING THE DRAWING, THE VALUES IN ZE WILL BE SCALED TO THE
C     RANGE Z1,Z2.
C     [22-NOV-74]
 
      DIMENSION   ZE(1)
 
      CALL PLTFR
      CALL VISNH
      CALL VISTS (Z1,ZE,Z2,N,M)
      CALL PLTEJ
      RETURN
 
      END
      SUBROUTINE  PLTUR (XA,X1,DX,X2,XB,YA,Y1,DY,Y2,YB,W,PL)
 
C     [UNIT RETICLE]
C     COVER THE PLOTTER PAGE WITH A NET OF FIDUCIAL MARKS INDICATING
C     UNIT INTERVALS OF DATA.
C     (XA,XB)  X-VALUES AT X-MARGINS
C     (X1,X2)  X-INTERVAL TO BE RETICLED
C     DX       X-DISTANCE BETWEEN CENTERS OF FIDUCIAL MARKS
C     (YA,YB)  Y-VALUES AT Y-MARGINS
C     (Y1,Y2)  Y-INTERVAL TO BE RETICLED
C     DY       Y-DISTANCE BETWEEN CENTERS OF FIDUCIAL MARKS
C     W        WIDTH OF FIDUCIAL MARK
C     PL       PEN MOVEMENT SUBROUTINE
C     DX AND DY MAY BE SIGNED OR MAY BE ABSOLUTE VALUES, LIKEWISE
C     THE X-, AND Y-INTERVALS MAY BE EITHER INCREASING OR DECREASING.
C     PLTUR ASSUMES THE UNIT SQUARE FOR ITS PAGE FORMAT, SO THAT
C     PL=PLTCA IS A SUITABLE ARGUMENT.
C     [05-JAN-75]
 
      EXTERNAL	  PL
 
      EX(X)=XS*(X-XA)
      WY(Y)=YS*(Y-YA)
 
      XS=1.0/(XB-XA)
      YS=1.0/(YB-YA)
      D=SIGN(DX,XB-XA)
      E=SIGN(DY,YB-YA)
      S=SIGN(1.0,D)
      T=SIGN(1.0,E)
      S1=S*X1
      S2=S*X2
      T1=T*Y1
      X=X2
      Y=Y2
   10 CALL PLTFM (EX(X),WY(Y),W,PL)
      X=X-D
      IF (((S*X).GE.S1).AND.((S*X).LE.S2)) GO TO 10
      D=-D
      X=X-D
      Y=Y-E
      IF ((T*Y).GE.T1) GO TO 10
      RETURN
 
      END
      SUBROUTINE  VISBO (X1,T1,B1,M,X0,T0,B0,N0,X,Y,N,I,PL)
 
C     [BOUNDS]
C     X(N)   ARRAY OF ARGUMENTS
C     Y(N)   ARRAY OF FUNCTION VALUES
C     I      DIRECTION OF PEN MOVEMENT
C     PL     PEN MOVEMENT SUBROUTINE
C     [29-MAY-74]
 
      LOGICAL	  L,PO,EQ,VV,VISSL
      DIMENSION   U(2),X(1),Y(1)
      DIMENSION   X0(1),T0(1),B0(1),X1(1),T1(1),B1(1)
      EQUIVALENCE (U1,U(1)),(U2,U(2))
      DATA	  EP/1.0E-4/
 
      II(J)=MAX0(MIN0(J+I,M),1)
      PO(X)=X.GT.EP
      EQ(X,Y)=ABS(X-Y).LE.(0.5E-4)
 
C === INITIALIZATION
 
      IF (N.LE.1) RETURN
      J=(N+1-I*(N-1))/2
      J1=((M+1)*(1-I))/2
      CALL PL (X(J),Y(J),.FALSE.)
      IF (N0.LE.1) GO TO 61
      S=FLOAT(I)
      ET= 1.0
      EB=-1.0

      L=.TRUE.
      K=(1-I)/2
      J0=(N0+1-I*(N0-1))/2
      Z=X(J)
      Z0=X0(J0)
      IF (EQ(Z,Z0)) GO TO 32
      IF (S*(Z-Z0)) 10,32,20
 
C --- IF THE FUNCTION IS DEFINED WHILE BOUNDS ARE NOT, THE FUNCTION
C --- IS VISIBLE AND MUST BE COPIED, ESTABLISHING NEW BOUNDS.
 
   10 J1=II(J1)
      CALL PL (X(J),Y(J),.TRUE.)
      X1(J1)=X(J)
      T1(J1)=Y(J)
      B1(J1)=Y(J)
      J=J+I
      IF (EQ(X(J),Z0)) GO TO 30
      IF (S*(X(J)-Z0)) 10,30,30
 
C --- IF BOUNDS, BUT NOT THE FUNCTION, ARE DEFINED, THEY PERSIST.
 
   20 J1=II(J1)
      X1(J1)=X0(J0)
      T1(J1)=T0(J0)
      B1(J1)=B0(J0)
      J0=J0+I
      IF (EQ(Z,X0(J0))) GO TO 30
      IF (S*(Z-X0(J0))) 30,30,20
 
C === MAIN LOOP
 
C --- AT A POINT WHERE EITHER THE FUNCTION OR THE BOUNDS ARE
C --- DEFINED, IT MAY BE NECESSARY TO OBTAIN THE OTHER BY LINEAR
C --- INTERPOLATION, UNLESS THEIR POINTS OF DEFINITION COINCIDE.
 
   30 IF ((J.LT.1).OR.(J.GT.N))    GO TO 50
      IF ((J0.LT.1).OR.(J0.GT.N0)) GO TO 60
      Z=X(J)
      Z0=X0(J0)
   32 EX=S*AMIN1(S*Z,S*Z0)
      WY=VISLI(EX,X,Y,MAX0(MIN0(J+K,N),2))
      TO=VISLI(EX,X0,T0,MAX0(MIN0(J0+K,N0),2))
      BO=VISLI(EX,X0,B0,MAX0(MIN0(J0+K,N0),2))
      IF (EQ(EX,Z0)) J0=J0+I
      IF (EQ(EX,Z))  J=J+I
 
C --- POSSIBLE INTERSECTIONS BETWEEN THE FUNCTION AND THE BOUNDS
C --- MUST BE RECORDED SO AS TO DESCRIBE THE NEW BOUNDS ACCURATELY.
C --- CARE IS NECESSARY TO AVOID TRIVIAL INTERSECTIONS, OR THOSE
C --- WHICH OCCUR AT ENDPOINTS.
 
      TE=AMAX1(WY,TO)
      BE=AMIN1(WY,BO)
      DT=WY-TO
      DB=WY-BO
      VT=ET+DT
      VB=EB+DB
      IF (L) GO TO 46
      JJ=0
      IF (SIGN(1.0,DT).EQ.SIGN(1.0,ET)) GO TO 41
      VT=DT-ET
      JJ=JJ+1
      U(JJ)=XX-ET*((EX-XX)/(DT-ET))
   41 IF (SIGN(1.0,DB).EQ.SIGN(1.0,EB)) GO TO 42
      VB=DB-EB
      JJ=JJ+1
      U(JJ)=XX-EB*((EX-XX)/(DB-EB))
   42 IF (JJ.EQ.0) GO TO 44
      DO 43 KK=1,JJ
      IF ((KK.EQ.1).AND.(JJ.EQ.1)) XI=U1
      IF ((KK.EQ.1).AND.(JJ.EQ.2)) XI=S*AMIN1(S*U1,S*U2)
      IF  (KK.EQ.2)		   XI=S*AMAX1(S*U1,S*U2)
      F=(XI-XX)/(EX-XX)
      YI=YY+F*(WY-YY)
      CALL PL (XI,YI,((KK.EQ.1).AND.VV))
      IF (EQ(XX,XI).OR.EQ(XI,EX))  GO TO 43
      IF ((KK.EQ.2).AND.EQ(U1,U2)) GO TO 43
      J1=II(J1)
      X1(J1)=XI
      T1(J1)=TT+F*(TO-TT)
      B1(J1)=BB+F*(BO-BB)
   43 CONTINUE
   44 IF ((J1.LT.2).OR.(J1.GT.M-1)) GO TO 46
      IF (.NOT.VISSL(EX,TE,X1,T1,J1+K)) GO TO 46
      IF (     VISSL(EX,BE,X1,B1,J1+K)) GO TO 48
   46 J1=II(J1)
   48 X1(J1)=EX
      T1(J1)=TE
      B1(J1)=BE
      VV=PO(VT).OR.PO(-VB)
      CALL PL (EX,WY,VV)
      L=.FALSE.
      ET=DT
      EB=DB
      XX=EX
      YY=WY
      TT=TO
      BB=BO
      GO TO 30
 
C === TERMINATION
 
C --- IF THE FUNCTION IS EXHAUSTED BEFORE THE BOUNDS, COPY THEM.
 
   50 IF ((J0.LT.1).OR.(J0.GT.N0)) GO TO 70
      J1=II(J1)
      X1(J1)=X0(J0)
      T1(J1)=T0(J0)
      B1(J1)=B0(J0)
      J0=J0+I
      GO TO 50
 
C --- IF THE BOUNDS ARE EXHAUSTED BEFORE THE FUNCTION, COPY THE
C     REMAINING PART OF THE FUNCTION, WHICH WILL BE VIRIBLE.
 
   60 CALL PL (EX,WY,.FALSE.)
   61 IF ((J.LT.1).OR.(J.GT.N)) GO TO 70
      CALL PL (X(J),Y(J),.TRUE.)
      J1=II(J1)
      X1(J1)=X(J)
      T1(J1)=Y(J)
      B1(J1)=Y(J)
      J=J+I
      GO TO 61
 
C --- COPY THE NEW BOUNDS OVER THE OLD ONES, SHIFTING THEM AS NECESSARY.
 
   70 N0=((M+1)*(1-I))/2+I*J1
      J1=(J1+1-I*(J1-1))/2
      DO 71 J0=1,N0
      X0(J0)=X1(J1)
      T0(J0)=T1(J1)
      B0(J0)=B1(J1)
   71 J1=J1+1
      RETURN
 
      END
      SUBROUTINE  VISDC (Z1,ZE,Z2,NZ,NX,MX,NY,MY,US,VS,L,PL)
 
C     [DIAGONAL CONTOURED SEQUENCE]
C     ZE(J,I) ARRAY OF FUNCTION VALUES
C     (Z1,Z2) SPAN OF Z VALUES
C     NX,NY   RANGES OF J AND I
C     MX,MY   MAXIMA ATTAINABLE BY J AND I
C     US,VS   TOTAL SHEARS IN U AND V DIRECTIONS
C     L       DIRECTION OF VIEW (1=WEST, -1=EAST)
C     PL      PEN MOVEMENT SUBROUTINE
C     [16-MAY-74]
 
      EXTERNAL	  PL
      DIMENSION   ZE(1)
      DIMENSION   A(501),B(501)
      DIMENSION   D(501),E(501),G(501)
      DIMENSION   U(501),V(501),W(501)
      DATA	  M,MK/1,501/
 
      IX(J,I)=(I-1)*MX+J
      SC(Z)=ZS*(Z-Z1)
 
      NA=0
      ND=0
      L=L*M
      MM=1
      EL=FLOAT(L)
      EM=FLOAT(M)
      TE=0.5*(EL+1.0)
      ZS=(1.0-VS)/(Z2-Z1)
      DZ=(Z2-Z1)/FLOAT(NZ-1)
      DUI=-(EL*US)/FLOAT(MY-1)
      DUJ=(1.0-US)/FLOAT(MX-1)
      DVI=VS/FLOAT(MY-1)
 
      I0=(NY+1-M*(NY-3))/2
      J0=(NX+1-L*(NX-1))/2
      K0=((MK+1)*(1-N))/2
   10 K=K0
      I=MAX0(MIN0(I0,NY+1),0)
      J=J0
      EU=TE*US+DUI*FLOAT(I-1)+DUJ*FLOAT(J-1)
      VE=      DVI*FLOAT(I-1)
   20 IF ((I.LT.1).OR.(I.GT.NY)) GO TO 22

      K=MAX0(MIN0(K+N,MK),1)

      U(K)=EU
      V(K)=VE+SC(ZE(IX(J,I)))
   22 I=I-M
      EU=EU-DUI
      VE=VE-EM*DVI
      IF ((I.LT.1).OR.(I.GT.NY)) GO TO 30
      K=MAX0(MIN0(K+N,MK),1)
	  U(K)=EU
      V(K)=VE+SC(ZE(IX(J,I)))
      J=J+L
      EU=EU+EL*DUJ
      IF ((J.GE.1).AND.(J.LE.NX)) GO TO 20
   30 ZI=Z1
      DO 60 IZ=1,NZ
      K=K0
	I=MAX0(MIN0(I0,NY+1),0)
      J=J0
	EU=TE*US+DUI*FLOAT(I-1)+DUJ*FLOAT(J-1)
      VE=      DVI*FLOAT(I-1)
   40 IF ((I.LT.1).OR.(I.GT.NY)) GO TO 42
      K=MAX0(MIN0(K+N,MK),1)
      W(K)=VE+SC(ZI)
   42 I=I-M
      EU=EU-DUI
      VE=VE-EM*DVI
      IF ((I.LT.1).OR.(I.GT.NY)) GO TO 50
      K=MAX0(MIN0(K+N,MK),1)
      W(K)=VE+SC(ZI)
      J=J+L
      EU=EU+EL*DUJ
      IF ((J.GE.1).AND.(J.LE.NX)) GO TO 40
   50 CALL VISRB (A,B,NA,MK,U,V,K,U,W,K,-1.0)
      CALL VISHH (D,E,G,ND,A,B,NA,MM,PL)
      MM=-MM
   60 ZI=ZI+DZ
      I0=I0+M
      IF ((I0.GE.0).AND.(I0.LE.NY+1)) GO TO 10
      J0=J0+L
      IF ((J0.GE.1).AND.(J0.LE.NX))   GO TO 10
      RETURN
 
      END
      SUBROUTINE  VISDO (Z1,S1,S2,Z2,NX,MX,NY,MY,US,VS,L,IS,PL)
 
C     [DOUBLE SURFACE]
C     S1,S2   ARRAYS CONTAINING THE TWO SURFACES
C     Z1,Z2   SPAN OF SURFACE VALUES
C     MX,MY   COMMON DIMENSION OF THE ARRAYS S1 AND S2
C     NX,NY   SECTIONS OF S1 AND S2 ACTUALLY USED
C     US,VS   TOTAL SHEARS IN U AND V DIRECTIONS
C     L       DIRECTION OF VIEW (1=WEST, -1=EAST)
C     IS      SEPARATION OPTION (1=YES , -1=NO)
C     PL      PEN MOVEMENT SUBROUTINE
C     [15-MAY-74]
 
      EXTERNAL	  PL
      LOGICAL	  P,Q
      DIMENSION   S1(1),S2(1)
      DIMENSION   X1(351),T1(351),B1(351)
      DIMENSION   X2(351),T2(351),B2(351)
      DIMENSION   D(351),E(351),G(351),H(351)
      DIMENSION   A(351),B(351)
      DIMENSION   U(201),V1(201),V2(201)
      DATA	  M,MK,MA/1,201,351/
 
      IX(J,I)=(I-1)*MX+J
      SC(Z)=ZS*(Z-Z1)
 
      N1=0
      N2=0
      N=L*M
      EF=1.0-VS
      EL=FLOAT(L)
      EM=FLOAT(M)
      ZS=EF/(Z2-Z1)
      TE=0.5*(EL+1.0)
      DUI=-(EL*US)/FLOAT(MY-1)
      DUJ=(1.0-US)/FLOAT(MX-1)
      DVI=VS/FLOAT(MY-1)
 
      I0=(NY+1-M*(NY-3))/2
      J0=(NX+1-L*(NX-1))/2
      K0=((MK+1)*(1-N))/2
   10 K=K0
      I=MAX0(MIN0(I0,NY+1),0)
      J=J0
      EU=TE*US+DUI*FLOAT(I-1)+DUJ*FLOAT(J-1)
      VE=      DVI*FLOAT(I-1)
   20 IF ((I.LT.1).OR.(I.GT.NY)) GO TO 22
      K=MAX0(MIN0(K+N,MK),1)
      U(K)=EU
      V1(K)=VE+SC(S1(IX(J,I)))
      V2(K)=VE+SC(S2(IX(J,I)))
   22 I=I-M
      EU=EU-DUI
      VE=VE-EM*DVI
      IF ((I.LT.1).OR.(I.GT.NY)) GO TO 30
      K=MAX0(MIN0(K+N,MK),1)
      U(K)=EU
      V1(K)=VE+SC(S1(IX(J,I)))
      V2(K)=VE+SC(S2(IX(J,I)))
      J=J+L
      EU=EU+EL*DUJ
      IF ((J.GE.1).AND.(J.LE.NX)) GO TO 20
 
   30 P=L.LT.0
      Q=L.GT.0
      IF (P) KK=MK-K+1
      IF (Q) CALL VISRB (A,B,NA,MA,U,V1,K,U,V2,K,1.0)
      IF (P) CALL VISRB (A,B,NA,MA,U(K),V1(K),KK,U(K),V2(K),KK,1.0)
      IF (Q) CALL VISRB (G,H,NG,MA,U,V1,K,U,V2,K,-1.0)
      IF (P) CALL VISRB (G,H,NG,MA,U(K),V1(K),KK,U(K),V2(K),KK,-1.0)
 
      IF (IS.LT.0) GO TO 40
      DO 36 II=1,NA
   36 B(II)=EF*B(II)+VS
      DO 38 II=1,NG
   38 H(II)=EF*H(II)
 
   40 CALL VISRB (D,E,ND,MA,A,B,NA,X1,T1,N1,1.0)
      CALL VISHH (X2,T2,B2,N2,D,E,ND,1,PL)
      CALL VISRB (D,E,ND,MA,G,H,NG,X2,B2,N2,-1.0)
      CALL VISHH (X1,T1,B1,N1,D,E,ND,-1,PL)
 
      I0=I0+M
      IF ((I0.GE.0).AND.(I0.LE.NY+1)) GO TO 10
      J0=J0+L
      IF ((J0.GE.1).AND.(J0.LE.NX))   GO TO 10
      RETURN
 
      END
      SUBROUTINE  VISDS (Z1,ZE,Z2,J1,J2,MX,I1,I2,MY,US,VS,L,M,PL)
 
C     [DIAGONAL SEQUENCE]
C     ZE(J,I) ARRAY OF FUNCTION VALUES
C     (Z1,Z2) SPAN OF Z VALUES
C     NX,NY   RANGES OF J AND I
C     MX,MY   MAXNMA ATTAINABLE BY J AND I
C     US,VS   TOTAL SHEARS IN U AND V DIRECTIONS
C     L       DIRECTION OF VIEW (1=WEST, -1=EAST)
C     PL      PEN MOVEMENT SUBROUTINE
C     [06-OCT-74]
 
      EXTERNAL	  PL
      DIMENSION   ZE(1),U(501),V(501)
      DATA	  MK/501/
 
      IX(J,I)=(I-1)*MX+J
      SC(Z)=ZS*(Z-Z1)
 
      NL=ISIGN(1,L)
      NM=ISIGN(1,M)
      N=NL*NM
      IL=I1-IABS(M)
      IU=I2+IABS(M)
      MM=1
      TE=0.5*FLOAT(NL+1)
      ZS=(1.0-VS)/(Z2-Z1)
      DUI=-(FLOAT(NL)*US)/FLOAT(MY-1)
      EUI=DUI*FLOAT(M)
      DUJ=(1.0-US)/FLOAT(MX-1)
      EUJ=DUJ*FLOAT(L)
      DVI=VS/FLOAT(MY-1)
      EVI=DVI*FLOAT(M)
      I0=(I2+I1-NM*(I2-I1))/2+M
      J0=(J2+J1-NL*(J2-J1))/2
      K0=((MK+1)*(1-N))/2
   10 K=K0
      I=MAX0(MIN0(I0,IU),IL)
      J=J0
      EU=TE*US+DUI*FLOAT(I-1)+DUJ*FLOAT(J-1)
      VE=      DVI*FLOAT(I-1)
   20 IF ((I.LT.I1).OR.(I.GT.I2)) GO TO 22
      K=MAX0(MIN0(K+N,MK),1)
      U(K)=EU
      V(K)=VE+SC(ZE(IX(J,I)))
   22 I=I-M
      EU=EU-EUI
      VE=VE-EVI
      IF ((I.LT.I1).OR.(I.GT.I2)) GO TO 30
      K=MAX0(MIN0(K+N,MK),1)
      U(K)=EU
      V(K)=VE+SC(ZE(IX(J,I)))
      J=J+L
      EU=EU+EUJ
      IF ((J.GE.J1).AND.(J.LE.J2)) GO TO 20
   30 KK=(K+1-N*(K-1))/2
      CALL VISHO (U(KK),V(KK),(MK+1-N*(MK-2*K+1))/2,MM,PL)
      MM=-MM
      I0=I0+M
      IF ((I0.GE.IL).AND.(I0.LE.IU)) GO TO 10
      J0=J0+L
      IF ((J0.GE.J1).AND.(J0.LE.J2)) GO TO 10
      RETURN
 
      END
      SUBROUTINE  VISES (Z1,ZE,Z2,X1,X2,NX,E1,E2,NE,L,M)
 
C     [ELLIPTICAL SEQUENCE]
C     (Z1,Z2)  RANGE OF ZE
C     (X1,X2)  RANGE OF XI
C     (E1,E2)  RANGE OF ETA
C     NX       NUMBER OF XI VALUES
C     NE       NUMBER OF ETA VALUES
C     L= 1  WESTERN VIEW,  L=-1  EASTERN VIEW
C     M= 1  SOUTHERN VIEW, M=-1  NORTHERN VIEW
C     [05-OCT-74]
 
      EXTERNAL	  PLTCA
      DIMENSION   ZE(1),U(501),V(501)
 
      IX(J,I)=(I-1)*NX+J
      KI(K)=MAX0(MIN0(K+N,MK),1)
      SC(Z)=ZS*(Z-Z1)+0.2
 
      N=L*M
      MM=1
      MK=501
      EL=FLOAT(L)
      EM=FLOAT(M)
      DX=(X2-X1)/FLOAT(NX-1)
      DE=(E2-E1)/FLOAT(NE-1)
      ZS=0.58/(Z2-Z1)
 
      I0=(NE+1-M*(NE-3))/2
      J0=((NX+1)-L*(NX-1))/2
   10 K=((MK+1)*(1-N))/2
      I=MAX0(MIN0(I0,NE+1),0)
      J=H0
      E=E1+DE*FLOAT(I-1)
      X=X1+DX*FLOAT(J-1)
   20 IF ((I.LT.1).OR.(I.GT.NE)) GO TO 22
      K=KI(K)
      U(K)=0.166*(COSH(X)*COS(E)+3.0)
      V(K)=SC(ZE(IX(J,I)))+0.075*SINH(X)*SIN(E)
   22 I=I-M
      E=E-EM*DE
      IF ((I.LT.1).OR.(I.GT.NE)) GO TO 30
      K=KI(K)
      U(K)=0.166*(COSH(X)*COS(E)+3.0)
      V(K)=SC(ZE(IX(J,I)))+0.075*SINH(X)*SIN(E)
      J=J+L
      X=X+EL*DX
      IF ((J.GE.1).AND.(J.LE.NX)) GO TO 20
   30 IF (N.GT.0)  CALL VISHO (U,V,K,MM,PLTCA)
      IF (N.LT.0)  CALL VISHO (U(K),V(K),MK-K+1,MM,PLTCA)
      MM=-MM
      I0=I0+M
      IF ((I0.GE.0).AND.(I0.LE.NE+1)) GO TO 10
      J0=J0+L
      IF ((J0.GE.1).AND.(J0.LE.NX))   GO TO 10
      RETURN
 
      END
      SUBROUTINE  VISHH (X0,T0,B0,N0,X,Y,N,I,PL)
 
C     [HALF HORIZON]
C     SOME OF THE HIDDEN LINE SUBROUTINES EMPLOY MORE THAN ONE HORIZON,
C     WHICH MEANS THAT THE ARRAYS CONTAINING THESE HORIZONS MUST APPEAR
C     AS EXPLICIT ARGUMENTS IN THE UPDATING SUBROUTINES. NEVERTHELESS,
C     THREE OF THE ARGUMENTS OF VISBO ARE NOTHING BUT WORKING ARRAYS
C     WHICH CAN STILL BE REMOVED FROM THE CALLING PROGRAMS IF THEY ARE
C     PLACED IN AN INTERMEDIATE SUBROUTINE SUCH AS THIS ONE.
C     X0(N0)  ARRAY OF ARGUMENTS FOR THE HORIZON
C     T0(N0)  ARRAY OF VALUES OF THE UPPER HORIZON
C     B0(N0)  ARRAY OF VALUES OF THE LOWER HORIZON
C     X(N)    ARRAY OF ARGUMENTS
C     Y(N)    ARRAY OF FUNCTION VALUES
C     I       PEN DIRECTION (1=FORWARD, -1=BACKWARD)
C     PL      PEN MOVEMENT SUBROUTINE
C     [29-MAY-74]
 
      EXTERNAL	  PL
      DIMENSION   X(1),Y(1)
      DIMENSION   X0(1),T0(1),B0(1)
      DIMENSION   X1(701),T1(701),B1(701)
      DATA	  M/701/
 
      CALL VISBO (X1,T1,B1,M,X0,T0,B0,N0,X,Y,N,I,PL)
      RETURN
 
      END
      SUBROUTINE  VISHO (X,Y,N,I,PL)
 
C     [HORIZONS]
C     SUBROUTINE WHICH UPDATES THE HORIZONS FOR THE HIDDEN LINE
C     SUBROUTINES.  THE ACTUAL WORK IS DONE BY VISBO, BUT VISHO
C     CONTAINS THE WORK ARRAYS NEEDED IN THE SIMPLEST APPLICATIONS,
C     AVOIDING THE NECESSITY TO DEFINE THEM SEPARATELY FOR EACH
C     INDIVIDUAL PROGRAM.
C     X(N)   ARRAY OF ARGUMENTS
C     Y(N)   ARPAY OF FUNCTION VALUES
C     I      DIRECTION OF PEN MOVEMENT
C     PL     PEN MOVEMENT SU@ROUTINE
C     [29-MAY-74]
 
      EXTERNAL	  PL
      DIMENSION   X(1),Y(1)
      DIMENSION   X0(701),T0(701),B0(701)
      DIMENSION   X1(701),T1(701),B1(701)
      COMMON/VIS/ N0
      DATA	  M/701/
 
      CALL VISBO (X1,T1,B1,M,X0,T0,B0,N0,X,Y,N,I,PL)
      RETURN
 
      END
      FUNCTION	  VISLI (Z,X,Y,I)
 
C     [LINEAR INTERPOLATION]
C     PERFORM THE BACKWARD LINEAR INTERPOLATIONS REQUIRED BY THE HORIZON
C     ROUTINES. VISLI RECEIVES SUCH INTENSE USAGE THAT IT DELIBERATELY
C     ESCHEWS A CHECK FOR A ZERO DENOMINATOR, WHICH STILL OCCASIONALLY
C     CREATES OVERFLOWS.
C     Z     POINT AT WHICH INTERPOLATION IS MADE
C     X     ARRAY IN WHICH Z IS INTERPOLATED
C     Y     ARRAY FROM WHICH TO TAKE THE INTERPOLATED VALUE
C     I     AN INDEX FOR WHICH X(I-1).LE.Z.LE.X(I)
C     [05-MAY-74]
 
      DIMENSION   X(1),Y(1)
 
      X1=X(I-1)
      X2=X(I)
      Y1=Y(I-1)
      Y2=Y(I)
      VISLI=Y1+(Y2-Y1)*((Z-X1)/(X2-X1))
      RETURN
 
      END
      SUBROUTINE  VISNH
 
C     [NULL HORIZON]
C     SETS UP THE NULL INITIAL HORIZON WHICH MOST OF THE HIDDEN LINE
C     SUBROUTINES REQUIRE.
C     [22-NOV-74]
 
      COMMON/VIS/ N0
 
      N0=0
      RETURN
 
      END
      SUBROUTINE  VISNP (PH,TH,JP,IT,NP,NT,O)
 
C     [NEAREST POINT]
C     DETERMINE THE COORDINATES OF THE LINE OF SIGHT JOINING THE
C     OBSERVER TO THE CENTER OF THE SPHERICAL COORDINATES.
C     (PH,TH)  SURFACE COORDINATES OF NEAREST POINT
C     (JP,IT)  INDICES OF NEAREST POINT
C     (NP,NT)  RANGE OF PHI, THETA INDICES
C     O(3,3)   ORTHOGONAL MATRIX DEFINING ROTATION
C     [16-JUN-74]
 
      DIMENSION   O(3,3)
 
      O31=O(3,1)
      O32=O(3,2)
      PH=ATAN2(O32,O31)
      IF (PH.LT.0.0) PH=6.28318+PH
      RH=SQRT(O31*O31+O32*O32)
      TH=ATAN2(RH,O(3,3))
      JP=1+IFIX(0.5+(FLOAT(NP)*PH)/6.28318)
      IT=MIN0(MAX0(IFIX(1.5+(TH*FLOAT(NT-1))/3.14159),1),NT)
      RETURN
 
      END
      SUBROUTINE  VISPS (Z1,ZE,Z2,R1,R2,NR,P1,P2,NP,L,M)
 
C     [POLAR SEQUENCE]
C     L= 1  WESTERN VIEW,  L=-1  EASTERN VIEW
C     M= 1  SOUTHERN VIEW, M=-1  NORTHERN VIEW
C     [05-OCT-74]
 
      EXTERNAL	  PLTCA
      DIMENSION   ZE(1),U(501),V(501)
 
      IX(J,I)=(I-1)*NR+J
      SC(Z)=0.167+ZS*(Z-Z1)
 
      N=L*M
      MM=1
      MK=501
      EL=FLOAT(L)
      EM=FLOAT(M)
      ZS=0.667/(Z2-Z1)
      DP=(P2-P1)/FLOAT(NP-1)
      DR=(R2-R1)/FLOAT(NR-1)
 
      I0=(NP+1-M*(NP-3))/2
      J0=((NR+1)-L*(NR-1))/2
      K0=((MK+1)*(1-N))/2
   10 K=K0
      I=MAX0(MIN0(I0,NP+1),0)
      J=J0
      P=P1+DP*FLOAT(I-1)
      R=R1+DR*FLOAT(J-1)
   20 IF ((I.GT.NP).OR.(I.LT.1)) GO TO 22
      K=MAX0(MIN0(K+N,MK),1)
      U(K)=0.5*(1.0+R*COS(P))
      V(K)=SC(ZE(IX(J,I)))+0.167*R*SIN(P)
   22 I=I-M
      P=P-EM*DP
      IF ((I.GT.NP).OR.(I.LT.1)) GO TO 30
      K=MAX0(MIN0(K+N,MK),1)
      U(K)=0.5*(1.0+R*COS(P))
      V(K)=SC(ZE(IX(J,I)))+0.167*R*SIN(P)
      J=J+L
      R=R+EL*DR
      IF ((J.LE.NR).AND.(J.GE.1)) GO TO 20
   30 IF (N.GT.0) CALL VISHO (U,V,K,MM,PLTCA)
      IF (N.LT.0) CALL VISHO (U(K),V(K),MK-K+1,MM,PLTCA)
      MM=-MM
      I0=I0+M
      IF ((I0.GE.0).AND.(I0.LE.NP+1)) GO TO 10
      J0=J0+L
      IF ((J0.GE.1).AND.(J0.LE.NR)) GO TO 10
      RETURN
 
      END
      SUBROUTINE  VISRB (X,Y,J,M,X1,Y1,N1,X2,Y2,N2,S)
 
C     [RESTRICTED BOUND]
C     THE BOUND IS RESTRICTED TO THE INTERVAL WHERE THE FIRST
C     FUNCTION IS DEFINED.
C     X(M),Y(M)     ARRAYS FOR THE BOUND OF TWO FUNCTIONS
C     J,M	    ACTUAL DIMENSION, MAXIMUM DIMENSION OF X,Y
C     X1(N1),Y1(N1) FIRST  FUNCTION
C     X2(N2),Y2(N2) SECOND FUNCTION
C     S 	    TYPE OF BOUND (S=1.0,UPPER; S=-1.0,LOWER)
C     [15-MAY-74]
 
      LOGICAL	  EQ,VISSL
      DIMENSION   X(1),Y(1),X1(1),Y1(1),X2(1),Y2(1)
 
      EQ(X,Y)=ABS(Y-X).LE.1.0E-5
 
      L=.TRUE.
      J=0
      J1=1
      J2=1
      Z1=X1(J1)
      Z2=X2(J2)
      IF (N1.LE.1)   RETURN
      IF (N2.LE.1)   GO TO 60
      IF (EQ(Z1,Z2)) GO TO 32
      IF (Z1-Z2)     10,32,20
 
   10 J=MIN0(J+1,M)
      X(J)=Z1
      Y(J)=Y1(J1)
      J1=J1+1
      Z1=X1(J1)
      IF (EQ(Z1,Z2)) GO TO 32
      IF (Z1-Z2)     10,32,32
 
   20 J2=J2+1
      Z2=X2(J2)
      IF (EQ(Z1,Z2)) GO TO 32
      IF (Z1-Z2)     32,32,20
 
   30 IF (J1.GT.N1) RETURN
      IF (J2.GT.N2) GO TO 60
      Z1=X1(J1)
      Z2=X2(J2)
   32 Z=AMIN1(Z1,Z2)
      W1=VISLI(Z,X1,Y1,MAX0(J1,2))
      W2=VISLI(Z,X2,Y2,MAX0(J2,2))
      IF (EQ(Z,Z1)) J1=J1+1
      IF (EQ(Z,Z2)) J2=J2+1
 
      W=S*AMAX1(S*W1,S*W2)
      D=W1-W2
      IF (L.OR.(J.LE.1)) GO TO 42
      IF ((EQ(D,0.0)).OR.(EQ(E,0.0))) GO TO 41
      IF (SIGN(1.0,D).EQ.SIGN(1.0,E)) GO TO 41
      J=MIN0(J+1,M)
      X(J)=U-E*((Z-U)/(D-E))
      Y(J)=WW+(X(J)-U)*((W1-WW)/(Z-U))
   41 IF (VISSL(Z,W,X,Y,J)) GO TO 43
   42 J=MIN0(J+1,M)
   43 X(J)=Z
      Y(J)=W
      L=.FALSE.
      E=D
      U=Z
      WW=W1
      GO TO 30
 
   60 IF (J1.GT.N1) RETURN
      J=MIN0(J+1,M)
      X(J)=X1(J1)
      Y(J)=Y1(J1)
      J1=J1+1
      GO TO 60
 
      END
      SUBROUTINE  VISRS (Z1,ZE,Z2,NX,MX,NY,MY,TH,PL)
 
C     [ROTATED SEQUENCE]
C     XE(J,I) ARRAY OF FUNCTION VALUES
C     (Z1,Z2) SPAN OF Z VALUES
C     NX,NY   RANGES OF J AND I
C     MX,MY   MAXIMA ATTAINABLE BY J AND I
C     TH      ANGLE OF ROTATION (DEGREES, CLOCKWISE).
C     L       DIRECTION OF VIEW (1=WEST, -1=EAST)
C     M       DIRECTION OF VIEW (1=SOUTH, -1=NORTH)
C     PL      PEN MOVEMENT SUBROUTINE
C     THE HORIZONTAL SCALE IS NOT ALWAYS CONSTANT, BUT IS ADJUSTED
C     SO THAT THE DRAWING WILL OCCUPY THE FULL BREADTH OF THE PAGE.
C     [19-DEC-74]
 
      EXTERNAL	  PL
      DIMENSION   ZE(1),U(501),V(501)
      DATA	  M,MK,VF/1,501,0.333/
 
      IX(J,I)=(I-1)*MX+J
      SC(Z)=ZS*(Z-Z1)
 
      IF (TH.LT.0.0) L=-1
      IF (TH.GE.0.0) L= 1
      N=L*M
      MM=1
      SI=SIND(TH)
      CO=COSD(TH)
      EL=FLOAT(L)
      EM=FLOAT(M)
      ZS=(1.0-VF)/(Z2-Z1)
      SF=1.0/(ABS(CO)+ABS(SI))
      U0=0.25*(EL+1.0)*(1.0-SF*(CO-SI))
      V0=0.15*(EL-1.0)*SF*SI
      DUI=-(SF*SI)/FLOAT(MY-1)
      DUJ= (SF*CO)/FLOAT(MX-1)
      DVI= (VF*SF*CO)/FLOAT(MY-1)
      DVJ= (VF*SF*SI)/FLOAT(MX-1)
 
      I0=(NY+1-M*(NY-3))/2
      J0=(NX+1-L*(NX-1))/2
      K0=((MK+1)*(1-N))/2
   10 K=K0
      I=MAX0(MIN0(I0,NY+1),0)
      J=J0
      EU=U0+DUI*FLOAT(I-1)+DUJ*FLOAT(J-1)
      VE=V0+DVI*FLOAT(I-1)+DVJ*FLOAT(J-1)
   20 IF ((I.LT.1).OR.(I.GT.NY)) GO TO 22
      K=MAX0(MIN0(K+N,MK),1)
      U(K)=EU
      V(K)=VE+SC(ZE(IX(J,I)))
   22 I=I-M
      EU=EU-EM*DUI
      VE=VE-EM*DVI
      IF ((I.LT.1).OR.(I.GT.NY)) GO TO 30
      K=MAX0(MIN0(K+N,MK),1)
      U(K)=EU
      V(K)=VE+SC(ZE(IX(J,I)))
      J=J+L
      EU=EU+EL*DUJ
      VE=VE+EL*DVJ
      IF ((J.GE.1).AND.(J.LE.NX)) GO TO 20
   30 IF (L.GT.0) CALL VISHO (U,V,K,MM,PL)
      IF (L.LT.0) CALL VISHO (U(K),V(K),MK-K+1,MM,PL)
      MM=-MM
      I0=I0+M
      IF ((I0.GE.0).AND.(I0.LE.NY+1)) GO TO 10
      J0=J0+L
      IF ((J0.GE.1).AND.(J0.LE.NX))   GO TO 10
      RETURN
 
      END
      LOGICAL FUNCTION	VISSL (EX,WY,X,Y,I)
 
C     [STRAIGHT LINE]
C     [15-MAY-74]
 
      DIMENSION   X(1),Y(1)
 
      X1=X(I-1)
      X2=X(I)
      Y1=Y(I-1)
      Y2=Y(I)
      D=(WY-Y1)*(X2-X1)-(EX-X1)*(Y2-Y1)
      VISSL=ABS(D).LT.(1.0E-10)
      RETURN
 
      END
      SUBROUTINE  VISSP (RHO,PHI,R,T,P,O)
 
C     [SPHERICAL PROJECTION]
C     DETERMINE THE PLANAR POLAR COORDINATES IN THE PLOTTER PAGE
C     OF A POINT PROJECTED FROM THE FUNCTIONAL SURFACE, WHICH IS
C     DEFINED IN SPHERICAL POLAR COORDINATES.
C     (RHO, PHI) PLANE POLAR COORDINATES
C     (R,T,P)	 SPHERICAL SURFACE COORDINATES
C     O(3,3)	 ORTHOGONAL MATRIX EXPRESSING THE FIGURE'S ROTATION
C     [20-MAY-74]
 
      DIMENSION   O(3,3)
 
      X=R*SIN(T)*COS(P)
      Y=R*SIN(T)*SIN(P)
      Z=R*COS(T)
      U=O(1,1)*X+O(1,2)*Y+O(1,3)*Z
      V=O(2,1)*X+O(2,2)*Y+O(2,3)*Z
      RHO=SQRT(U*U+V*V)
      PHI=ATAN2(V,U)/6.28318
      RETURN
 
      END
      SUBROUTINE  VISSS (FU,J1,J2,NP,I1,I2,NT,L,M,S,B,O,PR,PL)
 
C     [SPHERICAL SEQUENCE]
C     FU(NP,NT)  ARRAY OF FUNCTION VALUES
C     (J1,J2)	 INTERVAL OF PHI INDICES TO BE GRAPHED
C     (I1,I2)	 INTERVAL OF THETA INDICES TO BE GRAPHED
C     L= 1  WESTERN VIEW,  L=-1  EASTERN VIEW
C     M= 1  SOUTHERN VIEW, M=-1  NORTHERN VIEW
C     B= (ATAN2 CUT LINE DOES NOT FALL IN QUADRANT BEING GRAPHED)
C     O(3,3)	 ORTHOGONAL ROTATION MATRIX
C     PR	 PROJECTION SUBROUTINE
C     PL	 PEN MOVEMENT SUBROUTINE
C     [22-NOV-74]
 
      EXTERNAL	  PL
      LOGICAL	  B
      DIMENSION   FU(1),O(3,3)
      DIMENSION   AZ(501),RA(501),OR(501)
      DATA	  MK/501/
 
      TAN(X)=SIN(X)/COS(X)
      IX(J,I)=(I-1)*NP+MOD(NP+J-1,NP)+1
      ZP(X1,Y1,X2,Y2)=X1-Y1*((X2-X1)/(Y2-Y1))
      WH(T,P)=S*SIGN(1.0,1.57079-T)*(TAN(T)-TN*COS(P-PH))
 
      N=L*M
      NN=1
      EL=FLOAT(L)
      EM=FLOAT(M)
      DP=6.28/FLOAT(NP)
      DT=3.14/FLOAT(NT)
      CALL VISNP (PH,TH,JP,IT,NP,NT,O)
      TN=TAN(TH)
 
      I0=((I1+I2)+M*(I1-I2))/2+M
      J0=((J1+J2)+L*(J1-J2))/2
      K0=((MK+1)*(1-N))/2
   10 K=K0
      I=MAX0(MIN0(I0,I2+1),I1-1)
      J=J0
      T=TH+DT*FLOAT(I-IT)
      P=PH+DP*FLOAT(J-JP)
   20 IF ((I.GT.I2).OR.(I.LT.I1)) GO TO 22
      K=MAX0(MIN0(K+N,MK),1)
      CALL PR (RA(K),AZ(K),FU(IX(J,I)),T,P,O)
      OR(K)=WH(T,P)
   22 I=I-M
      T=T-EM*DT
      IF ((I.GT.I2).OR.(I.LT.I1)) GO TO 30
      K=MAX0(MIN0(K+N,MK),1)
      CALL PR (RA(K),AZ(K),FU(IX(J,I)),T,P,O)
      OR(K)=WH(T,P)
      J=J+L
      P=P+EL*DP
      IF ((J.LE.J2).AND.(J.GE.J1)) GO TO 20
 
   30 IF (N.GT.0) GO TO 32
      M1=MK-K+1
      DO 31 MM=1,M1
      AZ(MM)=AZ(K+MM-1)
      RA(MM)=RA(K+MM-1)
   31 OR(MM)=OR(K+MM-1)
      K=M1
   32 IF (K.LE.1) GO TO 50
      IF (B) GO TO 36
      DO 35 MM=1,K
   35 IF (AZ(MM).LT.0.0) AZ(MM)=AZ(MM)+1.0
   36 IF (AZ(1).LE.AZ(K)) GO TO 38
      DO 37 MM=1,K
      T1=AZ(MM)
      T2=RA(MM)
      T3=OR(MM)
      AZ(MM)=AZ(K-MM+1)
      RA(MM)=RA(K-MM+1)
      OR(MM)=OR(K-MM+1)
      AZ(K-MM+1)=T1
      RA(K-MM+1)=T2
   37 OR(K-MM+1)=T3
   38 K1=0
   40 K1=K1+1
      IF (K1.GE.K) GO TO 50
      IF (OR(K1).GT.0.0) GO TO 40
      K2=K1
      IF (K1.LE.1) GO TO 41
      K1=K1-1
      AZ(K1)=ZP(AZ(K1),OR(K1),AZ(K1+1),OR(K1+1))
      RA(K1)=ZP(RA(K1),OR(K1),RA(K1+1),OR(K1+1))
   41 IF (K2.GT.K) GO TO 43
      IF (OR(K2).GT.0.0) GO TO 42
      K2=K2+1
      GO TO 41
   42 IF (K2.GT.K) GO TO 43
      AZ(K2)=ZP(AZ(K2),OR(K2),AZ(K2-1),OR(K2-1))
      RA(K2)=ZP(RA(K2),OR(K2),RA(K2-1),OR(K2-1))
      K2=K2+1
   43 IF (K2-K1.LT.2) GO TO 44
      IF (AZ(K1).GE.AZ(K1+1)) AZ(K1)=AZ(K1+1)-0.0025
      IF (AZ(K2-1).LE.AZ(K2-2)) AZ(K2-1)=AZ(K2-2)+0.0025
   44 IF (K2-K1.GT.1) CALL VISHO (AZ(K1),RA(K1),K2-K1,NN,PL)
      K1=K2
      GO TO 40
 
   50 NN=-NN
      I0=I0+M
      IF ((I0.GE.I1-1).AND.(I0.LE.I2+1)) GO TO 10
      J0=J0+L
      IF ((J0.GE.J1).AND.(J0.LE.J2)) GO TO 10
      RETURN
 
      END
      SUBROUTINE  VISTR (Z1,S1,S2,S3,Z2,NX,MX,NY,MY,US,VS,VD,L,IS,PL)
 
C     [TRIPLE SURFACE]
C     US,VS   TOTAL SHEARS IN U AND V DIRECTIONS
C     VD      VERTICAL DISPLACEMENT BETWEEN FUNCION SEGMENTS
C     L       DIRECTION OF VIEW (1=WEST, -1=EAST)
C     IS      SEPARATION OPTION (1=YES,  -1=NO)
C     PL      PEN MOVEMENT SUBROUTINE
C     [16-MAY-74]
 
      EXTERNAL	  PL
      DIMENSION   S1(1),S2(1),S3(1)
      DIMENSION   G1(275),G2(275),G3(275),H1(275),H2(275),H3(275)
      DIMENSION   U1(275),U2(275),U3(275),F1(275),F2(275),F3(275)
      DIMENSION   X1(275),X2(275),X3(275)
      DIMENSION   B1(275),T1(275),B2(275),T2(275),B3(275),T3(275)
      DIMENSION   U(275),V1(275),V2(275),V3(275)
      DATA	  M,MK/1,275/
 
      IX(J,I)=(I-1)*MX+J
      SC(Z)=ZS*(Z-Z1)
 
      N1=0
      N2=0
      N3=0
      N=L*M
      DD=2.0*VD
      EF=1.0-2.0*VD
      EL=FLOAT(L)
      EM=FLOAT(M)
      TE=0.5*(EL+1.0)
      ZS=(1.0-VS)/(Z2-Z1)
      DUI=-(EL*US)/FLOAT(MY-1)
      DUJ=(1.0-US)/FLOAT(MX-1)
      DVI=VS/FLOAT(MY-1)
 
      I0=(NY+1-M*(NY-3))/2
      J0=(NX+1-L*(NX-1))/2
      K0=((MK+1)*(1-N))/2
   10 K=K0
      I=MAX0(MIN0(I0,NY+1),0)
      J=J0
      EU=TE*US+DUI*FLOAT(I-1)+DUJ*FLOAT(J-1)
      VE=      DVI*FLOAT(I-1)
   20 IF ((I.LT.1).OR.(I.GT.NY)) GO TO 22
      K=MAX0(MIN0(K+N,MK),1)
      U(K)=EU
      V1(K)=VE+SC(S1(IX(J,I)))
      V2(K)=VE+SC(S2(IX(J,I)))
      V3(K)=VE+SC(S3(IX(J,I)))
   22 I=I-M
      EU=EU-DUI
      VE=VE-EM*DVI
      IF ((I.LT.1).OR.(I.GT.NY)) GO TO 30
      K=MAX0(MIN0(K+N,MK),1)
      U(K)=EU
      V1(K)=VE+SC(S1(IX(J,I)))
      V2(K)=VE+SC(S2(IX(J,I)))
      V3(K)=VE+SC(S3(IX(J,I)))
      J=J+L
      EU=EU+EL*DUJ
      IF ((J.GE.1).AND.(J.LE.NX)) GO TO 20
 
   30 IF (L.GT.0) GO TO 32
      NK=MK-K+1
      DO 31 KK=1,NK
      II=K+KK-1
      U(KK)=U(II)
      V1(KK)=V1(II)
      V2(KK)=V2(II)
      V3(KK)=V3(II)
   31 II=II+1
      K=NK
 
   32 CALL VISRB (G1,H1,NG1,MK,U,V1,K,U,V2,K,1.0)
      CALL VISRB (U3,F3,NU3,MK,G1,H1,NG1,U,V3,K,1.0)
      CALL VISRB (G1,H1,NG1,MK,U,V2,K,U,V3,K,-1.0)
      CALL VISRB (G2,H2,NG2,MK,U,V3,K,U,V1,K,-1.0)
      CALL VISRB (G3,H3,NG3,MK,U,V1,K,U,V2,K,-1.0)
      CALL VISRB (U1,F1,NU1,MK,G1,H1,NG1,G2,H2,NG2,1.0)
      CALL VISRB (U2,F2,NU2,MK,U1,F1,NU1,G3,H3,NG3,1.0)
      CALL VISRB (U1,F1,NU1,MK,G3,H3,NG3,U,V3,K,-1.0)
 
      IF (IS.LT.0) GO TO 40
      DO 36 KK=1,NU1
   36 F1(KK)=EF*F1(KK)
      DO 37 KK=1,NU2
   37 F2(KK)=EF*F2(KK)+VD
      DO 38 KK=1,NU3
   38 F3(KK)=EF*F3(KK)+DD
 
   40 CALL VISRB (G1,H1,NG1,MK,U1,F1,NU1,X2,B2,N2,-1.0)
      CALL VISHH (X1,B1,T1,N1,G1,H1,NG1, 1,PL)
      CALL VISRB (G1,H1,NG1,MK,U2,F2,NU2,X1,T1,N1, 1.0)
      CALL VISRB (G2,H2,NG2,MK,G1,H1,NG1,X3,B3,N3,-1.0)
      CALL VISHH (X2,B2,T2,N2,G2,H2,NG2,-1,PL)
      CALL VISRB (G1,H1,NG1,MK,U3,F3,NU3,X2,T2,N2, 1.0)
      CALL VISHH (X3,B3,T3,N3,G1,H1,NG1, 1,PL)
 
      I0=I0+M
      IF ((I0.GE.0).AND.(I0.LE.NY+1)) GO TO 10
      J0=J0+L
      IF ((J0.GE.1).AND.(J0.LE.NX))   GO TO 10
      RETURN
 
      END
      SUBROUTINE  VISTS (Z1,ZE,Z2,N,M)
 
C     [TRIANGULAR SEQUENCE]
C     ZE(M,M)  ARRAY OF VALUES
C     Z1,Z2    RANGE OF VALUES
C     N        ACTUAL  LENGTH OF COLUMN
C     M        MAXIMUM LENGTH OF COLUMN
C     [05-OCT-74]
 
      EXTERNAL	  PLTCA
      DIMENSION   ZE(1),U(501),V(501)
 
      IX(J,I)=(I-1)*M+J
      AM(J,I)=ZS*(ZE(IX(J,I))-Z1)
 
      MK=501
      HZ=0.35
      ZS=(2.0*HZ)/(Z2-Z1)
      DU=1.0/FLOAT(N-1)
      DV=0.3/FLOAT(N-1)
      HU=0.5*DU
      VE=0.0
      L=N-1
      DO 30 I=1,L
      K=0
      EU=HU*FLOAT(I-1)
      DO 10 J=I,N
      K=MIN0(K+1,MK)
      U(K)=EU
      V(K)=VE+AM(J,I)
   10 EU=EU+DU
      CALL VISHO (U,V,K,1,PLTCA)
      K=0
      EU=HU*FLOAT(I-1)
      DO 20 J=I,L
      K=MIN0(K+1,MK)
      U(K)=EU
      V(K)=VE+AM(J,I)
      K=MIN0(K+1,MK)
      U(K)=EU+HU
      V(K)=VE+DV+AM(J+1,I+1)
   20 EU=EU+DU
      K=MIN0(K+1,MK)
      U(K)=EU
      V(K)=VE+AM(N,I)
      CALL VISHO (U,V,K,-1,PLTCA)
   30 VE=VE+DV
      RETURN
 
      END

