* * $Id: tricof.F,v 1.1.1.1 1996/04/01 15:03:14 mclareni Exp $ * * $Log: tricof.F,v $ * Revision 1.1.1.1 1996/04/01 15:03:14 mclareni * Mathlib gen * * #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