C
C     ..................................................................
C
C        SUBROUTINE FORIF
C
C        PURPOSE
C           FOURIER ANALYSIS OF A GIVEN PERIODIC FUNCTION IN THE
C           RANGE 0-2PI
C           COMPUTES THE COEFFICIENTS OF THE DESIRED NUMBER OF TERMS
C           IN THE FOURIER SERIES F(X)=A(0)+SUM(A(K)COS KX+B(K)SIN KX)
C           WHERE K=1,2,...,M TO APPROXIMATE THE COMPUTED VALUES OF A
C           GIVEN FUNCTION SUBPROGRAM
C
C        USAGE
C           CALL FORIF(FUN,N,M,A,B,IER)
C
C        DESCRIPTION OF PARAMETERS
C           FUN-NAME OF FUNCTION SUBPROGRAM TO BE USED FOR COMPUTING
C               DATA POINTS
C           N  -DEFINES THE INTERVAL SUCH THAT 2N+1 POINTS ARE TAKEN
C               OVER THE INTERVAL (0,2PI). THE SPACING IS THUS 2PI/2N+1
C           M  -THE MAXIMUM ORDER OF THE HARMONICS TO BE FITTED
C           A  -RESULTANT VECTOR OF FOURIER COSINE COEFFICIENTS OF
C               LENGTH M+1
C               A SUB 0, A SUB 1,..., A SUB M
C           B  -RESULTANT VECTOR OF FOURIER SINE COEFFICIENTS OF
C               LENGTH M+1
C               B SUB 0, B SUB 1,..., B SUB M
C           IER-RESULTANT ERROR CODE WHERE
C               IER=0  NO ERROR
C               IER=1  N NOT GREATER OR EQUAL TO M
C               IER=2  M LESS THAN 0
C
C        REMARKS
C           M MUST BE GREATER THAN OR EQUAL TO ZERO
C           N MUST BE GREATER THAN OR EQUAL TO M
C           THE FIRST ELEMENT IN VECTOR B IS ZERO IN ALL CASES
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           FUN-NAME OF USER FUNCTION SUBPROGRAM USED FOR COMPUTING
C               DATA POINTS
C           CALLING PROGRAM MUST HAVE FORTRAN EXTERNAL STATEMENT
C           CONTAINING NAMES OF FUNCTION SUBPROGRAMS LISTED IN CALL TO
C           FORIF
C
C        METHOD
C           USES RECURSIVE TECHNIQUE DESCRIBED IN A. RALSTON, H. WILF,
C           'MATHEMATICAL METHODS FOR DIGITAL COMPUTERS', JOHN WILEY
C           AND SONS, NEW YORK, 1960, CHAPTER 24. THE METHOD OF
C           INDEXING THROUGH THE PROCEDURE HAS BEEN MODIFIED TO
C           SIMPLIFY THE COMPUTATION.
C
C     ..................................................................
C
      SUBROUTINE FORIF(FUN,N,M,A,B,IER)
      DIMENSION A(1),B(1)
C
C        CHECK FOR PARAMETER ERRORS
C
      IER=0
   20 IF(M) 30,40,40
   30 IER=2
      RETURN
   40 IF(M-N) 60,60,50
   50 IER=1
      RETURN
C
C        COMPUTE AND PRESET CONSTANTS
C
   60 AN=N
      COEF=2.0/(2.0*AN+1.0)
      CONST=3.141593*COEF
      S1=SIN(CONST)
      C1=COS(CONST)
      C=1.0
      S=0.0
      J=1
      FUNZ=FUN(0.0)
   70 U2=0.0
      U1=0.0
      AI=2*N
C
C        FORM FOURIER COEFFICIENTS RECURSIVELY
C
   75 X=AI*CONST
      U0=FUN(X)+2.0*C*U1-U2
      U2=U1
      U1=U0
      AI=AI-1.0
      IF(AI) 80,80,75
   80 A(J)=COEF*(FUNZ+C*U1-U2)
      B(J)=COEF*S*U1
      IF(J-(M+1)) 90,100,100
   90 Q=C1*C-S1*S
      S=C1*S+S1*C
      C=Q
      J=J+1
      GO TO 70
  100 A(1)=A(1)*0.5
      RETURN
      END