      SUBROUTINE  VISES (Z1,ZE,Z2,X1,X2,NX,E1,E2,NE,L,M,PL)

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

      EXTERNAL    PL
      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=J0
      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,PL)
      IF (N.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.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     [10-MAY-75]

      EXTERNAL    PL
      LOGICAL     P(701)
      DIMENSION   X(1),Y(1),X0(1),T0(1),B0(1)
      DIMENSION   X1(701),T1(701),B1(701)
      DATA        M,P/701,701*T/

      CALL VISBO (X1,T1,B1,M,X0,T0,B0,N0,X,Y,P,N,I,PL)
      RETURN
      END


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


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

      EXTERNAL    PL
      LOGICAL     P(701)
      DIMENSION   X(1),Y(1)
      DIMENSION   X0(701),T0(701),B0(701)
      DIMENSION   X1(701),T1(701),B1(701)
      COMMON/VIS/ N0
      DATA        M,P/701,701*T/

      CALL VISBO (X1,T1,B1,M,X0,T0,B0,N0,X,Y,P,N,I,PL)
      RETURN
      END
      SUBROUTINE  VISIS (Z1,ZE,Z2,J1,J2,MX,I1,I2,MY,O,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     PL        PEN MOVEMENT SUBROUTINE, USUALLY PLTCA
C     [30-MAY-75]

      EXTERNAL    PL
      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*(ZE(IX)-Z1)
      U(K)=PR(1)
      V(K)=PR(2)
   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*(ZE(IX)-Z1)
      U(K)=PR(1)
      V(K)=PR(2)
      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


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


      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


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


      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,PL)

C     [POLAR SEQUENCE]
C     L= 1  WESTERN VIEW,  L=-1  EASTERN VIEW
C     M= 1  SOUTHERN VIEW, M=-1  NORTHERN VIEW
C     [17-MAY-75]

      EXTERNAL    PL
      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,PL)
      IF (N.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.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


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


      SUBROUTINE  VISRS (Z1,ZE,Z2,NX,MX,NY,MY,TH,PL)

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


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


      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,Q,B,S,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     Q         ACCORDING TO THE ORIENTATION OF EACH INFINITESIMAL
C               RECTANGLE TO THE LINE OF SIGHT AND SEARCH DIRECTION ITS
C               EDGES WILL FORM PART OF THE HORIZON OR NOT
C     B=        (ATAN2 CUT LINE DOES NOT FALL IN QUADRANT BEING GRAPHED)
C     S         S=1.0, GRAPH FU POSITIVE; S=-1.0, GRAPH FU NEGATIVE
C     O(3,3)    ORTHOGONAL ROTATION MATRIX
C     PR        PROJECTION SUBROUTINE
C     PL        PEN MOVEMENT SUBROUTINE
C     [15-MAY-75]

      EXTERNAL    PL
      LOGICAL     B,PE(701)
      DIMENSION   FU(1),O(3,3)
      DIMENSION   AZ(701),RA(701),OR(701)
      DATA        MK/701/

      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)=Q*SIGN(1.0,1.57079-T)*(TAN(T)-TN*COS(P-PH))

      NL=ISIGN(1,L)
      NM=ISIGN(1,M)
      N=NL*NM
      IL=I1-IABS(M)
      IU=I2+IABS(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)+NM*(I1-I2))/2+M
      J0=((J1+J2)+NL*(J1-J2))/2
      K0=((MK+1)*(1-N))/2
   10 K=K0
      I=MAX0(MIN0(I0,IU),IL)
      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)
      FV=FU(IX(J,I))
      CALL PR (RA(K),AZ(K),ABS(FV),T,P,O)
      PE(K)=S*FV.GT.0.0
      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)
      FV=FU(IX(J,I))
      CALL PR (RA(K),AZ(K),ABS(FV),T,P,O)
      PE(K)=S*FV.GT.0.0
      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)
      PE(MM)=PE(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=PE(MM)
      T4=OR(MM)
      AZ(MM)=AZ(K-MM+1)
      RA(MM)=RA(K-MM+1)
      PE(MM)=OR(K-MM+1)
      OR(MM)=OR(K-MM+1)
      AZ(K-MM+1)=T1
      RA(K-MM+1)=T2
      PE(K-MM+1)=T3
   37 OR(K-MM+1)=T4
   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 VISCH (AZ(K1),RA(K1),PE(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,PL)

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

      EXTERNAL    PL
      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,PL)
      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,PL)
   30 VE=VE+DV
      RETURN
      END

