10 REM  "CURWAVL/BAS"  FEEDBACK OF JET MAY 4, 1987
15 '   PLOTS JET CURRENT BELOW
20 DIM P1(64),P2(64),P(64),J1(64),M(64),X1(32),X2(32),A(66),B(32),C(10)
30 SCREEN 2 : KEY OFF :PI=3.141593
40 R%=32 : REM ROUND TRIP TIME IN TICKS
50 C1=1!:   'SATURATION JET CURRENT
60 G=2 : 'GAIN ON JET -RELATIVE
70 J%=10 :REM DELAY ON JET
80 S1=.00! : 'JET OFFSET, INITIAL
90 FOR F%=0 TO 10:READ C(F%):NEXT F%
100 DATA 0,-3.4,0,10,3.3,3,.8,1,.3,.4,.1 :REM IMITATES FLUTE
110  FOR F%=0 TO 36:P2(F%)=SIN(2*PI*F%/32):NEXT F%
120 CLS:LOCATE 25,1 : PRINT "PRESS SPACE BAR TO INTERRUPT"
130 FOR U%=37 TO 9999
140 T%=63 AND U%
150 V=G*M(63 AND (U%-J%))+S1
160 IF V>1 THEN W=1:GOTO 190  
170 IF V<-1 THEN W=0:GOTO 190  
180 W=1/2+3/4*V-V*V*V/4 :REM JET CHARACTERISTIC
190 J1(T%)=C1*W
200 S=0
210 FOR F%=0 TO 10:S=S+C(F%)*P2(63 AND (U%-R%+8-F%)):NEXT F%
220 P1(T%) =S/17.03
230 P2(T%)=P1(63 AND (U%-3))+.5*J1(T%)-.5*J1(63 AND (U%-3))
240 P(T%)=P1(T%)+P2(T%)
250 M3=P(T%)-J1(T%) : REM NEW DRIVING FUNCTION FOR MOUTH CURRENT
252 M(T%)=.99*M1 + .995*(M3-M2) : REM FILTERED MOUTH CURRENT
256 M1=M(T%) : REM SAVES PRESENT MOUTH CURRENT
258 M2=M3 : REM SAVES PRESENT MOUTH DRIVING FUNCTION
260 X%=X%+3:IF X%>599 THEN X%=X%-600:CLS
270 LOCATE 1,1:PRINT U%;"    ";
280 Y%=INT(80-4*P(T%)+.5)
290 PSET(X%,Y%):PSET(X%,150-12*J1(T%))
300       REM CAPTURE ZERO CROSSINGS U2% AND U1%
310 IF SGN(P(T%))= 1 AND SGN(P((U%-1) AND 63)) =-1 THEN U3%=U2%:U2%=U1%: U1%=U%
320 LOCATE 1,8
325 PRINT "                          ZERO CROSSINGS AT ";U2%;" AND ";U1%;
330 R$=INKEY$:IF R$=" " THEN GOSUB 350
340 NEXT U%
350 PRINT "DO YOU WANT TO: (1) VIEW DATA   (2) GET FOURIER ANALYSIS                            (3) CONTINUE CALCULATION    (4) PRINT OUT POINTS    (5) QUIT
360 INPUT CH
370 IF CH=1 THEN GOSUB 650: GOTO 350
380 IF CH=2 THEN GOSUB 540 :GOTO 350
390 IF CH=3 THEN LOCATE 1,1:PRINT STRING$(160," "):RETURN
395 IF CH=4 THEN GOSUB 410 : GOTO 350
400 IF CH=5 THEN END ELSE GOTO 350
410 LPRINT:LPRINT"TUBE DELAY=";R%;"   JET DELAY=";J%;"   GAIN  =";G;"   CURRENT=";C1;"   OFFSET=";S1
420 LPRINT
430 LPRINT "TIME  I/PIPE   I/JET   I/MOUTH"
440 A$="####  ##.###  ##.###  ##.###"
450 B$="     !    !    !    !    I    !    !    !    !    !"
460 FOR F%=40 TO 0 STEP -1
470 Z%= (U%-1-F%) AND 63
480 H=25+INT(5*P(Z%)+.5)
490 LPRINT USING A$;U%-1-F%;P(Z%);J1(Z%);M(Z%);
500 IF H<1 OR H>50 THEN LPRINT B$:GOTO 520  
510 LPRINT LEFT$(B$,H);"*";RIGHT$(B$,50-H)
520 NEXT F%
530 LPRINT:LPRINT:RETURN
540 REM    SUBR TO ESTIMATE  EXACT ZERO CROSSINGS
542 INPUT "PERIOD FACTOR (1 0R 2)";PF
544 IF PF=1 GOTO 550 ELSE IF PF=2 THEN U2%=U3%:GOTO 550
545 GOTO 542
550 D%=U2% :GOSUB 610
560 X2=X : D%=U1% : GOSUB 610
570 PE= U1%+X-U2%-X2 : PRINT "PERIOD =";PE;"TIME UNITS":PRINT
580 INPUT " (1) RETURN TO MENU   (2) CONTINUE FOURIER";CO
590 IF CO=1 THEN RETURN
600 IF CO=2 THEN GOTO 700  ELSE GOTO 580
610 Y2=P(D% AND 63): Y1=P((D%-1) AND 63): Y3=P((D%+1) AND 63)
620 A=(Y3+Y1)/2-Y2 : B=(Y3-Y1)/2 : C=Y2
630 X=(-B+SQR(B^2-4*A*C))/2/A
640 RETURN
650 CLS:PRINT "TIME    I/PIPE   I/JET   I/MOUTH        TIME   I/PIPE   I/JET   I/MOUTH":PRINT
655 B$="####   ##.###   ##.###   ##.###"
657 FOR F% = 40 TO 20 STEP -1
660 Z%= (U%-1-F%) AND 63
670 PRINT USING B$;U%-1-F%,P(Z%);J1(Z%);M(Z%);:PRINT "         ";
671 Y%=(U%+20-F%) AND 63:PRINT USING B$;U%+20-F%;P(Y%);J1(Y%);M(Y%)
684 NEXT F%
686 LOCATE 25,1:PRINT "PRESS SPACE BAR TO RETURN TO MENU";
687 IF INKEY$<>" " THEN GOTO 687
690 RETURN
700 REM ** CONVERTS TO 32 POINT EXACT CYCLE FOR FOURIER ANALYSIS
710 NP=U1%-U2% +2
720 REM   LENGTH OF PERIOD IN POINT SPACES IS PE FROM LINE 570
730 ' PUT IN Y COORDINATES FOR CURVE, STARTING WITH POINT 0
740 REM
750 FOR F%=0 TO NP : A(F%)=P((U2%+F%) AND 63) : NEXT F%
760  X1(0)=A(0)
770 FOR I=1 TO 31
780 DP=I*PE/32
790 DC=INT(DP+.5)
800 IF DC=0 THEN DP=DP+PE:GOTO 790
810 XD=DP-DC
820 X1(I)=(A(DC+1)+A(DC-1)-2*A(DC))/2*XD*XD+(A(DC+1)-A(DC-1))/2*XD+A(DC)
830 NEXT I
840 REM FOURIERA - NUMERICAL VALUES ARE COEFFICIENTS FOR EXPONENTIONAL NOTATION;                   MULTIPLY BY 2 (EXCEPT DC TERM) FOR R*COS + I*SIN
850 N=32
860 PI=3.141593 : L= 5
870 REM TIME FUNCTION X1(Z) WHERE Z RUNS FROM 0 TO N-1
880 SCREEN 2:WIDTH 80:CLS
890 FOR Z=0 TO 31:P=INT(Z/8):Q=Z-8*P
900 D=20*P+1: LOCATE Q +4,D:PRINT Z;: LOCATE Q+4,D+4 :PRINT X1(Z);
910 NEXT Z
920 PRINT
930 B=X1(0)
940 FOR Z=0 TO 31
950 IF ABS(X1(Z))>B THEN B=ABS(X1(Z))
960 NEXT Z
970  VIEW (0, 98)-(300,190):WINDOW (0,-B)-(32,B)
980 LINE (0,0)-(32,0)
990 FOR Z =0 TO 31
1000 PSET(Z,X1(Z)):IF X1(Z)=0 THEN PRESET (Z,X1(Z))
1010 NEXT Z
1030 LOCATE 24,40 :INPUT "READY";G$
1040 REM- SCALE TIME FUNCTION
1050 FOR Z=0 TO N-1
1060 X1(Z)=X1(Z)/N
1070 NEXT Z
1080 FOR I=0 TO 31:X2(I)=0: NEXT I
1090 REM - FFT IN PLACE ALGORITHM--
1100 PRINT "FFT CALCULATION IN PROGRESS"
1110 I1=N/2:I2=1:V=2*PI/N
1120 FOR I=1 TO L
1130 I3=0:I4=I1
1140 FOR K=1 TO I2
1150 X=INT(I3/I1)
1160 GOSUB 1660
1170 I5=Y
1180 Z1=COS(V*I5)
1190 Z2=-SIN(V*I5)
1200 FOR M=I3 TO I4-1
1210 A1=X1(M):A2=X2(M)
1220 B1=Z1*X1(M+I1)-Z2*X2(M+I1)
1230 B2=Z2*X1(M+I1)+Z1*X2(M+I1)
1240 X1(M)=A1+B1:X2(M)=A2+B2
1250 X1(M+I1)=A1-B1:X2(M+I1)=A2-B2
1260 NEXT M
1270 I3=I3+2*I1:I4=I4+2*I1
1280 NEXT K
1290 I1=I1/2:I2=2*I2
1300 NEXT I
1310 REM OUTPUT RESULTS
1320  REM --OUTPUT MAGNITUDE SPECTRUM PLOT
1330 B=0 : GOSUB 1800
1340 REM "CALCULATIONS IN PROGRESS-----"
1350 FOR Z=0 TO N/2
1360 X=Z
1370 GOSUB 1750
1380 IF X3>B THEN B=X3
1390 NEXT Z
1400 FOR Z= 0 TO N/2
1410 X=Z
1420 GOSUB 1750
1430 X4=X3/B*100 : IF Z=0 THEN X4=X4/2
1440 LINE (Z,0)-(.3+Z,X4),,BF
1450 NEXT Z
1460 LOCATE 25,1 : PRINT "PRESS ANY KEY WHEN READY";
1470 IF INKEY$="" THEN GOTO 1470
1480 REM-OUTPUT TABLE OF VALUES-
1490 U=0: Z=0 : VIEW: WINDOW: CLS
1500 PRINT"TUBE DELAY=";R%;"   JET DELAY=";J%;"   GAIN  =";G;"   CURRENT=";C1;"   OFFSET=";S1:PRINT
1510  PRINT"HARMONIC";TAB(14);"REAL";TAB(30);
1520 PRINT"IMAGINARY";TAB(50);"MAGNITUDE";TAB(65);"ANGLE"
1530 X=U
1540 GOSUB 1750
1550 IF X1(Y)=0 THEN AN=90*SGN(X2(Y)):GOTO 1570
1560 AN=57.296*ATN(X2(Y)/X1(Y)):IF X1(Y)<0 THEN AN=AN+180
1570 IF U=0 THEN PRINT U;TAB(10);X1(Y);TAB(30);X2(Y);TAB(50);X3;TAB(65);AN:GOTO 1580
1575 PRINT U;TAB(10);2*X1(Y);TAB(30);2*X2(Y);TAB(50);2*X3;TAB(65);AN
1580 U=U+1:Z=Z+1
1590 IF U>N/2 THEN GOTO 1610
1600 GOTO 1530
1610 PRINT: PRINT "TIME UNITS FROM";U2%;"TO";U1%;"    PERIOD = ";PE
1620 LOCATE 25,1:PRINT "   PRESS SPACE BAR TO RETURN TO MENU";
1630 IF INKEY$=" " GOTO 1640 ELSE GOTO 1630
1640 CLS:RETURN
1650 REM --SCRAMBLER SUBROUTINE--
1660 Y=0:N1=N
1670 FOR W=1 TO L
1680 N1=N1/2
1690  IF X<N1 THEN 1720
1700  Y=Y+2^(W-1)
1710 X=X-N1
1720 NEXT W
1730 RETURN
1740 REM--MAGNITUDE SUBROUTINE--
1750 GOSUB 1660
1760 X3=SQR(X1(Y)^2+X2(Y)^2)
1770 RETURN
1780 END
1790 REM ** MAGNITUDE PLOT ROUTINE **
1800 SCREEN 1 : WIDTH 80 : KEY OFF
1810 FOR I=0 TO 16:LOCATE 23,4*I+4:PRINT I; :NEXT I
1820 FOR I=0 TO 10:LOCATE 22-2*I,1 :PRINT 10*I;:NEXT I
1830 LINE (31,12)-(576,172),,B: VIEW(31,12)-(576,172)
1840 WINDOW(0,0)-(17,100)
1850 FOR I=0 TO 10
1860 LINE (0,10*I)-(17,10*I),,,&H8080 :NEXT I
1870 RETURN
1880 END