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


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

      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
