      PROGRAM FOUR1
C  THIS PROGRAMME COMPUTES A ONE-DIMENSIONAL FOURIER SERIES.
C  THE INPUT COEFFICIENTS CAN BE IN THE FORM EITHER A(h), B(h)
C  OR |F(h)|, PHI(h).
C  THE OUTPUT RANGE FOR X/A CAN BE EITHER (-1/2,1/2) OR (0,1).
      DIMENSION X(0:100),Y(0:100),IC(0:200),A(0:50),B(0:50),Z(0:200)
      INTEGER H,HMAX
      CHARACTER*4 NC,FX
      NC='n'
      FX='f(x)'
      OPEN(UNIT=9,FILE='LPT1')
      TPI=6.2831853
C  CONSTANT TO CONVERT DEGREES TO RADIANS
      DTR=0.0174533
      WRITE(6,'('' IF DATA IS IN FORM A &amp; B THEN INPUT 0 '')')
      WRITE(6,'('' IF IN THE FORM F &amp; PHI THEN INPUT 1 '')')
   13 READ(5,*)KEY
      IF(KEY.EQ.0.OR.KEY.EQ.1)GOTO 12
      GOTO 13
   12 WRITE(6,'('' INDICATE INTERVAL OF X/A AT WHICH FUNCTION IS'')')
      WRITE(6,'('' CALCULATED,1/N, BY INPUT OF EVEN INTEGER N <=200'')') 
      READ(5,*)N
      WRITE(6,'('' FOR OUTPUT OF X/A IN RANGE (-1/2,1/2) INPUT 0'')')
      WRITE(6,'('' FOR THE RANGE (0,1) INPUT 1'')')
    7 READ(5,*)IOUT
      IF(IOUT.EQ.0.OR.IOUT.EQ.1)GOTO 11
      GOTO 7
   11 WRITE(6,'('' INPUT h(+ve) FOLLOWED BY EITHER A(h),B(h) OR '')')
      WRITE(6,'('' F(h),PHI(h) [in degrees].   VALUES OF A, B'')')
      WRITE(6,'('' AND F MUST ALL BE < 1000.  FOR THE TERM OF'')')
      WRITE(6,'('' INDEX 0 INPUT  0, A [OR F], 0.  THE MAXIMUM'')')
      WRITE(6,'('' VALUE OF h IS 50.  TO TERMINATE DATA INPUT'')')
      WRITE(6,'('' H = 100. '')')
C  CLEAR A, B ARRAYS
      DO 50 I=0,50
      A(I)=0  
   50 B(I)=0
C  READ IN DATA. CONVERT TO A,B FORM IF NECESSARY. ALSO FIND HMAX.
      HMAX=0
    3 READ(5,*)H
      IF(H.EQ.100)GOTO 4
      READ(5,*)XX,YY
      IF(KEY.EQ.0)THEN
      A(H)=XX
      B(H)=YY
      GOTO 2
      ENDIF
      A(H)=XX*COS(DTR*YY)
      B(H)=XX*SIN(DTR*YY)
    2 IF(H.GT.HMAX)HMAX=H      
      GOTO 3
C  START THE SUMMATION.   SYMMETRY ABOUT THE ORIGIN IS USED BY DOING
C  SEPARATE SUMMATIONS FOR THE COSINE AND SINE TERMS FOR X/A FROM
C  0 TO 1/2 AND THEN COMBINING THEM.
    4 DO 5 I=0,N/2
      X(I)=A(0)
      Y(I)=0
      DO 6 J=1,HMAX
      ARG=TPI*FLOAT(J*I)/N
      X(I)=X(I)+2*A(J)*COS(ARG)
    6 Y(I)=Y(I)+2*B(J)*SIN(ARG)
    5 CONTINUE
      DO 8 I=0,N/2
      IF(IOUT.EQ.0)THEN
      IC(I)=I-N/2
      Z(I)=X(N/2-I)-Y(N/2-I)
      IC(N-I)=-IC(I)
      Z(N-1)=X(N/2-I)+Y(N/2-I)
      GOTO 8
      ENDIF
      IC(I)=I
      Z(I)=X(I)+Y(I)
      IC(N-I)=N-I
      Z(N-I)=X(I)-Y(I)
    8 CONTINUE
C  CALCULATION COMPLETE.   OUTPUT STAGE BEGINS.
      WRITE(9,100)N
  100 FORMAT(41H THE COORDINATES ARE X/A = n/N WHERE N = ,I3)
      WRITE(9,'('' '')')
      WRITE(9,150)NC,FX,NC,FX,NC,FX,NC,FX
  150 FORMAT(4(6X,A4,1X,A4))
      WRITE(9,200)(IC(I),Z(I),I=0,N)
  200 FORMAT(4(I8,F7.2))
      STOP
      END
    

