home *** CD-ROM | disk | FTP | other *** search
/ HAM Radio 1 / HamRadio.cdr / tech / eepub07 / phaslok.bas < prev    next >
Encoding:
BASIC Source File  |  1986-09-16  |  20.7 KB  |  895 lines
  1. 1 DEFSNG A-H: DEFINT I-L,N: DEFDBL M,O-Z: PI=3.141592653589794#
  2. 2 DIM D%(40),PR(40),PJ(40),ZR(40),ZJ(40),TPR(40),TPJ(40),TZR(40),TZJ(40)
  3. 3 DIM P(40),XR(40),XJ(40),ZC(40),PC(40)
  4. 4 DIM MAT1(40),MAT2(40),MAT3(40),ZM(40),ZT(40),PM(40),PT(40),RP(40),RM(40)
  5. 5 DIM RA(40),RT(40),RES(40),PHI(40),ROOT(80)
  6. 6 DIM XPASS!(2010),YPASS!(2010):ITERATION=0:NPP=1000:NCURVES=1
  7. 7 GOTO 19
  8. 8 CLS:PRINT"THIS PROGRAM WILL GENERATE MULTIPLE PLOTS"
  9. 9 INPUT"HOW MANY WILL YOU NEED";NCURVES
  10. 10 INPUT"HOW MANY POINTS PER PLOT";NPP
  11. 11 IF NPP*NCURVES>2000 THEN 9
  12. 12 ITERATION=0
  13. 13 GOTO 19
  14. 19 CLS
  15. 20 PRINT"MENU"
  16. 30 PRINT
  17. 40 PRINT"MAKE PROGRAM CHOICE BELOW"
  18. 50 PRINT
  19. 60 PRINT"TYPE A TO LOAD FREQUENCY VARIABLES FROM FILE"
  20. 70 PRINT
  21. 80 PRINT"TYPE B TO LOAD FREQUENCY VARIABLES FROM KEYBOARD"
  22. 90 PRINT
  23. 100 PRINT"TYPE C TO EDIT FREQUENCY VARIABLE LIST"
  24. 110 PRINT
  25. 120 PRINT"TYPE D TO RUN POLYNOMIAL FILTER SYNTHESIS"
  26. 130 PRINT
  27. 140 PRINT"TYPE E TO RUN OPEN & CLOSED LOOP FREQUENCY RESPONSE"
  28. 150 PRINT
  29. 160 PRINT"TYPE F TO DETERMINE CLOSED LOOP NOISE BANDWIDTH"
  30. 170 PRINT
  31. 180 PRINT"TYPE G TO COMPUTE CLOSED LOOP ROOT LOCUS"
  32. 181 PRINT
  33. 182 PRINT"TYPE H TO RUN CLOSED LOOP TRANSIENT RESPONSE"
  34. 183 PRINT
  35. 184 PRINT"TYPE  PLOT  FOR MULTIPLE PLOTTING PARAMETERS"
  36. 190 PRINT
  37. 192 PRINT"TYPE Q TO QUIT"
  38. 194 PRINT
  39. 200 INPUT"PROGRAM CHOICE";A$
  40. 210 PRINT
  41. 220 IF A$="A" THEN 300
  42. 230 IF A$="B" THEN 320
  43. 240 IF A$="C" THEN 340
  44. 250 IF A$="D" THEN 360
  45. 260 IF A$="E" THEN 380
  46. 270 IF A$="F" THEN 400
  47. 280 IF A$="G" THEN 420
  48. 281 IF A$="H" THEN 426
  49. 282 IF A$="Q" THEN 430
  50. 283 IF A$="PLOT" THEN 8
  51. 290 GOTO 19
  52. 300 GOSUB 1000
  53. 310 GOTO 19
  54. 320 GOSUB 5000
  55. 330 GOTO 19
  56. 340 GOSUB 10000
  57. 350 GOTO 19
  58. 360 GOSUB 15000
  59. 370 GOTO 19
  60. 380 GOSUB 20000
  61. 390 GOTO 19
  62. 400 GOSUB 25000
  63. 410 GOTO 19
  64. 420 GOSUB 30000
  65. 425 GOTO 19
  66. 426 GOSUB 50000
  67. 427 GOTO 19
  68. 430 PRINT
  69. 435 INPUT"QUIT (Y/N)";B$
  70. 440 IF B$="Y" THEN 460
  71. 450 GOTO 19
  72. 460 PRINT
  73. 470 INPUT"DO YOU WANT TO SAVE THE FREQUENCY VARIABLES IN A FILE (Y/N)";C$
  74. 480 IF C$="N" THEN 610
  75. 490 PRINT
  76. 500 INPUT"DECLARE FILENAME";D$
  77. 510 OPEN D$ FOR OUTPUT AS #1
  78. 515 WRITE #1,NZ%
  79. 520 IF NZ%=0 THEN 565
  80. 530 FOR I=1 TO NZ%
  81. 540 WRITE #1,ZR(I)
  82. 550 WRITE #1,ZJ(I)
  83. 560 NEXT I
  84. 565 WRITE #1,NP%
  85. 570 FOR I=1 TO NP%
  86. 580 WRITE #1,PR(I)
  87. 590 WRITE #1,PJ(I)
  88. 600 NEXT I
  89. 610 END
  90. 1000 CLS
  91. 1010 INPUT"DECLARE FILENAME";F$
  92. 1020 OPEN F$ FOR INPUT AS #1
  93. 1030 INPUT #1,NZ%
  94. 1040 IF NZ%=0 THEN 1090
  95. 1050 FOR I=1 TO NZ%
  96. 1060 INPUT #1,ZR(I)
  97. 1070 INPUT #1,ZJ(I)
  98. 1080 NEXT I
  99. 1090 INPUT #1,NP%
  100. 1100 FOR I=1 TO NP%
  101. 1110 INPUT #1,PR(I)
  102. 1120 INPUT #1,PJ(I)
  103. 1130 NEXT I
  104. 1135 CLOSE #1
  105. 1270 RETURN
  106. 5000 CLS
  107. 5010 PRINT"INPUT # OF ZEROS";
  108. 5020 INPUT NZ%
  109. 5030 PRINT"INPUT # OF POLES";
  110. 5040 INPUT NP%
  111. 5050 PRINT
  112. 5060 IF NZ%=0 THEN 5150
  113. 5070 FOR I=1 TO NZ%
  114. 5080 PRINT"ZERO #";I;"REAL PART";
  115. 5090 INPUT ZR(I)
  116. 5100 PRINT"ZERO #";I;"IMAGINARY PART";
  117. 5110 INPUT ZJ(I)
  118. 5130 NEXT I
  119. 5140 PRINT
  120. 5150 FOR I=1 TO NP%
  121. 5160 PRINT"POLE #";I;"REAL PART";
  122. 5170 INPUT PR(I)
  123. 5180 PRINT"POLE #";I;"IMAGINARY PART";
  124. 5190 INPUT PJ(I)
  125. 5200 NEXT I
  126. 5210 RETURN
  127. 10000 CLS
  128. 10010 PRINT"THERE ARE";NZ%;"ZEROES AND";NP%;"POLES"
  129. 10020 PRINT
  130. 10030 IF NZ%=0 THEN 10090
  131. 10040 FOR I=1 TO NZ%
  132. 10050 PRINT"ZERO #";I;"REAL PART IS";ZR(I)
  133. 10060 PRINT"ZERO #";I;"IMAGINARY PART IS";ZJ(I)
  134. 10070 NEXT I
  135. 10080 PRINT
  136. 10090 FOR I=1 TO NP%
  137. 10100 PRINT"POLE #";I;"REAL PART IS";PR(I)
  138. 10110 PRINT"POLE #";I;"IMAGINARY PART IS";PJ(I)
  139. 10120 NEXT I
  140. 10130 PRINT
  141. 10140 PRINT"EDIT ZEROES (Y/N)";
  142. 10150 INPUT B$
  143. 10160 PRINT
  144. 10170 IF B$="N" THEN 10750
  145. 10180 PRINT"ADD,DELETE,OR SUBSTITUTE (A,D,S)";
  146. 10190 INPUT C$
  147. 10200 PRINT
  148. 10210 IF C$="S" THEN 10630
  149. 10220 IF C$="D" THEN 10380
  150. 10230 IF C$="A" THEN 10240 ELSE 10180
  151. 10240 PRINT"HOW MANY";
  152. 10250 INPUT NC%
  153. 10260 ND%=NZ%+NC%
  154. 10270 PRINT
  155. 10280 FOR I=1 TO ND%
  156. 10290 IF I>NZ% THEN 10310
  157. 10300 GOTO 10350
  158. 10310 PRINT"ZER0 #";I;"REAL PART";
  159. 10320 INPUT ZR(I)
  160. 10330 PRINT"ZERO #";I;"IMAGINARY PART";
  161. 10340 INPUT ZJ(I)
  162. 10350 NEXT I
  163. 10360 NZ%=NZ%+NC%
  164. 10370 GOTO 10750
  165. 10380 PRINT"HOW MANY";
  166. 10390 INPUT NC%
  167. 10400 PRINT
  168. 10410 FOR J=1 TO NC%
  169. 10420 INPUT"DELETE ZERO #";D%(J)
  170. 10430 NEXT J
  171. 10440 M=0
  172. 10450 FOR J=1 TO NC%
  173. 10460 FOR I=1 TO NZ%
  174. 10470 IF I=D%(J) THEN 10510
  175. 10480 M=M+1
  176. 10490 TZR(M)=ZR(I)
  177. 10500 TZJ(M)=ZJ(I)
  178. 10510 NEXT I
  179. 10520 NEXT J
  180. 10530 NZ%=NZ%-NC%
  181. 10560 FOR I=1 TO NZ%
  182. 10570 ZR(I)=TZR(I)
  183. 10580 ZJ(I)=TZJ(I)
  184. 10590 NEXT I
  185. 10620 GOTO 10750
  186. 10630 FOR I=1 TO NZ%
  187. 10640 PRINT"ZERO #";I;"IS";ZR(I);"AND J";ZJ(I)
  188. 10650 PRINT
  189. 10660 PRINT"CHANGE (Y/N)";
  190. 10670 INPUT D$
  191. 10675 PRINT
  192. 10680 IF D$="N" THEN 10740
  193. 10700 PRINT"ZERO #";I;"REAL PART";
  194. 10710 INPUT ZR(I)
  195. 10720 PRINT"ZERO #";I;"IMAGINARY PART";
  196. 10730 INPUT ZJ(I)
  197. 10740 NEXT I
  198. 10750 PRINT
  199. 10760 PRINT"EDIT POLES (Y/N)";
  200. 10770 INPUT A$
  201. 10780 IF A$="N" THEN 11380
  202. 10790 PRINT
  203. 10800 REM START POLE EDIT
  204. 10810 PRINT"ADD,DELETE,OR SUBSTITUTE POLES (A,D,S)";
  205. 10820 INPUT A$
  206. 10830 PRINT
  207. 10840 IF A$="S" THEN 11260
  208. 10850 IF A$="D" THEN 11010
  209. 10860 IF A$="A" THEN 10870 ELSE 10810
  210. 10870 PRINT"HOW MANY";
  211. 10880 INPUT NC%
  212. 10890 ND%=NP%+NC%
  213. 10900 PRINT
  214. 10910 FOR I=1 TO ND%
  215. 10920 IF I>NP% THEN 10940
  216. 10930 GOTO 10980
  217. 10940 PRINT"POLE #";I;"REAL PART";
  218. 10950 INPUT PR(I)
  219. 10960 PRINT"POLE #";I;"IMAGINARY PART";
  220. 10970 INPUT PJ(I)
  221. 10980 NEXT I
  222. 10990 NP%=NP%+NC%
  223. 11000 GOTO 11380
  224. 11010 PRINT"HOW MANY";
  225. 11020 INPUT NC%
  226. 11030 PRINT
  227. 11040 FOR J=1 TO NC%
  228. 11050 INPUT"DELETE POLE #";D%(J)
  229. 11060 NEXT J
  230. 11070 M=O
  231. 11080 FOR J=1 TO NC%
  232. 11090 FOR I=1 TO NP%
  233. 11100 IF I=D%(J) THEN 11140
  234. 11110 M=M+1
  235. 11120 TPR(M)=PR(I)
  236. 11130 TPJ(M)=PJ(I)
  237. 11140 NEXT I
  238. 11150 NEXT J
  239. 11160 NP%=NP%-NC%
  240. 11190 FOR I=1 TO NP%
  241. 11200 PR(I)=TPR(I)
  242. 11210 PJ(I)=TPJ(I)
  243. 11220 NEXT I
  244. 11250 GOTO 11380
  245. 11260 FOR I=1 TO NP%
  246. 11270 PRINT"POLE #";I;"IS";PR(I);"AND J";PJ(I)
  247. 11280 PRINT
  248. 11290 PRINT"CHANGE (Y/N)";
  249. 11300 INPUT C$
  250. 11305 PRINT
  251. 11310 IF C$="N" THEN 11370
  252. 11330 PRINT"POLE #";I;"REAL PART";
  253. 11340 INPUT PR(I)
  254. 11350 PRINT"POLE #";I;"IMAGINARY PART";
  255. 11360 INPUT PJ(I)
  256. 11370 NEXT I
  257. 11380 RETURN
  258. 15000 CLS
  259. 15010 PRINT"REQUIRED ATTENUATION IN dB";
  260. 15020 INPUT Z21
  261. 15040 PRINT
  262. 15050 PRINT"INPUT CUTTOFF FREQUENCY RATIO OF REQUIRED ATTENUATION";
  263. 15060 INPUT W
  264. 15070 IF W<=1# THEN 15050
  265. 15080 PRINT
  266. 15090 PRINT"PERMISSABLE PASSBAND RIPPLE IN dB , 0dB FOR BUTTERWORTH";
  267. 15100 INPUT R
  268. 15120 K#=10#^(Z21/10#)
  269. 15125 IF R=0 THEN 15440
  270. 15130 RES=(10#^(R/10#))-1#
  271. 15140 RE=SQR(RES)
  272. 15150 TW=SQR((K#-1#)/RES)
  273. 15160 FOR I=1 TO 20
  274. 15170 KA#=((W+SQR(W^2#-1#))^I+(W+SQR(W^2#-1#))^(1#/I))/2#
  275. 15180 IF KA#=>TW THEN 15200
  276. 15190 NEXT I
  277. 15200 PRINT
  278. 15210 PRINT"N=";I
  279. 15230 PRINT
  280. 15240 PRINT"RESULTS SATISFACTORY (Y/N)";
  281. 15250 INPUT A$
  282. 15260 IF A$="Y" THEN 15270 ELSE 15000
  283. 15270 X=(SQR((1#/RES)+1#)+1#/RE)^(1#/I)
  284. 15280 Y=1#/X
  285. 15290 SINH=(X-Y)/2#
  286. 15300 COSH#=(X+Y)/2#
  287. 15310 N%=I-1
  288. 15320 PRINT
  289. 15330 PRINT"FREQUENCY SCALE FACTOR";
  290. 15340 INPUT KF#
  291. 15350 PRINT
  292. 15370 PRINT
  293. 15380 FOR J=0 TO N%
  294. 15390 PR(J+1+NP%)=KF#*(-SINH*(SIN(((2#*J+1#)*PI)/(2#*I))))
  295. 15400 PJ(J+1+NP%)=KF#*(COSH#*(COS(((2#*J+1#)*PI)/(2#*I))))
  296. 15420 NEXT J
  297. 15430 GOTO 15610
  298. 15440 FOR I=1 TO 20
  299. 15450 KA#=1+W^(2#*I)
  300. 15460 IF KA#=>K# THEN 15480
  301. 15470 NEXT I
  302. 15480 PRINT
  303. 15490 PRINT"N=";I
  304. 15500 PRINT
  305. 15510 INPUT"RESULTS SATISFACTORY";A$
  306. 15520 IF A$="Y" THEN 15530 ELSE 15000
  307. 15530 PRINT
  308. 15540 INPUT"FREQUENCY SCALE FACTOR";KF#
  309. 15550 N%=I
  310. 15560 FOR J=1 TO N%
  311. 15570 P(J)=PI/2#+(PI/(2#*N%))*(2#*J-1#)
  312. 15580 PR(J+NP%)=KF#*(COS(P(J)))
  313. 15590 PJ(J+NP%)=KF#*(SIN(P(J)))
  314. 15600 NEXT J
  315. 15610 N=I
  316. 15620 J=0
  317. 15622 FOR L=1 TO N
  318. 15624 TPR(L)=PR(NP%+L)
  319. 15626 TPJ(L)=PJ(NP%+L)
  320. 15628 NEXT L
  321. 15630 FOR L=1 TO N
  322. 15640 J=J+1
  323. 15670 IF ABS(TPJ(L))=>ABS((10^-6)*TPR(L)) THEN 15680 ELSE 15750
  324. 15680 IF TPJ(L)<0 THEN 15770
  325. 15690 PR(NP%+J)=TPR(L)
  326. 15700 PJ(NP%+J)=TPJ(L)
  327. 15710 PR(NP%+J+1)=TPR(L)
  328. 15720 PJ(NP%+J+1)=-TPJ(L)
  329. 15730 J=J+1
  330. 15740 GOTO 15770
  331. 15750 PR(NP%+J)=TPR(L)
  332. 15760 PJ(NP%+J)=0#
  333. 15770 NEXT L
  334. 15780 NP%=NP%+I
  335. 15790 RETURN
  336. 20000 CLS
  337. 20010 KN#=1#
  338. 20020 IF NZ%=0 THEN 20080
  339. 20030 FOR I=1 TO NZ%
  340. 20040 ZT=ABS(ZR(I))+ABS(ZJ(I))
  341. 20050 IF ZT=0 THEN 20070
  342. 20060 KN#=KN#*SQR(ZR(I)^2+ZJ(I)^2)
  343. 20070 NEXT I
  344. 20080 KD#=1#
  345. 20090 FOR I=1 TO NP%
  346. 20100 PT=ABS(PR(I))+ABS(PJ(I))
  347. 20110 IF PT=0 THEN 20130
  348. 20120 KD#=KD#*SQR(PR(I)^2+PJ(I)^2)
  349. 20130 NEXT I
  350. 20140 PRINT
  351. 20150 PRINT"INPUT LOOP GAIN";
  352. 20160 INPUT U
  353. 20170 K#=(KD#/KN#)*U
  354. 20172 PRINT
  355. 20174 INPUT"PLOT (Y/N)";P$
  356. 20180 PRINT
  357. 20190 PRINT "START FREQ.";
  358. 20200 INPUT WS
  359. 20210 PRINT"NUMBER OF DECADES";
  360. 20220 INPUT WF
  361. 20222 IF P$="Y" THEN 20252
  362. 20230 PRINT"DECADE INCREMENT";
  363. 20240 INPUT WD
  364. 20250 N%=WF/WD
  365. 20251 GOTO 20260
  366. 20252 N%=500:WD=WF/500
  367. 20260 PRINT
  368. 20270 PRINT
  369. 20272 IF P$="Y" THEN 20290
  370. 20280 PRINT USING"\            \";"W","H(W)-dB","THETA(W)","H(W)'-dB","THETA(W)'"
  371. 20290 FOR I=0 TO N%
  372. 20300 W=WS*(10^(I*WD))
  373. 20310 PM=1
  374. 20320 PT=0
  375. 20330 TD=0
  376. 20340 FOR J=1 TO NP%
  377. 20350 PM=PM/(SQR(PR(J)^2+(W-PJ(J))^2))
  378. 20360 IF PR(J)=0 THEN 20390
  379. 20370 PT=PT+ATN((W-PJ(J))/PR(J))
  380. 20380 GOTO 20400
  381. 20390 PT=PT-PI/2
  382. 20400 NEXT J
  383. 20410 IF NZ%=0 THEN 20490
  384. 20420 FOR H=1 TO NZ%
  385. 20430 PM=PM*SQR(ZR(H)^2+(W-ZJ(H))^2)
  386. 20440 IF ZR(H)=0 THEN 20470
  387. 20450 PT=PT-ATN((W-ZJ(H))/ZR(H))
  388. 20460 GOTO 20480
  389. 20470 PT=PT+PI/2
  390. 20480 NEXT H
  391. 20490 PM=PM*K#
  392. 20500 R=PM*COS(PT)
  393. 20510 S=PM*SIN(PT)
  394. 20520 R=R+1
  395. 20530 MC#=PM/(SQR(R*R+S*S))
  396. 20540 X=R/(SQR(R*R+S*S))
  397. 20550 IF 1-X*X=0 THEN 20570
  398. 20560 PC=PT-SGN(S)*(PI/2-ATN(X/SQR(1-X*X)))
  399. 20570 XM=8.68589*LOG(PM)
  400. 20580 XC=8.68589*LOG(MC#)
  401. 20582 IF P$="Y" THEN 20592
  402. 20590 PRINT USING"##.###^^^^    ";W,XM,PT,XC,PC
  403. 20591 GOTO 20600
  404. 20592 XPASS!(I)=W:YPASS!(I)=XM
  405. 20600 NEXT I
  406. 20601 IF P$="Y" THEN 20602 ELSE 20610
  407. 20602 CLS
  408. 20603 NCURVES!=1
  409. 20604 NPOINTS!=500
  410. 20605 CALL PLOT((NCURVES!),(NPOINTS!),XPASS!(),YPASS!())
  411. 20610 PRINT
  412. 20620 PRINT"NEW GAIN VALUE (Y/N)";
  413. 20630 INPUT A$
  414. 20640 PRINT
  415. 20650 IF A$="Y" THEN 20150
  416. 20660 RETURN
  417. 25000 CLS
  418. 25010 KN#=1#
  419. 25020 IF NZ%=0 THEN 25080
  420. 25030 FOR I=1 TO NZ%
  421. 25040 ZT=ABS(ZR(I))+ABS(ZJ(I))
  422. 25050 IF ZT=0 THEN 25070
  423. 25060 KN#=KN#*SQR(ZR(I)^2+ZJ(I)^2)
  424. 25070 NEXT I
  425. 25080 KD#=1#
  426. 25090 FOR I=1 TO NP%
  427. 25100 PT=ABS(PR(I))+ABS(PJ(I))
  428. 25110 IF PT=0 THEN 25130
  429. 25120 KD#=KD#*SQR(PR(I)^2+PJ(I)^2)
  430. 25130 NEXT I
  431. 25140 PRINT
  432. 25150 PRINT"INPUT LOOP GAIN";
  433. 25160 INPUT U
  434. 25170 K#=(KD#/KN#)*U
  435. 25180 PRINT
  436. 25190 PRINT"INPUT APPROXIMATE -3dB RADIAN FREQUENCY";
  437. 25200 INPUT WC
  438. 25210 PRINT
  439. 25220 XN=0
  440. 25230 FOR A%=1 TO 3
  441. 25240 IF A%=1 THEN 25270
  442. 25250 WD=WC*(10^(A%-3))
  443. 25260 GOTO 25540
  444. 25270 FOR B%=1 TO 100
  445. 25280 WD=WC/100#
  446. 25290 W=WD*B%
  447. 25300 PM=1
  448. 25310 PT=0
  449. 25320 FOR J=1 TO NP%
  450. 25330 PM=PM/(SQR(PR(J)^2+(W-PJ(J))^2))
  451. 25340 IF PR(J)=0 THEN 25370
  452. 25350 PT=PT+ATN((W-PJ(J))/PR(J))
  453. 25360 GOTO 25380
  454. 25370 PT=PT-PI/2
  455. 25380 NEXT J
  456. 25390 IF NZ%=0 THEN 25470
  457. 25400 FOR H=1 TO NZ%
  458. 25410 PM=PM*SQR(ZR(H)^2+(W-ZJ(H))^2)
  459. 25420 IF ZR(H)=0 THEN 25450
  460. 25430 PT=PT-ATN((W-ZJ(H))/ZR(H))
  461. 25440 GOTO 25460
  462. 25450 PT=PT+PI/2
  463. 25460 NEXT H
  464. 25470 PM=PM*K#
  465. 25472 IF PM=>1# AND ABS(PT)>PI THEN 25474 ELSE 25480
  466. 25474 PRINT
  467. 25476 PRINT"THE SYSTEM IS UNSTABLE"
  468. 25478 GOTO 25810
  469. 25480 R=PM*COS(PT)
  470. 25490 S=PM*SIN(PT)
  471. 25500 R=R+1
  472. 25510 XN=XN+((PM*PM)/(R*R+S*S))*WD
  473. 25520 NEXT B%
  474. 25530 IF A%=1 THEN 25790
  475. 25540 FOR C%=1 TO 90
  476. 25550 W=WD*(10+C%)
  477. 25560 PM=1
  478. 25570 PT=0
  479. 25580 FOR J=1 TO NP%
  480. 25590 PM=PM/(SQR(PR(J)^2+(W-PJ(J))^2))
  481. 25600 IF PR(J)=0 THEN 25630
  482. 25610 PT=PT+ATN((W-PJ(J))/PR(J))
  483. 25620 GOTO 25640
  484. 25630 PT=PT-PI/2
  485. 25640 NEXT J
  486. 25650 IF NZ%=0 THEN 25730
  487. 25660 FOR H=1 TO NZ%
  488. 25670 PM=PM*SQR(ZR(H)^2+(W-ZJ(H))^2)
  489. 25680 IF ZR(H)=0 THEN 25710
  490. 25690 PT=PT-ATN((W-ZJ(H))/ZR(H))
  491. 25700 GOTO 25720
  492. 25710 PT=PT+PI/2
  493. 25720 NEXT H
  494. 25730 PM=PM*K#
  495. 25740 R=PM*COS(PT)
  496. 25750 S=PM*SIN(PT)
  497. 25760 R=R+1
  498. 25770 XN=XN+((PM*PM)/(R*R+S*S))*WD
  499. 25780 NEXT C%
  500. 25790 NEXT A%
  501. 25800 PRINT "THE LOOP NOISE BANDWIDTH IS";XN;"RAD/SEC"
  502. 25810 PRINT
  503. 25820 PRINT"NEW LOOP GAIN (Y/N)";
  504. 25830 INPUT A$
  505. 25840 PRINT
  506. 25850 IF A$="N" THEN 25860 ELSE 25150
  507. 25860 RETURN
  508. 30000 CLS
  509. 30010 FOR I=1 TO NZ%
  510. 30020 XR(I)=ZR(I)
  511. 30030 XJ(I)=ZJ(I)
  512. 30040 NEXT I
  513. 30050 N=NZ%:NX%=NZ%:I=NZ%
  514. 30060 GOSUB 40000
  515. 30070 FOR I=1 TO NZ%
  516. 30080 ZR(I)=XR(I)
  517. 30090 ZJ(I)=XJ(I)
  518. 30100 NEXT I
  519. 30110 FOR I=1 TO NP%
  520. 30120 XR(I)=PR(I)
  521. 30130 XJ(I)=PJ(I)
  522. 30140 NEXT I
  523. 30150 N=NP%:NX%=NP%:I=NP%
  524. 30160 GOSUB 40000
  525. 30170 FOR I=1 TO NP%
  526. 30180 PR(I)=XR(I)
  527. 30190 PJ(I)=XJ(I)
  528. 30200 NEXT I
  529. 30210 GOSUB 42000
  530. 30220 GOSUB 48000
  531. 30221 IF ITAG=4 THEN 30239
  532. 30222 IF ITAG=1 THEN 30225
  533. 30223 IF ITAG=2 THEN 30228
  534. 30224 IF ITAG=3 THEN 30231
  535. 30225 PRINT
  536. 30226 PRINT"DIVIDE BY ZERO IN POLYNOMIAL ROOT SUBROUTINE"
  537. 30227 GOTO 30233
  538. 30228 PRINT
  539. 30229 PRINT"POLYNOMIAL ROOT SUBROUTINE IS DIVERGING"
  540. 30230 GOTO 30233
  541. 30231 PRINT
  542. 30232 PRINT"POLYNOMIAL ROOT SUBROUTINE ITERATION LIMIT EXCEEDED"
  543. 30233 PRINT
  544. 30234 PRINT"PRESS ANY KEY WHEN READY"
  545. 30235 A$=INKEY$:IF A$="" THEN 30235
  546. 30236 GOTO 30280
  547. 30239 J=0
  548. 30240 FOR I=1 TO 2*NP% STEP 2
  549. 30242 J=J+1
  550. 30250 PR(J)=ROOT(I):PJ(J)=ROOT(I+1)
  551. 30260 NEXT I
  552. 30280 RETURN
  553. 40000 J=0
  554. 40010 FOR L=1 TO N
  555. 40020 TPR(L)=XR(NX%+L)
  556. 40030 TPJ(L)=XJ(NX%+L)
  557. 40040 NEXT L
  558. 40050 FOR L=1 TO N
  559. 40060 J=J+1
  560. 40070 IF ABS(TPJ(L))=>ABS((10^-6)*TPR(L)) THEN 40080 ELSE 40150
  561. 40080 IF TPJ(L)<0 THEN 40170
  562. 40090 XR(NX%+J)=TPR(L)
  563. 40100 XJ(NX%+J)=TPJ(L)
  564. 40110 XR(NX%+J+1)=TPR(L)
  565. 40120 XJ(NX%+J+1)=-TPJ(L)
  566. 40130 J=J+1
  567. 40140 GOTO 40170
  568. 40150 XR(NX%+J)=TPR(L)
  569. 40160 XJ(NX%+J)=0#
  570. 40170 NEXT L
  571. 40180 NX%=NX%+I
  572. 40190 RETURN
  573. 42000 KN#=1#
  574. 42010 IF NZ%=0 THEN 42070
  575. 42020 FOR I=1 TO NZ%
  576. 42030 ZT=ABS(ZR(I))+ABS(ZJ(I))
  577. 42040 IF ZT=0 THEN 42060
  578. 42050 KN#=KN#*SQR(ZR(I)^2+ZJ(I)^2)
  579. 42060 NEXT I
  580. 42070 KD#=1#
  581. 42080 FOR I=1 TO NP%
  582. 42090 PT=ABS(PR(I))+ABS(PJ(I))
  583. 42100 IF PT=0 THEN 42120
  584. 42110 KD#=KD#*SQR(PR(I)^2+PJ(I)^2)
  585. 42120 NEXT I
  586. 42130 CLS
  587. 42140 INPUT"LOOP GAIN";K#
  588. 42150 K#=(KD#/KN#)*K#
  589. 42160 IF NZ%=0 THEN 42290
  590. 42170 IF NZ%=1 THEN 42310
  591. 42180 IF NZ%=2 THEN 42340
  592. 42190 NX%=NZ%
  593. 42200 FOR I=1 TO NX%
  594. 42210 XR(I)=ZR(I)
  595. 42220 XJ(I)=ZJ(I)
  596. 42230 NEXT I
  597. 42240 GOSUB 44000
  598. 42250 FOR I=1 TO NX%+1
  599. 42260 ZC(I)=K#*MAT3(I)
  600. 42270 NEXT I
  601. 42280 GOTO 42370
  602. 42290 ZC(1)=K#
  603. 42300 GOTO 42370
  604. 42310 ZC(1)=-ZR(1)*K#
  605. 42320 ZC(2)=K#
  606. 42330 GOTO 42370
  607. 42340 ZC(1)=K#*(ZR(1)*ZR(2)+ABS(ZJ(1)*ZJ(2)))
  608. 42350 ZC(2)=-(ZR(1)+ZR(2))*K#
  609. 42360 ZC(3)=K#
  610. 42370 IF NP%=1 THEN 42490
  611. 42380 IF NP%=2 THEN 42520
  612. 42390 NX%=NP%
  613. 42400 FOR I=1 TO NX%
  614. 42410 XR(I)=PR(I)
  615. 42420 XJ(I)=PJ(I)
  616. 42430 NEXT I
  617. 42440 GOSUB 44000
  618. 42450 FOR I=1 TO NX%+1
  619. 42460 PC(I)=MAT3(I)
  620. 42470 NEXT I
  621. 42480 GOTO 42550
  622. 42490 PC(1)=-PR(I)
  623. 42500 PC(2)=1#
  624. 42510 GOTO 42550
  625. 42520 PC(1)=PR(1)*PR(2)+ABS(PJ(1)*PJ(2))
  626. 42530 PC(2)=-(PR(1)+PR(2))
  627. 42540 PC(3)=1#
  628. 42550 FOR I=1 TO NZ%+1
  629. 42560 PC(I)=PC(I)+ZC(I)
  630. 42570 NEXT I
  631. 42580 PRINT
  632. 42590 PRINT"THE NUMERATOR COEFFICIENTS ARE :"
  633. 42600 PRINT
  634. 42610 FOR I=0 TO NZ%
  635. 42620 PRINT"S";I;" = ";ZC(I+1)
  636. 42630 NEXT I
  637. 42640 PRINT
  638. 42650 PRINT"THE DENOMINATOR COEFFICIENTS ARE :"
  639. 42660 PRINT
  640. 42670 FOR I=0 TO NP%
  641. 42680 PRINT"S";I;" = ";PC(I+1)
  642. 42690 NEXT I
  643. 42700 PRINT
  644. 42710 PRINT"PRESS ANY KEY WHEN READY"
  645. 42720 A$=INKEY$:IF A$="" THEN 42720
  646. 42730 RETURN
  647. 44000 I=1
  648. 44010 IF XJ(I)=0 THEN 44020 ELSE 44070
  649. 44020 MAT1(1)=-XR(I)
  650. 44030 MAT1(2)=1#
  651. 44040 NDA=1
  652. 44050 I=I+1
  653. 44060 GOTO 44120
  654. 44070 MAT1(1)=XR(I)^2#+XJ(I)^2#
  655. 44080 MAT1(2)=-2#*XR(I)
  656. 44090 MAT1(3)=1#
  657. 44100 NDA=2
  658. 44110 I=I+2
  659. 44120 IF XJ(I)=0 THEN 44130 ELSE 44180
  660. 44130 MAT2(1)=-XR(I)
  661. 44140 MAT2(2)=1#
  662. 44150 NDB=1
  663. 44160 I=I+1
  664. 44170 GOTO 44230
  665. 44180 MAT2(1)=XR(I)^2#+XJ(I)^2#
  666. 44190 MAT2(2)=-2#*XR(I)
  667. 44200 MAT2(3)=1#
  668. 44210 NDB=2
  669. 44220 I=I+2
  670. 44230 GOSUB 46000
  671. 44240 FOR J=1 TO NDC+1
  672. 44250 MAT1(J)=MAT3(J)
  673. 44260 NEXT J
  674. 44270 NDA=NDC
  675. 44280 IF NDC=NX% THEN 44300
  676. 44290 GOTO 44120
  677. 44300 RETURN
  678. 46000 NDC = NDA + NDB
  679. 46010 QNDAP1% = NDA + 1
  680. 46020 QNDCP1% = NDC + 1
  681. 46030  FOR QI% = 1 TO QNDCP1%
  682. 46040   MAT3(QI%) = 0#
  683. 46050  NEXT QI%
  684. 46060 QNDBP1% = NDB + 1
  685. 46070  FOR QI% = 1 TO QNDAP1%
  686. 46080   FOR QJ% = 1 TO QNDBP1%
  687. 46090    QIPJ% = QI% + QJ% - 1
  688. 46100   MAT3(QIPJ%) = MAT3(QIPJ%) + MAT1(QI%)*MAT2(QJ%)
  689. 46110  NEXT QJ%
  690. 46120 NEXT QI%
  691. 46130 RETURN
  692. 48000 FOR I=NP% TO 0 STEP -1
  693. 48010 MAT1(NP%-I+1)=PC(I+1)
  694. 48020 NEXT I
  695. 48030 N=NP%
  696. 48031 SF=(MAT1(N+1))^(1/N)
  697. 48040 P1=1
  698. 48050 Q1=1
  699. 48060 ITAG=1000
  700. 48070 EPS#=.0001#
  701. 48080 QIITAG = ITAG
  702. 48090 ITAG = 0
  703. 48100 QD = MAT1(1)
  704. 48110  FOR QI% = 1 TO N
  705. 48120   MAT1(QI%) = MAT1(QI% + 1)/QD
  706. 48130  NEXT QI%
  707. 48140 IF N > 1 THEN 48220
  708. 48150 ITAG = ITAG + 1
  709. 48160 QL = ITAG + ITAG
  710. 48170 ROOT(QL) = 0!
  711. 48180 ROOT(QL-1) = -MAT1(1)
  712. 48190 N = ITAG
  713. 48200 ITAG = 4
  714. 48210 RETURN
  715. 48220 IF N > 2 THEN 48260
  716. 48230 QP = MAT1(1)
  717. 48240 QQ = MAT1(2)
  718. 48250 GOTO 48680
  719. 48260 QP = P1
  720. 48270 QQ = Q1
  721. 48280 QM = 1
  722. 48290 MAT2(1) = MAT1(1) - QP
  723. 48300 MAT2(2) = MAT1(2) - QP*MAT2(1) - QQ
  724. 48310 QL = N - 1
  725. 48320 MAT3(1) = MAT2(1) - QP
  726. 48330 MAT3(2) = MAT2(2) - QP*MAT3(1) - QQ
  727. 48340 IF QL = 2 THEN QL = 3
  728. 48350  FOR QJ% = 3 TO QL
  729. 48360   MAT2(QJ%) = MAT1(QJ%) - QP*MAT2(QJ%-1) - QQ*MAT2(QJ%-2)
  730. 48370   MAT3(QJ%) = MAT2(QJ%) - QP*MAT3(QJ%-1) - QQ*MAT3(QJ%-2)
  731. 48380  NEXT QJ%
  732. 48390 QL = N - 1
  733. 48400 QCBARL = MAT3(QL) - MAT2(QL)
  734. 48410 QDEN = -QCBARL
  735. 48420 IF N <> 3 THEN QDEN = QDEN*MAT3(N-3)
  736. 48430 QDEN = QDEN + MAT3(N-2)*MAT3(N-2)
  737. 48440 IF QDEN <> 0 THEN 48480
  738. 48450 N = ITAG
  739. 48460 ITAG = 1
  740. 48470 RETURN
  741. 48480 MAT2(N) = MAT1(N) - QP*MAT2(N-1) - QQ*MAT2(N-2)
  742. 48490 QDELTP = -MAT2(N)
  743. 48500 IF N <> 3 THEN QDELTP = QDELTP*MAT3(N-3)
  744. 48510 QDELTP = (MAT2(N-1)*MAT3(N-2) + QDELTP)/QDEN
  745. 48520 QDELTQ = (MAT2(N)*MAT3(N-2) - MAT2(N-1)*QCBARL)/QDEN
  746. 48530 QP = QP + QDELTP
  747. 48540 QQ = QQ + QDELTQ
  748. 48550 QSUM = ABS(QDELTP) + ABS(QDELTQ)
  749. 48560 IF QM = 1 THEN QSUM1 = QSUM
  750. 48570 IF QM <> 5 OR QSUM <= QSUM1 THEN 48610
  751. 48580 N = ITAG
  752. 48590 ITAG = 2
  753. 48600 RETURN
  754. 48610 IF QSUM <= EPS# THEN 48680
  755. 48620 IF QM < QIITAG THEN 48660
  756. 48630 N = ITAG
  757. 48640 ITAG = 3
  758. 48650 RETURN
  759. 48660 QM = QM + 1
  760. 48670 GOTO 48290
  761. 48680 QD = -QP*.5
  762. 48690 ITAG = ITAG + 2
  763. 48700 QF = QQ - QD*QD
  764. 48710 QE = SQR(ABS(QF))
  765. 48720 QL = ITAG + ITAG
  766. 48730 IF QF > 0! THEN 48790
  767. 48740 ROOT(QL) = 0!
  768. 48750 ROOT(QL-1) = QD-QE
  769. 48760 ROOT(QL-2) = 0!
  770. 48770 ROOT(QL-3) = QD + QE
  771. 48780 GOTO 48830
  772. 48790 ROOT(QL) = -QE
  773. 48800 ROOT(QL-1) = QD
  774. 48810 ROOT(QL-2) = QE
  775. 48820 ROOT(QL-3) = QD
  776. 48830 N = N-2
  777. 48840 IF N > 0 THEN 48880
  778. 48850 N = ITAG
  779. 48860 ITAG = 4
  780. 48870 RETURN
  781. 48880  FOR QI% = 1 TO N
  782. 48890   MAT1(QI%) = MAT2(QI%)
  783. 48900  NEXT QI%
  784. 48910 GOTO 48140
  785. 48920 RETURN
  786. 50000 CLS
  787. 50010 KN#=1#
  788. 50020 IF NZ%=0 THEN 50080
  789. 50030 FOR I=1 TO NZ%
  790. 50040 ZT=ABS(ZR(I))+ABS(ZJ(I))
  791. 50050 IF ZT=0 THEN 50070
  792. 50060 KN#=KN#*SQR(ZR(I)^2#+ZJ(I)^2#)
  793. 50070 NEXT I
  794. 50080 KD#=1#
  795. 50090 FOR I=1 TO NP%
  796. 50100 PT=ABS(PR(I))+ABS(PJ(I))
  797. 50110 IF PT=0 THEN 50130
  798. 50120 KD#=KD#*SQR(PR(I)^2#+PJ(I)^2#)
  799. 50130 NEXT I
  800. 50140 K#=KD#/KN#
  801. 50150 PRINT
  802. 50160 FOR I=1 TO NP%
  803. 50170 P=1#
  804. 50180 T=0
  805. 50190 FOR J=1 TO NP%
  806. 50200 PA=PR(I)-PR(J)
  807. 50210 PB=PJ(I)-PJ(J)
  808. 50220 PD=ABS(PA)+ABS(PB)
  809. 50230 IF PD=0 THEN 50330
  810. 50240 IF PB=0 THEN 50290
  811. 50250 PM=SQR(PA^2#+PB^2#)
  812. 50260 X=PA/PM
  813. 50270 PT=(PI/2#-ATN(X/SQR(1-X*X)))*SGN(PB)
  814. 50280 GOTO 50310
  815. 50290 PM=ABS(PA)
  816. 50300 PT=PI/2#+PI/2#*SGN(-PA)
  817. 50310 T=T+PT
  818. 50320 P=P*PM
  819. 50330 NEXT J
  820. 50340 RM(I)=P
  821. 50350 RP(I)=T
  822. 50360 NEXT I
  823. 50370 PRINT
  824. 50380 IF NZ%=0 THEN 50620
  825. 50390 FOR I=1 TO NP%
  826. 50400 P=1#
  827. 50410 T=0
  828. 50420 FOR J=1 TO NZ%
  829. 50430 ZA=PR(I)-ZR(J)
  830. 50440 ZB=PJ(I)-ZJ(J)
  831. 50450 ZD=ABS(ZA)+ABS(ZB)
  832. 50460 IF ZD=0 THEN 50560
  833. 50470 IF ZB=0 THEN 50520
  834. 50480 ZM=SQR(ZA^2#+ZB^2#)
  835. 50490 X=ZA/ZM
  836. 50500 ZT=(PI/2#-ATN(X/SQR(1-X*X)))*SGN(ZB)
  837. 50510 GOTO 50540
  838. 50520 ZM=ABS(ZA)
  839. 50530 ZT=PI/2#+PI/2#*SGN(-ZA)
  840. 50540 T=T+ZT
  841. 50550 P=P*ZM
  842. 50560 NEXT J
  843. 50570 RT(I)=T
  844. 50580 RA(I)=P
  845. 50590 NEXT I
  846. 50600 PRINT
  847. 50610 GOTO 50700
  848. 50620 PRINT USING"\             \";"POLE REAL";"POLE IMAG.";"RESIDUE MAG.";"RESIDUE ARG."
  849. 50630 PRINT
  850. 50640 FOR I=1 TO NP%
  851. 50650 RES(I)=K#/RM(I)
  852. 50660 PHI(I)=-RP(I)
  853. 50670 PRINT USING"##.####^^^^    ";PR(I),PJ(I),RES(I),PHI(I)
  854. 50680 NEXT I
  855. 50690 GOTO 50760
  856. 50700 PRINT USING"\             \";"POLE REAL";"POLE IMAG.";"RESIDUE MAG.";"RESIDUE ARG."
  857. 50710 FOR I=1 TO NP%
  858. 50720 RES(I)=K#*(RA(I)/RM(I))
  859. 50730 PHI(I)=RT(I)-RP(I)
  860. 50740 PRINT USING"##.####^^^^    ";PR(I),PJ(I),RES(I),PHI(I)
  861. 50750 NEXT I
  862. 50760 PRINT
  863. 50770 PRINT"PRESS ANY KEY TO CONTINUE"
  864. 50780 A$=INKEY$: IF A$="" THEN 50780
  865. 50782 IF NCURVES>1 AND ITERATION>0 THEN 50860
  866. 50790 PRINT
  867. 50800 PRINT"INPUT START TIME";
  868. 50810 INPUT TS
  869. 50820 PRINT"STOP TIME";
  870. 50830 INPUT TF
  871. 50850 TD=(TF-TS)/(NPP-1)
  872. 50860 CLS
  873. 50862 PRINT"COMPUTING TIME RESPONSE OF CURVE #";ITERATION+1
  874. 50890 FOR J= 0 TO (NPP-1)
  875. 50900 T=TS+J*TD
  876. 50910 RT=0
  877. 50920 FOR I=1 TO NP%
  878. 50930 IF PJ(I)=0 THEN 50960
  879. 50940 RR=2*RES(I)*EXP(PR(I)*T)*COS(PJ(I)*T+PHI(I))
  880. 50950 GOTO 50980
  881. 50960 RR=RES(I)*EXP(PR(I)*T)*COS(PHI(I))
  882. 50970 GOTO 50990
  883. 50980 I=I+1
  884. 50990 RT=RT+RR
  885. 51000 NEXT I
  886. 51010 XPASS!(J+1+NPP*ITERATION)=T:YPASS!(J+1+NPP*ITERATION)=RT
  887. 51020 NEXT J
  888. 51022 ITERATION=ITERATION+1
  889. 51024 IF NCURVES>ITERATION THEN 51050
  890. 51030 NCURVES!=NCURVES:NPP!=NPP
  891. 51040 CALL PLOT((NCURVES!),(NPP!),XPASS!(),YPASS!())
  892. 51047 NCURVES=1:ITERATION=0:NPP=1000:REM REINITIALIZE
  893. 51050 PRINT
  894. 51060 RETURN
  895.