      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,S,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     S         S=1, GRAPH FU POSITIVE; S=-1, GRAPH FU NEGATIVE
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     [15-MAY-75]

      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,S,O,PR,PL)
      CALL VISSS (FU,J1,J2,NP,I1,I2,NT,1, 1,S2,B,S,O,PR,PL)
      CALL VISSS (FU,J1,J2,NP,I3,I4,NT,1, 1,S2,B,S,O,PR,PL)
      CALL VISSS (FU,J1,J2,NP,I3,I4,NT,1,-1,S1,B,S,O,PR,PL)
      CALL VISNH
      CALL VISSS (FU,J3,J4,NP,I1,I2,NT,-1,-1,S1,C,S,O,PR,PL)
      CALL VISSS (FU,J3,J4,NP,I1,I2,NT,-1, 1,S2,C,S,O,PR,PL)
      CALL VISSS (FU,J3,J4,NP,I3,I4,NT,-1, 1,S2,C,S,O,PR,PL)
      CALL VISSS (FU,J3,J4,NP,I3,I4,NT,-1,-1,S1,C,S,O,PR,PL)
      RETURN
      END
      SUBROUTINE  PLTSW (Z1,ZE,Z2,NX,NY,PL)

C     [SOUTHWEST VIEW]
C     ZE(NX,NY) ARRAY OF FUNCTION VALUES
C     Z1,Z2     RANGE OF FUNCTION VALUES
C     PL        PEN MOVEMENT SUBROUTINE
C     [10-MAY-75]

      EXTERNAL    PL
      DIMENSION   ZE(1)

      CALL VISNH
      CALL VISDS (Z1,ZE,Z2,1,NX,NX,1,NY,NY,0.2,0.2,1,1,PL)
      RETURN
      END

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

      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

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

      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


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


      SUBROUTINE  PLTTR (X,Y,P)

C     [TRIANGULAR ]
C     TAKING X,Y AS TWO OF THREE HOMOGENEOUS COORDINATES WITH
C     X+Y+Z=1, CALCULATE THE PLANAR EQUIVALENT FOR GRAPHING PURPOSES
C     [07-JUN-75]

      LOGICAL     P

      CALL PLTCA (0.500*(1.0-X+Y),0.866*(1.0-X-Y),(P.AND.(X+Y.LE.1.0)))
      RETURN
      END


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


      SUBROUTINE  PLTTV (Z1,ZE,Z2,N,M,PL)

C     [TRIANGULAR VIEW]
C     PROGRAM TO PRODUCE A PERSPECTIVE DRAWING OF A FUNCTION DEFINED
C     OVER THE UNIT TRIANGLE IN TERMS OF HOMOGENEOUS COORDINATES. THE
C     FUNCTION VALUES OCCUPY THE LOWER 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.  PL IS A PEN MOVEMENT SUBROUTINE, PERHAPS PLTCA.
C     [11-MAY-75]

      EXTERNAL    PL
      DIMENSION   ZE(1)

      CALL VISNH
      CALL VISTS (Z1,ZE,Z2,N,M,PL)
      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  PVIDS (Z1,ZE,Z2,J1,J2,MX,I1,I2,MY,US,VS,L,M,S,PL)

C     [DIAGONAL SEQUENCE]
C     ZE(MX,MY) ARRAY OF FUNCTION VALUES
C     (Z1,Z2)   SPAN OF Z VALUES
C     J1,J2     RANGE OF X (HORIZONTAL COORDINATE)
C     I1,I2     RANGE OF Y (DEPTH COORDINATE)
C     US,VS     TOTAL SHEARS IN U AND V DIRECTIONS (MAXIMUM=1.0)
C     L         DIRECTION OF VIEW AND INCREMENT (+:WEST,  -:EAST)
C     M         DIRECTION OF VIEW AND INCREMENT (+:SOUTH, -:NORTH)
C     S         =1.0, PLOT POSITIVE PART; =-1.0, PLOT NEGATIVE PART
C     PL        PEN MOVEMENT SUBROUTINE, USUALLY PLTCA
C     [19-MAY-75]

      EXTERNAL    PL
      LOGICAL     P(501)
      DIMENSION   ZE(1),U(501),V(501)
      DATA        MK/501/
      NL=ISIGN(1,L)
      NM=ISIGN(1,M)
      N=NL*NM
      IL=I1-IABS(M)
      IU=I2+IABS(M)
      MM=M*MX
      NN=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
      IX=J+(I-1)*MX
      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+ZS*(ABS(ZE(IX))-Z1)
      P(K)=(S*ZE(IX).GE.0.0)
   22 I=I-M
      IX=IX-MM
      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+ZS*(ABS(ZE(IX))-Z1)
      P(K)=(S*ZE(IX).GE.0.0)
      J=J+L
      IX=IX+L
      EU=EU+EUJ
      IF ((J.GE.J1).AND.(J.LE.J2)) GO TO 20
   30 KK=(K+1-N*(K-1))/2
      CALL VISCH (U(KK),V(KK),P(KK),(MK+1-N*(MK-2*K+1))/2,NN,PL)
      NN=-NN
      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


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


      SUBROUTINE  PVIIS (Z1,ZE,Z2,J1,J2,MX,I1,I2,MY,O,S,PL)

