*
* $Id$
*
* $Log$
* Revision 1.1  1996/04/01 15:03:13  mclareni
* Initial revision
*
*
#include "sys/CERNLIB_machine.h"
#include "_gen/pilot.h"
      SUBROUTINE TRICOF(F,NF,A,NE,B,NO,IOP)
C
C     CERN LIBRARY NO E 206
C
C     REVISED VERSION JUNE 1971.
C
C     PURPOSE = TO COMPUTE THE COEFFICIENTS IN A TRIGONOMETRIC EXPANSION
C               FOR A FUNCTION GIVEN IN EQUIDISTANT POINTS
C
C     PARAMETERS
C
C     F       = AN ARRAY USED FOR STORING THE FUNCTION VALUES F(X)
C     NF      = THE NUMBER OF FUNCTION VALUES IN THE ARRAY F.NF MUST
C               HAVE THE STRUCTURE , NF = 2*N+1 , THE GENERAL CASE
C                                    NF = N+1   , THE EVEN CASE
C                                    NF = N-1   , THE ODD CASE
C     A       = AN ARRAY USED FOR RETURNING THE COEFFICIENTS OF THE COS-
C               INE TERMS
C     NE      = THE NUMBER OF COEFFICIENTS IN THE ARRAY A , NE.LE.(N+1)
C     B       = AN ARRAY USED FOR RETURNING THE COEFFICIENTS OF THE SINE
C               TERMS
C     NO      = THE NUMBER OF COEFFICIENTS IN THE ARRAY B , NO.LE.(N-1)
C     IOP     = OPTION NUMBER , IOP = 1 , THE GENERAL CASE
C                               IOP = 2 , THE EVEN CASE
C                               IOP = 3 , THE ODD CASE
C
      DIMENSION F(NF) , A(NE) , B(NO)
C
      REAL KSI0 , KSI1 , KSIK
      DATA ZERO,HALF,ONE,TWO,PI/0.,.5,1.,2.,3.1415926535898/
C
C     COMPUTE THE NUMBER N (SEE EXPLANATION OF PARAMETERS)
C
C*UL 1000 N=0
      N=0
      IF (IOP.EQ.1) N=(NF-1)/2
      IF (IOP.EQ.2) N=NF-1
      IF (IOP.EQ.3) N=NF+1
      MESSAG=1
      IF (N.EQ.0) GO TO 1280
C
C     STOP IF IOP DOES NOT HAVE A CORRECT VALUE
C
      IF (IOP.GT.1) GO TO 1030
      MESSAG=2
      IF ((2*N-NF+1).NE.0) GO TO 1280
C
C     STOP IF NF DOES NOT HAVE THE CORRECT STRUCTURE IN THE GENERAL CASE
C
C     SPLIT THE FUNCTION F(X) IN AN EVEN AND ODD PART
C
      M=N+1
      DO 1020 J=1,N
      COF1=HALF*(F(M+J)+F(M-J))
      COF2=HALF*(F(M+J)-F(M-J))
      F(M+J)=COF2
      F(M-J)=COF1
 1020 CONTINUE
C
C     REWRITE N IN POWERS OF 2 I.E. N=NBASE*2**NEXP
C
 1030 NBASE=N
      NEXP =0
 1040 NINR =NBASE/2
      IF ((NBASE-2*NINR).NE.0) GO TO 1050
      NBASE=NINR
      NEXP =NEXP+1
      GO TO 1040
C
C     DO SOME INITIAL CALCULATIONS
C
 1050 REALN=NBASE
      ARG  =HALF*PI/REALN
      KSI0 =COS(ARG)
      ETA0 =SIN(ARG)
C
C     START CALCULATION OF COEFFICIENTS
C
      IF (IOP.EQ.3) GO TO 1160
C
C     **********     EVEN COEFFICIENT CALCULATION     **********
C
C     CHECK IF NE IS IN THE CORRECT RANGE
C
      MESSAG=3
      IF ((NE.LE.2).OR.(NE.GT.(N+1))) GO TO 1280
C
C     COMPUTE THE BASIC COEFFICIENTS A(K) , K=1(1)(NBASE+1)
C
C     START CALCULATION OF A(1)
C
      NN    =NBASE-1
      NPOINT=1
      NINCRE=2**NEXP
      NLOCAL=NINCRE+1
      BASEIN=ONE/REALN
      A(1)  =HALF*(F(1)+F(N+1))
      IF (NN.EQ.0) GO TO 1065
      DO 1060 J=1,NN
      A(1)  =A(1)+F(NLOCAL)
      NLOCAL=NLOCAL+NINCRE
 1060 CONTINUE
 1065 A(1)  =TWO*BASEIN*A(1)
C
C     START CALCULATION OF A(K) , K=2(1)(NBASE+1)
C
      KSI1=KSI0
      KSIK=KSI1
      ETA1=ETA0
      ETAK=ETA1
      CONST=HALF*F(N+1)
