Web pdp-10.trailing-edge.com

Trailing-Edge - PDP-10 Archives - decuslib20-03 - decus/20-0082/calcom.src
There is 1 other file named calcom.src in the archive. Click here to see a list.
```      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     [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

```