Trailing-Edge
-
PDP-10 Archives
-
decuslib10-09
-
43,50466/csmp3.f4
There are no other files named csmp3.f4 in the archive.
SUBROUTINE CSM8(V,H,IOFSET,NSTP)
C PRINT CONTROLLER
INTEGER OU
REAL IPLOT(0/100)
DIMENSION V(1),H(1),KI(4)
DIMENSION C(76)
DIMENSION SYM(4)
COMMON REALS(395),INTS(547)
EQUIVALENCE (INTS(386),KEY7),(INTS(388),KEY9)
EQUIVALENCE (INTS(535),K1),(INTS(536),K2)
EQUIVALENCE (INTS(537),K3),(INTS(538),K4)
EQUIVALENCE (INTS(539),NK)
EQUIVALENCE (REALS(2),C(1)),(REALS(77),T)
EQUIVALENCE (REALS(392),VDEL),(REALS(394),VMIN)
EQUIVALENCE (KI(1),K1)
DATA FBLANK,FDASH,FI,FPLUS/1H ,1H-,1HI,1H+/
DATA SYM/'A','B','C','D'/
COMMON /ODEVIM/OU
GO TO (160,9,140,9),KEY7+2
9 NMAX=100
NZ=.5+50.*(-VMIN)/VDEL
IF(OU.NE.5.AND.OU.NE.8.AND.OU.NE.19)GO TO 70
NMAX=0
DO 8 K=0,50
8 IPLOT(K)=FBLANK
DO 40 J=1,NK
N=.5+50.*(C(KI(J))-VMIN)/VDEL
IF(NZ.GE.0.AND.NZ.LE.50)IPLOT(NZ)=FI
IF(N.GE.0.AND.N.LE.50)IPLOT(N)=SYM(J)
40 NMAX=MAX0(NMAX,MIN0(N,50),MIN0(NZ,50))
GO TO 100
70 NZ=.5+100.*(-VMIN)/VDEL
DO 92 K=0,100
92 IPLOT(K)=FBLANK
DO 90 J=1,NK
N=.5+100.*(C(KI(J))-VMIN)/VDEL
IF(NZ.GE.0.AND.NZ.LE.100.)IPLOT(NZ)=FI
90 IF(N.GE.0.AND.N.LE.100.)IPLOT(N)=SYM(J)
100 WRITE(OU,130)T,(IPLOT(I),I=0,NMAX)
130 FORMAT(1H ,G9.3,1X,101A1)
RETURN
C PRINT ONLY
140 IPLOT(1)=C(K1)
IPLOT(2)=C(K2)
IPLOT(3)=C(K3)
IPLOT(4)=C(K4)
WRITE(OU,150) T,(IPLOT(I),I=1,NK)
150 FORMAT(1H ,G11.5,4G13.6)
RETURN
C GRAPHING
160 DO 162 I=1,NK
ISS=NSTP+(I-1)*IOFSET
V(NSTP+(I-1)*IOFSET)=C(KI(I))
162 H(NSTP+(I-1)*IOFSET)=T
END
SUBROUTINE CSM8A(IOFSET)
C PRINT SPECIFICATIONS
INTEGER PRINT(4),TEST8,OU
COMMON REALS(395),INTS(547)
EQUIVALENCE (INTS(384),KEY5),(INTS(385),KEY6)
EQUIVALENCE (INTS(386),KEY7),(INTS(388),KEY9)
EQUIVALENCE (INTS(532),TEST8),(INTS(535),PRINT(1)),(INTS(539),NK)
EQUIVALENCE (REALS(78),DT),(REALS(391),TSAMP)
EQUIVALENCE (REALS(392),VDEL),(REALS(394),VMIN)
EQUIVALENCE (REALS(80),TTOT)
C
COMMON /ODEVIM/OU
COMMON /NOSTOP/ITHROU
C
IF (ITHROU.EQ.1)GO TO 270
WRITE(30,10)
10 FORMAT(/10X,14HOUTPUT CONTROL/)
C PRINT INTERVAL SPECIFICATION
20 WRITE(30,30)
30 FORMAT(16H PRINT INTERVAL=$)
TSAMP=FINPUT(0,IERR)
IF (IERR.NE.0) GO TO 20
IF (FINPUT(-1,IERR).NE.0.0) GO TO 20
IOFSET=INT(TTOT/TSAMP)+1
IF (TSAMP.GE.DT) GO TO 60
WRITE(30,50)
50 FORMAT(56H PRINT INTERVAL CANNOT BE LESS THAN INTEGRATION INTERVAL
1)
GO TO 20
C
C PRINT VARIABLES SPECIFICATION
60 CONTINUE
IF(KEY7.EQ.1)GO TO 160
C PRINT AND PLOT AND GRAPH
100 FORMAT(18H NOT A LEGAL BLOCK/)
110 WRITE(30,120)
120 FORMAT(25H MINIMUM, MAXIMUM VALUES=$)
VMIN=FINPUT(0,IERR)
IF (IERR.NE.0) GO TO 110
VMAX=FINPUT(0,IERR)
IF (IERR.NE.0) GO TO 110
IF (FINPUT(-1,IERR).NE.0.0) GO TO 110
VDEL=VMAX-VMIN
IF (VDEL.GT.0.0) GO TO 160
C LET GRAPH PROG CALC MAX,MIN(OPTIONAL)
IF(VDEL.EQ.0.AND.KEY7.LT.0)GO TO 160
WRITE(30,150)
150 FORMAT(36H MAXIMUM CANNOT BE LESS THAN MINIMUM/)
GO TO 110
160 WRITE(30,170)
170 FORMAT(36H BLOCK A, BLOCK B, BLOCK C, BLOCK D=$)
NK=0
DO 190 I=1,4
PRINT(I)=KINPUT(0,IERR)
IF (IERR) 200,180,160
180 IF (PRINT(I).LT.1.OR.PRINT(I).GT.75) GO TO 220
NK=NK+1
190 CONTINUE
IF (FINPUT(-1,IERR).NE.0.0) GO TO 160
200 IF (NK) 160,160,320
220 WRITE(30,100)
GO TO 160
C PRINT ONLY HEADING
250 IF (TEST8.EQ.2.AND.KEY7.NE.1) GO TO 270
BLOCK='BLOCK'
WRITE(OU,260) (BLOCK,PRINT(I),I=1,NK)
260 FORMAT(1H1/5X,5HTIME ,4(5X,A5,I3))
GO TO 320
C
270 GO TO (320,280,250,280),KEY7+2
C PRINT AND PLOT HEADING
280 IF (OU.NE.5.AND.OU.NE.8.AND.OU.NE.19) GO TO 300
WRITE(30,290)VMIN,VMAX
290 FORMAT(1H1,4HTIME,2X,G13.6,36X,G13.6/,11X,'+',10('----+'))
GO TO 320
300 WRITE(OU,310)VMIN,VMAX
310 FORMAT(1H1,1X,4HTIME,2X,G13.6,86X,G13.6,/,11X,'+',20('----+'))
C TEST8=1 UNTIL FIRST TIME THROUGH CSM8A
C TEST8=2 AFTER FIRST TIME THROUGH CSM8A
320 TEST8=2
RETURN
END
SUBROUTINE CSM10(IOFSET)
C ALLOCATOR FOR GRAPHING(DEALLOCATES OTHERWISE)
DIMENSION S(1)
COMMON REALS(395),INTS(547)
EQUIVALENCE (INTS(386),KEY7),(INTS(539),NK),(INTS(529),TEST5)
IF(KEY7.LT.0)GO TO 5
CALL ALLCOR(1,IERR,I1,S)
I2=I1
GO TO 10
5 MAX=NK*IOFSET*2
CALL ALLCOR(MAX,IERR,I1,S)
I2=I1+IOFSET*NK
IF(IERR.GE.0)GO TO 10
WRITE(30,8)
8 FORMAT(' NOT ENOUGH CORE TO GRAPH, TRY INCREASING THE PRINT',/,
1 ' INTERVAL OR DECREASING THE TOTAL TIME',/)
TEST5=4
RETURN
10 CONTINUE
CALL CSM10A(S(I1),S(I2),IOFSET)
RETURN
END
SUBROUTINE CSM10A(V,H,IOFSET)
C CONTROLS THE COMPUTATION AND OUTPUT
DIMENSION V(1),H(1)
INTEGER TEST5,ORDER(76)
LOGICAL RSAC
DIMENSION INTG(25),C(76),PAR1(75),Y(25),DYDT(25),YK(25)
COMMON REALS(395),INTS(547)
EQUIVALENCE (INTS(388),KEY9)
EQUIVALENCE (INTS(396),INTG(1)),(INTS(449),ORDER(1))
EQUIVALENCE (INTS(529),TEST5),(INTS(534),NLIST)
EQUIVALENCE (INTS(542),NEQ),(INTS(546),IR)
EQUIVALENCE (REALS(2),C(1)),(REALS(77),T),(REALS(78),DT)
EQUIVALENCE (REALS(79),DTS2),(REALS(80),TTOT)
EQUIVALENCE (REALS(81),PAR1(1)),(REALS(341),Y(1))
EQUIVALENCE (REALS(366),DYDT(1)),(REALS(391),TSAMP)
EQUIVALENCE (INTS(539),NK),(REALS(392),VDEL),(REALS(394),VMIN)
EQUIVALENCE (INTS(386),KEY7)
INTEGER OU
COMMON /ODEVIM/OU
NP=IOFSET
C NORMAL SETUP
DO 10 NEXT=2,NLIST
I=ORDER(NEXT)
10 C(I)=PAR1(I)
T=0.0
TZERO=0.0
DO 20 NEXT=1,NEQ
I=INTG(NEXT)
20 Y(NEXT)=C(I)
IR=7243
EPSLN=DTS2/(TSAMP*2.0)
TEST5=1
N=1
NSTP=0
NN=T/TSAMP+1.0
CALL CSM11
C
VMN=VMIN
VMX=VDEL+VMN
HMN=0
HMX=TTOT
C
C START EXECUTION
30 IF (RSAC(0)) GO TO 110
40 GO TO (50,80,100,110,110,110),TEST5
50 NSTP=NSTP+1
CALL CSM8(V,H,IOFSET,NSTP)
C FIRST HALF-STEP
60 TEST5=2
DO 70 NEXT=1,NEQ
YK(NEXT)=Y(NEXT)
70 Y(NEXT)=Y(NEXT)+DTS2*DYDT(NEXT)
AXX=N
TNEXT=AXX*DT+TZERO
T=TNEXT-DTS2
CALL CSM11
GO TO 40
C SECOND HALF STEP
80 TEST5=3
DO 90 NEXT=1,NEQ
90 Y(NEXT)=YK(NEXT)+DT*DYDT(NEXT)
T=TNEXT
N=N+1
CALL CSM11
GO TO 30
C TIME TO PRINT
100 M=T/TSAMP+EPSLN
IF (M.LT.NN) GO TO 120
110 NSTP=NSTP+1
CALL CSM8(V,H,IOFSET,NSTP)
NN=M+1
C IS RUN FINISHED
120 IF (TEST5.GT.3) GO TO 150
130 IF (RSAC(0)) GO TO 140
IF (T-TTOT+DTS2) 60,150,150
140 TEST5=5
150 IF(OU.NE.5.AND.KEY7.GE.0)WRITE(OU,160)
160 FORMAT(//1H1)
C
IF(KEY7.GE.0)GO TO 920
C GRAPHING
C VS - VERT. SCAL. PARAM.
VS=-1
C HS - HORZ. SCAL. PARAM.
HS=-1
C IV - Y AXIS PARAM.(LINEAR)
IV=0
C IH - X AXIS PARAM.(LINEAR)
IH=0
C IOPT - OPTION PARAM.
C IOPT1 - PLOT DATA
IOPT1=1
C IOPT2 - PLOT AXIS WITH LABELS
IOPT2=20
C IOPT3 - NO BOX
IOPT3=0
C IOPT4 - PLOT ALL DATA
IOPT4=0
C IOPT5 - FULL SCREEN PLOT
IOPT5=10000
IOPT=IOPT1+IOPT2+IOPT3+IOPT4+IOPT5
C IERASE - ERASE SCREEN BEFORE PLOTTING
IERASE=1
C ERR - PLOT MEDIUM
ERR=20
K=0
C CONTINUOUS PLOT
M=0
DO 200 I=1,NP*NK,NP
K=K+1
CALL GPLOT(V(I),H(I),NP,VS,HS,IV,IH,IOPT,M,ERR,IERASE,
1VMX,VMN,HMX,HMN)
IF(I.NE.1)GO TO 85
IERASE=0
IOPT=IOPT-20
C SET 5 OR 6 LABELS ON THE CURVE SO THEY DONT OVER LAP
85 L=0
DO 87 J=0,NP-K-1,NP/6
L=L+1
V(L)=V(I+J+K)
H(L)=H(I+J+K)
87 CONTINUE
CALL GPLOT(V,H,L,VS,HS,IV,IH,IOPT,K+1,ERR,IERASE,
1VMX,VMN,HMX,HMN)
200 CONTINUE
CALL ANMODE
CALL HOME
BNK=' '
TYPE 910,(BLK,INTS(535+KK-1),KK=1,NK)
910 FORMAT(' PRINT SYMBOLS:',/,A1,'"SQUARE" - BLOCK ',I2,A1,
1' "X" - BLOCK ',I2,A1,' "+" - BLOCK ',I2,A1,
1' "TRIANGLE" - BLOCK ',I2,/)
920 RETURN
END
SUBROUTINE CSM11
C DOES THE COMPUTATION REQUIRED
C TO EVALUATE THE DERIVATIVE VECTOR
C FOR ONE-HALF TIME STEP
INTEGER TEST5,ORDER(76)
LOGICAL RSAC
DIMENSION INTG(25),C(76),F(3,11),Y(25),DYDT(25)
DIMENSION MTRX1(75),MTRX2(75),MTRX3(75),MTRX4(75),MTRX5(75)
DIMENSION PAR1(75),PAR2(75),PAR3(75)
COMMON REALS(395),INTS(547)
EQUIVALENCE (INTS(1),MTRX1(1)),(INTS(76),MTRX2(1))
EQUIVALENCE (INTS(151),MTRX3(1)),(INTS(226),MTRX4(1))
EQUIVALENCE (INTS(301),MTRX5(1))
EQUIVALENCE (INTS(396),INTG(1)),(INTS(449),ORDER(1))
EQUIVALENCE (INTS(529),TEST5),(INTS(534),NLIST)
EQUIVALENCE (INTS(540),NCON),(INTS(542),NEQ),(INTS(546),IR)
EQUIVALENCE (REALS(2),C(1)),(REALS(78),DT),(REALS(79),DTS2)
EQUIVALENCE (REALS(81),PAR1(1)),(REALS(156),PAR2(1))
EQUIVALENCE (REALS(231),PAR3(1)),(REALS(306),F(1,1))
EQUIVALENCE (REALS(341),Y(1)),(REALS(366),DYDT(1))
C
DO 10 I=1,NEQ
J=INTG(I)
10 C(J)=Y(I)
NEXT=NCON
20 I=ORDER(NEXT)
P1=PAR1(I)
P2=PAR2(I)
P3=PAR3(I)
J=MTRX2(I)
K=MTRX3(I)
L=MTRX4(I)
IF (J.GE.0.AND.J.LE.76) CJ=C(J)
IF (K.GE.0.AND.K.LE.76) CK=C(K)
IF (L.GE.0.AND.L.LE.76) CL=C(L)
M=MTRX1(I)
C MODIFIED FOR BLOCKS A,C,E 25 APR 74.
IF (M.LE.10) GO TO (1750,30,2750,40,3750,80,110,120,130,140),M
M=M-10
IF (M.LE.10) GO TO (650,150,180,190,210,220,230,240,270,290),M
M=M-10
GO TO (340,350,360,370,380,390,410,510,520),M
C SPECIAL LOADABLE FUNCTIONS A,C,E ADDED 25 APR 74.
C 'A' (DEFAULT IN CSMP4.F4 IS LOG(CJ) BASE E)
1750 CI=BLOCKA(CJ,CK,CL,P1,P2,P3,$750)
GO TO 600
C 'C' (DEFAULT IN CSMP4.F4 IS COS(CJ+3.14159265*P1) )
C (CJ IN RADIANS)
2750 CI=BLOCKC(CJ,CK,CL,P1,P2,P3,$750)
GO TO 600
C 'E' (DEFAULT IN CSMP4.F4 IS EXP(CJ) )
3750 CI=BLOCKE(CJ,CK,CL,P1,P2,P3,$750)
GO TO 600
C B - BANG-BANG
30 CI=SIGN(1.0,CJ)
GO TO 600
C D - DEAD SPACE
40 IF (CJ) 50,200,60
50 DIFF=CJ-P2
IF (DIFF) 70,200,200
60 DIFF=CJ-P1
IF (DIFF) 200,200,70
70 CI=DIFF
GO TO 600
C F - FUNCTION GENERATOR
80 NF=MTRX5(I)
P3=P1-P2
IF (P3.LE.0.0) GO TO 750
P1=10.0*(CJ-P2)/P3
IF (P1.GT.0.0) GO TO 90
CI=F(NF,1)
GO TO 600
90 NSECT=P1
IF (NSECT.LT.10) GO TO 100
CI=F(NF,11)
GO TO 600
100 P2=NSECT
P3=P1-P2
P1=F(NF,NSECT+1)
P2=F(NF,NSECT+2)
CI=P1+P3*(P2-P1)
GO TO 600
C G - GAIN
110 CI=P1*CJ
GO TO 600
C H - HALF POWER (SQUARE ROOT)
120 IF (CJ.LT.0.0) GO TO 750
CI=SQRT(CJ)
GO TO 600
C I - INTEGRATOR (MAXIMUM 25 ELEMENTS)
130 M=MTRX5(I)
DYDT(M)=CJ+P2*CK+P3*CL
GO TO 650
C J - JITTER (RANDOM NUMBER GENERATOR BETWEEN + AND - 1)
140 IR=259*IR
CI=FLOAT(IR)/131072.0
GO TO 600
C K - CONSTANT
C L - LIMITER
150 IF (CJ.LT.P1) GO TO 160
CI=P1
GO TO 600
160 IF (CJ.GT.P2) GO TO 280
170 CI=P2
GO TO 600
C M - MAGNITUDE
180 CI=ABS(CJ)
GO TO 600
C N - NEGATIVE CLIPPER
190 IF (CJ.GT.0.0) GO TO 280
200 CI=0.0
GO TO 600
C O - OFFSET
210 CI=CJ+P1
GO TO 600
C P - POSITIVE CLIPPER
220 IF (CJ) 280,200,200
C Q - QUIT
230 IF (CJ-CK) 650,650,850
C R - RELAY
240 IF (CJ.LT.0.0) GO TO 260
250 CI=CK
GO TO 600
260 CI=CL
GO TO 600
C S - SWITCH
270 M=P1
IF (RSAC(M)) GO TO 250
280 CI=CJ
GO TO 600
C T -TIME PULSE GENERATOR
290 IF (TEST5-2) 300,200,330
300 MTRX5(I)=0
310 IF (CJ.LT.0.0) GO TO 200
MTRX5(I)=1
320 PAR2(I)=-P1+DTS2+DT
CI=1.0
GO TO 600
330 IF (MTRX5(I).EQ.0) GO TO 310
IF (P2.GE.0.0) GO TO 320
PAR2(I)=P2+DT
GO TO 200
C U - UNIT DELAY
340 IF (TEST5.NE.1) C(I)=P2
PAR2(I)=CJ
GO TO 650
C V - VACUOUS (USED IN CONJUNCTION WITH WYE ELEMENT)
350 IF (TEST5.EQ.1) MTRX5(I)=NEXT
GO TO 650
C W - WEIGHTED SUMMER
360 CI=CJ*P1+CK*P2+CL*P3
GO TO 600
C X - MULTIPLIER
370 CI=CJ*CK
GO TO 600
C Y - WYE(USED IN CONJUNCTION WITH VACUOUS ELEMENT)
380 IF (ABS(1.0-CK/CJ).LE.P1) GO TO 280
IF (RSAC(0)) GO TO 800
C(K)=(1.0-P2)*CJ+P2*CK
NEXT=MTRX5(K)
GO TO 20
C Z - ZERO ORDER HOLD
390 IF (TEST5.NE.1) GO TO 400
PAR2(I)=C(I)
P2=C(I)
400 IF (CK.LE.0.0) GO TO 170
PAR2(I)=CJ
GO TO 280
C + - SUMMER
410 IF (J) 420,430,440
420 J=-J
CI=-C(J)
GO TO 450
430 CI=0.0
GO TO 450
440 CI=CJ
450 IF (K) 460,480,470
460 K=-K
CI=CI-C(K)
GO TO 480
470 CI=CI+CK
480 IF (L) 490,600,500
490 L=-L
CI=CI-C(L)
GO TO 600
500 CI=CI+CL
GO TO 600
C - - SIGN INVERTER
510 CI=-CJ
GO TO 600
C / - DIVIDER
520 IF (CK.EQ.0.0) GO TO 750
CI=CJ/CK
C 1 - SPECIAL ELEMENT NUMBER 1
C 2 - SPECIAL ELEMENT NUMBER 2
C 3 - SPECIAL ELEMENT NUMBER 3
C 4 - SPECIAL ELEMENT NUMBER 4
C 5 - SPECIAL ELEMENT NUMBER 5
C HAVE ALL BEEN DELETED
600 C(I)=CI
650 IF (NEXT-NLIST) 700,900,750
700 NEXT=NEXT+1
GO TO 20
C PROCESSING ERROR
750 TEST5=4
RETURN
C RUN TERMINATED BY SWITCH 0
800 TEST5=5
RETURN
C RUN TERMINATED BY QUIT ELEMENT
850 TEST5=6
900 RETURN
END
SUBROUTINE CSM13
C BLOCK OUTPUT INTERROGATION
DIMENSION C(75)
COMMON REALS(395),INTS(547)
EQUIVALENCE (REALS(2),C(1))
WRITE(30,10)
10 FORMAT (/10X,28H OUTPUT INTERROGATION OPTION/)
20 WRITE(30,30)
30 FORMAT(7H BLOCK=$)
I=KINPUT(0,IERR)
IF (IERR) 90,40,50
40 IF (FINPUT(-1,IERR).NE.0.0) GO TO 50
IF (I) 50,90,70
50 WRITE(30,60)
60 FORMAT(5H WHAT)
GO TO 20
70 IF (I.GT.75) GO TO 50
WRITE(30,80) I,C(I)
80 FORMAT(16H OUTPUT OF BLOCK,I3,4H IS ,G15.8)
GO TO 20
90 RETURN
END
C PLOT ROUTINE FOR H-P PLOTTER AND TEK TERMINAL
C FEATURES AUTOMATIC SELECTION OF DATA
C POINTS TO BE PLOTTED TO REPRODUCE THE
C CURVE TO DESIRED ACCURACY.
C***MODIFIED TO USE USER SUPPLIED MAX AND MIN FOR MULTIPLE PLOTS
C
SUBROUTINE GPLOT(V,H,N,VS,HS,IV,IH,IOPT,M,ERR,IERASE,
1VMAX,VMIN,HMAX,HMIN)
C
C V,H LINEAR ARRAYS FOR HORIZ AND VERT DATA
C N LENGTH OF V AND H
C VS,HS 0. AUTO SCALING
C <0. AUTO SCALING, AXIS SPAN OR # CYCLES
C RETURNED IN VS AND/OR HS
C >0. SCALING SPECIFIED, VS AND/OR HS ARE
C EQUAL TO THE DESIRED SPAN
C IV,IH 0 LINEAR
C 1 LOG
C 2 POLAR
C IOPT IOPT=IOPT1+IOPT2+IOPT3+IOPT4
C IOPT1 0 NOPLOT IOPT3 0 NO BOX
C 1 PLOT DATA 100 LARGE BOX
C 2 SAVE DATA 200 SMALL BOX
C
C IOPT2 0 NO AXES IOPT4 0 PLOT ALL
C 10 AXES ONLY 1000 OMIT LAST POINT
C 20 AXES+LABELS
C IOPT5 0 REGULAR PLOT
C 10000 FULL SCREEN PLOT
C (TEK ONLY)
C M 0 LINES
C 1-6 SYMBOLS
C 'L' WHERE L IS ANY CHARACTER TO BE PLOTTED
C AT THE DATA POINTS
C ERR ERROR (FINE=5.,MEDIUM=20.,COARSE=50.)
C IERASE 0 OVERPLOT
C 1 ERASE
C
DIMENSION V(1),H(1),ST(2)
LOGICAL PLTDATA,LASTPT,BOX,AXES,LABELS,ERASE,SBOX,BIGPLOT,SAVE
C
2 FORMAT(' PLTL')
3 FORMAT(' PLTT'/)
4 FORMAT(1X,I4,1X,I4)
5 FORMAT(' X AXIS SCALE IS',G,'PER DIVISION')
6 FORMAT(' Y AXIS SCALE IS',G,'PER DIVISION')
7 FORMAT(1X,I4,1X,I4,'^')
8 FORMAT(' X AXIS LOG ORIGIN AT',G)
9 FORMAT(' Y AXIS LOG ORIGIN AT',G)
10 FORMAT(' CAN''T GRAPH NEGATIVES OR ZERO ON LOG SCALE')
11 FORMAT(7X,'NUMBER OF CYCLES IS',I2)
12 FORMAT(' X AXIS MINIMUM IS',G)
13 FORMAT(' Y AXIS MINIMUM IS',G/)
14 FORMAT('0IMPROPER VALUE FOR GPLOT SUBROUTINE ARGUMENT!!'//)
15 FORMAT(G10.4)
16 FORMAT(I3)
17 FORMAT('0VERTICAL DATA RANGE TOO SMALL !! '//)
18 FORMAT('0HORIZONTAL DATA RANGE TOO SMALL !! '//)
C
C CHECK ARGUMENTS
C
20 IF(N.LE.0)GO TO 90
IF(((IV/3).NE.0).OR.((IH/3).NE.0))GO TO 90
IF(IV.EQ.2)IH=2
IF(IH.EQ.2)IV=2
BIGPLOT=.FALSE.
SAVE=.FALSE.
LASTPT=.TRUE.
ERASE=.TRUE.
BOX=.TRUE.
SBOX=.TRUE.
AXES=.TRUE.
PLTDATA=.TRUE.
LABELS=.TRUE.
IX=MOD(IOPT,10)
IF(IX.EQ.0)PLTDATA=.FALSE.
IF(IX.EQ.2)SAVE=.TRUE.
IX=MOD(IOPT/10,10)
IF(IX.LT.2)LABELS=.FALSE.
IF(IX.LT.1)AXES=.FALSE.
IX=MOD(IOPT/100,10)
IF(IX.LT.2)SBOX=.FALSE.
IF(IX.LT.1)BOX=.FALSE.
IX=MOD(IOPT/1000,10)
IF(IX.EQ.1)LASTPT=.FALSE.
IX=MOD(IOPT/10000,10)
IF(IX.EQ.1)BIGPLOT=.TRUE.
IF(IERASE.EQ.0)ERASE=.FALSE.
IF(.NOT.PLTDATA.AND.(IOPT.NE.0))GOTO 1040
C
IF(VMAX.NE.VMIN)GO TO 9123
CALL MINIMU(V,N,VMIN)
CALL MAXIMU(V,N,VMAX)
9123 IF(HMAX.NE.HMIN)GO TO 9124
CALL MINIMU(H,N,HMIN)
CALL MAXIMU(H,N,HMAX)
9124 CALL PLOTTER(TYPE,ITYPE)
IF(SAVE)CALL SAVSET(1)
25 IF(VS.LE.0.)GO TO 100
C
C V SCALING SPECIFIED
C
30 VSPAN=VS
VSCALE=VSPAN/10.
IF(IV-1)50,35,50
35 NVCYCL=VS
VSCALE=7999./VS
GOTO 300
50 IF(HS.LE.0.)GO TO 200
C
C H SCALING SPECIFIED
C
60 HSPAN=HS
HSCALE=HSPAN/15.
IF(IH-1)1000,70,80
70 NHCYCL=HS
HSCALE=7999./HS
GOTO 400
80 SPAN=HS
GOTO 260
C
C POLAR AUTO SCALING
C
240 IND=1
SPAN=VSPAN
SCALE=VSCALE
GO TO 180
260 IND=0
VSPAN=AMAX1(SPAN,VSPAN)
VSCALE=VSPAN/10.
HSPAN=VSPAN*1.5
HSCALE=VSCALE
GO TO 1000
C
C ERROR RETURN
C
90 TYPE 14
RETURN
92 TYPE 17
RETURN
94 TYPE 18
RETURN
C
C LINEAR AUTO SCALING
C
100 IF(IV.EQ.1)GO TO 300
RANGE=VMAX-VMIN
IF(RANGE.LT.1.E-35)GOTO 92
IND=0
IF(IV.EQ.0)GO TO 160
RANGE=2.*AMAX1(VMAX,ABS(VMIN))
160 NN=INT(ALOG10(RANGE))
IF(RANGE.LT.1.)NN=NN-1
P=10.**NN
R = RANGE / P
180 VSPAN=10.*P
IF(R.LE.7.51)VSPAN=7.5*P
IF(R.LE.5.01)VSPAN=5.*P
IF(R.LE.2.51)VSPAN=2.5*P
IF(R.LE.2.01)VSPAN=2.*P
IF(R.LE.1.01)VSPAN=P
VSCALE=VSPAN/10.
IF(IND.EQ.1)GO TO 260
GO TO 50
200 IF(IH.EQ.1)GO TO 400
RANGE=HMAX-HMIN
IF(RANGE.LT.1.E-35)GOTO 94
IF(IH.EQ.0)GO TO 220
RANGE=2.*AMAX1(HMAX,ABS(HMIN))
220 NN=INT(ALOG10(RANGE))
IF(RANGE.LT.1.)NN=NN-1
P=10.**NN
R=RANGE/P
IF(IH.EQ.2)GO TO 240
HSPAN=10.*P
IF(R.LE.7.51)HSPAN=7.5*P
IF(R.LE.5.01)HSPAN=5.*P
IF(R.LE.3.76)HSPAN=3.75*P
IF(R.LE.3.01)HSPAN=3.*P
IF(R.LE.1.51)HSPAN=1.5*P
HSCALE=HSPAN/15.
GO TO 1000
C
C AUTO SCALING - LOG
C
300 IF((VMAX.GT.0.).AND.(VMIN.LE.0.))GO TO 390
TT=ALOG10(ABS(VMIN))
NVMIN=INT(TT+((SIGN(1.002,TT)-1.)/2.))
IF(VS.GT.0.)GOTO 50
TT=ALOG10(ABS(VMAX))
NVMAX=INT(TT+((SIGN(.998,TT)+1.)/2.))
NVCYCL=IABS(NVMAX-NVMIN)
VSCALE=7999./NVCYCL
GO TO 50
390 TYPE 10
RETURN
400 IF((HMAX.GT.0.).AND.(HMIN.LE.0.))GO TO 390
TT=ALOG10(ABS(HMIN))
NHMIN=INT(TT+((SIGN(1.002,TT)-1.)/2.))
IF(HS.GT.0.)GOTO 1000
TT=ALOG10(ABS(HMAX))
NHMAX=INT(TT+((SIGN(.998,TT)+1.)/2.))
NHCYCL=IABS(NHMAX-NHMIN)
HSCALE=7999./NHCYCL
GO TO 1000
C
C SET PLOT WINDOWS
C
1000 IF(ITYPE.NE.2)GO TO 1002
IF(.NOT.ERASE)GO TO 990
CALL NEWPAG
990 IF(IH.EQ.2)GO TO 1003
CALL VWINDO(0.,9999.,0.,9999.)
CALL SWINDO(100,900,50,600)
IF(BIGPLOT)CALL SWINDO(0,1000,50,666)
GO TO 1002
1003 CALL VWINDO(0.,7333.,0.,9999.)
CALL SWINDO(250,770,50,700)
IF(BIGPLOT)CALL SWINDO(140,858,0,780)
C
C SAVE DATA
C
1002 IF(.NOT.PLTDATA)GOTO 1004
HHMIN=10.**NHMIN
VVMIN=10.**NVMIN
CALL CHGDEV('DSK14',4)
CALL DEFINE FILE(4,0,NR,'DATA.TMP',0,0)
WRITE(4)(H(I),I=1,N)
WRITE(4)(V(I),I=1,N)
END FILE (4)
IF(ITYPE.NE.2)TYPE 2
C
C SCALE DATA
C
1004 VSHIFT=0.
HSHIFT=0.
IF(IV.EQ.0)VSHIFT=AMOD(VMIN,VSCALE)
IF(IH.EQ.0)HSHIFT=AMOD(HMIN,HSCALE)
VSHIFT=VSHIFT+VSCALE*(.5-SIGN(.5,VSHIFT))
HSHIFT=HSHIFT+HSCALE*(.5-SIGN(.5,HSHIFT))
IF(.NOT.PLTDATA)GOTO 692
1005 DO 1007 I=1,N
IF(IV.EQ.0)V(I)=((V(I)-VMIN+VSHIFT)*7999./VSPAN)
IF(IV.EQ.1)V(I)=(VSCALE*ALOG10(ABS(V(I)/VVMIN)))
IF(IV.EQ.2)V(I)=V(I)*7999./VSPAN+3999.
V(I)=V(I)+999.
C
IF(IH.EQ.0)H(I)=((H(I)-HMIN+HSHIFT)*7999./HSPAN)
IF(IH.EQ.1)H(I)=(HSCALE*ALOG10(ABS(H(I)/HHMIN)))
IF(IH.EQ.2)H(I)=H(I)*7999./HSPAN+2666.
H(I)=H(I)+1333.
1007 CONTINUE
C
IF(M.EQ.0)GO TO 600
C
C PLOT SYMBOLS
C
DO 1009 I=1,N
JX=H(I)
JY=V(I)
IF((I.EQ.N).AND..NOT.LASTPT)GOTO 690
CALL SYMPLT(JX,JY,M,1.,ITYPE)
1009 CONTINUE
GO TO 690
C
C PLOT FIRST POINT
C
600 JX=H(1)
JY=V(1)
JXL=JX
JYL=JY
IS=3
CALL PLOUT(JX,JY,ITYPE,1)
C
C PLOT THIS POINT?
C
610 IF(IS.GT.N)GO TO 655
DO 650 I=IS,N
SLOPE=(V(I)-JYL)/(H(I)-JXL)
ASLOPE=ABS(SLOPE)
IF=I-1
DO 630 J=IS-1,IF
IF(ASLOPE.LE.1.)DERR=ABS(V(J)-SLOPE*(H(J)-JXL)-JYL)
IF(ASLOPE.GT.1.)DERR=ABS(H(J)-(V(J)-JYL)/SLOPE-JXL)
IF(DERR.GE.ERR)GO TO 660
IF((IF-IS).GT.50)GO TO 660
630 CONTINUE
650 CONTINUE
C
C FINISH
C
655 IF(.NOT.LASTPT)GOTO 690
JX=H(N)
JY=V(N)
CALL PLOUT(JX,JY,ITYPE,3)
GOTO 690
C
C SET-UP POINT TO BE PLOTTED
C
660 JX=H(I-1)
JY=V(I-1)
IS=I+1
GO TO 1006
C
C RESTORE DATA
C
690 CALL DEFINE FILE(4,0,NR,'DATA.TMP',0,0)
READ(4)(H(I),I=1,N)
READ(4)(V(I),I=1,N)
CALL RELEAS(4)
CALL DELETE('DATA.TMP')
CALL RSTDEV
692 IF(VS.GE.0.)GOTO 695
VS=VSPAN
IF(IV.EQ.1)VS=NVCYCL
695 IF(HS.GE.0.)GOTO 1010
HS=HSPAN
IF(IH.EQ.1)HS=NHCYCL
GO TO 1010
C
C CHECK LENGTH OF LINE AND PLOT
C
1006 VECT=0.
IF(ITYPE.EQ.1)
1 VECT=SQRT((15.*(JX-JXL)/7999.)**2.+(10.*(JY-JYL)/7999.)**2.)
NDIV=(VECT/3.01+1.0)
JXD=(JX-JXL)/NDIV
JYD=(JY-JYL)/NDIV
DO 1008 IC=1,NDIV
JXL=JXL+JXD
JYL=JYL+JYD
CALL PLOUT(JXL,JYL,ITYPE,3)
1008 CONTINUE
GO TO 610
C
C OUTPUT SCALE INFO
C
1010 HMINSH=HMIN-HSHIFT
VMINSH=VMIN-VSHIFT
IF(IV.NE.2)GOTO 1015
HMINSH=-5.*HSCALE
VMINSH=HMINSH
1015 IF(ITYPE.EQ.2)GO TO 1040
TYPE 3
IF(IH.EQ.1)GO TO 1050
1020 TYPE 5,HSCALE
IF(IV.EQ.1)GO TO 1060
1030 TYPE 6,VSCALE
1035 TYPE 12,HMINSH
TYPE 13,VMINSH
1040 IF(AXES)GOTO 2000
IF(BOX)GOTO 3000
1045 IF(ITYPE.EQ.2)GO TO 2500
GOTO9000
1050 TYPE 8,HHMIN
TYPE 11,NHCYCL
IF(IV.NE.1)GO TO 1030
1060 TYPE 9,VVMIN
TYPE 11,NVCYCL
GO TO 1035
C
C HORIZONTAL AXIS PRINT-OUT
C
2000 NDIV=1
IF(ITYPE.EQ.1)NDIV=4
JX=1333
IF(IV.NE.1)JY=(-VMIN+VSHIFT)*7999./VSPAN+999
IF(JY.LT.999)JY=999
IF(JY.GT.8999)JY=8999
IF(IV.EQ.1)JY=999
IF(IH.EQ.2)JY=4999
IF(ITYPE.EQ.1)TYPE 2
CALL PLOUT(JX,JY,ITYPE,1)
INC=8000/NDIV
IF(IH.EQ.2)INC=5333/NDIV
DO 2010 I=1,NDIV
JX=JX+INC
2010 CALL PLOUT(JX,JY,ITYPE,3)
JYZ=JY+125
JYZZ=JY-125
JXST=1333
JXFIN=JX
JYZERO=JY
JX=1333
IF(IH.EQ.1)GOTO 2100
C
C LINEAR HORIZONTAL TICS
C
DO 2020 I=1,16
CALL PLOUT(JX,JYZ,ITYPE,1)
CALL PLOUT(JX,JYZZ,ITYPE,3)
JX=JX+533
IF((I.GT.10).AND.(IH.EQ.2))GOTO 2030
2020 CONTINUE
2030 JX=JX-266
2035 CALL PLOUT(JX+60,JYZ,ITYPE,1)
CALL PLOUT(JX-60,JYZZ,ITYPE,3)
CALL PLOUT(JX-60,JYZ,ITYPE,1)
CALL PLOUT(JX+60,JYZZ,ITYPE,3)
C
C VERTICAL AXIS PRINT-OUT
C
2050 JY=999
IF(IH.NE.1)JX=(-HMIN+HSHIFT)*7999./HSPAN+1333
IF(JX.LT.1333)JX=1333
IF(JX.GT.9333.)JX=9333
IF(IH.EQ.1)JX=1333
IF(IV.EQ.2)JX=3999
CALL PLOUT(JX,JY,ITYPE,1)
NDIV=1
IF(ITYPE.EQ.1)NDIV=4
INC=8000/NDIV
DO 2060 I=1,NDIV
JY=JY+INC
2060 CALL PLOUT(JX,JY,ITYPE,3)
JYST=999
JYFIN=8999
JXZERO=JX
JY=999
JXZ=JX+90
JXZZ=JX-90
IF(IV.EQ.1)GOTO 2300
C
C LINEAR VERTICAL TICS
C
DO 2080 I=1,11
CALL PLOUT(JXZ,JY,ITYPE,1)
CALL PLOUT(JXZZ,JY,ITYPE,3)
2080 JY=JY+800
JY=JY-400
2085 CALL PLOUT(JX-60,JY+125,ITYPE,1)
CALL PLOUT(JX,JY,ITYPE,3)
CALL PLOUT(JX+60,JY+125,ITYPE,1)
CALL PLOUT(JX,JY,ITYPE,3)
CALL PLOUT(JX,JY-125,ITYPE,3)
GOTO 2450
C
C LOG HORIZONTAL TICS
C
2100 NC=NHCYCL
JYQ=JY+250
JYQQ=JY-250
IT=(NC+4)/5
IF(IT.LE.3)GOTO 2110
NC=2-(NC-(NC/2)*2)
2110 NSEG=7999/NC
SEG=FLOAT(NSEG)
DO 2140 I=1,NC
CALL PLOUT(JX,JYQ,ITYPE,1)
CALL PLOUT(JX,JYQQ,ITYPE,3)
IF(IT.GT.3)GOTO 2130
DO 2120 J=2,10-IT,IT
JXP=JX+INT(SEG*ALOG10(FLOAT(J)))
CALL PLOUT(JXP,JYZ,ITYPE,1)
CALL PLOUT(JXP,JYZZ,ITYPE,3)
2120 CONTINUE
2130 JX=JX+NSEG
2140 CONTINUE
CALL PLOUT(JX,JYQ,ITYPE,1)
CALL PLOUT(JX,JYQQ,ITYPE,3)
JX=JX+266
GOTO 2035
C
C LOG VERTICAL TICS
C
2300 NC=NVCYCL
JXQ=JX+180
JXQQ=JX-180
IT=(NC+4)/5
IF(IT.LE.2)GOTO 2310
NC=2-(NC-(NC/2)*2)
2310 NSEG=7999/NC
SEG=(FLOAT(NSEG))
DO 2340 I=1,NC
CALL PLOUT(JXQ,JY,ITYPE,1)
CALL PLOUT(JXQQ,JY,ITYPE,3)
IF(IT.GT.2)GOTO 2330
DO 2320 J=2,10-IT,IT
JYP=JY+INT(SEG*ALOG10(FLOAT(J)))
CALL PLOUT(JXZ,JYP,ITYPE,1)
CALL PLOUT(JXZZ,JYP,ITYPE,3)
2320 CONTINUE
2330 JY=JY+NSEG
2340 CONTINUE
CALL PLOUT(JXQ,JY,ITYPE,1)
CALL PLOUT(JXQQ,JY,ITYPE,3)
JY=JY+400
GOTO 2085
2450 IF(LABELS)GOTO 2600
IF(BOX)GOTO 3000
IF(ITYPE.EQ.2)GO TO 2500
2460 TYPE 3
GOTO9000
2500 CALL HOME
CALL ANMODE
GOTO9000
C
C LABEL AXES
C
2600 IF(.NOT.SAVE)GOTO2605
IF(BOX)GOTO3000
GOTO9000
2605 IF(ITYPE.EQ.1)GOTO 2700
C
C LABELS FOR TEK TERMINAL
C
C X AXIS LABELS
C
IF(IH.EQ.1)GOTO 2610
CALL PLOUT(JXST-266,JYZERO-600,ITYPE,1)
CALL ANMODE
CALL LOUT(HMINSH,6,3,0)
CALL PLOUT(JXFIN-266,JYZERO-600,ITYPE,1)
CALL ANMODE
Q=HMINSH+HSCALE*(15-5*IH/2)
CALL LOUT(Q,6,3,0)
GO TO 2650
2610 CALL PLOUT(JXST,JYZERO-800,ITYPE,1)
CALL BAKSP
CALL BAKSP
CALL LOUT(10,0)
CALL TOUTPT(11)
CALL LOUT(NHMIN,0)
CALL PLOUT(JXFIN,JYZERO-800,ITYPE,1)
CALL BAKSP
CALL BAKSP
CALL LOUT(10,0)
CALL TOUTPT(11)
NQ=NHMIN+NHCYCL
CALL LOUT(NQ,0)
C
C Y AXIS LABELS
C
2650 IF(IV.EQ.1)GOTO 2660
CALL PLOUT(JXZERO-1333,JYST+200,ITYPE,1)
CALL ANMODE
CALL LOUT(VMINSH,6,3,0)
CALL PLOUT(JXZERO-1333,JYFIN-400,ITYPE,1)
CALL ANMODE
Q=VMINSH+VSCALE*10.
CALL LOUT(Q,6,3,0)
IF(BOX)GOTO 3000
GOTO 2500
2660 CALL PLOUT(JXZERO-666,JYST,ITYPE,1)
CALL BAKSP
CALL BAKSP
CALL LOUT(10,0)
CALL TOUTPT(11)
CALL LOUT(NVMIN,0)
CALL PLOUT(JXZERO-666,JYFIN-600,ITYPE,1)
CALL BAKSP
CALL BAKSP
CALL LOUT(10,0)
CALL TOUTPT(11)
NQ=NVMIN+NVCYCL
CALL LOUT(NQ,0)
IF(BOX)GOTO 3000
GOTO 2500
C
C LABELS FOR HP PLOTTER
C
C X AXIS LABELS
C
2700 IF(IH.EQ.1)GOTO 2720
ENCODE(10,15,ST)HMINSH
CALL TITLE(ST,2,JXST-600,JYZERO-400,2.,0.)
Q=HMINSH+HSCALE*(15-5*IH/2)
ENCODE(10,15,ST)Q
CALL TITLE(ST,2,JXFIN-600,JYZERO-400,2.,0.)
GOTO 2750
2720 CALL TITLE(' 10',1,JXST-600,JYZERO-400,2.,0.)
ENCODE(5,16,ST)NHMIN
CALL LJUST(ST,1)
CALL TITLE(ST,1,JXST,JYZERO-300,2.,0.)
CALL TITLE(' 10',1,JXFIN-600,JYZERO-400,2.,0.)
NQ=NHMIN+NHCYCL
ENCODE(5,16,ST)NQ
CALL LJUST(ST,1)
CALL TITLE(ST,1,JXFIN,JYZERO-300,2.,0.)
C
C Y AXIS LABELS
C
2750 IF(IV.EQ.1)GOTO 2770
ENCODE(10,15,ST)VMINSH
CALL TITLE(ST,2,JXZERO-1333,JYST+100,2.,0.)
Q=VMINSH+VSCALE*10.
ENCODE(10,15,ST)Q
CALL TITLE(ST,2,JXZERO-1333,JYFIN-200,2.,0.)
IF(BOX)GOTO 3000
GOTO 2460
2770 CALL TITLE(' 10',1,JXST-960,JYST,2.,0.)
ENCODE(5,16,ST)NVMIN
CALL LJUST(ST,1)
CALL TITLE(ST,1,JXST-360,JYST+100,2.,0.)
CALL TITLE(' 10',1,JXST-960,JYFIN-200,2.,0.)
NQ=NVMIN+NVCYCL
ENCODE(5,16,ST)NQ
CALL LJUST(ST,1)
CALL TITLE(ST,1,JXST-360,JYFIN-100,2.,0.)
IF(BOX)GOTO 3000
GOTO 2460
C
C DRAW BOX
C
3000 IF(SBOX)GOTO 3040
JXST=0
JXFIN=9999
IF(IH.EQ.2)JXFIN=7333
JYST=0
JYFIN=9999
GOTO 3050
3040 JXST=1333
JYST=999
JXFIN=9333
IF(IH.EQ.2)JXFIN=6666
JYFIN=8999
3050 CALL PLOUT(JXST,JYST,ITYPE,1)
INCX=(JXFIN-JXST)/NDIV
INCY=(JYFIN-JYST)/NDIV
JX=JXST
JY=JYST
DO 3060 I=1,NDIV
JX=JX+INCX
3060 CALL PLOUT(JX,JY,ITYPE,3)
DO 3070 I=1,NDIV
JY=JY+INCY
3070 CALL PLOUT(JX,JY,ITYPE,3)
DO 3080 I=1,NDIV
JX=JX-INCX
3080 CALL PLOUT(JX,JY,ITYPE,3)
DO 3090 I=1,NDIV
JY=JY-INCY
3090 CALL PLOUT(JX,JY,ITYPE,3)
IF(ITYPE.EQ.1)GOTO 2460
GOTO 2500
9000 IF(SAVE)CALL SAVSET(-1)
RETURN
END
C
C AUXILLARY SUBROUTINE
C
SUBROUTINE LJUST(ST,N)
C
C ST ASCII STRING
C N LENGTH OF ST
C
C ROUTINE LEFT-JUSTIFIES CONTENTS OF ST
C
DIMENSION ST(1)
C
NUM=5*N
DO 10 I=1,NUM
CALL GETCHR(ST,I,CHR)
IF(CHR.NE.' ')GOTO 20
10 CONTINUE
RETURN
20 IF(I.EQ.1)RETURN
DO 30 J=I,NUM
CALL GETCHR(ST,J,CHR)
30 CALL PUTCHR(ST,J-I+1,CHR)
DO 40 J=NUM-I+2,NUM
40 CALL PUTCHR(ST,J,' ')
RETURN
END