C     ------------------------------------------------------------------


      FUNCTION    CARG (Z)

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

      COMPLEX     Z

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


C     ------------------------------------------------------------------


      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

   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


C     ------------------------------------------------------------------
C     ------------------------------------------------------------------


      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


C     ------------------------------------------------------------------


      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


C     ------------------------------------------------------------------
C     ------------------------------------------------------------------


      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




C     ------------------------------------------------------------------




      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 ORTHOGRAPHIC RELIEF.
C     Z0        CONTOUR LEVEL SOUGHT
C     (XF,YF)   LIGHTING DIRECTION FOR ORTHOGRAPHIC RELIEF
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     [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


C     ------------------------------------------------------------------
      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     FU        CABS OR CARG ACCORDING TO CONTOURS DESIRED
C     PL        PEN MOVEMENT SUBROUTINE, USUALLY PLTCA
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
C     ------------------------------------------------------------------


      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


C     ------------------------------------------------------------------



      SUBROUTINE  PLT00

C     [INITIALIZATION]
C     INITIALIZING SUBROUTINE TO START OFF A SERIES OF GRAPHS. CALLS
C     <PLOTS>, THEN MOVES THE PEN AT MOST 11" TO THE RIGHT TO INSURE
C     ITS PROPER POSITIONING.
C     [18-FEB-73]

      CALL PLOTS (I)
      CALL PLOT  (0.0,-11.0,-3)
      RETURN
      END


C     ------------------------------------------------------------------
      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

C     ------------------------------------------------------------------


      SUBROUTINE  PLTBH (X,Y,P)

C     [BOTTOM HALF]
C     SCALE THE CARTESIAN COORDINATES X,Y SO AS TO PLACE A GRAPH IN
C     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


C     ------------------------------------------------------------------


      SUBROUTINE  PLTBO

C     [BORDER]
C     SET UP AN 8-1/2" X 11" FRAME WITH AN INNER FRAME 1" INSIDE OF IT,
C     THEN LOCATE 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


C     ------------------------------------------------------------------
C     ------------------------------------------------------------------


      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     [10-MAY-75]

      DATA        SX,SY/5.50,4.25/

      CALL PLOTS (I)
      CALL PLOT  (0.0,0.0,2)
      CALL PLOT  (0.0,0.0,3)
      CALL PLOT  (-SY-0.02,SX,-3)
      RETURN
      END


C     ------------------------------------------------------------------
      SUBROUTINE  PLTBV (Z1,ZE,Z2,NX,NY,PL)

C     [BIRDSEYE VIEW]
C     [17-MAY-75]

      EXTERNAL    PL
      DIMENSION   ZE(1)
      DATA        HX,HY/4.50,3.25/

      F(I,J)=(R*ZS)/(Z0-ZE(I+NX*(J-1)))

      K=1
      R=5.0
      Z0=1.3*Z2
      ZS=0.125*(Z2-Z1)
      DX=(1.75*HX)/FLOAT(NX-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(NX-1)
      DY=-(1.75*HY)/FLOAT(NY-1)
      X= 0.875*HX
      Y= 0.875*HY
      K=-1
      DO 40 I=1,NX
      J1=((NY+1)-K*(NY-1))/2
      J2=((NY+1)+K*(NY-1))/2
      EF=F(NX-I+1,J1)
      CALL PLTMS (EF*X,EF*Y,.FALSE.)
      DO 30 J=J1,J2,K
      EF=F(NX-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 PLTMC (EX,WY,P)
      RETURN
      END


C     ------------------------------------------------------------------

      SUBROUTINE  PLTCI (X,Y,R,PL)

C     [CIRCLE]
C     DRAW A CIRCLE ON THE PLOTTER WITH CENTER (X,Y), RADIUS R. PL,
C     NORMALLY PLTCA, MOVES THE PEN.
C     [09-MAY-75]

      N=60
      DT=6.28318/FLOAT(N)
      TH=DT
      CALL PL (X+R,Y,.FALSE.)
      DO 10 I=1,N
      CALL PL (X+R*COS(TH),Y+R*SIN(TH),.TRUE.)
   10 TH=TH+DT
      RETURN
      END


C     ------------------------------------------------------------------

      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 SPECIFIED
C     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


C     ------------------------------------------------------------------


      SUBROUTINE  PLTEU (O,E1,E2,E3)

C     [EULER ANGLES]
C     O=R1*R2*R3.  R1 AND R3 ARE COUNTERCLOCKWISE ROTATIONS IN THE X-Y
C     PLANE THROUGH ANGLES E1 AND E3 RESPECTIVELY; R2 IS A ROTATION IN
C     THE Y-Z PLANE THROUGH THE COUNTERCLOCKWISE ANGLE E2.
C     O(3,3)      MATRIX IN WHICH ROTATION IS STORED
C     E1,E2,E3    EULER ANGELS (DEGREES) OF ROTATION

C     [30-MAY-75]

      DIMENSION   O(3,3)

      C1=COSD(E1)
      C2=COSD(E2)
      C3=COSD(E3)
      S1=SIND(E1)
      S2=SIND(E2)
      S3=SIND(E3)
      O(1,1)= C1*C3-S1*C2*S3
      O(1,2)=-C1*S3-S1*C2*C3
      O(1,3)= S1*S2
      O(2,1)= S1*C3+C1*C2*S3
      O(2,2)=-S1*S3+C1*C2*C3
      O(2,3)=-C1*S2
      O(3,1)= S2*S3
      O(3,2)= S2*C3
      O(3,3)= C2
      RETURN
      END
      SUBROUTINE  PLTEV (Z1,ZE,Z2,NX,NE,PL)

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     PL        PEN MOVEMENT SUBROUTINE, NORMALLY PLTCA
C     [17-MAY-75]

      EXTERNAL    PL
      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 VISNH
      CALL VISES (Z1,ZE(3*NN+1),Z2,X1,X2,NX,Q3,Q4,NQ,-1, 1,PL)
      CALL VISES (Z1,ZE        ,Z2,X1,X2,NX,Q0,S1,NQ, 1, 1,PL)
      CALL VISNH
      CALL VISES (Z1,ZE(2*NN+1),Z2,X1,X2,NX,Q2,S3,NQ,-1,-1,PL)
      CALL VISES (Z1,ZE(  NN+1),Z2,X1,X2,NX,Q1,Q2,NQ, 1,-1,PL)
      RETURN
      END
      SUBROUTINE  PLTFI (Y1,WY,Y2,N,PL)

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     PL     PEN MOVEMENT SUBROUTINE
C     [17-MAY-75]

      EXTERNAL    PL
      DIMENSION   WY(1)

      IF (N.LT.2) RETURN
      CALL PLTIG (1.0,Y1,1,PL)
      CALL PLTIG (FLOAT(N),Y2,2,PL)
      CALL PLTIG (1.0,WY(1),3,PL)
      DO 10 I=2,N
   10 CALL PLTIG (FLOAT(I),WY(I),4,PL)
      RETURN
      END


C     ------------------------------------------------------------------


      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 CENTER
C     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

C     ------------------------------------------------------------------

      SUBROUTINE  PLTGA (X1,X,X2,Y1,Y,Y2,N,PL)
C     [GRAPH ARRAY]
C     PLOT A GRAPH BY CONNECTING THE DATA POINTS BY STRAIGHT LINES. THE
C     POINTS DEFINING THE GRAPH ARE TAKEN FROM TWO ARRAYS, ONE HOLDING
C     THE X-VALUES AND THE OTHER CONTAINING THE Y-VALUES. THE RESPECTIVE
C     SCALES ARE INDICATED BY THE VALUES TO BE ASSIGNED TO THE MARGINS
C     OF THE GRAPH. ORDINARILY THE MARGINS WOULD BE GIVEN ROUNDED VALUES
C     SLIGHTLY LARGER THAN  THE EXTREME DATA VALUES.  HOWEVER, THE GRAPH
C     MAY BE CENTERED IN VARIOUS WAYS BY ASSIGNING ONE OR MORE MARGINS
C     CONSIDERABLY LARGER VALUES.  LIKEWISE EXCERPTS FROM THE GRAPH MAY
C     BE CHOSEN BY GIVING 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     PL     PEN MOVEMENT SUBROUTINE, USUALLY PLTCA
C     [07-JUN-75]
      DIMENSION   X(1),Y(1)
      EX(X)=(X-X1)*SCX-HX
      WY(Y)=(Y-Y1)*SCY-HY
      IF (N.LT.2) RETURN
      SCX=1.0/(X2-X1)
      SCY=1.0/(Y2-Y1)
      CALL PL (EX(X(1)),WY(Y(1)),.FALSE.)
      DO 10 I=2,N
   10 CALL PL (EX(X(I)),WY(Y(I)),.TRUE.)
      RETURN
      END
      SUBROUTINE  PLTHP (NR,NT,PL)

C     [HYPERBOLIC POLAR]
C     POLAR COORDINATE GRID WITH RADIAL HYPERBOLIC TANGENT DISTORTION.
C     NR   NUMBER OF CIRCLES OF CONSTANT RADIUS
C     NT   NUMBER OF RADII OF CONSTANT ANGLE
C     PL   PEN MOVEMENT SUBROUTINE, NORMALLY PLTCA
C     [17-MAY-75]

      EXTERNAL    PL
      DATA        HX,HY/4.50,3.25/
      DATA        RR,UU/3.25,3.00/

      CALL PLTBO
      CALL PLTIG (-HX,-HY,1,PL)
      CALL PLTIG ( HX, HY,2,PL)
      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),PL)
   11 CALL PLTCI (0.0,0.0,RR,PL)

      T=0.0
      SS=(0.25*RR)/UU
      DO 20 I=1,NT,2
      CALL PLTIG (RR*COS(T),RR*SIN(T),3,PL)
      CALL PLTIG (SS*COS(T),SS*SIN(T),4,PL)
      T=T+DT
      CALL PLTIG (SS*COS(T),SS*SIN(T),3,PL)
      CALL PLTIG (RR*COS(T),RR*SIN(T),4,PL)
   20 T=T+DT
      RETURN
      END
      SUBROUTINE  PLTIG (X,Y,L,PL)

C     [INCREMENTAL GRAPH]
C     PLOT A GRAPH POINT BY POINT. THE ORIGIN OF THE GRAPH IS THE LOWER
C     LEFT HAND CORNER OF A 1.0 X 1.0 SQUARE.  THE ACTUAL POSITION ON A
C     LETTER SIZED PLOTTER SHEET MUST BE CALCULATED BY THE PEN MOVEMENT
C     SUBROUTINE PL. THE 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     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     [06-JUN-75]
      EXTERNAL    PL
      DATA        N,DT/101,0.01/
      EX(X)=(X-X1)*SX
      WY(Y)=(Y-Y1)*SY
      GO TO (10,20,5,40,5,5,5),L
    5 SX=1.0/(X2-X1)
      SY=1.0/(Y2-Y1)
      GO TO (7,7,30,7,50,60,70),L
    6 X0=X
      Y0=Y
    7 RETURN
   10 X1=X
      Y1=Y
      GO TO 7
   20 X2=X
      Y2=Y
      GO TO 7
   30 CALL PL (EX(X),WY(Y),.FALSE.)
      GO TO 6
   40 CALL PL (EX(X),WY(Y),.TRUE.)
      GO TO 6
   50 TE=0.0
      CALL PL (TE,WY(Y),.FALSE.)
      DO 51 I=1,N
      CALL PL (TE,WY(Y),.TRUE.)
   51 TE=TE+DT
      TE=0.0
      CALL PL (EX(X),TE,.FALSE.)
      DO 52 I=1,N
      CALL PL (EX(X),TE,.TRUE.)
   52 TE=TE+DT
      GO TO 7
   60 CALL PLTFM (EX(X),WY(Y),0.010,PL)
      CALL PL (EX(X0),WY(Y0),.FALSE.)
      GO TO 7
   70 CALL PLTFM (EX(X),WY(Y),0.015,PL)
      CALL PL (EX(X0),WY(Y0),.FALSE.)
      GO TO 7
      END
C     ------------------------------------------------------------------


      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


C     ------------------------------------------------------------------


      SUBROUTINE  PLTIV (Z1,ZE,Z2,NX,NY,RO,TI,PL)

C     [INCLINED VIEW]
C     THE SURFACE MAY BE TILTED IN THE DIRECTION OF THE OBSERVER AND
C     THEN ROTATED ABOUT A VERTICAL AXIS BEFORE GENERATING A HIDDEN
C     LINE VIEW.  TILT IS ZERO WHEN SEEN DIRECTLY OVERHEAD, 90 DEGREES
C     WHEN SEEN DIRECTLY FROM THE GROUND.  POSITIVE TILT IS TOWARD THE
C     OBSERVER, NEGATIVE TILT AWAY FROM HIM.  THE ANGLE OF ROTATION IS
C     ZERO WHEN THE Y-AXIS RUNS DIRECTLY AWAY FROM THE OBSERVER, AND
C     IS POSITIVE WHEN THE POSITIVE X-AXIS MOVES TOWARD HIM.
C     ZE(NX,NY) ARRAY OF FUNCTION VALUES
C     Z1,Z2     RANGE OF FUNCTION VALUES
C     RO,TI     ANGLES OF ROTATION, TILT (DEGREES)
C     PL        PEN MOVEMENT SUBROUTINE, USUALLY PLTCA
C     [30-MAY-75]

      EXTERNAL    PL
      DIMENSION   ZE(1),O(3,3)

      CALL PLTEU (O,RO,TI,0.0)
      CALL VISNH
      CALL VISIS (Z1,ZE,Z2,1,NX,NX,1,NY,NY,O,PL)
      RETURN
      END


C     ------------------------------------------------------------------
      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

