      PROGRAM SIMP1
C  THIS EVALUATES THE INTEGRALS FOR FINDING A(h) AND B(h) BY
C  SIMPSON'S RULE.   BECAUSE OF THE PERIODIC NATURE OF f(x) THE
C  SUMMATION IS DONE OVER THE POINTS 0 TO N-1 (NOT N) AND THE 
C  POINT 0 IS GIVEN WEIGHT 2 RATHER THAN 1 AS USUAL.
      REAL F(0:199),A(0:50),B(0:50)
      INTEGER H,HMAX
      CHARACTER*4 HC,AC,BC,ANS*1
      HC='h'
      AC='A(h)'
      BC='B(h)'
      OPEN(UNIT=9,FILE='LPT1')
      TPI=6.2831853
      WRITE(6,'('' INPUT THE INTERVAL, 1/N, FOR QUADRATURE AS AN'')')
      WRITE(6,'('' EVEN INTEGER N <= 200 '')')
      READ(5,*)NN
      WRITE(6,'('' INPUT N VALUES OF f(x) FOR x = 0 TO 1-1/N '')')
      WRITE(6,'('' IN STEPS OF 1/N. IF YOU MAKE AN ERROR THERE  '')')
      WRITE(6,'('' WILL BE AN OPPORTUNITY TO CORRECT IT LATER. '')')
      DO 1 I=0,NN-1
      WRITE(6,100)I
  100 FORMAT(11H INPUT f{X[,I3,1H])
      READ(5,*)F(I)
    1 CONTINUE
    3 WRITE(6,'('' DO YOU WANT TO CHECK THE VALUES OF f(x)?  Y/N '')')
      READ(5,25)ANS
   25 FORMAT(A1)
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 2
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')GOTO 4
      GOTO 3
    4 DO 5 I=0,NN-1
      WRITE(6,200)I,NN,F(I)
  200 FORMAT(4H X =,I4,1H/,I3,7H f(x) =,F6.2)
    7 WRITE(6,'('' DO YOU WANT TO CHANGE THIS? Y/N '')')
      READ(5,25)ANS
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 5
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')GOTO 6
      GOTO 7
    6 WRITE(6,'('' INPUT REPLACEMENT f(x) '')')
      READ(5,*)F(I)
    5 CONTINUE
    2 WRITE(6,'('' INPUT HMAX FOR THE CALCULATION <=50] '')')
      READ(5,*)HMAX
C  THE VALUES OF A(h) AND B(h) ARE NOW TO BE CALCULATED
      DO 10 H=0,HMAX
      AA=0
      BB=0
      W=4
C  SUMMATIONS FOR A(h) AND B(h)
      DO 12 J=0,NN-1
      W=6-W
      ARG=TPI*H*J/NN
      AA=AA+W*F(J)*COS(ARG)
      BB=BB+W*F(J)*SIN(ARG)
   12 CONTINUE
      A(H)=AA/3.0/NN
      B(H)=BB/3.0/NN
   10 CONTINUE
C  OUTPUT RESULTS
      WRITE(9,150)HC,AC,BC,HC,AC,BC,HC,AC,BC,HC,AC,BC
  150 FORMAT(4(4X,A4,A4,3X,A4))
      WRITE(9,250)(I,A(I),B(I),I=0,HMAX)
  250 FORMAT(4(I5,2F7.2))
      STOP
      END


