home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
HAM Radio 1
/
HamRadio.cdr
/
tech
/
eepub07
/
phaslok.bas
< prev
next >
Wrap
BASIC Source File
|
1986-09-15
|
21KB
|
895 lines
1 DEFSNG A-H: DEFINT I-L,N: DEFDBL M,O-Z: PI=3.141592653589794#
2 DIM D%(40),PR(40),PJ(40),ZR(40),ZJ(40),TPR(40),TPJ(40),TZR(40),TZJ(40)
3 DIM P(40),XR(40),XJ(40),ZC(40),PC(40)
4 DIM MAT1(40),MAT2(40),MAT3(40),ZM(40),ZT(40),PM(40),PT(40),RP(40),RM(40)
5 DIM RA(40),RT(40),RES(40),PHI(40),ROOT(80)
6 DIM XPASS!(2010),YPASS!(2010):ITERATION=0:NPP=1000:NCURVES=1
7 GOTO 19
8 CLS:PRINT"THIS PROGRAM WILL GENERATE MULTIPLE PLOTS"
9 INPUT"HOW MANY WILL YOU NEED";NCURVES
10 INPUT"HOW MANY POINTS PER PLOT";NPP
11 IF NPP*NCURVES>2000 THEN 9
12 ITERATION=0
13 GOTO 19
19 CLS
20 PRINT"MENU"
30 PRINT
40 PRINT"MAKE PROGRAM CHOICE BELOW"
50 PRINT
60 PRINT"TYPE A TO LOAD FREQUENCY VARIABLES FROM FILE"
70 PRINT
80 PRINT"TYPE B TO LOAD FREQUENCY VARIABLES FROM KEYBOARD"
90 PRINT
100 PRINT"TYPE C TO EDIT FREQUENCY VARIABLE LIST"
110 PRINT
120 PRINT"TYPE D TO RUN POLYNOMIAL FILTER SYNTHESIS"
130 PRINT
140 PRINT"TYPE E TO RUN OPEN & CLOSED LOOP FREQUENCY RESPONSE"
150 PRINT
160 PRINT"TYPE F TO DETERMINE CLOSED LOOP NOISE BANDWIDTH"
170 PRINT
180 PRINT"TYPE G TO COMPUTE CLOSED LOOP ROOT LOCUS"
181 PRINT
182 PRINT"TYPE H TO RUN CLOSED LOOP TRANSIENT RESPONSE"
183 PRINT
184 PRINT"TYPE PLOT FOR MULTIPLE PLOTTING PARAMETERS"
190 PRINT
192 PRINT"TYPE Q TO QUIT"
194 PRINT
200 INPUT"PROGRAM CHOICE";A$
210 PRINT
220 IF A$="A" THEN 300
230 IF A$="B" THEN 320
240 IF A$="C" THEN 340
250 IF A$="D" THEN 360
260 IF A$="E" THEN 380
270 IF A$="F" THEN 400
280 IF A$="G" THEN 420
281 IF A$="H" THEN 426
282 IF A$="Q" THEN 430
283 IF A$="PLOT" THEN 8
290 GOTO 19
300 GOSUB 1000
310 GOTO 19
320 GOSUB 5000
330 GOTO 19
340 GOSUB 10000
350 GOTO 19
360 GOSUB 15000
370 GOTO 19
380 GOSUB 20000
390 GOTO 19
400 GOSUB 25000
410 GOTO 19
420 GOSUB 30000
425 GOTO 19
426 GOSUB 50000
427 GOTO 19
430 PRINT
435 INPUT"QUIT (Y/N)";B$
440 IF B$="Y" THEN 460
450 GOTO 19
460 PRINT
470 INPUT"DO YOU WANT TO SAVE THE FREQUENCY VARIABLES IN A FILE (Y/N)";C$
480 IF C$="N" THEN 610
490 PRINT
500 INPUT"DECLARE FILENAME";D$
510 OPEN D$ FOR OUTPUT AS #1
515 WRITE #1,NZ%
520 IF NZ%=0 THEN 565
530 FOR I=1 TO NZ%
540 WRITE #1,ZR(I)
550 WRITE #1,ZJ(I)
560 NEXT I
565 WRITE #1,NP%
570 FOR I=1 TO NP%
580 WRITE #1,PR(I)
590 WRITE #1,PJ(I)
600 NEXT I
610 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#=1#
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
20320 PT=0
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 XPASS!(I)=W:YPASS!(I)=XM
20600 NEXT I
20601 IF P$="Y" THEN 20602 ELSE 20610
20602 CLS
20603 NCURVES!=1
20604 NPOINTS!=500
20605 CALL PLOT((NCURVES!),(NPOINTS!),XPASS!(),YPASS!())
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#=1#
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
25310 PT=0
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
25570 PT=0
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 FOR I=1 TO NZ%
30020 XR(I)=ZR(I)
30030 XJ(I)=ZJ(I)
30040 NEXT I
30050 N=NZ%:NX%=NZ%:I=NZ%
30060 GOSUB 40000
30070 FOR I=1 TO NZ%
30080 ZR(I)=XR(I)
30090 ZJ(I)=XJ(I)
30100 NEXT I
30110 FOR I=1 TO NP%
30120 XR(I)=PR(I)
30130 XJ(I)=PJ(I)
30140 NEXT I
30150 N=NP%:NX%=NP%:I=NP%
30160 GOSUB 40000
30170 FOR I=1 TO NP%
30180 PR(I)=XR(I)
30190 PJ(I)=XJ(I)
30200 NEXT I
30210 GOSUB 42000
30220 GOSUB 48000
30221 IF ITAG=4 THEN 30239
30222 IF ITAG=1 THEN 30225
30223 IF ITAG=2 THEN 30228
30224 IF ITAG=3 THEN 30231
30225 PRINT
30226 PRINT"DIVIDE BY ZERO IN POLYNOMIAL ROOT SUBROUTINE"
30227 GOTO 30233
30228 PRINT
30229 PRINT"POLYNOMIAL ROOT SUBROUTINE IS DIVERGING"
30230 GOTO 30233
30231 PRINT
30232 PRINT"POLYNOMIAL ROOT SUBROUTINE ITERATION LIMIT EXCEEDED"
30233 PRINT
30234 PRINT"PRESS ANY KEY WHEN READY"
30235 A$=INKEY$:IF A$="" THEN 30235
30236 GOTO 30280
30239 J=0
30240 FOR I=1 TO 2*NP% STEP 2
30242 J=J+1
30250 PR(J)=ROOT(I):PJ(J)=ROOT(I+1)
30260 NEXT I
30280 RETURN
40000 J=0
40010 FOR L=1 TO N
40020 TPR(L)=XR(NX%+L)
40030 TPJ(L)=XJ(NX%+L)
40040 NEXT L
40050 FOR L=1 TO N
40060 J=J+1
40070 IF ABS(TPJ(L))=>ABS((10^-6)*TPR(L)) THEN 40080 ELSE 40150
40080 IF TPJ(L)<0 THEN 40170
40090 XR(NX%+J)=TPR(L)
40100 XJ(NX%+J)=TPJ(L)
40110 XR(NX%+J+1)=TPR(L)
40120 XJ(NX%+J+1)=-TPJ(L)
40130 J=J+1
40140 GOTO 40170
40150 XR(NX%+J)=TPR(L)
40160 XJ(NX%+J)=0#
40170 NEXT L
40180 NX%=NX%+I
40190 RETURN
42000 KN#=1#
42010 IF NZ%=0 THEN 42070
42020 FOR I=1 TO NZ%
42030 ZT=ABS(ZR(I))+ABS(ZJ(I))
42040 IF ZT=0 THEN 42060
42050 KN#=KN#*SQR(ZR(I)^2+ZJ(I)^2)
42060 NEXT I
42070 KD#=1#
42080 FOR I=1 TO NP%
42090 PT=ABS(PR(I))+ABS(PJ(I))
42100 IF PT=0 THEN 42120
42110 KD#=KD#*SQR(PR(I)^2+PJ(I)^2)
42120 NEXT I
42130 CLS
42140 INPUT"LOOP GAIN";K#
42150 K#=(KD#/KN#)*K#
42160 IF NZ%=0 THEN 42290
42170 IF NZ%=1 THEN 42310
42180 IF NZ%=2 THEN 42340
42190 NX%=NZ%
42200 FOR I=1 TO NX%
42210 XR(I)=ZR(I)
42220 XJ(I)=ZJ(I)
42230 NEXT I
42240 GOSUB 44000
42250 FOR I=1 TO NX%+1
42260 ZC(I)=K#*MAT3(I)
42270 NEXT I
42280 GOTO 42370
42290 ZC(1)=K#
42300 GOTO 42370
42310 ZC(1)=-ZR(1)*K#
42320 ZC(2)=K#
42330 GOTO 42370
42340 ZC(1)=K#*(ZR(1)*ZR(2)+ABS(ZJ(1)*ZJ(2)))
42350 ZC(2)=-(ZR(1)+ZR(2))*K#
42360 ZC(3)=K#
42370 IF NP%=1 THEN 42490
42380 IF NP%=2 THEN 42520
42390 NX%=NP%
42400 FOR I=1 TO NX%
42410 XR(I)=PR(I)
42420 XJ(I)=PJ(I)
42430 NEXT I
42440 GOSUB 44000
42450 FOR I=1 TO NX%+1
42460 PC(I)=MAT3(I)
42470 NEXT I
42480 GOTO 42550
42490 PC(1)=-PR(I)
42500 PC(2)=1#
42510 GOTO 42550
42520 PC(1)=PR(1)*PR(2)+ABS(PJ(1)*PJ(2))
42530 PC(2)=-(PR(1)+PR(2))
42540 PC(3)=1#
42550 FOR I=1 TO NZ%+1
42560 PC(I)=PC(I)+ZC(I)
42570 NEXT I
42580 PRINT
42590 PRINT"THE NUMERATOR COEFFICIENTS ARE :"
42600 PRINT
42610 FOR I=0 TO NZ%
42620 PRINT"S";I;" = ";ZC(I+1)
42630 NEXT I
42640 PRINT
42650 PRINT"THE DENOMINATOR COEFFICIENTS ARE :"
42660 PRINT
42670 FOR I=0 TO NP%
42680 PRINT"S";I;" = ";PC(I+1)
42690 NEXT I
42700 PRINT
42710 PRINT"PRESS ANY KEY WHEN READY"
42720 A$=INKEY$:IF A$="" THEN 42720
42730 RETURN
44000 I=1
44010 IF XJ(I)=0 THEN 44020 ELSE 44070
44020 MAT1(1)=-XR(I)
44030 MAT1(2)=1#
44040 NDA=1
44050 I=I+1
44060 GOTO 44120
44070 MAT1(1)=XR(I)^2#+XJ(I)^2#
44080 MAT1(2)=-2#*XR(I)
44090 MAT1(3)=1#
44100 NDA=2
44110 I=I+2
44120 IF XJ(I)=0 THEN 44130 ELSE 44180
44130 MAT2(1)=-XR(I)
44140 MAT2(2)=1#
44150 NDB=1
44160 I=I+1
44170 GOTO 44230
44180 MAT2(1)=XR(I)^2#+XJ(I)^2#
44190 MAT2(2)=-2#*XR(I)
44200 MAT2(3)=1#
44210 NDB=2
44220 I=I+2
44230 GOSUB 46000
44240 FOR J=1 TO NDC+1
44250 MAT1(J)=MAT3(J)
44260 NEXT J
44270 NDA=NDC
44280 IF NDC=NX% THEN 44300
44290 GOTO 44120
44300 RETURN
46000 NDC = NDA + NDB
46010 QNDAP1% = NDA + 1
46020 QNDCP1% = NDC + 1
46030 FOR QI% = 1 TO QNDCP1%
46040 MAT3(QI%) = 0#
46050 NEXT QI%
46060 QNDBP1% = NDB + 1
46070 FOR QI% = 1 TO QNDAP1%
46080 FOR QJ% = 1 TO QNDBP1%
46090 QIPJ% = QI% + QJ% - 1
46100 MAT3(QIPJ%) = MAT3(QIPJ%) + MAT1(QI%)*MAT2(QJ%)
46110 NEXT QJ%
46120 NEXT QI%
46130 RETURN
48000 FOR I=NP% TO 0 STEP -1
48010 MAT1(NP%-I+1)=PC(I+1)
48020 NEXT I
48030 N=NP%
48031 SF=(MAT1(N+1))^(1/N)
48040 P1=1
48050 Q1=1
48060 ITAG=1000
48070 EPS#=.0001#
48080 QIITAG = ITAG
48090 ITAG = 0
48100 QD = MAT1(1)
48110 FOR QI% = 1 TO N
48120 MAT1(QI%) = MAT1(QI% + 1)/QD
48130 NEXT QI%
48140 IF N > 1 THEN 48220
48150 ITAG = ITAG + 1
48160 QL = ITAG + ITAG
48170 ROOT(QL) = 0!
48180 ROOT(QL-1) = -MAT1(1)
48190 N = ITAG
48200 ITAG = 4
48210 RETURN
48220 IF N > 2 THEN 48260
48230 QP = MAT1(1)
48240 QQ = MAT1(2)
48250 GOTO 48680
48260 QP = P1
48270 QQ = Q1
48280 QM = 1
48290 MAT2(1) = MAT1(1) - QP
48300 MAT2(2) = MAT1(2) - QP*MAT2(1) - QQ
48310 QL = N - 1
48320 MAT3(1) = MAT2(1) - QP
48330 MAT3(2) = MAT2(2) - QP*MAT3(1) - QQ
48340 IF QL = 2 THEN QL = 3
48350 FOR QJ% = 3 TO QL
48360 MAT2(QJ%) = MAT1(QJ%) - QP*MAT2(QJ%-1) - QQ*MAT2(QJ%-2)
48370 MAT3(QJ%) = MAT2(QJ%) - QP*MAT3(QJ%-1) - QQ*MAT3(QJ%-2)
48380 NEXT QJ%
48390 QL = N - 1
48400 QCBARL = MAT3(QL) - MAT2(QL)
48410 QDEN = -QCBARL
48420 IF N <> 3 THEN QDEN = QDEN*MAT3(N-3)
48430 QDEN = QDEN + MAT3(N-2)*MAT3(N-2)
48440 IF QDEN <> 0 THEN 48480
48450 N = ITAG
48460 ITAG = 1
48470 RETURN
48480 MAT2(N) = MAT1(N) - QP*MAT2(N-1) - QQ*MAT2(N-2)
48490 QDELTP = -MAT2(N)
48500 IF N <> 3 THEN QDELTP = QDELTP*MAT3(N-3)
48510 QDELTP = (MAT2(N-1)*MAT3(N-2) + QDELTP)/QDEN
48520 QDELTQ = (MAT2(N)*MAT3(N-2) - MAT2(N-1)*QCBARL)/QDEN
48530 QP = QP + QDELTP
48540 QQ = QQ + QDELTQ
48550 QSUM = ABS(QDELTP) + ABS(QDELTQ)
48560 IF QM = 1 THEN QSUM1 = QSUM
48570 IF QM <> 5 OR QSUM <= QSUM1 THEN 48610
48580 N = ITAG
48590 ITAG = 2
48600 RETURN
48610 IF QSUM <= EPS# THEN 48680
48620 IF QM < QIITAG THEN 48660
48630 N = ITAG
48640 ITAG = 3
48650 RETURN
48660 QM = QM + 1
48670 GOTO 48290
48680 QD = -QP*.5
48690 ITAG = ITAG + 2
48700 QF = QQ - QD*QD
48710 QE = SQR(ABS(QF))
48720 QL = ITAG + ITAG
48730 IF QF > 0! THEN 48790
48740 ROOT(QL) = 0!
48750 ROOT(QL-1) = QD-QE
48760 ROOT(QL-2) = 0!
48770 ROOT(QL-3) = QD + QE
48780 GOTO 48830
48790 ROOT(QL) = -QE
48800 ROOT(QL-1) = QD
48810 ROOT(QL-2) = QE
48820 ROOT(QL-3) = QD
48830 N = N-2
48840 IF N > 0 THEN 48880
48850 N = ITAG
48860 ITAG = 4
48870 RETURN
48880 FOR QI% = 1 TO N
48890 MAT1(QI%) = MAT2(QI%)
48900 NEXT QI%
48910 GOTO 48140
48920 RETURN
50000 CLS
50010 KN#=1#
50020 IF NZ%=0 THEN 50080
50030 FOR I=1 TO NZ%
50040 ZT=ABS(ZR(I))+ABS(ZJ(I))
50050 IF ZT=0 THEN 50070
50060 KN#=KN#*SQR(ZR(I)^2#+ZJ(I)^2#)
50070 NEXT I
50080 KD#=1#
50090 FOR I=1 TO NP%
50100 PT=ABS(PR(I))+ABS(PJ(I))
50110 IF PT=0 THEN 50130
50120 KD#=KD#*SQR(PR(I)^2#+PJ(I)^2#)
50130 NEXT I
50140 K#=KD#/KN#
50150 PRINT
50160 FOR I=1 TO NP%
50170 P=1#
50180 T=0
50190 FOR J=1 TO NP%
50200 PA=PR(I)-PR(J)
50210 PB=PJ(I)-PJ(J)
50220 PD=ABS(PA)+ABS(PB)
50230 IF PD=0 THEN 50330
50240 IF PB=0 THEN 50290
50250 PM=SQR(PA^2#+PB^2#)
50260 X=PA/PM
50270 PT=(PI/2#-ATN(X/SQR(1-X*X)))*SGN(PB)
50280 GOTO 50310
50290 PM=ABS(PA)
50300 PT=PI/2#+PI/2#*SGN(-PA)
50310 T=T+PT
50320 P=P*PM
50330 NEXT J
50340 RM(I)=P
50350 RP(I)=T
50360 NEXT I
50370 PRINT
50380 IF NZ%=0 THEN 50620
50390 FOR I=1 TO NP%
50400 P=1#
50410 T=0
50420 FOR J=1 TO NZ%
50430 ZA=PR(I)-ZR(J)
50440 ZB=PJ(I)-ZJ(J)
50450 ZD=ABS(ZA)+ABS(ZB)
50460 IF ZD=0 THEN 50560
50470 IF ZB=0 THEN 50520
50480 ZM=SQR(ZA^2#+ZB^2#)
50490 X=ZA/ZM
50500 ZT=(PI/2#-ATN(X/SQR(1-X*X)))*SGN(ZB)
50510 GOTO 50540
50520 ZM=ABS(ZA)
50530 ZT=PI/2#+PI/2#*SGN(-ZA)
50540 T=T+ZT
50550 P=P*ZM
50560 NEXT J
50570 RT(I)=T
50580 RA(I)=P
50590 NEXT I
50600 PRINT
50610 GOTO 50700
50620 PRINT USING"\ \";"POLE REAL";"POLE IMAG.";"RESIDUE MAG.";"RESIDUE ARG."
50630 PRINT
50640 FOR I=1 TO NP%
50650 RES(I)=K#/RM(I)
50660 PHI(I)=-RP(I)
50670 PRINT USING"##.####^^^^ ";PR(I),PJ(I),RES(I),PHI(I)
50680 NEXT I
50690 GOTO 50760
50700 PRINT USING"\ \";"POLE REAL";"POLE IMAG.";"RESIDUE MAG.";"RESIDUE ARG."
50710 FOR I=1 TO NP%
50720 RES(I)=K#*(RA(I)/RM(I))
50730 PHI(I)=RT(I)-RP(I)
50740 PRINT USING"##.####^^^^ ";PR(I),PJ(I),RES(I),PHI(I)
50750 NEXT I
50760 PRINT
50770 PRINT"PRESS ANY KEY TO CONTINUE"
50780 A$=INKEY$: IF A$="" THEN 50780
50782 IF NCURVES>1 AND ITERATION>0 THEN 50860
50790 PRINT
50800 PRINT"INPUT START TIME";
50810 INPUT TS
50820 PRINT"STOP TIME";
50830 INPUT TF
50850 TD=(TF-TS)/(NPP-1)
50860 CLS
50862 PRINT"COMPUTING TIME RESPONSE OF CURVE #";ITERATION+1
50890 FOR J= 0 TO (NPP-1)
50900 T=TS+J*TD
50910 RT=0
50920 FOR I=1 TO NP%
50930 IF PJ(I)=0 THEN 50960
50940 RR=2*RES(I)*EXP(PR(I)*T)*COS(PJ(I)*T+PHI(I))
50950 GOTO 50980
50960 RR=RES(I)*EXP(PR(I)*T)*COS(PHI(I))
50970 GOTO 50990
50980 I=I+1
50990 RT=RT+RR
51000 NEXT I
51010 XPASS!(J+1+NPP*ITERATION)=T:YPASS!(J+1+NPP*ITERATION)=RT
51020 NEXT J
51022 ITERATION=ITERATION+1
51024 IF NCURVES>ITERATION THEN 51050
51030 NCURVES!=NCURVES:NPP!=NPP
51040 CALL PLOT((NCURVES!),(NPP!),XPASS!(),YPASS!())
51047 NCURVES=1:ITERATION=0:NPP=1000:REM REINITIALIZE
51050 PRINT
51060 RETURN