C     [INCLINED SEQUENCE]
C     ZE(MX,MY) ARRAY OF FUNCTION VALUES
C     (Z1,Z2)   SPAN OF Z VALUES
C     J1,J2     RANGE OF X (HORIZONTAL COORDINATE)
C     I1,I2     RANGE OF Y (DEPTH COORDINATE)
C     O(3,3)    ORTHOGONAL MATRIX DEFINING INCLINATION
C     S         =1.0, PLOT POSITIVE; =-1.0, PLOT NEGATIVE
C     PL        PEN MOVEMENT SUBROUTINE, USUALLY PLTCA
C     [31-MAY-75]

      EXTERNAL    PL
      LOGICAL     P(501)
      DIMENSION   O(3,3),ZE(1),U(501),V(501)
      DATA        MK/501/

      PR(I)=0.667*(O(1,I)*EU+O(2,I)*VE+O(3,I)*UU)+0.5

      L=-IFIX(SIGN(1.0,O(2,1)))
      M= IFIX(SIGN(1.0,O(1,1)))
      NL=ISIGN(1,L)
      NM=ISIGN(1,M)
      N=NL*NM
      IL=I1-IABS(M)
      IU=I2+IABS(M)
      MM=M*MX
      NN=1
      TE=0.5*FLOAT(NL+1)
      ZS=1.0/(Z2-Z1)
      DUJ=1.0/FLOAT(MX-1)
      EUJ=DUJ*FLOAT(L)
      DVI=1.0/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
      IX=J+(I-1)*MX
      EU=DUJ*FLOAT(J-1)-0.5
      VE=DVI*FLOAT(I-1)-0.5
   20 IF ((I.LT.I1).OR.(I.GT.I2)) GO TO 22
      K=MAX0(MIN0(K+N,MK),1)
      UU=ZS*(ABS(ZE(IX))-Z1)
      U(K)=PR(1)
      V(K)=PR(2)
      P(K)=(S*ZE(IX).GE.0.0)
   22 I=I-M
      IX=IX-MM
      VE=VE-EVI
      IF ((I.LT.I1).OR.(I.GT.I2)) GO TO 30
      K=MAX0(MIN0(K+N,MK),1)
      UU=ZS*(ABS(ZE(IX))-Z1)
      U(K)=PR(1)
      V(K)=PR(2)
      P(K)=(S*ZE(IX).GE.0.0)
      J=J+L
      IX=IX+L
      EU=EU+EUJ
      IF ((J.GE.J1).AND.(J.LE.J2)) GO TO 20
   30 KK=(K+1-N*(K-1))/2
      CALL VISCH (U(KK),V(KK),P(KK),(MK+1-N*(MK-2*K+1))/2,NN,PL)
      NN=-NN
      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  PVIIV (Z1,ZE,Z2,NX,NY,RO,TI,S,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 AND TILT, IN DEGREES
C     S         =1.0, GRAPH POSITIVE; =-1.0, GRAPH NEGATIVE
C     PL        PEN MOVEMENT SUBROUTINE, USUALLY PLTCA
C     [31-MAY-75]

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

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


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


      SUBROUTINE  PVISE (Z1,ZE,Z2,NX,NY,S,PL)

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(NX,NY) ARRAY OF FUNCTION VALUES
C     Z1,Z2     RANGE OF FUNCTION VALUES
C      S        =1.0, PLOT POSITIVE PART; =-1.0, PLOT NEGATIVE PART
C     PL        PEN MOVEMENT SUBROUTINE, PERHAPS PLTCA
C     [19-MAY-75]

      EXTERNAL    PL
      DIMENSION   ZE(1)

      CALL VISNH
      CALL PVIDS (Z1,ZE,Z2,1,NX,NX,1,NY,NY,0.2,0.2,-1,1,S,PL)
      RETURN
      END
      SUBROUTINE  PVISW (Z1,ZE,Z2,NX,NY,S,PL)

C     [SOUTHWEST VIEW]
C     ZE(NX,NY) ARRAY OF FUNCTION VALUES
C     Z1,Z2     RANGE OF FUNCTION VALUES
C      S        =1.0, PLOT POSITIVE PART; =-1.0, PLOT NEGATIVE PART
C     PL        PEN MOVEMENT SUBROUTINE
C     [20-MAY-75]

      EXTERNAL    PL
      DIMENSION   ZE(1)

      CALL VISNH
      CALL PVIDS (Z1,ZE,Z2,1,NX,NX,1,NY,NY,0.2,0.2,1,1,S,PL)
      RETURN
      END


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


      SUBROUTINE  PVITV (Z1,ZE,Z2,N,M,S,PL)

C     [TRIANGULAR VIEW]
C     PROGRAM TO PRODUCE A PERSPECTIVE DRAWING OF A FUNCTION DEFINED
C     OVER 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.  PL IS A PEN MOVEMENT SUBROUTINE, PERHAPS PLTCA.
C     [18-MAY-75]

      EXTERNAL    PL
      DIMENSION   ZE(1)

      CALL VISNH
      CALL PVITS (Z1,ZE,Z2,N,M,S,PL)
      RETURN
      END
      SUBROUTINE  PVITS (Z1,ZE,Z2,N,M,S,PL)

C     [TRIANGULAR SEQUENCE]
C     ZE(M,M)  ARRAY OF VALUES
C     Z1,Z2    RANGE OF VALUES
C     S        S=1.0, PLOT POSITIVE; S=-1.0, PLOT NEGATIVE
C     PL       PEN MOVEMENT SUBROUTINE, PERHAPS PLTCA
C     [18-MAY-75]

      EXTERNAL    PL
      LOGICAL     P(501)
      DIMENSION   ZE(1),U(501),V(501)

      IX(J,I)=(I-1)*M+J
      SC(Z)=ZS*(Z-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)
      EF=ZE(IX(J,I))
      U(K)=EU
      V(K)=VE+SC(ABS(EF))
      P(K)=S*EF.GE.0.0
   10 EU=EU+DU
      CALL VISCH (U,V,P,K,1,PL)
      K=0
      EU=HU*FLOAT(I-1)
      DO 20 J=I,L
      K=MIN0(K+1,MK)
      EF=ZE(IX(J,I))
      U(K)=EU
      V(K)=VE+SC(ABS(EF))
      P(K)=S*EF.GE.0.0
      K=MIN0(K+1,MK)
      EF=ZE(IX(J+1,I+1))
      U(K)=EU+HU
      V(K)=VE+DV+SC(ABS(EF))
      P(K)=S*EF.GE.0.0
   20 EU=EU+DU
      K=MIN0(K+1,MK)
      EF=ZE(IX(N,I))
      U(K)=EU
      V(K)=VE+SC(ABS(EF))
      P(K)=S*EF.GE.0.0
      CALL VISCH (U,V,P,K,-1,PL)
   30 VE=VE+DV
      RETURN
      END
      SUBROUTINE  VISBO (X1,T1,B1,M,X0,T0,B0,N0,X,Y,P,N,I,PL)

C     [BOUNDS]
C     X(N)   ARRAY OF ARGUMENTS
C     Y(N)   ARRAY OF FUNCTION VALUES
C     P(N)   ARRAY OF VISIBILITY FLAGS
C     I      DIRECTION OF PEN MOVEMENT (1=FORWARD, -1=BACKWARD)
C     PL     PEN MOVEMENT SUBROUTINE
C     [10-MAY-75]

      LOGICAL     L,P(1),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),P(J))
      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)).AND.P(J)
      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 VISIBLE.

   60 CALL PL (EX,WY,.FALSE.)
   61 IF ((J.LT.1).OR.(J.GT.N)) GO TO 70
      CALL PL (X(J),Y(J),P(J))
      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  VISCH (X,Y,P,N,I,PL)

