      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 *****