PROGRAM analyze
*     determine the Fourier coefficients a_k and b_k
      IMPLICIT NONE
      REAL*8 delta, period
      INTEGER N, nmax, IC1_, IC2_, IR_
      DATA IC1_/1/
      DATA IC2_/2/
      CALL GWopen(IR_, 0)
      CALL parameters(N,nmax,delta,period)
      CALL screen(nmax,period,IC1_,IC2_)
      CALL coefficents(N,nmax,delta,period,IC1_,IC2_)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE parameters(N,nmax,delta,period)
      IMPLICIT NONE
      REAL*8 delta, period
      INTEGER N, nmax
      WRITE(*,'(A,$)') 'number of data points N (even) = '
      READ(*,*) N
      WRITE(*,'(A,$)') 'sampling time dt = '
      READ(*,*) delta
      period = N*delta                      ! assumed period
*     maximum value of mode corresponding to Nyquist frequency
      nmax = N/2
      END
C
      SUBROUTINE screen(nmax,period,IC1_,IC2_)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 pi, ymax, period
      REAL ticksize, GWaspect
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER nmax, IR_, IC1_, IC2_
      ymax = 2
      ticksize = ymax/50
      CALL GWvport(IR_, 0.0, 0.5, GWaspect(1), 1.0)
      WRITE(Stmp1_, '(A,F7.5)') 
     +      '     a_k        frequency interval =', 2*pi/period
      CALL GWindow(IR_, -1.0, REAL(-ymax), REAL(nmax+1), REAL(ymax))
      CALL GWputtxt(IR_, 0.0, REAL(ymax*0.7D0), Stmp1_)
      CALL plotaxis(nmax,ticksize)
      CALL GWsetpen(IR_, 13, -1, -1, -1)
      CALL GWsavevp(IR_, IC1_)
      CALL GWreset(IR_)
      CALL GWvport(IR_, 0.0, 0.0, GWaspect(1), 0.5)
      WRITE(Stmp1_, '(A)') '     b_k'
      CALL GWindow(IR_, -1.0, REAL(-ymax), REAL(nmax+1), REAL(ymax))
      CALL GWputtxt(IR_, 0.0, REAL(ymax*0.7D0), Stmp1_)
      CALL plotaxis(nmax,ticksize)
      CALL GWsetpen(IR_, 13, -1, -1, -1)
      CALL GWsavevp(IR_, IC2_)
      END
C
      SUBROUTINE plotaxis(nmax,ticksize)
      IMPLICIT NONE
      REAL ticksize
      INTEGER k, nmax, IR_
      CALL GWline(IR_, 0.0, 0.0, REAL(nmax), 0.0)
      DO k = 1, nmax
         CALL GWline(IR_, REAL(k), -ticksize, REAL(k), ticksize)
      END DO
      END
C
      SUBROUTINE coefficents(N,nmax,delta,period,IC1_,IC2_)
      IMPLICIT NONE
      REAL*8 ak, bk, pi, t, wk, delta, period, f
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER IR_, i, k, N, nmax, IC1_, IC2_
      DO k = 0, nmax
         ak = 0
         bk = 0
         wk = 2*pi*k/period
*        rectangular approximation
         DO i = 0, N - 1
            t = i*delta
            ak = ak + f(t)*cos(wk*t)
            bk = bk + f(t)*sin(wk*t)
         END DO
         ak = 2*ak/N
         bk = 2*bk/N
         CALL GWselvp(IR_, IC1_)
         CALL GWline(IR_, REAL(k), 0.0, REAL(k), REAL(ak))
         CALL GWselvp(IR_, IC2_)
         CALL GWline(IR_, REAL(k), 0.0, REAL(k), REAL(bk))
      END DO
      END
C
      FUNCTION f(t)
      IMPLICIT NONE
      REAL*8 f, pi, w0, t
      PARAMETER(pi = 3.141592653589793D0)
      w0 = 0.1D0*pi
      f = sin(w0*t)                         ! simple example
      END