home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
World of Shareware - Software Farm 2
/
wosw_2.zip
/
wosw_2
/
QBAS
/
CHAOSEXE.ZIP
/
XFFT.TRU
< prev
next >
Wrap
Text File
|
1980-01-01
|
9KB
|
275 lines
1000 !PROGRAM PEND/FFT.TRU
1001 CLEAR
1002 PRINT" ***PENDULUM - FAST FOURIER TRANSFORM OF VARIABLES***"
1003 PRINT
1004 PRINT"THIS PROGRAM PROVIDES A PHASE DIAGRAM, TIME SERIES, AND FFT OF THE"
1005 PRINT"ANGLE OR ANG. VELOCITY OF THE PENDULUM. A BLINKING CURSOR INDICATES"
1006 PRINT"THAT THE PROGRAM IS READY FOR THE NEXT STEP. PRESS ANY KEY TO CONTNUE."
1007 PRINT"THE HANNING OPTION IS USED TO SMOOTH THE EFFECT OF THE ABRUPT WINDOW"
1008 PRINT"AND IS RECOMMENDED IN MOST CASES."
1009 PRINT
1010 LIBRARY "SGLIB.TRC"
1020 !
1030 DIM thetadata(5000),thetadotdata(5000),xreal(0 to 5000),ximag(0 to 10000)
1040 DIM tpoint(0 to 5000),power(2048),frequency(2048)
1050 DECLARE DEF accel
1060 DECLARE DEF bitr
1070 INPUT prompt"Max frequency (try 0.5): ":maxfreq
1080 INPUT prompt"Input driving force strength: ":g
1090 INPUT prompt"Input damping term:":q
1100 INPUT prompt"Input initial position: ": xint
1110 INPUT prompt"Input initial velocity: ": vint
1120 INPUT Prompt"Input min.time:":tmin
1130 INPUT prompt"No. of FFT points(..256,512,1024,2048,4096) : ":number
1140 PRINT"Desired power spectrum quantity"
1150 PRINT" 1)Power spectrum of angular velocity"
1160 PRINT" 2)Power spectrum of angle"
1170 INPUT prompt"Choose 1 or 2 :":ps
1180
1190 LET del=.5/maxfreq
1200 LET tmax=number*del+tmin
1210 LET w=0.6666667
1220 LET eps=1.0e-6
1230 LET tstep=0.5
1240
1250 LET t=0
1260 LET x=xint
1270 LET v=vint
1280 LET sumvel=0
1290 LET count=0
1300 LET p=1
1310 CALL SETTEXT("PENDULUM PHASE DIAGRAM","ANGLE","ANGULAR VELOCITY")
1320 CALL RESERVELEGEND
1330 PRINT"computing data"
1340 !
1350 FOR i=1 to 10000000
1360 CALL rk4(x,v,tstep,xnew,vnew,t,w,g,q) ! Take a step of size tstep
1370 LET tshalf=tstep/2
1380 CALL rk4(x,v,tshalf,xnh,vnh,t,w,g,q) !Take two half steps
1390 CALL rk4(xnh,vnh,tshalf,xn,vn,t+tshalf,w,g,q)
1400 LET d1=abs(xn-xnew)
1410 LET d2=abs(vn-vnew)
1420 LET delta=max(d1,d2)
1430 IF delta<eps then
1440 IF t>tmin then
1450 LET tnew=t + tstep
1460 LET w1=mod(-t,del)
1470 LET w2=mod(tnew,del)
1480 IF w1<tstep then
1490 IF w2<tstep then
1500 LET ts=w1
1510 CALL rk4(x,v,ts,xp,vp,t,w,g,q)
1520 IF abs(xp)>pi then LET xp=xp-2*pi*abs(x)/x
1530 LET thetadata(p)=xp
1540 LET thetadotdata(p)=vp
1550 LET p=p+1
1560 LET sumvel= sumvel + v
1570 LET count = count + 1
1580 END IF
1590 END IF
1600 END IF
1610 LET x=xnew
1620 LET v=vnew
1630 LET t=t+tstep
1640 LET tstep=tstep*.95*(eps/delta)^.25
1650 IF abs(x)>pi then
1660 LET x=x-2*pi*abs(x)/x
1670 END IF
1680 ELSE
1690 LET tstep=tstep*.95*(eps/delta)^.2
1700 END IF
1710 IF t>tmax then LET i=10000001
1720 NEXT i
1730 LET n=p-1
1740 LET meanvel=sumvel/count
1750 CLEAR
1760 CALL setgraphtype("xy")
1770 CALL datagraph(thetadata,thetadotdata,1,0,"white")
1780 CALL ADDLEGEND("G="&STR$(G)&" Q="&str$(Q),0,1,"WHITE")
1790 CALL DRAWLEGEND
1800 GET KEY keyvariable
1810 !
1820 !
1830 !PREPARATION OF THE FFT DATA
1840 CLEAR
1850 INPUT prompt" HANNING OPTION Y/N? ": hanning$
1860 LET tgamma=log2(n)
1870 IF abs(int(tgamma)-tgamma)=0 then
1880 LET gamma=tgamma
1890 GOTO 1920
1900 END IF
1910 LET gamma=int(tgamma)+1
1920 PRINT "gamma= ";gamma
1930 LET newn=2^gamma
1940 LET nu=gamma
1950 FOR i=n+1 to newn
1960 LET xreal(i)=0
1970 NEXT i
1980 LET n=newn
1990 PRINT"n=";n
2000 CLEAR
2010 IF ps=1 then LET title$="ANGULAR VELOCITY"
2020 IF ps=2 then LET title$="ANGLE"
2030 CALL settext("TIME SERIES","TIME",title$)
2040 CALL setxscale(tmin,tmax)
2050 FOR k=0 to n-1
2060
2070 IF ps = 1 then
2080 LET xreal(k)=thetadotdata(k+1)
2090 ELSE IF ps = 2 then
2100 LET xreal(k)=thetadata(k+1)
2110 END IF
2120 IF hanning$="y" then LET xreal(k) = xreal(k)*(.5-.5*cos(2*pi*k/(n-1)))
2130 LET ximag(k)=0
2140 LET tpoint(k)=tmin+k*del
2150 NEXT k
2160 CALL setgraphtype("")
2170 CALL datagraph(tpoint,xreal,1,0,"white")
2180 GET KEY keyvariable
2190 FOR i= 1 to 100
2200 NEXT i
2210
2220 !FFT ALGORITHM
2230 CLEAR
2240 PRINT "Calculating FFT"
2250 LET n2=n/2
2260 LET nu1=nu-1
2270 LET k=0
2280 FOR l=1 to nu
2290 DO while k<(n-1)
2300 FOR i=1 to n2
2310 LET argument=k/2^nu1
2320 LET garbage=int(argument)
2330 LET p=bitr(garbage,nu)
2340 LET arg =2*pi*p/n
2350 LET c=cos(arg)
2360 LET s=sin(arg)
2370 LET k1=k+1
2380 LET k1n2=k1+n2
2390 LET treal=xreal(k1n2)*c+ximag(k1n2)*s
2400 LET timag=ximag(k1n2)*c-xreal(k1n2)*s
2410 LET xreal(k1n2)=xreal(k1)-treal
2420 LET ximag(k1n2)=ximag(k1)-timag
2430 LET xreal(k1)=xreal(k1)+treal
2440 LET ximag(k1)=ximag(k1)+timag
2450 LET k=k+1
2460 NEXT i
2470 LET k=k+n2
2480 LOOP
2490 LET k=0
2500 LET nu1=nu1-1
2510 LET n2=int(n2/2)
2520 NEXT l
2530
2540 FOR k=1 to n
2550 LET i=bitr(k-1,nu)+1
2560 IF i<=k then GOTO 2630
2570 LET treal=xreal(k)
2580 LET timag=ximag(k)
2590 LET xreal(k)=xreal(i)
2600 LET ximag(k)=ximag(i)
2610 LET xreal(i)=treal
2620 LET ximag(i)=timag
2630 NEXT k
2640
2650 !GRAPHING THE FFT
2660 CLEAR
2670 INPUT prompt"Plot 1)power spectrum, or 2)log power spectrum: ":pps
2680 INPUT prompt"Frequency variable - 1)linear, or 2)log: ":freqvar
2690 LET maxfreq=.5/del
2700 LET minfreq=1/(number*del)
2710
2720 CLEAR
2730 !Y-AXIS
2740 IF pps = 1 then
2750 LET TITLE$="POWER SPECTRUM"
2760 LET YAXIS$="POWER"
2770 ELSE
2780 LET TITLE$="LOG POWER SPECTRUM"
2790 LET YAXIS$="LOG POWER"
2800 END IF
2810 !X-AXIS
2820 IF freqvar=2 then
2830 LET XAXIS$="LOG FREQUENCY"
2840 ELSE
2850 LET XAXIS$="FREQUENCY"
2860 END IF
2870
2880 !DRAW AXES
2890 CLEAR
2900 CALL setxscale(minfreq,maxfreq)
2910 CALL setyscale(1e-6,.99)
2920 CALL SETTEXT(TITLE$,XAXIS$,YAXIS$)
2930 CALL RESERVELEGEND
2940
2950 !PLOT POINTS
2960 FOR i=1 to n/2
2970 LET frequency(i)=i/(n*del)
2980 LET power(i)=(((xreal(i))^2+(ximag(i))^2))/(n^2)
2990
3000 NEXT i
3010 !PLOT TEXT
3020 CALL setaxes(0)
3030 IF pps=1 then
3040 IF freqvar=1 then CALL setgraphtype("xy")
3050 IF freqvar=2 then CALL setgraphtype("logx")
3060 END IF
3070 IF pps=2 then
3080 IF freqvar=1 then CALL setgraphtype("logy")
3090 IF freqvar=2 then CALL setgraphtype("logxy")
3100 END IF
3102 IF NUMBER=4096 THEN
3110 CALL datagraph(frequency,power,1,1,"white")
3112 ELSE
3114 CALL DATAGRAPH(FREQUENCY,POWER,1,0,"WHITE")
3116 END IF
3120 CALL ADDLEGEND("N="&STR$(N)&" MAX FREQ ="&STR$(MAXFREQ)&" DEL F="&STR$(MINFREQ),0,1,"WHITE")
3130 CALL ADDLEGEND(" G="&STR$(G)&" Q="&STR$(Q),0,1,"WHITE")
3140 CALL drawlegend
3150 IF hanning$="y" then
3160 ! CALL ADDLEGEND(" HANNING",0,1,"WHITE")
3170 END IF
3180 GET KEY keyvariable
3190 INPUT PROMPT "Another with Hanning? y/n: ":hann$
3200 IF hann$= "y" THEN GOTO 1840
3210 INPUT PROMPT "FFT of different quantity (y/n)? :": diffquant$
3220 IF diffquant$ = "y" THEN
3230 IF ps = 2 THEN LET ps = 1
3240 IF ps = 1 THEN LET ps = 2
3250 GOTO 2000
3260 END IF
3270 INPUT PROMPT "Different presentation of same FFT? (y/n): ":diffplot$
3280 IF diffplot$ ="y" THEN GOTO 2660
3290 END
3300 !
3310 SUB rk4(x,v,tstep,xnew,vnew,t,w,g,q)
3320 DECLARE DEF accel
3330 LET xk1=tstep*v
3340 LET vk1=tstep*accel(x,v,t,w,g,q)
3350 LET xk2=tstep*(v+vk1/2)
3360 LET vk2=tstep*accel(x+xk1/2,v+vk1/2,t+tstep/2,w,g,q)
3370 LET xk3=tstep*(v+vk2/2)
3380 LET vk3=tstep*accel(x+xk2/2,v+vk2/2,t+tstep/2,w,g,q)
3390 LET xk4=tstep*(v+vk3)
3400 LET vk4=tstep*accel(x+xk3,v+vk3,t+tstep,w,g,q)
3410 LET vnew=v+(vk1+2*vk2+2*vk3+vk4)/6
3420 LET xnew=x+(xk1+2*xk2+2*xk3+xk4)/6
3430 END SUB
3440 !
3450 !
3460 DEF accel(x,v,t,w,g,q)
3470 LET accel=-sin(x)-(1/q)*v+g*cos(w*t)
3480 END DEF
3490
3500 !BIT REVERSER FUNCTION
3510 DEF bitr(j,nu)
3520 LET j1=j
3530 LET ibitr=0
3540 FOR i=1 to nu
3550 LET j2 = int(j1/2)
3560 LET ibitr=ibitr*2+(j1-2*j2)
3570 LET j1=j2
3580 NEXT i
3590 LET bitr=ibitr
3600 END DEF