      PROGRAM synthesize
      IMPLICIT NONE
      INTEGER IR_
      CALL GWopen(IR_, 0)
      CALL plotf(0.0,0.5,0.5,1.0)
      CALL plotf(0.0,0.5,0.0,0.5)
      CALL plotf(0.5,1.0,0.5,1.0)
      CALL plotf(0.5,1.0,0.0,0.5)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE plotf(xmin,xmax,ymin,ymax)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL pi, xmin, xmax, ymin, ymax, GWaspect
      PARAMETER(pi = 3.141592)
      INTEGER IR_, N
      CALL GWvport(IR_, xmin*GWaspect(1), ymin, xmax*GWaspect(1), ymax)
      CALL GWindow(IR_, -4.0, -2.0, 4.0, 2.0)
      CALL GWsetpen(IR_, 0, -1, -1, -1)
      CALL GWline(IR_, -pi, 0.0, pi, 0.0)
      CALL GWline(IR_, 0.0, -1.5, 0.0, 1.5)
      CALL GWinput(IR_, 'number of modes = ', Stmp1_)
      READ(Stmp1_,*) N
      CALL GWsetpen(IR_, -1, -1, -1, -1)
      CALL GWsetpen(IR_, 13, -1, -1, -1)
      WRITE(Stmp1_,*) 'number of modes =', N
      CALL GWputtxt(IR_, -pi, 1.5, Stmp1_);
      CALL fourier(N)
      END
C
      SUBROUTINE fourier(N)
*     compute Fourier series and plot function
      IMPLICIT NONE
      REAL*8 a(0:1000), b(1000), dt, f, pi, t
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER IR_, N, k, nplot
      CALL coefficients(N,a,b)
      nplot = 100
      t = -pi
      dt = pi/100
      DO WHILE(t .LE. pi)
         f = a(0)/2
         DO k = 1, N
            IF(a(k) .NE. 0) f = f + a(k)*cos(k*t)
            IF(b(k) .NE. 0) f = f + b(k)*sin(k*t)
         END DO
         CALL GWline2(IR_, REAL(t), REAL(f))
         t = t + dt
      END DO
      END
C
      SUBROUTINE coefficients(N,a,b)
*     generate Fourier coefficients for special case 
      IMPLICIT NONE
      REAL*8 pi, a(0:1000), b(1000)
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER N, k
      a(0) = 0
      DO k = 1, N
         a(k) = 0
         IF(mod(k,2) .NE. 0) THEN
            b(k) = 2/(k*pi)
         ELSE
            b(k) = 0
         END IF
      END DO
      END
