Trailing-Edge
-
PDP-10 Archives
-
decuslib20-02
-
decus/20-0075/spplt.bas
There is 1 other file named spplt.bas in the archive. Click here to see a list.
00002 ' ******************************************************************
00003 ' * THIS IS A HYBRID ORBITAL PLOTTING PROGRAM ORIGINATED BY *
00004 ' * STEPHEN L. HOLMGREN, LAWRENCE UNIVERSITY CLASS OF 1973, IN *
00005 ' * THE FORM OF CHAIN-LINKED PROGRAMS "PLOTV5" AND HYPLT7" WRITTEN *
00006 ' * IN BASIC-PLUS FOR USE ON A PDP-11/20 WITH A TSP-212 PLOTTER. *
00007 ' * THIS VERSION WAS ADAPTED FOR THE U OF O PDP-10 BY JAMES S. *
00008 ' * EVANS, VISITING FACULTY MEMBER FROM LAWRENCE UNIVERSITY. *
00009 ' ******************************************************************
00070 ' **********DEFINE CONSTANTS AND FUNCTIONS*****
00072 C0 = 0.25/SQR(2.*3.14159) ' PART OF NORMALIZATION CONSTANT
00074 A0 = 0.52915 ' BOHR RADIUS IN ANGSTROMS
00076 D0 = 0.05 ' CONVERGENCE CRITERION FOR ROOT POSITIONS = 5%
00078 S0 = 1.0E-20 ' DIVISION BY ZERO PREVENTATIVE
00080 B0 = 8.72664E-2 ' 5 DEGREES = PI/36 RADIANS
00082 DEF FNA(S)=C1*E*(S2*(2.-Z1*S)+S1*A*C*Z1*S) ' PSI(R,THETA)
00084 DEF FNB(S)=C2*E*(S2*(-2.+.5*Z1*S)+S1*A*C*(1.-.5*Z1*S)) ' D/DR(PSI)
00086 DEF FNE(R)=EXP(-Z1*R/2)
00088 DIM P(3), A1(38)
00090 DIM R1(37,3), R2(37,3), R3(37,3), R4(37,3)
00100 ' **********ENTRY FOR NEW CALCULATIONS*****
00102 LET C=S2=1.0
00104 LET A1(1)=0.0
00110 PRINT "WHAT EFFECTIVE NUCLEAR CHARGE";
00112 INPUT Z0
00120 PRINT "WHAT AMOUNT OF P CHARACTER";
00122 INPUT A
00124 IF A<20 THEN 128
00126 S2=0.0 ' DEFAULT TO PURE "P" IF A>20
00128 A=SQR(A)
00130 PRINT "HOW MANY CONTOURS";
00132 INPUT I5
00134 IF I5<=3 THEN 140
00136 I5=3
00138 PRINT "MAXIMUM OF 3 CONTOURS WILL BE USED"
00140 ' **********CALCULATE NORMALIZATION FACTORS*****
00142 Z1=Z0/A0
00144 Z9=1/SQR(Z0) ' MODIFY OUTER LIMITS FOR ROOT SEARCHES
00146 C1=(C0*Z1^1.5)/SQR(S2+A*A)
00148 C2=C1*Z1
00150 ' **********BEGIN ANGLE-STEP LOOP*****
00160 FOR J=1 TO 37
00162 LET J1=J2=1 ' J1,J2=2 WHEN MAX,MIN FOUND
00164 R0=0.0
00180 IF J=1 THEN 190
00182 IF J<=11 THEN 400
00190 IF A>1 THEN 200
00192 J1=2
00194 M1=0.0 ' MAX PSI OCCURS AT R=0 IF 0<A<=1
00200 ' **********BEGIN LOOP FOR FINDING MAX,MIN*****
00210 FOR I=1 TO 100
00212 E=FNE(R0)
00220 IF J1+J2=4 THEN 332
00230 ON J1 GO TO 240,280 ' HAS MAX BEEN FOUND?
00240 S1=1.0 ' POSITIVE LOBE SWITCH
00249 ' POS. SLOPE MEANS MAX FURTHER RIGHT
00250 IF FNB(R0)>=0. THEN 280
00260 M1=R0-0.025 ' MAX PSI SET HALFWAY BETWEEN PRESENT & PAST R0
00262 J1=2
00280 ON J2 GO TO 290,320 ' HAS MIN BEEN FOUND?
00290 S1=-1.0 ' NEGATIVE LOBE SWITCH
00300 IF FNB(R0)<=0. THEN 320
00310 M2=R0-0.025 ' MIN PSI SET HALFWAY BETWEEN PRESENT & PAST R0
00312 J2=2
00320 R0=R0+0.05 ' PRECISION OF 0.025 ANGSTROM FOR MAX,MIN
00330 NEXT I ' **********END LOOP FOR FINDING MAX,MIN*****
00332 IF J<>1 THEN 400
00333 S1=+1.0
00334 E=FNE(M1)
00335 PRINT "APPROX. PSI MAX =";FNA(M1);"AT R =";M1
00336 PRINT "CHOOSE YOUR OWN CONTOURS";
00338 INPUT F$
00340 IF F$="YES" THEN 349
00342 P(1)=0.03
00344 P(2)=0.06
00346 P(3)=0.12
00348 GO TO 354
00349 PRINT "PLEASE ENTER SMALLEST FIRST"
00350 FOR I=1 TO I5
00351 INPUT P(I)
00352 NEXT I
00353 PRINT
00354 PRINT "CONTOURS USED"
00356 FOR I=1 TO I5
00357 PRINT P(I);
00358 NEXT I
00359 PRINT
00400 C = COS(A1(J)) ' CALCULATE COSINE AT ANGLE DETERMINED BY J
00410 ' CALCULATE PSIMAX,MIN AND NODE POSITION FOR EACH ANGLE
00420 S1=+1.0
00422 E=FNE(M1)
00424 P2=FNA(M1) ' MAX PSI (IN POSITIVE LOBE)
00430 S1=-1.0
00432 E=FNE(M2)
00434 P1=FNA(M2) ' MIN PSI (IN NEGATIVE LOBE)
00440 LET M5=P4=0.0
00441 M6=M2
00442 R8=0.2
00443 C$=" NODE "
00444 IF S2=1. THEN 447
00445 R9=0.0 ' NODE ALWAYS AT 0 FOR PURE "P"
00446 GO TO 448
00447 GOSUB 5001
00448 N0=R9 ' LOCATION OF NODE IN PSI AT CURRENT ANGLE
00450 ' MAKE INITIAL GUESSES FOR EACH ROOT
00452 LET R1(0,1)=(M1+6.0*Z9)/2.
00454 LET R1(0,2)=R1(0,3)=M1+0.75
00460 FOR I=1 TO 3
00462 R2(0,I)=M1/2 +S0 ' PREVENT DIVISION BY ZERO
00464 R3(0,I)=(M2+N0)/2.
00466 R4(0,I)=(M2+6.0*Z9)/2.
00468 NEXT I
00470 ' **********BEGIN LOOP OVER CONTOURS*****
00472 FOR I=1 TO I5
00474 P4 = -ABS( P(I) ) ' SELECT CONTOUR VALUE
00489 ' TEST FOR EXISTENCE OF ROOTS 3&4
00490 IF P4<P1 THEN 590
00500 ' .....ROOT 3 DETERMINED FIRST.....
00502 R8=R3(J-1,I) ' TAKE ROOT FROM PREVIOUS ANGLE AS INITIAL GUESS
00503 ' ABANDON ROOTS 3&4 IF NOT FOUND AT LAST ANGLE
00504 IF R8<=-3. THEN 590
00506 S1=-1.0
00508 M5=N0
00510 M6=M2
00512 C$=" ROOT 3 "
00514 GOSUB 5001
00516 R3(J,I)=R9
00520 ' .....ROOT 4 SHOULD EXIST IF ROOT 3 WAS FOUND.....
00522 R8=R4(J-1,I) ' INITIAL GUESS
00523 ' ABANDON SEARCH IF NOT FOUND AT LAST ANGLE
00524 IF R8<=-3. THEN 590
00526 S1=-1.0
00528 M5=M2
00530 M6=7.0*Z9
00532 C$=" ROOT 4 "
00534 GOSUB 5001
00536 R4(J,I)=R9
00538 GO TO 600
00590 LET R3(J,I)=R4(J,I)=-3. ' FLAG FOR NON-EXISTENCE
00600 ' ROOTS 1&2 ARE TRICKIER!!
00601 ' ROOTS 1&2 ARE REDUNDANT FOR THETA>90 DEG
00602 IF J>19 THEN 750
00604 P4=+ABS( P(I) ) ' SELECT CONTOUR VALUE
00609 ' TEST FOR EXISTENCE OF ROOTS 1&2
00610 IF P4>=P2 THEN 750
00620 ' ROOT 2 CAN BE ON EITHER SIDE OF ORIGIN!! MUST DETERMINE WHICH
00622 E=1.0
00624 Q=P4-FNA(0.0)
00626 IF Q>0. THEN 640
00628 IF Q<0. THEN 660
00630 R2(J,I)=0.0 ' ROOT 2 AT ORIGIN IF Q=0
00632 GO TO 680
00640 ' .....ROOT 2 ON "RIGHT SIDE" OF ORIGIN.....
00642 R8=ABS( R2(J-1,I) ) ' USE PREVIOUS ROOT AS INITIAL GUESS
00646 S1=+1.0
00648 M5=0.0
00650 M6=M1
00652 C$=" ROOT 2+ "
00654 GO SUB 5001
00656 R2(J,I)=R9
00658 GO TO 680
00660 ' .....ROOT 2 ON "LEFT SIDE" OF ORIGIN.....
00662 R8=ABS( R2(J-1,I) ) ' USE PREVIOUS ROOT AS INITIAL GUESS
00666 S1=-1.0
00668 M5=0.0
00670 M6=N0
00672 C$=" ROOT 2- "
00674 GOSUB 5001
00676 R2(J,I)=-R9 ' NEGATIVE SIGN FLAG FOR LEFT OF ORIGIN
00680 ' .....ROOT 1.....
00682 R8=R1(J-1,I) ' USE PREVIOUS ROOT AS INITIAL GUESS
00686 S1=+1.0
00688 M5=M1*(1+SGN(M1))/2 ' ROOT 1 ALWAYS ON "RIGHT SIDE"
00690 M6=6.0
00692 C$=" ROOT 1 "
00694 GOSUB 5001
00696 R1(J,I)=R9
00698 GO TO 760
00750 LET R2(J,I)=R1(J,I)=-3.
00760 NEXT I ' **********END LOOP OVER CONTOURS*****
00770 A1(J+1)=A1(J)+B0 ' INCREMENT ANGLE BY 5 DEG
00780 NEXT J ' **********END LOOP OVER ANGLES*****
00902 R1(0,0)=0 ' SORT INDICATOR
00904 GO TO 2900 ' *****EXIT TO GRAPHICAL OUTPUT*****
01000 '**********MODIFIED SOFTWARE FOR TSP-12 PLOTTER CONTROLLER*****
01010 ' VARIABLES: L,P$,P0-P3,P6,P8,P9,X,X0-X4,Y,Y0-Y4
01020 DIM P$(70)
01100 ' ENTRY POINT: SWITCH OUTPUT FROM COMPUTER TO PLOTTER
01110 P$(0)=CHR$(16)
01120 P9=1
01130 RETURN
01200 ' ENTRY POINT: SWITCH OUTPUT FROM COMPUTER BACK TO TTY
01210 IF P9=1 THEN 1230
01220 GOSUB 1700
01230 PRINT
01240 PRINT CHR$(0);
01250 RETURN
01300 ' ENTRY POINT: MOVE PEN TO POINT X,Y; PEN UP IF L<>0
01304 LET P8=P6=0
01310 X2=INT(510*(X-X0)/(X1-X0))+INT(S)*P8
01312 IF X2>=0 THEN 1316
01314 X2=0
01316 IF X2<=510 THEN 1320
01318 X2=510
01320 Y2=INT(510*(Y-Y0)/(Y1-Y0))+INT(S)*P6
01322 IF Y2>=0 THEN 1326
01324 Y2=0
01326 IF Y2<=510 THEN 1330
01328 Y2=510
01330 IF L=0 THEN 1350
01340 L=64
01350 X3=X2-X3
01360 Y3=Y2-Y3
01370 P0=INT(ABS(X3)/50+ABS(Y3)/50)+1
01380 FOR P1=1 TO P0
01390 X4=X2-(X3/P0)*(P0-P1)
01400 Y4=Y2-(Y3/P0)*(P0-P1)
01410 P3=64+INT(X4/8)
01420 GOSUB 1620
01430 P3=64+INT(Y4/8)
01440 GOSUB 1620
01450 P3=INT(Y4)-INT(Y4/8)*8
01460 P3=P3+8*(INT(X4)-INT(X4/8)*8)+L
01470 GOSUB 1550
01480 IF P9<70 THEN 1500
01490 GOSUB 1700
01500 NEXT P1
01510 L=0
01520 X3=X2
01530 Y3=Y2
01540 RETURN
01550 IF P3>13 THEN 1600
01560 IF P3=0 THEN 1580
01570 IF P3<9 THEN 1640
01580 P3=P3+8
01590 GO TO 1640
01600 IF P3<127 THEN 1640
01610 P3=126
01620 IF P3<127 THEN 1640
01630 P3=63
01640 P$(P9)=CHR$(P3)
01650 P9=P9+1
01660 RETURN
01700 PRINT
01710 FOR P2=0 TO P9-1
01720 PRINT P$(P2);
01730 NEXT P2
01740 P9=1
01750 RETURN
02900 ' **********ENTRY FROM CALCULATION OF ROOTS*****
02920 ' FUNCTION TO ADJUST SUBSCRIPTS FOR ROTATIONS
02930 DEF FNC(S)
02932 B=M+S
02934 IF B>=0 THEN 2940
02936 B=B+72
02938 GO TO 2934
02940 IF B<=72 THEN 2946
02942 B=B-72
02944 GO TO 2940
02946 FNC=B
02948 FNEND
02950 S=1 ' INITIALIZE SCALE,TRANS,ROTATE PARAMETERS
02952 LET T1=T2=M=0
02970 DIM S(72), C(72), X(150), Y(150)
02980 FOR I=0 TO 72
02982 B = B0*I ' 5 DEG INCREMENTS
02984 S(I)=SIN(B)
02986 C(I)=COS(B)
02988 NEXT I ' STORE SINES & COSINES NEEDED
02990 PRINT "COMPARE OR MAX SCALE";
02992 INPUT F2$
03000 X8=-R4(1,1) ' *****DETERMINE RANGE OF DATA*****
03002 X9=R1(1,1)
03004 Y8=-1.2*S(14)*R4(14,1)
03006 Y9=-Y8
03008 IF F2$<>"MAX SCALE" THEN 3028
03009 ' IN "MAX SCALE" MODE, OUTER CONTOUR WILL
03010 IF Y8>=X8 THEN 3016
03012 LET X0=Y0=Y8 ' FILL THE PAGE (UNLESS "SCALE" IS EXERCISED).
03014 GO TO 3018 ' PLOT SPAN IS SAME ON BOTH AXES.
03016 LET X0=Y0=X8
03018 IF Y9<=X9 THEN 3024
03020 LET X1=Y1=Y9
03022 GO TO 3032
03024 LET X1=Y1=X9
03026 GO TO 3032
03028 LET X0=Y0=-6 ' IN "COMPARE" MODE, BOUNDARIES OF ALL GRAPHS
03030 LET X1=Y1=6 ' ARE +/- 6 ANGSTROMS
03032 PRINT "OPTIONS: SCALE,TRANS,ROTATE,PLOT,AXIS,MORE,FINI"
03034 PRINT
03036 GO TO 4020 ' SEEK ROUTING OPTION
03040 ' **********CALCULATE & PLOT POSITIVE LOBES*****
03042 GOSUB 1100 ' TURNS ON THE PLOTTER
03044 FOR I=1 TO I5
03046 N=-1
03048 FOR K=1 TO 37
03050 R=R2(K,I)
03052 K2=K-1 ' ADJUST SUBSCRIPTS: HERE 0 MEANS 0 DEG
03059 ' TEST FOR FLAGS
03060 IF R<=-2 THEN 3120
03070 N=N+2
03072 ' EACH POINT AND ITS REFLECTION IN THE AXIS OF SYMMETRY
03074 ' ARE STORED IN CONSECUTIVE ORDER IN THE ARRAYS X&Y.
03089 ' IF R2 IS "LEFT" OF ORIGIN, THEN USE
03090 IF R<0 THEN 3096
03092 K3=K2 ' DIFFERENT ANGLE.....QUITE TRICKY!!
03094 GO TO 3100
03096 K3=72-K2
03100 GOSUB 5400 ' ASSIGN VALUES TO X,Y
03120 NEXT K
03130 FOR K=37 TO 1 STEP -1 ' NOW CONTINUE WITH ROOT 1
03140 K2=K-1
03142 R=R1(K,I)
03149 ' TEST FOR FLAGS
03150 IF R<0 THEN 3164
03158 N=N+2
03160 K3=K2
03162 GOSUB 5400 ' ASSIGN VALUES TO X,Y
03164 NEXT K
03166 X(0)=X(1) ' SET LAST POINT = FIRST POINT
03168 Y(0)=Y(1)
03170 GOSUB 5500 ' PLOT POS. LOBE OF THIS CONTOUR
03172 NEXT I ' **********END LOOP OVER POSITIVE LOBES*****
03180 L=1 ' PREPARE FOR NEGATIVE LOBE
03182 GOSUB 1300
03184 IF R1(0,0)=10 THEN 3200
03186 GOSUB 5600 ' PERFORM CONTINUITY CHECK ONLY ONCE PER ARRAY
03200 FOR I=1 TO I5 ' **********CALCULATE & PLOT NEG. LOBES*****
03202 N=-1
03210 FOR K=1 TO 37
03212 R=-R3(K,I) ' WE ARE TO THE LEFT OF ORIGIN
03214 K2=K-1 ' ADJUST SUBSCRIPTS: HERE 0 MEANS 0 DEG
03219 ' TEST FOR FLAGS
03220 IF R>0 THEN 3260
03230 N=N+2
03232 K3=K2
03240 GOSUB 5400 ' ASSIGN VALUES TO X,Y
03260 NEXT K
03270 FOR K=37 TO 1 STEP -1 ' NOW CONTINUE WITH ROOT 4
03272 K2=K-1
03280 R=-R4(K,I)
03281 ' TEST FOR FLAGS
03282 IF R>0 THEN 3310
03300 N=N+2
03302 K3=K2
03304 GOSUB 5400 ' ASSIGN VALUES TO X,Y
03310 NEXT K
03320 X(0)=X(1) ' SET LAST POINT = FIRST POINT
03322 Y(0)=Y(1)
03330 GOSUB 5500 ' PLOT NEG. LOBE OF THIS CONTOUR
03340 NEXT I ' *****END LOOP OVER NEGATIVE LOBES*****
04000 ' **********ROUTING OPTIONS & ASSOCIATED CALCULATIONS*****
04010 GOSUB 1200 ' TURN OFF PLOTTER
04012 LET S=1 ' RE-INITIALIZE SCALE,TRANS,ROTATE PARAMETERS
04014 LET T1=T2=M=0
04020 PRINT "YOUR CHOICE"; ' SEEK ROUTING OPTION
04022 INPUT C$
04030 IF C$<>"SCALE" THEN 4070
04040 PRINT "REDUCE SCALE BY FACTOR OF";
04042 INPUT S
04050 IF S>0 THEN 4020
04060 PRINT "POSITIVE SCALE FACTOR ONLY"
04062 S=1
04064 GO TO 4020 ' END "SCALE" SECTION
04070 IF C$<>"TRANS" THEN 4110
04080 X=ABS(S*X0-X8) ' COMPUTE ALLOWABLE TRANSLATIONS
04082 Y=ABS(S*Y0-Y8)
04084 PRINT "MAXIMUM TRANSLATION = ";X;" IN X OR Y DIRECTION"
04090 PRINT "X-TRANS =";
04092 INPUT T1
04094 IF ABS(T1)>X THEN 4104
04096 PRINT "Y-TRANS =";
04098 INPUT T2
04102 IF ABS(T2)<=Y THEN 4020
04104 PRINT "TRY AGAIN"
04106 GO TO 4084 ' END "TRANS" SECTION
04109 ' EXIT TO PLOT INSTRUCTIONS
04110 IF C$="PLOT" THEN 3040
04120 IF C$<>"ROTATE" THEN 4150
04130 PRINT "DEGREES OF ROTATION =";
04132 INPUT D
04134 IF D>=0 THEN 4140
04136 D=360+D
04140 M=INT(D/5 + .5) ' ROUNDED TO NEAREST 5 DEG
04142 GO TO 4020 ' END "ROTATE" SECTION
04150 IF C$<>"AXIS" THEN 4230
04154 IF (S-1)*T1*T2<>0 THEN 4170
04158 PRINT "WISH TO SCALE OR TRANS FIRST";
04160 INPUT M$
04162 IF M$="YES" THEN 4020
04170 S5=.02*ABS(Y0)
04172 Y=(Y8+T2)/S + S5
04174 X=(-1+T1)/S
04180 L=1
04182 GOSUB 1100
04184 GOSUB 1300
04190 Y=Y-S5
04192 GOSUB 1300
04194 X=0+T1/S
04196 GOSUB 1300
04198 Y=Y+S5
04200 GOSUB 1300
04202 Y=Y-S5
04204 GOSUB 1300
04206 X=(1+T1)/S
04208 GOSUB 1300
04210 Y=Y+S5
04212 GOSUB 1300
04214 GOSUB 1200
04220 GO TO 4020 ' END "AXIS" SECTION
04230 L=1 ' FOR "MORE" OR "FINI", PULL PEN ASIDE
04232 X=X1
04234 GOSUB 1100
04236 GOSUB 1300
04238 GOSUB 1200 ' TURN PLOTTER OFF, RETURN TO TTY
04240 IF C$<>"MORE" THEN 4250
04242 GO TO 100 ' SEEK NEW CONTOUR PARAMETERS
04250 IF C$<>"FINI" THEN 4020
04252 STOP
05001 ' **********ROOT-FINDING SUBROUTINES*****
05002 ' R8 = INITIAL GUESS R9 = FINAL ROOT
05003 ' M6 = MAX. OF INTERVAL M5 = MIN. OF INTERVAL
05004 ' **********BEGIN NEWTON-RAPHSON METHOD*****
05006 T=-1 ' T = COUNTER FOR 3 TRIES ALLOWED THIS METHOD
05008 Q=0 ' Q = ON-OFF SWITCH FOR SECOND METHOD
05010 R9=R8
05012 D1=2.0
05020 FOR K=1 TO 30 ' MAX 30 ITERATIONS PER TRY
05022 E=FNE(R9)
05030 C3 = (FNA(R9)-P4)/(FNB(R9)+S0) ' DIV. BY 0 PROTECTION
05031 ' DIVERGENCE LIMIT 10 ANGSTROMS
05032 IF ABS(C3)>10. THEN 5050
05040 R8=R9-C3
05042 D1=ABS(C3/(R9+S0)) ' DIV. BY 0 PROTECTION
05044 R9=R8
05045 ' CONVERGENCE CRITERION IS 100*D0%
05046 IF D1<D0 THEN 5060
05048 GO TO 5070
05050 R9=-1.0
05060 K=30
05070 NEXT K
05072 IF R9<M6 THEN 5076
05074 GO TO 5078
05075 ' RETURN IF ROOT IN PROPER RANGE
05076 IF R9>M5 THEN 5094
05078 T=T+2
05080 IF T>3 THEN 5090
05082 R9=(M6-M5)*T/4. + M5 ' TRY ANOTHER STARTING VALUE WITHIN INTERVAL
05084 GO TO 5012
05089 ' SWITCH TO SECOND METHOD AFTER 3 TRIES
05090 IF Q=0 THEN 5200
05092 R9=-2.
05094 RETURN ' **********END NEWTON-RAPHSON METHOD*****
05200 ' **********BEGIN INTERVAL-HALVING METHOD*****
05202 ' R6,R5 = UPPER,LOWER LIMITS FOR AXIS CROSSING
05220 Q=1
05230 R6=M6
05232 R5=M5
05234 R8=(R6+R5)/2.
05240 E=FNE(R8)
05242 D1=FNA(R8)-P4
05244 E=FNE(R5)
05246 IF D1*(FNA(R5)-P4)>0 THEN 5270
05250 R9=(R8+R5)/2. ' CROSSING BETWEEN R5&R8
05252 R6=R8
05260 GO TO 5280
05270 R9=(R8+R6)/2. ' CROSSING BETWEEN R8&R6
05272 R5=R8
05280 D1=(R9-R8)/R8
05289 ' 0.1% PRECISION DEMANDED
05290 IF ABS(D1)<0.001 THEN 5300
05292 R8=R9
05294 GO TO 5240
05300 R9=R8
05310 E=FNE(R6)
05312 D1=FNA(R6)-P4
05314 E=FNE(R5)
05316 IF D1*(FNA(R5)-P4)<0 THEN 5330
05320 R9=-4 ' FLAG FOR FAILURE TO FIND ROOT
05330 RETURN ' **********END INTERVAL-HALVING METHOD*****
05400 ' **********TO ASSIGN VALUES TO X,Y*****
05402 X(N)=R*C(FNC(K3))
05404 X(N+1)=R*C(FNC(-K3))
05406 Y(N)=R*S(FNC(K3))
05408 Y(N+1)=R*S(FNC(-K3))
05410 RETURN ' ****END SUBROUTINE*****
05500 ' **********TO PLOT ONE CONTOUR OF ONE LOBE AT A TIME*****
05502 L=1 ' PUT PLOTTER AT FIRST POINT
05504 X=(X(1)+T1)/S
05506 Y=(Y(1)+T2)/S
05510 GOSUB 1300
05522 I1=1
05524 I2=N
05526 I3=2
05540 FOR K=I1 TO I2 STEP I3
05550 X=(X(K)+T1)/S
05552 Y=(Y(K)+T2)/S
05554 GOSUB 1300
05560 NEXT K
05562 IF I3=-2 THEN 5580
05570 I1=N+1
05572 I2=0
05574 I3=-2
05576 GO TO 5540
05580 RETURN ' ******END SUBROUTINE*****
05600 ' **********TO SORT POINTS & ELIMINATE DISCONTINUITIES*****
05610 FOR I=1 TO I5
05612 ON I GO TO 5614, 5618, 5622
05614 N1=16
05616 GO TO 5630
05618 N1=13
05620 GO TO 5630
05622 N1=10
05630 FOR J=N1 TO 37
05632 IF R4(J,I)<=R3(J,I) THEN 5660
05634 IF R4(J,I)>R4(J-1,I) THEN 5660
05650 IF R3(J,I)>=R3(J-1,I) THEN 5670
05652 IF J+5>=37 THEN 5658
05654 J2=J+5
05656 GO TO 5660
05658 J2=J
05660 FOR K=J TO J2 ' ELIMINATE BAD ROOTS AND ALL SUBSEQUENT ONES
05662 LET R3(K,I)=R4(K,I)=-5
05664 NEXT K
05670 NEXT J
05680 NEXT I
05682 R1(0,0)=10 ' SET INDICATOR WHEN SORT IS DONE
05690 RETURN
09999 END
IIIAAAA>>AAA> AA@@@@IIIAFIII1 @@@@``III6~~FIII1 *U*1