FUNCTION CARG (Z) 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 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 ORTHOGRAPHIC RELIEF. C Z0 CONTOUR LEVEL SOUGHT C (XF,YF) LIGHTING DIRECTION FOR ORTHOGRAPHIC RELIEF C (IA,IB) X-INTERVAL TO BE CONTOURED C (JA,JB) Y-INTERVAL TO BE CONTOURED C ZE(NX,NY) ARRAY OF FUNCTION VALUES C PL PEN MOVEMENT SUBROUTINE C [06-JAN-75] LOGICAL*1 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*1 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 , 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*1 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*1 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*1 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*1 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*1 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*1 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*1 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 WHOSE 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*1 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*1 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*1 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*1 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*1 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*1 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*1 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*1 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*1 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 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 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*1 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*1 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 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. 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*1 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 VISIBLE. 60 CALL PL (EX,WY,.FALSE.) 61 IF ((J.LT.1).OR.(J.GT.N)) GO TO 70 CALL PL (X(J),Y(J),.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 N=L*M MM=1 EL=FLOAT(L) EM=FLOAT(M) TE=0.5*(EL+1.0) ZS=(1.0-VS)/(Z2-Z1) DZ=(Z2-Z1)/FLOAT(NZ-1) DUI=-(EL*US)/FLOAT(MY-1) DUJ=(1.0-US)/FLOAT(MX-1) DVI=VS/FLOAT(MY-1) I0=(NY+1-M*(NY-3))/2 J0=(NX+1-L*(NX-1))/2 K0=((MK+1)*(1-N))/2 10 K=K0 I=MAX0(MIN0(I0,NY+1),0) J=J0 EU=TE*US+DUI*FLOAT(I-1)+DUJ*FLOAT(J-1) VE= DVI*FLOAT(I-1) 20 IF ((I.LT.1).OR.(I.GT.NY)) GO TO 22 K=MAX0(MIN0(K+N,MK),1) U(K)=EU V(K)=VE+SC(ZE(IX(J,I))) 22 I=I-M EU=EU-DUI VE=VE-EM*DVI IF ((I.LT.1).OR.(I.GT.NY)) GO TO 30 K=MAX0(MIN0(K+N,MK),1) U(K)=EU V(K)=VE+SC(ZE(IX(J,I))) J=J+L EU=EU+EL*DUJ IF ((J.GE.1).AND.(J.LE.NX)) GO TO 20 30 ZI=Z1 DO 60 IZ=1,NZ K=K0 I=MAX0(MIN0(I0,NY+1),0) J=J0 EU=TE*US+DUI*FLOAT(I-1)+DUJ*FLOAT(J-1) VE= DVI*FLOAT(I-1) 40 IF ((I.LT.1).OR.(I.GT.NY)) GO TO 42 K=MAX0(MIN0(K+N,MK),1) W(K)=VE+SC(ZI) 42 I=I-M EU=EU-DUI VE=VE-EM*DVI IF ((I.LT.1).OR.(I.GT.NY)) GO TO 50 K=MAX0(MIN0(K+N,MK),1) W(K)=VE+SC(ZI) J=J+L EU=EU+EL*DUJ IF ((J.GE.1).AND.(J.LE.NX)) GO TO 40 50 CALL VISRB (A,B,NA,MK,U,V,K,U,W,K,-1.0) CALL VISHH (D,E,G,ND,A,B,NA,MM,PL) MM=-MM 60 ZI=ZI+DZ I0=I0+M IF ((I0.GE.0).AND.(I0.LE.NY+1)) GO TO 10 J0=J0+L IF ((J0.GE.1).AND.(J0.LE.NX)) GO TO 10 RETURN END SUBROUTINE VISDO (Z1,S1,S2,Z2,NX,MX,NY,MY,US,VS,L,IS,PL) C [DOUBLE SURFACE] C S1,S2 ARRAYS CONTAINING THE TWO SURFACES C Z1,Z2 SPAN OF SURFACE VALUES C MX,MY COMMON DIMENSION OF THE ARRAYS S1 AND S2 C NX,NY SECTIONS OF S1 AND S2 ACTUALLY USED C US,VS TOTAL SHEARS IN U AND V DIRECTIONS C L DIRECTION OF VIEW (1=WEST, -1=EAST) C IS SEPARATION OPTION (1=YES , -1=NO) C PL PEN MOVEMENT SUBROUTINE C [15-MAY-74] EXTERNAL PL LOGICAL*1 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=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,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) ARRAY OF FUNCTION VALUES C I DIRECTION OF PEN MOVEMENT C PL PEN MOVEMENT SUBROUTINE 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*1 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 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 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*1 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