C     [COLORED HORIZON]
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)   ARRAY OF FUNCTION VALUES
C     I      DIRECTION OF PEN MOVEMENT
C     PL     PEN MOVEMENT SUBROUTINE
C     [10-MAY-75]

      EXTERNAL    PL
      LOGICAL     P(1)
      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,P,N,I,PL)
      RETURN
      END

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

      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
      N=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


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


      SUBROUTINE  VISDS (Z1,ZE,Z2,J1,J2,MX,I1,I2,MY,US,VS,L,M,PL)

C     [DIAGONAL SEQUENCE]
C     ZE(MX,MY) ARRAY OF FUNCTION VALUES
C     (Z1,Z2)   SPAN OF Z VALUES
C     J1,J2     RANGE OF X (HORIZONTAL COORDINATE)
C     I1,I2     RANGE OF Y (DEPTH COORDINATE)
C     US,VS     TOTAL SHEARS IN U AND V DIRECTIONS (MAXIMUM=1.0)
C     L         DIRECTION OF VIEW AND INCREMENT (+:WEST,  -:EAST)
C     M         DIRECTION OF VIEW AND INCREMENT (+:SOUTH, -:NORTH)
C     PL        PEN MOVEMENT SUBROUTINE, USUALLY PLTCA
C     [10-MAY-75]

      EXTERNAL    PL
      DIMENSION   ZE(1),U(501),V(501)
      DATA        MK/501/

      NL=ISIGN(1,L)
      NM=ISIGN(1,M)
      N=NL*NM
      IL=I1-IABS(M)
      IU=I2+IABS(M)
      MM=M*MX
      NN=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
      IX=J+(I-1)*MX
      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+ZS*(ZE(IX)-Z1)
   22 I=I-M
      IX=IX-MM
      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+ZS*(ZE(IX)-Z1)
      J=J+L
      IX=IX+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,NN,PL)
      NN=-NN
      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