C
      NUPPER=NBASE
      IF (NBASE.GT.(NE-1)) NUPPER=NE-1
C
      DO 1090 K=1,NUPPER
      COF1=TWO*(TWO*KSIK**2-ONE)
      A2  =ZERO
      A1  =A2
      A0  =CONST
      NLOCAL=N+1-NINCRE
      DO 1070 J=1,NBASE
      A2=A1
      A1=A0
      A0=F(NLOCAL)+COF1*A1-A2
      NLOCAL=NLOCAL-NINCRE
 1070 CONTINUE
C
C*UL 1080 A(K+1)=BASEIN*(A0-A2)
      A(K+1)=BASEIN*(A0-A2)
      COF1  =KSIK
      COF2  =ETAK
      KSIK =KSI1*COF1-ETA1*COF2
      ETAK =ETA1*COF1+KSI1*COF2
 1090 CONTINUE
C
C     CALCULATION OF THE BASIC EVEN COEFFICIENTS FINISHED
C
      IF (NEXP.EQ.0) GO TO 1145
C
C     CONTINUE CALCULATION OF EVEN COEFFICIENTS
C
      NUMCOF=NBASE
      DO 1140 NSTEP=1,NEXP
      NINCRE=2**(NEXP-NSTEP)
      NPOINT=NINCRE+1
      NINCRE=2*NINCRE
      NLOCAL=NPOINT
      NUMBER=2*NUMCOF+1
C
C     COMPUTE CONSTANT TERM IN MID-POINT APPROXIMATION I.E. K=1
C
      SUM=ZERO
      DO 1100 J=1,NUMCOF
      SUM=SUM+F(NLOCAL)
      NLOCAL=NLOCAL+NINCRE
 1100 CONTINUE
C
      SUM =TWO*BASEIN*SUM
      COF1=A(1)
      A(1)=HALF*(COF1+SUM)
      IF(NUMBER.LE.NE) A(NUMBER)=HALF*(COF1-SUM)
C
      IF (NUMCOF.EQ.1) GO TO 1134
C
C     COMPUTE MID-POINT APPROXIMATION FOR K=2(1)NUMCOF
C
C*UL 1105 NN  =NUMCOF-1
      NN  =NUMCOF-1
      NUPPER=NN
      IF (NUPPER.GT.(NE-1)) NUPPER=NE-1
      KSIK=KSI1
      ETAK=ETA1
      DO 1130 K=1,NUPPER
      COF1=TWO*(TWO*KSIK**2-ONE)
      A2=ZERO
      A1=A2
      NLOCAL=N+2-NPOINT
      A0=F(NLOCAL)
      DO 1110 J=1,NN
      A2=A1
      A1=A0
      NLOCAL=NLOCAL-NINCRE
      A0=F(NLOCAL)+COF1*A1-A2
 1110 CONTINUE
C
C*UL 1120  SUM=TWO*BASEIN*(A0-A1)*KSIK
      SUM=TWO*BASEIN*(A0-A1)*KSIK
      COF1=A(K+1)
      A(K+1)=HALF*(COF1+SUM)
      INDEX=NUMBER-K
      IF (INDEX.GT.NE) GO TO 1025
      A(INDEX)=HALF*(COF1-SUM)
C
 1025 COF1=KSIK
      COF2=ETAK
      KSIK=KSI1*COF1-ETA1*COF2
      ETAK=ETA1*COF1+KSI1*COF2
C
 1130 CONTINUE
 1134 INDEX=NUMCOF+1
      IF (INDEX.GT.NE) GO TO 1136
      A(INDEX)=HALF*A(INDEX)
C
C     CALCULATIONS OF MID-POINT APPROXIMATIONS FINISHED
C
C     DO CHANGES RELATED TO HALVING OF THE INTERVAL
C
 1136 ARG =HALF*ARG
      COF1=ETA1
      ETA1=SIN(ARG)
      KSI1=HALF*COF1/ETA1
      BASEIN=HALF*BASEIN
      NUMCOF=2*NUMCOF
C
 1140 CONTINUE
 1145 IF (NEXP.EQ.0) NUMBER=NBASE+1
      IF (NUMBER.GT.NE) GO TO 1150
      A(NUMBER)=HALF*A(NUMBER)
C
C     CALULATION OF EVEN COEFFICIENTS FINISHED
C
 1150 IF (IOP.EQ.2) RETURN
C
C     RETURN TO CALLING PROGRAM IF F(X) WAS AN EVEN FUNCTION
C     IF IOP=1 CHANGE SIGN OF EACH SECOND COEFFICIENTS
C
      NINR=NE/2
      IF (NINR.EQ.0) GO TO 1166
      DO 1164 K=1,NINR
      A(2*K)=-A(2*K)
 1164 CONTINUE
