HAM Radio 1/HamRadio.cdr/tech/design4/smploop.bas

BASIC Source File | 1979-12-31 | 15KB | 625 lines

10 CLS
20 DEFSNG A-H: DEFINT I-N: DEFDBL O-Z: PI=3.141592653589794#
30 DIM D%(20),PR(20),PJ(20),ZR(20),ZJ(20),TPR(20),TPJ(20),TZR(20),TZJ(20)
40 DIM C#(2048),P(20),CXPLOT(1026),CYPLOT(1026)
50 CLS
60 INPUT"SAMPLING FREQUENCY (IN HZ)";SAMPLEF
70 PRINT
80 INPUT"SAMPLING EFFICIENCY (0>ETA<1)";SAMPLEN
90 SSAMPLE=SAMPLEF*(LOG(1-SAMPLEN))
100 WSAMPLE=2*PI*SAMPLEF
110 CLS
120 PRINT"MENU"
130 PRINT
140 PRINT"MAKE PROGRAM CHOICE BELOW"
150 PRINT
160 PRINT"TYPE A TO LOAD FREQUENCY VARIABLES FROM FILE"
170 PRINT
180 PRINT"TYPE B TO LOAD FREQUENCY VARIABLES FROM KEYBOARD"
190 PRINT
200 PRINT"TYPE C TO EDIT FREQUENCY VARIABLE LIST"
210 PRINT
220 PRINT"TYPE D TO RUN POLYNOMIAL FILTER SYNTHESIS"
230 PRINT
240 PRINT"TYPE E TO RUN OPEN & CLOSED LOOP FREQUENCY RESPONSE"
250 PRINT
260 PRINT"TYPE F TO DETERMINE CLOSED LOOP NOISE BANDWIDTH"
270 PRINT
280 PRINT"TYPE G TO RUN CLOSED LOOP TRANSIENT RESPONSE"
290 PRINT
300 PRINT"TYPE Q TO QUIT"
310 PRINT
320 INPUT"PROGRAM CHOICE";A$
330 PRINT
340 IF A$="A" THEN 430
350 IF A$="B" THEN 450
360 IF A$="C" THEN 470
370 IF A$="D" THEN 490
380 IF A$="E" THEN 510
390 IF A$="F" THEN 530
400 IF A$="G" THEN 550
410 IF A$="Q" THEN 570
420 GOTO 110
430 GOSUB 1000
440 GOTO 110
450 GOSUB 5000
460 GOTO 110
470 GOSUB 10000
480 GOTO 110
490 GOSUB 15000
500 GOTO 570
510 GOSUB 20000
520 GOTO 570
530 GOSUB 25000
540 GOTO 570
550 GOSUB 30000
570 PRINT
580 INPUT"QUIT (Y/N)";B$
590 IF B$="Y" THEN 610
600 GOTO 50
610 PRINT
620 INPUT"DO YOU WANT TO SAVE THE FREQUENCY VARIABLES IN A FILE (Y/N)";C$
630 IF C$="N" THEN 780
640 PRINT
650 INPUT"DECLARE FILENAME";D$
660 OPEN D$ FOR OUTPUT AS #1
670 WRITE #1,NZ%
680 IF NZ%=0 THEN 730
690 FOR I=1 TO NZ%
700 WRITE #1,ZR(I)
710 WRITE #1,ZJ(I)
720 NEXT I
730 WRITE #1,NP%
740 FOR I=1 TO NP%
750 WRITE #1,PR(I)
760 WRITE #1,PJ(I)
770 NEXT I
780 END
1000 CLS
1010 INPUT"DECLARE FILENAME";F$
1020 OPEN F$ FOR INPUT AS #1
1030 INPUT #1,NZ%
1040 IF NZ%=0 THEN 1090
1050 FOR I=1 TO NZ%
1060 INPUT #1,ZR(I)
1070 INPUT #1,ZJ(I)
1080 NEXT I
1090 INPUT #1,NP%
1100 FOR I=1 TO NP%
1110 INPUT #1,PR(I)
1120 INPUT #1,PJ(I)
1130 NEXT I
1135 CLOSE #1
1270 RETURN
5000 CLS
5010 PRINT"INPUT # OF ZEROS";
5020 INPUT NZ%
5030 PRINT"INPUT # OF POLES";
5040 INPUT NP%
5050 PRINT
5060 IF NZ%=0 THEN 5150
5070 FOR I=1 TO NZ%
5080 PRINT"ZERO #";I;"REAL PART";
5090 INPUT ZR(I)
5100 PRINT"ZERO #";I;"IMAGINARY PART";
5110 INPUT ZJ(I)
5130 NEXT I
5140 PRINT
5150 FOR I=1 TO NP%
5160 PRINT"POLE #";I;"REAL PART";
5170 INPUT PR(I)
5180 PRINT"POLE #";I;"IMAGINARY PART";
5190 INPUT PJ(I)
5200 NEXT I
5210 RETURN
10000 CLS
10010 PRINT"THERE ARE";NZ%;"ZEROES AND";NP%;"POLES"
10020 PRINT
10030 IF NZ%=0 THEN 10090
10040 FOR I=1 TO NZ%
10050 PRINT"ZERO #";I;"REAL PART IS";ZR(I)
10060 PRINT"ZERO #";I;"IMAGINARY PART IS";ZJ(I)
10070 NEXT I
10080 PRINT
10090 FOR I=1 TO NP%
10100 PRINT"POLE #";I;"REAL PART IS";PR(I)
10110 PRINT"POLE #";I;"IMAGINARY PART IS";PJ(I)
10120 NEXT I
10130 PRINT
10140 PRINT"EDIT ZEROES (Y/N)";
10150 INPUT B$
10160 PRINT
10170 IF B$="N" THEN 10750
10180 PRINT"ADD,DELETE,OR SUBSTITUTE (A,D,S)";
10190 INPUT C$
10200 PRINT
10210 IF C$="S" THEN 10630
10220 IF C$="D" THEN 10380
10230 IF C$="A" THEN 10240 ELSE 10180
10240 PRINT"HOW MANY";
10250 INPUT NC%
10260 ND%=NZ%+NC%
10270 PRINT
10280 FOR I=1 TO ND%
10290 IF I>NZ% THEN 10310
10300 GOTO 10350
10310 PRINT"ZER0 #";I;"REAL PART";
10320 INPUT ZR(I)
10330 PRINT"ZERO #";I;"IMAGINARY PART";
10340 INPUT ZJ(I)
10350 NEXT I
10360 NZ%=NZ%+NC%
10370 GOTO 10750
10380 PRINT"HOW MANY";
10390 INPUT NC%
10400 PRINT
10410 FOR J=1 TO NC%
10420 INPUT"DELETE ZERO #";D%(J)
10430 NEXT J
10440 M=0
10450 FOR J=1 TO NC%
10460 FOR I=1 TO NZ%
10470 IF I=D%(J) THEN 10510
10480 M=M+1
10490 TZR(M)=ZR(I)
10500 TZJ(M)=ZJ(I)
10510 NEXT I
10520 NEXT J
10530 NZ%=NZ%-NC%
10560 FOR I=1 TO NZ%
10570 ZR(I)=TZR(I)
10580 ZJ(I)=TZJ(I)
10590 NEXT I
10620 GOTO 10750
10630 FOR I=1 TO NZ%
10640 PRINT"ZERO #";I;"IS";ZR(I);"AND J";ZJ(I)
10650 PRINT
10660 PRINT"CHANGE (Y/N)";
10670 INPUT D$
10675 PRINT
10680 IF D$="N" THEN 10740
10700 PRINT"ZERO #";I;"REAL PART";
10710 INPUT ZR(I)
10720 PRINT"ZERO #";I;"IMAGINARY PART";
10730 INPUT ZJ(I)
10740 NEXT I
10750 PRINT
10760 PRINT"EDIT POLES (Y/N)";
10770 INPUT A$
10780 IF A$="N" THEN 11380
10790 PRINT
10800 REM START POLE EDIT
10810 PRINT"ADD,DELETE,OR SUBSTITUTE POLES (A,D,S)";
10820 INPUT A$
10830 PRINT
10840 IF A$="S" THEN 11260
10850 IF A$="D" THEN 11010
10860 IF A$="A" THEN 10870 ELSE 10810
10870 PRINT"HOW MANY";
10880 INPUT NC%
10890 ND%=NP%+NC%
10900 PRINT
10910 FOR I=1 TO ND%
10920 IF I>NP% THEN 10940
10930 GOTO 10980
10940 PRINT"POLE #";I;"REAL PART";
10950 INPUT PR(I)
10960 PRINT"POLE #";I;"IMAGINARY PART";
10970 INPUT PJ(I)
10980 NEXT I
10990 NP%=NP%+NC%
11000 GOTO 11380
11010 PRINT"HOW MANY";
11020 INPUT NC%
11030 PRINT
11040 FOR J=1 TO NC%
11050 INPUT"DELETE POLE #";D%(J)
11060 NEXT J
11070 M=O
11080 FOR J=1 TO NC%
11090 FOR I=1 TO NP%
11100 IF I=D%(J) THEN 11140
11110 M=M+1
11120 TPR(M)=PR(I)
11130 TPJ(M)=PJ(I)
11140 NEXT I
11150 NEXT J
11160 NP%=NP%-NC%
11190 FOR I=1 TO NP%
11200 PR(I)=TPR(I)
11210 PJ(I)=TPJ(I)
11220 NEXT I
11250 GOTO 11380
11260 FOR I=1 TO NP%
11270 PRINT"POLE #";I;"IS";PR(I);"AND J";PJ(I)
11280 PRINT
11290 PRINT"CHANGE (Y/N)";
11300 INPUT C$
11305 PRINT
11310 IF C$="N" THEN 11370
11330 PRINT"POLE #";I;"REAL PART";
11340 INPUT PR(I)
11350 PRINT"POLE #";I;"IMAGINARY PART";
11360 INPUT PJ(I)
11370 NEXT I
11380 RETURN
15000 CLS
15010 PRINT"REQUIRED ATTENUATION IN dB";
15020 INPUT Z21
15040 PRINT
15050 PRINT"INPUT CUTTOFF FREQUENCY RATIO OF REQUIRED ATTENUATION";
15060 INPUT W
15070 IF W<=1# THEN 15050
15080 PRINT
15090 PRINT"PERMISSABLE PASSBAND RIPPLE IN dB , 0dB FOR BUTTERWORTH";
15100 INPUT R
15120 K#=10#^(Z21/10#)
15125 IF R=0 THEN 15440
15130 RES=(10#^(R/10#))-1#
15140 RE=SQR(RES)
15150 TW=SQR((K#-1#)/RES)
15160 FOR I=1 TO 20
15170 KA#=((W+SQR(W^2#-1#))^I+(W+SQR(W^2#-1#))^(1#/I))/2#
15180 IF KA#=>TW THEN 15200
15190 NEXT I
15200 PRINT
15210 PRINT"N=";I
15230 PRINT
15240 PRINT"RESULTS SATISFACTORY (Y/N)";
15250 INPUT A$
15260 IF A$="Y" THEN 15270 ELSE 15000
15270 X=(SQR((1#/RES)+1#)+1#/RE)^(1#/I)
15280 Y=1#/X
15290 SINH=(X-Y)/2#
15300 COSH#=(X+Y)/2#
15310 N%=I-1
15320 PRINT
15330 PRINT"FREQUENCY SCALE FACTOR";
15340 INPUT KF#
15350 PRINT
15370 PRINT
15380 FOR J=0 TO N%
15390 PR(J+1+NP%)=KF#*(-SINH*(SIN(((2#*J+1#)*PI)/(2#*I))))
15400 PJ(J+1+NP%)=KF#*(COSH#*(COS(((2#*J+1#)*PI)/(2#*I))))
15420 NEXT J
15430 GOTO 15610
15440 FOR I=1 TO 20
15450 KA#=1+W^(2#*I)
15460 IF KA#=>K# THEN 15480
15470 NEXT I
15480 PRINT
15490 PRINT"N=";I
15500 PRINT
15510 INPUT"RESULTS SATISFACTORY";A$
15520 IF A$="Y" THEN 15530 ELSE 15000
15530 PRINT
15540 INPUT"FREQUENCY SCALE FACTOR";KF#
15550 N%=I
15560 FOR J=1 TO N%
15570 P(J)=PI/2#+(PI/(2#*N%))*(2#*J-1#)
15580 PR(J+NP%)=KF#*(COS(P(J)))
15590 PJ(J+NP%)=KF#*(SIN(P(J)))
15600 NEXT J
15610 N=I
15620 J=0
15622 FOR L=1 TO N
15624 TPR(L)=PR(NP%+L)
15626 TPJ(L)=PJ(NP%+L)
15628 NEXT L
15630 FOR L=1 TO N
15640 J=J+1
15670 IF ABS(TPJ(L))=>ABS((10^-6)*TPR(L)) THEN 15680 ELSE 15750
15680 IF TPJ(L)<0 THEN 15770
15690 PR(NP%+J)=TPR(L)
15700 PJ(NP%+J)=TPJ(L)
15710 PR(NP%+J+1)=TPR(L)
15720 PJ(NP%+J+1)=-TPJ(L)
15730 J=J+1
15740 GOTO 15770
15750 PR(NP%+J)=TPR(L)
15760 PJ(NP%+J)=0#
15770 NEXT L
15780 NP%=NP%+I
15790 RETURN
20000 CLS
20010 KN#=1#
20020 IF NZ%=0 THEN 20080
20030 FOR I=1 TO NZ%
20040 ZT=ABS(ZR(I))+ABS(ZJ(I))
20050 IF ZT=0 THEN 20070
20060 KN#=KN#*SQR(ZR(I)^2+ZJ(I)^2)
20070 NEXT I
20080 KD#=ABS(SSAMPLE)
20090 FOR I=1 TO NP%
20100 PT=ABS(PR(I))+ABS(PJ(I))
20110 IF PT=0 THEN 20130
20120 KD#=KD#*SQR(PR(I)^2+PJ(I)^2)
20130 NEXT I
20140 PRINT
20150 PRINT"INPUT LOOP GAIN";
20160 INPUT U
20170 K#=(KD#/KN#)*U
20172 PRINT
20174 INPUT"PLOT (Y/N)";P$
20180 PRINT
20190 PRINT "START FREQ.";
20200 INPUT WS
20210 PRINT"NUMBER OF DECADES";
20220 INPUT WF
20222 IF P$="Y" THEN 20252
20230 PRINT"DECADE INCREMENT";
20240 INPUT WD
20250 N%=WF/WD
20251 GOTO 20260
20252 N%=500:WD=WF/500
20260 PRINT
20270 PRINT
20272 IF P$="Y" THEN 20290
20280 PRINT USING"\            \";"W","H(W)-dB","THETA(W)","H(W)'-dB","THETA(W)'"
20290 FOR I=0 TO N%
20300 W=WS*(10^(I*WD))
20310 PM=(1/(SQR(W^2+SSAMPLE^2)))*ABS((SIN(PI*(W/WSAMPLE)))/((PI*W)/WSAMPLE))
20320 PT=-PI*(W/WSAMPLE)+ATN(W/SSAMPLE)
20330 TD=0
20340 FOR J=1 TO NP%
20350 PM=PM/(SQR(PR(J)^2+(W-PJ(J))^2))
20360 IF PR(J)=0 THEN 20390
20370 PT=PT+ATN((W-PJ(J))/PR(J))
20380 GOTO 20400
20390 PT=PT-PI/2
20400 NEXT J
20410 IF NZ%=0 THEN 20490
20420 FOR H=1 TO NZ%
20430 PM=PM*SQR(ZR(H)^2+(W-ZJ(H))^2)
20440 IF ZR(H)=0 THEN 20470
20450 PT=PT-ATN((W-ZJ(H))/ZR(H))
20460 GOTO 20480
20470 PT=PT+PI/2
20480 NEXT H
20490 PM=PM*K#
20500 R=PM*COS(PT)
20510 S=PM*SIN(PT)
20520 R=R+1
20530 MC#=PM/(SQR(R*R+S*S))
20540 X=R/(SQR(R*R+S*S))
20550 IF 1-X*X=0 THEN 20570
20560 PC=PT-SGN(S)*(PI/2-ATN(X/SQR(1-X*X)))
20570 XM=8.68589*LOG(PM)
20580 XC=8.68589*LOG(MC#)
20582 IF P$="Y" THEN 20592
20590 PRINT USING"##.###^^^^    ";W,XM,PT,XC,PC
20591 GOTO 20600
20592 CXPLOT(I)=W:CYPLOT(I)=XC
20600 NEXT I
20601 IF P$="Y" THEN 20602 ELSE 20610
20602 CLS
20603 NCURVES!=1
20604 NPOINTS!=500
20605 CALL PLOT((NCURVES!),(NPOINTS!),CXPLOT(),CYPLOT())
20610 PRINT
20620 PRINT"NEW GAIN VALUE (Y/N)";
20630 INPUT A$
20640 PRINT
20650 IF A$="Y" THEN 20150
20660 RETURN
25000 CLS
25010 KN#=1#
25020 IF NZ%=0 THEN 25080
25030 FOR I=1 TO NZ%
25040 ZT=ABS(ZR(I))+ABS(ZJ(I))
25050 IF ZT=0 THEN 25070
25060 KN#=KN#*SQR(ZR(I)^2+ZJ(I)^2)
25070 NEXT I
25080 KD#=ABS(SSAMPLE)
25090 FOR I=1 TO NP%
25100 PT=ABS(PR(I))+ABS(PJ(I))
25110 IF PT=0 THEN 25130
25120 KD#=KD#*SQR(PR(I)^2+PJ(I)^2)
25130 NEXT I
25140 PRINT
25150 PRINT"INPUT LOOP GAIN";
25160 INPUT U
25170 K#=(KD#/KN#)*U
25180 PRINT
25190 PRINT"INPUT APPROXIMATE -3dB RADIAN FREQUENCY";
25200 INPUT WC
25210 PRINT
25220 XN=0
25230 FOR A%=1 TO 3
25240 IF A%=1 THEN 25270
25250 WD=WC*(10^(A%-3))
25260 GOTO 25540
25270 FOR B%=1 TO 100
25280 WD=WC/100#
25290 W=WD*B%
25300 PM=(1/(SQR(W^2+SSAMPLE^2)))*ABS((SIN(PI*(W/WSAMPLE)))/((PI*W)/WSAMPLE))
25310 PT=ATN(W/SSAMPLE)-PI*(W/WSAMPLE)
25320 FOR J=1 TO NP%
25330 PM=PM/(SQR(PR(J)^2+(W-PJ(J))^2))
25340 IF PR(J)=0 THEN 25370
25350 PT=PT+ATN((W-PJ(J))/PR(J))
25360 GOTO 25380
25370 PT=PT-PI/2
25380 NEXT J
25390 IF NZ%=0 THEN 25470
25400 FOR H=1 TO NZ%
25410 PM=PM*SQR(ZR(H)^2+(W-ZJ(H))^2)
25420 IF ZR(H)=0 THEN 25450
25430 PT=PT-ATN((W-ZJ(H))/ZR(H))
25440 GOTO 25460
25450 PT=PT+PI/2
25460 NEXT H
25470 PM=PM*K#
25472 IF PM=>1# AND ABS(PT)>PI THEN 25474 ELSE 25480
25474 PRINT
25476 PRINT"THE SYSTEM IS UNSTABLE"
25478 GOTO 25810
25480 R=PM*COS(PT)
25490 S=PM*SIN(PT)
25500 R=R+1
25510 XN=XN+((PM*PM)/(R*R+S*S))*WD
25520 NEXT B%
25530 IF A%=1 THEN 25790
25540 FOR C%=1 TO 90
25550 W=WD*(10+C%)
25560 PM=(1/(SQR(W^2+SSAMPLE^2)))*ABS((SIN(PI*(W/WSAMPLE)))/((PI*W)/WSAMPLE))
25570 PT=ATN(W/SSAMPLE)-PI*(W/WSAMPLE)
25580 FOR J=1 TO NP%
25590 PM=PM/(SQR(PR(J)^2+(W-PJ(J))^2))
25600 IF PR(J)=0 THEN 25630
25610 PT=PT+ATN((W-PJ(J))/PR(J))
25620 GOTO 25640
25630 PT=PT-PI/2
25640 NEXT J
25650 IF NZ%=0 THEN 25730
25660 FOR H=1 TO NZ%
25670 PM=PM*SQR(ZR(H)^2+(W-ZJ(H))^2)
25680 IF ZR(H)=0 THEN 25710
25690 PT=PT-ATN((W-ZJ(H))/ZR(H))
25700 GOTO 25720
25710 PT=PT+PI/2
25720 NEXT H
25730 PM=PM*K#
25740 R=PM*COS(PT)
25750 S=PM*SIN(PT)
25760 R=R+1
25770 XN=XN+((PM*PM)/(R*R+S*S))*WD
25780 NEXT C%
25790 NEXT A%
25800 PRINT "THE LOOP NOISE BANDWIDTH IS";XN;"RAD/SEC"
25810 PRINT
25820 PRINT"NEW LOOP GAIN (Y/N)";
25830 INPUT A$
25840 PRINT
25850 IF A$="N" THEN 25860 ELSE 25150
25860 RETURN
30000 CLS
30010 KN#=1
30020 IF NZ%=0 THEN 30080
30030 FOR I=1 TO NZ%
30040 ZT=ABS(ZR(I))+ABS(ZJ(I))
30050 IF ZT=0 THEN 30070
30060 KN#=KN#*SQR(ZR(I)^2+ZJ(I)^2)
30070 NEXT I
30080 KD#=ABS(SSAMPLE)
30090 FOR I=1 TO NP%
30100 PT=ABS(PR(I))+ABS(PJ(I))
30110 IF PT=0 THEN 30130
30120 KD#=KD#*SQR(PR(I)^2+PJ(I)^2)
30130 NEXT I
30140 PRINT
30150 PRINT"INPUT LOOP GAIN";
30160 INPUT U
30170 K#=(KD#/KN#)*U
30180 PRINT
30190 INPUT"AUTOSCALE TIME BASE (Y/N)";Q$
30200 IF Q$="N" THEN 30590
30210 NO%=0
30220 FOR I=1 TO NP%
30230 IF PR(I)=0 AND PJ(I)=0 THEN 30240 ELSE 30250
30240 NO%=NO%+1
30250 NEXT I
30260 IF NO%=0 THEN 30570
30270 WS=U^(1/NO%)
30280 WL=WS/1024
30290 WU=0
30300 FOR I=20 TO 0 STEP -1
30310 W=WL*(2#^I)+WU
30320 PM=(1/(SQR(W^2+SSAMPLE^2)))*ABS((SIN(PI*(W/WSAMPLE)))/((PI*W)/WSAMPLE))
30330 PT=ATN(W/SSAMPLE)-PI*(W/WSAMPLE)
30340 TD=0
30350 FOR J=1 TO NP%
30360 PM=PM/(SQR(PR(J)^2+(W-PJ(J))^2))
30370 IF PR(J)=0 THEN 30400
30380 PT=PT+ATN((W-PJ(J))/PR(J))
30390 GOTO 30410
30400 PT=PT-PI/2
30410 NEXT J
30420 IF NZ%=0 THEN 30500
30430 FOR H=1 TO NZ%
30440 PM=PM*SQR(ZR(H)^2+(W-ZJ(H))^2)
30450 IF ZR(H)=0 THEN 30480
30460 PT=PT-ATN((W-ZJ(H))/ZR(H))
30470 GOTO 30490
30480 PT=PT+PI/2
30490 NEXT H
30500 PM=PM*K#
30510 IF PM<1# THEN 30530
30520 WU=WU+W
30530 NEXT I
30540 WS=.025#*PI*WU
30550 T=(2#*PI)/WS
30560 GOTO 30620
30570 PRINT
30580 PRINT"CANNOT AUTOSCALE TIME BASE"
30590 PRINT
30600 INPUT"FULL SCALE TIME BASE (IN SECONDS)";T
30610 WS=(2#*PI)/T
30620 PRINT
30630 INPUT"STEP OR IMPULSE RESPONSE (S/I)";A$
30640 PRINT
30650 C#(0)=1#
30660 C#(1)=0#
30670 FOR I=1 TO 1023
30680 W=WS*I
30690 PM=(1/(SQR(W^2+SSAMPLE^2)))*ABS((SIN(PI*(W/WSAMPLE)))/((PI*W)/WSAMPLE))
30700 PT=ATN(W/SSAMPLE)-PI*(W/WSAMPLE)
30710 TD=0
30720 FOR J=1 TO NP%
30730 PM=PM/(SQR(PR(J)^2+(W-PJ(J))^2))
30740 IF PR(J)=0 THEN 30770
30750 PT=PT+ATN((W-PJ(J))/PR(J))
30760 GOTO 30780
30770 PT=PT-PI/2
30780 NEXT J
30790 IF NZ%=0 THEN 30870
30800 FOR H=1 TO NZ%
30810 PM=PM*SQR(ZR(H)^2+(W-ZJ(H))^2)
30820 IF ZR(H)=0 THEN 30850
30830 PT=PT-ATN((W-ZJ(H))/ZR(H))
30840 GOTO 30860
30850 PT=PT+PI/2
30860 NEXT H
30870 PM=PM*K#
30880 IF PM =>1# AND ABS(PT)>PI THEN 30890 ELSE 30920
30890 PRINT
30900 PRINT"THE SYSTEM IS UNSTABLE"
30910 GOTO 31230
30920 R=PM*COS(PT)
30930 S=PM*SIN(PT)
30940 R=R+1
30950 MC#=PM/(SQR(R*R+S*S))
30955 MC#=MC#*(1/(1+(W/(WS*316))^6))
30960 X=R/(SQR(R*R+S*S))
30970 IF 1-X*X=0 THEN 30990
30980 PC=PT-SGN(S)*(PI/2-ATN(X/SQR(1-X*X)))
30990 PC=PC-(80*PI*I)/1023
31000 C#(2*I)=MC#*COS(PC)
31010 C#(2*I+1)=MC#*SIN(PC)
31020 NEXT I
31030 PRINT
31040 N%=10
31050 NORM%=2
31060 SCALE=1#
31070 CALL RIFT(C#(0),N%,NORM%,SCALE)
31080 PRINT
31090 PRINT
31100 IF A$="I" THEN 31140
31110 FOR I=1 TO 1023
31120 C#(I)=C#(I)+C#(I-1)-C#(0)
31130 NEXT I
31132 FOR I=1 TO 1023
31134 C#(I)=C#(I)/C#(1023)
31136 NEXT I
31140 FOR I=1 TO 1024
31150 CYPLOT(I)=C#(I-1):CXPLOT(I)=(I-1)*(T/1023)
31160 NEXT I
31162 CLS
31170 NCURVES!=1
31180 NPOINTS!=1024
31190 CALL PLOT((NCURVES!),(NPOINTS!),CXPLOT(),CYPLOT())
31200 INPUT"NEW LOOP GAIN";A$
31210 IF A$="N" THEN 31230
31220 GOTO 30000
31230 RETURN