PROGRAM FOURIE C C FOURIER SYNTHESIS OF WAVEFORMS C C ----------------------------------------------- C AUTHOR: ROBERT EHRLICH (GEORGE MASON U 1973) C Revised: Peter Signell (Mich. St. U. 3/86) C Revised: Gary Herzenstiel (Mich. St. U. 3/86) C Revised: E.J.D.Kales (Mich. St. U. 11/87) C Revised: E.J.D.Kales (Mich. St. U. 11/97) C Revised: Peter Signell (Mich. St. U. 3/02) C C Does a Fourier synthesis of four functions: adds C sine and cosine waves, each multiplied by an a C coefficient, and plots the resultant function. The C four functions are SQUARE WAVE, TRIANGLE, SPIKE, and C PROFILE. Calculates and plots average deviation when C the number of waves is 100 and type is not PROFILE. C C SUBROUTINES: C GETINF - get values for R,L,C,DT,N,V0,W C POINTS - loop over 51 points, put in graph array, C plot. C WAVES - loop over number of waves to be added C into the sum C QUERY - ask if user wants to run the program C again C C F O U R I E R C : C --------------------------------- C : : : C GETINF POINTS QUERY C (N1,TYPE) (N1,TYPE) (AGAIN) C : C WAVES C (N1,TYPE,J,PI,X,DX,FN,DS) C C ***** DECLARATIONS ***** REAL N INTEGER N1,TYPE CHARACTER Z*101,AGAIN C OPEN(UNIT=3,FILE='m352p1f.out') OPEN(UNIT=3,FILE='(T:84)m352p1f.out') CC OPEN(UNIT=3,FILE='(T:120)m352p1f.out') C C ***** M A I N B O D Y ***** C 10 WRITE(3,'(A1)') 12 CALL GETINF ( N1,TYPE) CALL POINTS (N1,TYPE ) CALL QUERY ( AGAIN) IF(AGAIN.EQ.'Y') GOTO 10 STOP END C C ***** G E T I N F ***** C SUBROUTINE GETINF( N1,TYPE) C C get values for N, & TYPE. N1 is the integer value of C N. C C ***** DECLARATIONS ***** REAL N INTEGER N1,TYPE C C type of calculation = 0,1,2,3 C 0 = SQUARE WAVE 1 = TRIANGLE C 2 = SPIKE 3 = PROFILE C N = number of time steps C PRINT *,'Please supply N & TYPE values:' WRITE(*,20) ' N=' READ *,N WRITE(*,20) 'TYPE=' READ *,TYPE 20 FORMAT(A5) N1=N WRITE(*,34) N1,TYPE WRITE(3,34) N1,TYPE 34 FORMAT('m352p1.for Fourier Synthesis of Waveforms'// +' Number of waves=',I4/' Wave Type=',I2/) RETURN END C C ***** P O I N T S ***** C SUBROUTINE POINTS(N1,TYPE ) C C loops over 51 points C C ***** DECLARATIONS ***** REAL PI,X,DX,FN,DS(100) INTEGER N1,TYPE,J,M CHARACTER Z*101 C C there are 51 points per cycle, so DX = 2*PI/50 C PI=3.14159265 X=-PI DX=2.0*PI/50.0 WRITE(3,40) 40 FORMAT(' X F(X)') WRITE(3,41) 41 FORMAT(19X,41H........................................., +24H........................) CC note: use 60H for "large" C DO 48,J=1,51 C initialize Z and set boundaries: Z='. ' CC small or large: C +' .' +//' .' CC + //' .' IF(TYPE.NE.3) Z(33:33)='|' CC IF(TYPE.NE.3) Z(51:51)='|' C FN=0.0 C the sum of N waves initially set to zero C CALL WAVES(N1,TYPE,J,PI,X,DX,FN,DS ) C C compute M in range M=1..80(or 100) to correspond C to value of resultant amplitude, store in FN. C Note that when FN=0, M=41 (51) C CC small or large: M=16.0*(FN+2.06251) CC M=25.0*(FN+2.06) C FOR PROFILE, MAKE HORIZ SCALE SAME AS IN ORIGINAL CURVE IF (TYPE.EQ.3) M = 30.0*FN + 32 C Z(M:M)='*' C store '*' at position M in the graph C C display graph one line at a time WRITE(3,45) X,FN,Z(1:80) CC small or large: 45 FORMAT(F8.4,F9.4,2X,A65) CC 45 FORMAT(F8.4,F9.4,2X,A101) C X=X+DX C increment by DX, calculate and plot FN at next point 48 CONTINUE WRITE(3,41) C C give deviation plot only when N=100 C IF(N1.EQ.100) THEN WRITE(3,'(A1//)') 12 WRITE(3,51)TYPE 51 FORMAT( + 'Average Deviations for N=2,4,...100 Waves of Type',I2// + ' N DEV') WRITE(3,41) DO 55,J=2,100,2 M=60.0*DS(J)/DS(2)+1.0 Z='. ' + //' ' Z(M:M)='*' WRITE(3,53) J,DS(J),Z(1:64),'.' CC small or large: 53 FORMAT(I4,F8.4,7X,A64,A1) CC WRITE(3,53) J,DS(J),Z,'.' CC 53 FORMAT(I4,F8.4,7X,A101,A1) 55 CONTINUE WRITE(3,41) ENDIF RETURN END C C ***** W A V E S ***** C SUBROUTINE WAVES(N1,TYPE,J,PI,X,DX,FN,DS ) C C loop over the number of waves to be added in the sum C C ***** DECLARATIONS ***** INTEGER N1,TYPE,J,K REAL PI,X,DX,FN,DS(100) REAL A(48),B(48),XK,AK,BK,F,C,N C DATA A/ .1322, .0138, .0816,-.0360, .0122, + -.0379,-.0046,-.0147, .0089, .0036, + .0084, .0050,-.0027,-.0015,-.0055, + .0013,-.0042, .0050,-.0017, .0030, + -.0041, .0010,-.0047,-.0002,-.0019, + .0020, .0001, .0017, .0002,-.0011, + -.0016,-.0021,-.0001,-.0008, .0023, + .0002, .0015,-.0016, .0003,-.0016, + .0001,-.0002, .0013, .0004, .0007, + -.0003,-.0006,-.0009/ DATA B/-.1534,-.0468, .0777,-.0440, .0593, + -.0544, .0387,-.0457, .0348,-.0322, + .0381,-.0229, .0306,-.0210, .0205, + -.0246, .0133,-.0202, .0127,-.0130, + .0142,-.0103, .0099,-.0121, .0066, + -.0117, .0075,-.0081, .0088,-.0081, + .0079,-.0090, .0065,-.0093, .0081, + -.0081, .0093,-.0074, .0089,-.0092, + .0076,-.0098, .0074,-.0088, .0091, + -.0077, .0084,-.0088/ C DO 60,K=1,N1 IF(J.EQ.1) DS(K)=0.0 XK=K C C AK and BK are the Fourier coefficients C AK=0.0 BK=0.0 C IF(TYPE.EQ.0) THEN C compute coefficients for the square wave. C F is the theoretical value of the function, C X. C IF(J.EQ.1.OR.J.EQ.51) THEN F=ABS(X)/X ELSE F=0.0 ENDIF C C factor C makes BK zero for even values of K. C C=K-2*(K/2) BK=4.0*C/(PI*XK) C ELSEIF(TYPE.EQ.1) THEN C compute coefficients for triangle. C factor C makes AK zero for even values of K. C F=1.0-2.0*ABS(X)/PI C=K-2*(K/2) AK=8.0*C/(PI*XK)**2 C ELSEIF(TYPE.EQ.2) THEN C compute coefficients for spike. C N=N1 AK=1.0/N IF(J.EQ.26) THEN F=1.0 ELSE F=0.0 ENDIF C C (FLAG = 3) ELSEIF(K.LE.48) THEN C compute coefficients for profile. C *NOTE* can't handle more than 48 waves in C this case. C AK=A(K) BK=B(K) C ENDIF C C add this wave to the sum C FN=FN+AK*COS(XK*X)+BK*SIN(XK*X) C C compute average deviation from theoretical C function. C IF(N1.EQ.100) DS(K)=DS(K)+ABS(F-FN)/51.0 60 CONTINUE RETURN END C C ***** Q U E R Y ***** C SUBROUTINE QUERY( AGAIN) CHARACTER AGAIN PRINT * WRITE(*,62) 62 FORMAT( +'Would you like to do this again using new values? +(Y/N)? ') READ( *,'(A1)')AGAIN RETURN END C ***** P R O G R A M E N D *****