C
C     **********     ODD COEFFICIENT CALCULATION     **********
C
C     CHECK IF NO IS IN THE CORRECT RANGE
C
      MESSAG=4
      IF ((NO.LE.1).OR.(NO.GT.(N-1))) GO TO 1280
C
C     COMPUTE THE BASIC COEFFICIENTS B(K) , K=1(1)NBASE
C
 1166 ARG=HALF*PI/REALN
 1160 IF (IOP.EQ.1) NMAX=2*N+1
      IF (IOP.EQ.3) NMAX=N
      NINCRE=2**NEXP
      NPOINT=NMAX-NINCRE
      NLOCAL=NPOINT
      BASEIN=ONE/REALN
      B(1)=ZERO
      IF (NBASE.EQ.1) GO TO 1200
      KSI1=TWO*KSI0**2-ONE
      KSIK=KSI1
      ETA1=TWO*KSI0*ETA0
      ETAK=ETA1
      NN =NBASE-1
      NNN=NN-1
      IF (NN.GT.NO) NN=NO
      DO 1190 K=1,NN
      COF1=TWO*KSIK
      A2  =ZERO
      A1  =A2
      A0  =F(NPOINT)
      NLOCAL=NPOINT-NINCRE
      IF (NNN.EQ.0) GO TO 1180
      DO 1170 J=1,NNN
      A2=A1
      A1=A0
      A0=F(NLOCAL)+COF1*A1-A2
      NLOCAL=NLOCAL-NINCRE
 1170 CONTINUE
C
 1180 B(K)=TWO*BASEIN*A0*ETAK
      COF1=KSIK
      COF2=ETAK
      KSIK=KSI1*COF1-ETA1*COF2
      ETAK=ETA1*COF1+KSI1*COF2
 1190 CONTINUE
C
C     CALCULATION OF THE BASIC ODD COEFFICIENTS FINISHED
C
 1200 IF (NEXP.EQ.0) GO TO 1260
C
C     CONTINUE CALCULATION OF ODD COEFFICIENTS
C
      KSI1=KSI0
      ETA1=ETA0
C
      NUMCOF=NBASE
      DO 1250 NSTEP=1,NEXP
      KSIK=KSI1
      ETAK=ETA1
      NINCRE=2**(NEXP-NSTEP)
      NPOINT=NMAX-NINCRE
      NINCRE=2*NINCRE
      NUMBER=2*NUMCOF
      IF (NUMCOF .LE. NO)  B(NUMCOF) = ZERO
C
C     COMPUTE MID-POINT APPROXIMATIONS FOR K=1(1)NUMCOF
C
      NN  =NUMCOF-1
      NUPPER=NUMCOF
      IF (NUPPER.GT.NO) NUPPER=NO
      DO 1240 K=1,NUPPER
      COF1=TWO*(TWO*KSIK**2-ONE)
      A2  =ZERO
      A1  =A2
      NLOCAL=NPOINT
      A0  =F(NLOCAL)
      IF (NN.EQ.0) GO TO 1220
      DO 1210 J=1,NN
      A2=A1
      A1=A0
      NLOCAL=NLOCAL-NINCRE
      A0=F(NLOCAL)+COF1*A1-A2
 1210 CONTINUE
C
 1220 SUM=TWO*BASEIN*(A0+A1)*ETAK
      COF1=B(K)
      B(K)=HALF*(COF1+SUM)
      IF (K.EQ.NUMCOF) GO TO 1230
      INDEX=NUMBER-K
      IF (INDEX.GT.NO) GO TO 1230
      B(INDEX)=-HALF*(COF1-SUM)
C
 1230 COF1=KSIK
      COF2=ETAK
      KSIK=KSI1*COF1-ETA1*COF2
      ETAK=ETA1*COF1+KSI1*COF2
C
 1240 CONTINUE
C
C     CALCULATION OF MID-POINT APPROXIMATION FINISHED
C
C     DO CHANGES RELATED TO HALVING OF INTERVAL
C
      ARG =HALF*ARG
      COF1=ETA1
      ETA1=SIN(ARG)
      KSI1=HALF*COF1/ETA1
      BASEIN=HALF*BASEIN
      NUMCOF=2*NUMCOF
C
 1250 CONTINUE
C
C     CALCULATION OF ODD COEFFICIENTS FINISHED
C
 1260 IF (IOP.EQ.3) RETURN
C
C     IF IOP=1 RECOMPUTE FUNCTION VALUES
C
      DO 1270 J=1,N
      COF2=F(M+J)
      COF1=F(M-J)
      F(M+J)=COF1+COF2
      F(M-J)=COF1-COF2
 1270 CONTINUE
C
      RETURN
C
C     PRINT ERROR MESSAGE.
C
 1280 WRITE(6,1290)  MESSAG
      STOP
C
 1290 FORMAT(5X,'SUBROUTINE TRICOF ERROR NO.',I2)
C
      END
