C ALGORITHM 649, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 13, NO. 1, P. 97. C C*********************************************************************** C* * C* ******** PACKAGE FOURCO ******** * C* * C* * C* ..... PACKAGE FOR COMPUTING SINE AND COSINE FOURIER ..... * C* ..... COEFFICIENTS OF A FUNCTION OVER AN INTERVAL ..... * C* * C* DOUBLE PRECISION VERSION - SEPTEMBER 1986 * C* * C* * C***** AUTHORS : GIULIO GIUNTA AND ALMERICO MURLI * C* DIPARTIMENTO DI MATEMATICA E APPLICAZIONI * C* UNIVERSITA' DI NAPOLI * C* NAPLES - ITALY * C* * C*********************************************************************** C C***** GENERAL PURPOSE - THIS PACKAGE EVALUATES THE FIRST ..NCOEF.. C TRIGONOMETRICAL FOURIER COEFFICIENTS C C C(J)=2/(B-A)*INTEGRAL BETWEEN A,B(FUN(X)* C COS(2*PI*(J-1)*(X-A)/(B-A))) C C S(J)=2/(B-A)*INTEGRAL BETWEEN A,B(FUN(X)* C SIN(2*PI*(J-1)*(X-A)/(B-A))) C C J=1,2,....,NCOEF C C (HOPEFULLY) TO A UNIFORM ABSOLUTE ACCURACY ..TOL.. C C THE PACKAGE HAS TWO TOP-ENDS (FCOEFS,FCOETD) THAT C THE USER CAN SELECT DEPENDING ON HIS INPUT PROBLEM. C C --FCOEFS. TO USE THIS TOP-END THE USER MUST PROVIDE A C FUNCTION SUBPROGRAM FOR EVALUATING ..FUN.. AT ANY C POINT OF THE INTEGRATION INTERVAL (A,B). C C --FCOETD. TO USE THIS TOP-END THE USER MUST PROVIDE C SAMPLED VALUES OF FUN AT EQUALLY SPACED POINTS OF C THE INTEGRATION INTERVAL (A,B). C C ALL THAT THE USER MUST KNOW IN ORDER TO USE THE C PACKAGE IS THE CALLING SEQUENCE OF THE TOP-END C SUITABLE FOR HIS PROBLEM C C C***** STRUCTURE AND PORTABILITY - THE PACKAGE COMPRISES A BASIC C ROUTINE FCOEBA,TWO TOP-ENDS,A SET OF FFT ROUTINES FOR C REAL DATA,A SET OF ROUTINES FOR NUMERICAL DIFFEREN- C TIATION AND TWO AUXILIARY ROUTINES. C ALL IS WRITTEN IN FORTRAN 77. ALL MACHINE DEPENDENT C CONSTANTS ARE SPECIFIED BY D1MACH FORTRAN FUNCTION C USED IN PORT LIBRARY. ALL THE CALLS TO INTRINSIC C FORTRAN FUNCTIONS USE GENERIC NAMES. MATHEMATICAL C CONSTANTS ARE COMPUTED TO FULL MACHINE ACCURACY IN C THE CODE. C C***** FORTRAN SUPPLIED FUNCTIONS - DBLE,SIN,COS,ABS,MIN,MAX,SQRT,ATAN C MOD. C***** OTHER REQUIRED FUNCTIONS - D1MACH C C C********************************************************************** C********************************************************************** C STANDARD TOP END FCOEFS * C********************************************************************** SUBROUTINE FCOEFS(A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,C,S,IERR,ESTER,MACT) C C********************************************************************** C C TOP-END FCOEFS ALLOWS TO COMPUTE TRIGONOMETRICAL FOURIER COEFFI- C CIENTS WHEN THE INTEGRAND FUNCTION ..FUN.. IS PROVIDED AS A C FUNCTION SUBPROGRAM C C C INPUT PARAMETERS - A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK,XLAMB. C C A,B -DOUBLE PRECISION- LOWER AND UPPER LIMIT OF C INTEGRATION. C C FUN -DOUBLE PRECISION FUNCTION OF ONE DOUBLE C PRECISION ARGUMENT - THE INTEGRAND C FUNCTION. C (EXTERNAL USER SUPPLIED). FUN MUST BE A C SMOOTH REAL VALUED FUNCTION ON THE CLOSED C INTERVAL (A,B). THE ACTUAL NAME FOR C FUN NEEDS TO BE DECLARED EXTERNAL IN THE C CALLING PROGRAM. C C NCOEF -INTEGER- THE REQUIRED NUMBER OF FOURIER COEF C C TOL -DOUBLE PRECISION - THE REQUIRED ABSOLUTE C ACCURACY OF THE COEFFICIENTS. NOTE THAT C IF A VALUE OF TOL .LT. THE MACHINE EPSILON C IS PROVIDED,THE ROUTINE TRIES TO DO THE C BEST IT CAN (SEE NOTE 1). C HOWEVER ..ESTER.. GREATER THAN THE C MACHINE EPSILON AND ..IERR.. NE 0 SHOULD C BE EXPECTED IN OUTPUT. C C MTOP -INTEGER- THE MAXIMUM NUMBER OF FUNCTION C EVALUATIONS IN FCOEBA IS MTOP+3. MTOP SHOULD C BE .GT. 16 AND AN INTEGER POWER OF 2. C IF MTOP IS ERRONEOUSLY SET .LT. 16 IT IS C SET TO 16 BY DEFAULT (TAKE CARE TO THE SIZE C THE WORKSPACE ARRAY). IF IT IS NOT SET TO C A POWER OF 2 THE ACTUAL VALUE USED IS THE C LARGEST POWER OF 2 LESS THAN MTOP. C IF MTOP IS SET NEGATIVE THEN ONE C OF THE FEATURES OF THE ABNORMAL TERMINATION C CRITERION IS SWITCHED OFF(SEE NOTE 1). C C IDER -INTEGER- ALLOWS THE CHOICE OF DIFFERENTIATIO C ROUTINE. C =0 USE THE STANDARD DIFFERENT.ROUTINE(CHPEN C 18 VALUES OF ..FUN.. IN THE INTERVAL C (A-(B-A)*7/64.,B+(B-A)*7/64.) ARE USED. C =1 USE A MORE ACCURATE AND EXPANSIVE DIFF. C ROUTINE(ENDACE).42 VALUES OF ..FUN.. IN C INTERVAL (A-(B-A)*19/64.,B+(B-A)*19/64.) C ARE USED. CHOOSE THIS OPTION IF THE C ACCURACY REQUIREMENTS ARE STRINGENT. C =2 USE A ONE-SIDED DIFFERENT.ROUTINE. C 14 VALUES OF ..FUN.. IN THE INTERVAL C (A,B) ARE USED. C =-P THE USER PROVIDES THE VALUE OF THE C FIRST P DERIVATIVES OF ..FUN.. AT ..A.. C AND ..B.. (SEE ..XLAMB.. BELOW). C P MUST BE LESS OR EQUAL TO 14,OTHERWISE C ONLY THE FIRST 14 DERIVATIVES ARE USED. C C IFIRST -INTEGER- NORMALLY MUST BE SET TO ZERO. C IF SET NONZERO IT CAUSES THE RESTART OF THE C PROCEDURE USING RESULTS FROM THE PREVIOUS C CALL TO THE TOP-END. IN THIS CASE ..IDER, C WORK,XLAMB,MACT,IERR.. FROM THE PREVIOUS C RUN MUST BE UNCHANGED (SEE ..IERR.. BELOW). C C WORK -DOUBLE PRECISION ARRAY(4*MTOP+20)- WORK SPAC C ARRAY. C C XLAMB -DOUBLE PRECISION(15)- C NORMALLY NOT USED. IF ..IDER.. IS SET C NEGATIVE THEN XLAMB MUST CONTAIN C (I) (I) C XLAMB(I)=FUN (B) - FUN (A) C I=1,2,...,ABS(IDER) C C I.E. THE DIFFERENCES OF THE FIRST ABS(IDER) C DERIVATIVES OF ..FUN.. AT ..B.. AND ..A.. C C C OUTPUT PARAMETERS - C,S,IERR,ESTER,MACT C C C -DOUBLE PRECISION(NCOEF)- ARRAY CONTAINING THE C FIRST NCOEF COSINE FOURIER COEFFICIENTS C C S -DOUBLE PRECISION(NCOEF)- ARRAY CONTAINING THE C FIRST NCOEF SINE FOURIER COEFFICIENTS C C /// WARNING : C(J),S(J) ARE THE (J-1)-TH /// C /// TRIGONOMETRIC FOURIER COEFFICIENTS /// C C IERR -INTEGER- ERROR INDICATOR C C IERR=-1 INPUT ERROR. ..NCOEF.. LE 0 OR REST C INVOKED WITH IERR FROM THE PREVIOUS C RUN NOT POSITIVE. C IERR=0 NORMAL TERMINATION.TOLERANCE SATISF C IERR.GT.0 PHYSICAL LIMIT TERMINATION. C (POSSIBILITY OF RESTART) C IERR=-2 ABNORMAL TERMINATION C C /// SEE NOTE 1. /// C C ESTER -DOUBLE PRECISION- THE ESTIMATED ABSOLUTE ACCU C NEGATIVE WHEN ..IERR.. NE 0. C C MACT -INTEGER- THE TOTAL NUMBER OF FUNCTION C EVALUATIONS ACTUALLY USED. C C*********************************************************************** C C NOTE 1. C C AFTER STAGE S (AFTER EXPENDING 2**S +3 FUNCTION VALUES) THE C ROUTINE FORMS AN ACCURACY ESTIMATE E(S) C C (1) IF E(S) .LT. TOL THE CALCULATION TERMINATES WITH IERR=0 C (THIS IS NORMAL TERMINATION) C C (2) IF (1) IS NOT SATISFIED, BUT 2**S = MTOP C THE CALCULATION TERMINATES WITH IERR=1 C (THIS IS PHYSICAL LIMIT TERMINATION) C C (3) IF (1) AND (2) ARE NOT SATISFIED BUT BOTH C C (A) E(S) .LT. ABNC1 * MAX(E(J):J=1,2,..,S-1) C C AND C C (B) E(S) .GT. ABNC2 * E(S-1) C C THE CALCULATION TERMINATES WITH IERR=-2 C (THIS IS AN ABNORMAL TERMINATION) C ABNC1=0.01 AND ABNC2=SQRT(0.1) HAVE BEEN CHOSEN AS A RESULT C OF EXPERIMENTATION. C C THE USER CAN ELIMINATE THE ABNORMAL FEATURE (3) BY USING C -ABS(MTOP) IN PLACE OF MTOP AS AN INPUT ARGUMENT.(THIS C REPLACES ABNC1 BY 0 INTERNALLY) C THE USER MAY ENSURE A RESULT BASED ON MTOP FUNCTION VALUES C BY BOTH USING -ABS(MTOP) IN PLACE OF MTOP AND SETTING C TOL=0.0D0. C C COMMENT. UNDER NORMAL CONDITIONS IT IS USUALLY THE CASE THAT C IF IERR.GE.1 ONE MIGHT OBTAIN BETTER ACCURACY BY RAISING THE C VALUE OF MTOP AND RE-RUNNING WITH IFIRST=1 C IF IERR=-2 ONLY MARGINALLY BETTER ACCURACY CAN BE OBTAINED AT C SIGNIFICANT EXPENSE SINCE THE PRESENT ACCURACY IS PROBABLY C GOVERNED BY NOISE LEVEL IN THE FUNCTION VALUES.NOTE THAT C IN THE LATTER CASE RESTART IS NOT ALLOWED. C********************************************************************** C FIRST DECLARATION STATEMENT OF FCOEFS C********************************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER NCOEF,IERR,MTOP,MACT,IDER,IFIRST INTEGER IPM1,IPM2 DIMENSION C(NCOEF),S(NCOEF),WORK(*) DIMENSION XLAMB(*),EREST1(14),EREST2(14) EXTERNAL FUN DATA EPS/0.D0/ C C********************************************************************** C FIRST EXECUTABLE STATEMENT OF FCOEFS C********************************************************************** C IF (EPS.EQ.0.D0) EPS=D1MACH(4) IF(IFIRST.EQ.0)GOTO 20 IF(IERR.GT.0) GO TO 10 IERR = -1 RETURN C 10 CONTINUE IPM1=IERR IF(IDER.LT.0) INITM=MACT-3 IF(IDER.EQ.0) INITM=MACT-21 IF(IDER.EQ.1) INITM=MACT-45 IF(IDER.EQ.2) INITM=MACT-17 DO 15 I=1,INITM+2 WORK(INITM+ABS(MTOP)+5-I) = WORK(2*INITM+5-I) 15 CONTINUE GO TO 130 C 20 CONTINUE INITM=8 IF (NCOEF.GT.0) GO TO 30 IERR = -1 RETURN C 30 CONTINUE IF (IDER.GE.0) GOTO 40 IPM1 = MIN(1-IDER,15) IPM2 = IPM1-1 GO TO 110 C 40 CONTINUE C * * C * HBASE IS THE STEPSIZE FOR THE DIFFERENTIATION ROUTINES * C * * HBASE = (B-A)/64.D0 C C CHOOSE THE DIFFERENTIATION ROUTINE C IF (IDER.NE.1) GO TO 70 CALL ENDACE(B,14,HBASE,WORK,EREST1,FUN) CALL ENDACE(A,14,HBASE,XLAMB,EREST2,FUN) DO 50 K=1,14 IF (EREST1(K).LT.0.D0.OR.EREST2(K).LT.0.D0) GO TO 60 50 CONTINUE C IPM1 = 15 GO TO 90 C 60 CONTINUE IPM1 = K GO TO 90 C 70 CONTINUE IF (IDER.NE.2) GO TO 80 CALL ONESID(-1,B,HBASE,WORK,FUN) CALL ONESID(1,A,HBASE,XLAMB,FUN) IPM1 = 6 GO TO 90 C 80 CONTINUE IPM1 = 6 CALL CHPEND(B,HBASE,WORK,FUN) CALL CHPEND(A,HBASE,XLAMB,FUN) 90 CONTINUE C IPM2 = IPM1 - 1 DO 100 K = 1,IPM2 XLAMB(K) = WORK(K) - XLAMB(K) 100 CONTINUE C C CHECK FOR PERIODIC DERIVATIVES C 110 CONTINUE DO 120 K = 1,IPM2 I = IPM2 - K + 1 IF(ABS(XLAMB(I)).GT.EPS * 1.0D4)GO TO 130 IPM1 = IPM1 - 1 120 CONTINUE C 130 CONTINUE C CALL FCOEBA(A,B,FUN,NCOEF,TOL,MTOP,IPM1,INITM,IFIRST,WORK(1), + WORK(ABS(MTOP)+3),WORK(2*ABS(MTOP)+5),XLAMB,C,S, + IERR,ESTER,MACT) C IF(IFIRST.NE.0) MACT = MACT-1 IF(IDER.EQ.0) MACT = MACT+18 IF(IDER.EQ.1) MACT = MACT+42 IF(IDER.EQ.2) MACT = MACT+14 RETURN C C C*********************************************************************** C C ******* END OF SUBROUTINE FCOEFS ******* C C*********************************************************************** END SUBROUTINE FCOETD(A,B,NCOEF,MTOP,FUNVAL,WORK,C,S,IERR + ,ESTER) C**************************************************************** C********************************************************************** C TOP END FCOETD * C********************************************************************** C C TOP-END FCOETD ALLOWS TO COMPUTE TRIGONOMETRICAL FOURIER C COEFFICIENTS WHEN THE INTEGRAND FUNCTION FUN IS PROVIDED C BY MEANS OF A SET OF FUNCTION VALUES AT EQUALLY SPACED POINTS C OF THE INTEGRATION INTERVAL (A,B). C C C INPUT PARAMETERS - A,B,NCOEF,MTOP,FUNVAL,WORK. C C A,B -DOUBLE PRECISION- THE BOUNDS OF INTEGRATIO C INTERVAL C C NCOEF -INTEGER- THE REQUIRED NUMBER OF FOURIER COE C C MTOP -INTEGER- THE NUMBER OF FUNCTION VALUES C STORED IN THE ARRAY FUNVAL. IT MUST BE C MTOP = 2**K + 1 , K .GE.4 AND LE 16 C IF MTOP IS ERRONEOUSLY SET IERR=-1 IS C RETURNED. C C FUNVAL -DOUBLE PRECISION(MTOP)- ARRAY CONTAINING TH C FUNCTION VALUES. C THE VALUE OF FUN(A+(I-1)*(B-A)/(MTOP-1)) IS C STORED IN FUNVAL(I) , I=1,2,..,MTOP. C C WORK -DOUBLE PRECISION(4*MTOP+20)- WORK SPACE ARR C C C OUTPUT PARAMETERS - C,S,IERR,ESTER C C C -DOUBLE PRECISION(NCOEF)- ARRAY CONTAINING T C FIRST NCOEF COSINE FOURIER COEFFICIENTS C C S -DOUBLE PRECISION(NCOEF)- ARRAY CONTAINING T C FIRST NCOEF SINE FOURIER COEFFICIENTS C C /// WARNING : C(J),S(J) ARE THE (J-1)-TH /// C /// TRIGONOMETRIC FOURIER COEFFICIENTS /// C C IERR -INTEGER- ERROR INDICATOR C C IERR=-1 INPUT ERROR.NCOEF.LT.0 OR MTOP C NOT EQUAL TO 2**K +1 , K.GE.4 C IERR.GE.0 NORMAL TERMINATION. C C ESTER - DOUBLE PRECISION- THE ESTIMATED ABSOLUTE A C C C*********************************************************************** C C*********************************************************************** C FIRST DECLARATION STATEMENT OF FCOETD C*********************************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER NCOEF,MTOP,IERR INTEGER IPOW2,IXVAL,IPM1,IPM2 DIMENSION FUNVAL(MTOP),WORK(*),C(NCOEF),S(NCOEF) DIMENSION XLAMB(6) C EXTERNAL DUMFUN C DATA EPS/0.D0/ C*********************************************************************** C FIRST EXECUTABLE STATEMENT OF FCOETD C*********************************************************************** IF (EPS.EQ.0.D0) EPS=D1MACH(4) C C **** TEST ON INPUT ARGUMENTS **** C IF (NCOEF.GT.0.AND.MTOP.GT.16) GO TO 10 IERR = -1 RETURN C 10 CONTINUE IPOW2 = 16 DO 20 K=1,12 IF (MTOP.EQ.IPOW2+1) GO TO 30 IPOW2 = IPOW2*2 20 CONTINUE C IERR = -1 RETURN C 30 CONTINUE C C **** COMPUTE APPROXIMATED DERIVATIVES **** C NMTOP = MTOP -1 HBASE = (B-A) / NMTOP IXVAL = MTOP CALL ONESMP(-1,FUNVAL,IXVAL,HBASE,WORK) IXVAL = 1 CALL ONESMP(1,FUNVAL,IXVAL,HBASE,XLAMB) IPM1 = 6 IPM2 = 5 DO 40 K=1,IPM2 XLAMB(K) = WORK(K)-XLAMB(K) 40 CONTINUE C C C CHECK FOR PERIODIC DERIVATIVES C DO 50 K=1,IPM2 I = IPM2-K+1 IF (ABS(XLAMB(I)).GT.EPS*1.0D4) GO TO 60 IPM1 = IPM1-1 50 CONTINUE C 60 CONTINUE IFIRST = 0 INITM = -NMTOP DO 70 K=1,MTOP WORK(K) = FUNVAL(K) 70 CONTINUE C TOL=0.D0 CALL FCOEBA(A,B,DUMFUN,NCOEF,TOL,NMTOP,IPM1,INITM,IFIRST + ,WORK(1),WORK(ABS(MTOP)+3),WORK(2*ABS(MTOP)+5) + ,XLAMB,C,S,IERR,ESTER,MACT) C ESTER=-ESTER RETURN C*********************************************************************** C C ******* END OF SUBROUTINE FCOETD ******* C C*********************************************************************** END DOUBLE PRECISION FUNCTION DUMFUN(X) C THIS IS A DUMMY FUNCTION, JUST TO FILL THE CALLING SEQUENCE C OF FCOEBA IN FCOETD. IN FACT FCOETD USES SAMPLED FUNCTION C VALUES STORED IN FUNVAL. DOUBLE PRECISION X DUMFUN = X RETURN END SUBROUTINE FCOEBA(A,B,FUN,NCOEF,TOL,MTOP,IPM1,INITM,IFIRST, + WORK1,WORK2,WORK3,XLAMB,C,S,IERR,ESTER,MACT) C*********************************************************************** C C C ******** BEGINNING OF SUBROUTINE **** FCOEBA ************ C C C*********************************************************************** C C C*********************************************************************** C C INPUT PARAMETERS - A,B,FUN,NCOEF,TOL,MTOP,IPM1,INITM,IFIRST,WORK1, C WORK2,WORK3,XLAMB. C C A,B,FUN,NCOEF,TOL,MTOP,IFIRST AS IN FCOEFS C C IPM1 -INTEGER- THE NUMBER OF INPUT DERIVATIVES +1 C C INITM -INTEGER- THE NUMBER OF FUNCTION EVALUATIONS C USED FOR FIRST APPROXIMATION. C IT IS SET NEGATIVE IF FUN IS SAMPLED (TOP-END C FCOETD) C C WORK1,WORK2,WORK3 -DOUBLE PRECISION ARRAYS- WORKSPACE C C XLAMB -DOUBLE PRECISION(IPM1)- THE ARRAY OF LENGTH C IPM1 CONTAINING THE FOLLOWING C QUANTITIES C C (I) (I) C XLAMB(I)=FUN (B) - FUN (A) C C I=1, 2, ..., IPM1 -1 C C OUTPUT PARAMETERS - C,S,IERR,ESTER,MACT AS IN FCOEFS C C C*********************************************************************** C FIRST DECLARATION STATEMENT OF FCOEBA C*********************************************************************** C IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL RNDLEV INTEGER NCOEF,IPM1,IERR,MTOP,MACT,INITM,IFIRST INTEGER MAXM,IPM2,NHALF,N,NQUART DIMENSION XLAMB(*),C(NCOEF),S(NCOEF),WORK1(*),WORK2(*),WORK3(*) DATA EPS/0.D0/,TWOPI/0.D0/ C * * C******************************************************************* C * FIRST EXECUTABLE STATEMENT OF FCOEBA * C******************************************************************* C **** INITIALIZATION **** C -------------- C IF (EPS.EQ.0.D0) EPS=D1MACH(4) IF (TWOPI.EQ.0.D0) TWOPI=8.D0*ATAN(1.D0) C XL = B - A MACT = 0 C * * C * THE ABNORMAL TERMINATION CRITERION IS BASED ON THE VALUES * C * OF COSTANTS ABNC1,ABNC2 AND IT IS SWITCHED OFF WHEN MTOP * C * IS NEGATIVE * C * * ABNC1 = 0.01D0 ABNC2 = SQRT(0.1D0) IF (MTOP.LT.0) ABNC1 = 0.D0 C RNDLEV = .FALSE. MAXM = MAX(16,ABS(MTOP)) IERR = 0 ESTMAX = 0.D0 N = ABS(INITM) NHALF = N/2 NQUART = N/4 IPM2 = IPM1 -1 IF (INITM.LT.0) GO TO 20 IF (IFIRST.NE.0) GO TO 130 DO 10 I = 1,MAXM+2 WORK1(I) = 0.0 WORK2(I) = 0.0 10 CONTINUE C FA = FUN(A) FB = FUN(B) C C * * C * MODIFICATION OF ARRAY XLAMB * C * * 20 CONTINUE IF (IPM2.EQ.0) GO TO 40 C XLK = XL** (IPM2) C DO 30 I = 1,IPM2 K = IPM1 - I XLAMB(K+1) = XLAMB(K)*XLK XLK = XLK/XL 30 CONTINUE C 40 CONTINUE IF (INITM.GT.0) XLAMB(1) = FB - FA IF (INITM.LT.0) XLAMB(1) = WORK1(N+1) - WORK1(1) IF (ABS(XLAMB(1)).LE.EPS*1.0D4 .AND. IPM1.EQ.1) IPM1 = 0 C C********************************************************************** C C **** FIRST APPROXIMATION **** C ------------------- C DO 60 K = 1,NHALF - 1 ARGNOR = DBLE(K)/DBLE(N) CALL HPMIN1(ARGNOR,XLAMB,IPM1,HP,HPS) IF (INITM.GT.0) GO TO 50 WORK1(K+1) = WORK1(K+1) - HP WORK1(N-K+1) = WORK1(N-K+1) - HPS GO TO 60 50 CONTINUE APP = ARGNOR*XL WORK1(K+1) = FUN(APP+A) - HP WORK1(N-K+1) = FUN(B-APP) - HPS 60 CONTINUE C CALL HPMIN1(0.0D0,XLAMB,IPM1,HP,HPS) CALL HPMIN1(0.5D0,XLAMB,IPM1,HPHAL,DUMMY) IF (INITM.GT.0) GO TO 70 WORK1(1) = (WORK1(1)+WORK1(N+1)-(HP+HPS)) * 0.5D0 WORK1(NHALF+1) = WORK1(NHALF+1) - HPHAL WORK1(N+1) = 0.D0 GO TO 80 70 CONTINUE WORK1(1) = 0.5D0* (FA+FB-(HP+HPS)) WORK1(NHALF+1) = FUN(A+XL/2.0D0) - HPHAL C 80 CONTINUE C C DFT -N POINTS ON WORK1 - IN PLACE. C CALL RFFTI(N,WORK3) CALL RFFTF(N,WORK1,WORK3) C OVN = 1./DBLE(N) WORK1(N+2) = 0.D0 DO 90 K=N+1,3,-1 WORK1(K) = OVN*WORK1(K-1) IF(MOD(K,2).EQ.0) WORK1(K) = -WORK1(K) 90 CONTINUE WORK1(2) = 0.D0 WORK1(1) = OVN*WORK1(1) C********************************************************************** C **** MAIN LOOP **** C ----------- C C (REPEAT UNTILL THE STOPPING CRITERION IS SATISFIED) C 100 CONTINUE C WORK2(1) = WORK1(1) WORK2(2) = WORK1(2) WORK2(N+1) = WORK1(N+1) WORK2(N+2) = WORK1(N+2) DO 110 K=3,N WORK2(K) = 2.0D0*WORK1(K) 110 CONTINUE C CALL TRASF(WORK2,N,1) C TRIGG1=0.D0 TRIGG2=0.D0 DO 120 K=1,N+2,2 TRIGG1=TRIGG1+WORK2(K) TRIGG2=TRIGG2+WORK2(K+1) 120 CONTINUE C 130 CONTINUE NQUART = NHALF NHALF = N N = 2*N C * * C * ESTIMATE THE ACCURACY OF THE APPROXIMATION * C * * IF (INITM.GT.0) GO TO 140 GIBB1 = 0.D0 GIBB2 = 0.D0 GO TO 150 C 140 CONTINUE ARGNOR = 1./DBLE(N) AN = ARGNOR*XL CALL HPMIN1(ARGNOR,XLAMB,IPM1,HP,HPS) APP1 = FUN(AN+A) - HP APP2 = FUN(B-AN) - HPS GIBB1 = ABS(APP1-TRIGG1) GIBB2 = ABS(APP2-TRIGG2) 150 CONTINUE EST = SQRT((WORK1(NHALF-1)*2.D0)**2+ (WORK1(NHALF)*2.D0)**2+ + WORK1(NHALF+1)**2) ESTER = MAX(2.D0*GIBB1,2.D0*GIBB2,EST,EPS) C * * C * THIS IS THE STOPPING CRITERION. IF SATISFIED JUMP OUT OF C * THE MAIN LOOP * C * * C * IF THE ACCURACY REQUIREMENT IS ACHIEVED STOP THE UPDATING * C * * IF (ESTER.LE.TOL) GO TO 240 C * * C * * C * IF IT IS IMPOSSIBLE TO ACHIEVE THE USER REQUIRED ACCURACY * C * WITH THE COST MAXM THEN FOURIER COEFFICIENTS OF BEST ACCURACY * C * COMPUTABLE BY THE PROCEDURE ARE RETURNED. * C * * IF (N.GT.MAXM) GO TO 230 C IF (RNDLEV) GO TO 160 IF (NHALF.EQ.4) ESTMAX = ESTER IF (ESTMAX.LT.ESTER) ESTMAX = ESTER IF (NHALF.GE.8 .AND. ESTER.LT.ABNC1*ESTMAX) RNDLEV = .TRUE. GO TO 170 C * * C * IF THE ROUNDOFF LEVEL HAS BEEN ACHIEVED AND THE PRESENT * C * ESTIMATED ERROR IS GT THE PREVIOUS ONE TIMES THE CONSTANT * C * ABCN2 THEN STOP THE UPDATING STAGE * C * * 160 CONTINUE C IF (ESTER.GE.ABNC2*PREVES) GO TO 230 C 170 CONTINUE C * * C * COMPUTE NEW APPROXIMATION * C * * PREVES = ESTER C WORK2(1) = APP1 WORK2(NHALF) = APP2 C C DO 180 K = 2,NQUART ARGNOR = DBLE(2*K-1)/DBLE(N) AN = ARGNOR*XL CALL HPMIN1(ARGNOR,XLAMB,IPM1,HP,HPS) WORK2(K) = FUN(AN+A) - HP WORK2(NHALF+1-K) = FUN(B-AN) - HPS 180 CONTINUE C * * C * CALLING FFT + TRASF A MIDPOINT TRAPEZOIDAL RULE IS COMPUTED * C * * C DFT -NHALF POINTS ON WORK2 - IN PLACE C TEMP=WORK2(NHALF) DO 190 K=NHALF-1,1,-1 WORK2(K+1)=WORK2(K) 190 CONTINUE C WORK2(1)=TEMP WORK2(NHALF+1)=0.D0 WORK2(NHALF+2)=0.D0 CALL RFFTI(NHALF,WORK3) CALL RFFTF(NHALF,WORK2,WORK3) C OVNH = 1./DBLE(NHALF) WORK2(NHALF+2) = 0.D0 DO 200 K=NHALF+1,3,-1 WORK2(K) = OVNH*WORK2(K-1) IF(MOD(K,2).EQ.0) WORK2(K) = -WORK2(K) 200 CONTINUE WORK2(2) = 0.D0 WORK2(1) = OVNH*WORK2(1) C CALL TRASF(WORK2,NHALF,0) C DO 210 K = 1,NHALF,2 IK=N+3-K WORK1(IK) = 0.5* (WORK2(K+1)-WORK1(K+1)) WORK1(IK-1) = 0.5* (WORK1(K)-WORK2(K)) 210 CONTINUE C DO 220 K = 1,NHALF+2 WORK1(K) = 0.5* (WORK1(K)+WORK2(K)) 220 CONTINUE C GO TO 100 C C **** END OF MAIN LOOP **** C*********************************************************************** C **** COMPUTE FINAL RESULTS **** C --------------------- C 230 CONTINUE IERR = IPM1 IF (RNDLEV .AND. ESTER.GE.ABNC2*PREVES) IERR = -2 ESTER = -ESTER C 240 CONTINUE C MACT = NHALF+3 C(1) = WORK1(1) S(1) = 0.0 IF (NCOEF.EQ.1) RETURN C * * C * PERIODIC FUNCTIONS AND NON PERIODIC FUNCTIONS WITH * C * PERIODIC DERIVATIVES ARE TREATED AS A SPECIAL CASE * C * * IF (IPM1.GT.0) GO TO 270 C C DO 260 K = 2,NCOEF IF (K.GT.NQUART) GO TO 250 C(K) = 2.0D0*WORK1(2*K-1) S(K) = 2.0D0*WORK1(2*K) GO TO 260 C 250 CONTINUE C(K) = 0.0D0 IF (K.EQ.NQUART+1) C(K) = WORK1(NHALF+1) S(K) = 0.0D0 260 CONTINUE C IF (IERR.EQ.0.AND.ESTER.LT.0.D0) IERR = 1 C RETURN C 270 CONTINUE IF (IPM1.NE.1) GO TO 300 C C DO 290 K = 2,NCOEF IF (K.GT.NQUART) GO TO 280 C(K) = 2.0D0*WORK1(2*K-1) S(K) = -2.0D0*XLAMB(1)/ (TWOPI*DBLE(K-1)) + 2.0D0*WORK1(2*K) GO TO 290 C 280 CONTINUE C(K) = 0.0D0 IF (K.EQ.NQUART+1) C(K) = WORK1(NHALF+1) S(K) = -2.0D0*XLAMB(1)/ (TWOPI*DBLE(K-1)) 290 CONTINUE C C RETURN C 300 CONTINUE C * * C * THIS IS THE GENERAL CASE * C * * C DO 340 K = 2,NCOEF C OV2PK = 1.D0 / (TWOPI*DBLE(K-1)) OV2PK2 = OV2PK*OV2PK CFACT = OV2PK2 SFACT = -OV2PK CSUM = XLAMB(2) * CFACT SSUM = XLAMB(1) * SFACT IF (IPM1.EQ.2) GO TO 320 DO 310 KK=3,IPM1,2 CFACT = -CFACT * OV2PK2 SFACT = -SFACT * OV2PK2 IF (KK+1.LE.IPM1) CSUM = CSUM + XLAMB(KK+1) * CFACT SSUM = SSUM + XLAMB(KK) * SFACT 310 CONTINUE C 320 CONTINUE IF (K.GT.NQUART) GO TO 330 C(K) = 2.0D0*(CSUM + WORK1(2*K-1)) S(K) = 2.0D0*(SSUM + WORK1(2*K)) GO TO 340 C 330 CONTINUE C(K) = 2.0D0*CSUM S(K) = 2.0D0*SSUM IF (K.EQ.NQUART+1) C(K) = C(K) + WORK1(NHALF+1) C 340 CONTINUE C C C C RETURN C*********************************************************************** C C ************** END OF SUBROUTINE FCOEBA *************** C C*********************************************************************** C END SUBROUTINE HPMIN1(ARGNOR,XLAMB,IPM1,HP,HPS) C*********************************************************************** C C BEGINNING OF SUBROUTINE ****** HPMIN1 ****** C C C*********************************************************************** C C C HPMIN1 EVALUATES THE POLYNOMIAL C C HPMIN1(X)=SUM(I=1,IPM1)(XLAMB(I)*BERNOULLI(I,X)/FACTORIAL(I)) C C WHERE BERNOULLI(I,X) IS THE I-TH BERNOULLI POLYNOMIAL C C AT THE POINTS ..ARGNOR.. AND 1-ARGNOR. C IT RETURNS SUCH THE VALUES IN HP AND HPS RESPECTIVELY C C THE DEGREE OF POLYNOMIAL HPMIN1 IS IPM1. MAXIMUM VALUE FOR C IPM1 IS 15. C C C INPUT - ARGNOR,XLAMB,IPM1 C ARGNOR - DOUBLE PRECISION - (0,1)-NORMALIZED ARGUMENT C XLAMB -DOUBLE PRECISION(IPM1)-ARRAY OF FUNCTION DERIVATIVES C SEE INITIAL COMMENTS OF FCOEBA. C IPM1 - INTEGER- THE DEGREE OF POLYNOMIAL HPMIN1 C C OUTPUT - HP,HPS C HP - DOUBLE PRECISION - THE VALUE OF HPMIN1 AT POINT ARGNO C HPS - DOUBLE PRECISION - THE VALUE OF HPMIN1 AT POINT 1-AR C C*********************************************************************** C*********************************************************************** C FIRST DECLARATION STATEMENT OF HPMIN1 C*********************************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER IPM1 C DIMENSION XLAMB(*) C DIMENSION BERPOL(15) LOGICAL ODD C*********************************************************************** C FIRST EXECUTABLE STATEMENT OF HPMIN1 C*********************************************************************** IF (IPM1.GT.0) GO TO 10 HP = 0.0D0 HPS = 0.0D0 RETURN C 10 CONTINUE BERPOL(1) = ARGNOR - 0.5D0 IF (IPM1.EQ.1) GO TO 20 Y = ARGNOR* (ARGNOR-1.0D0) Z = Y* (ARGNOR-0.5D0) BERPOL(2) = Y + (1.0D0/6.0D0) IF (IPM1.EQ.2) GO TO 20 BERPOL(3) = Z IF (IPM1.EQ.3) GO TO 20 YTO2 = Y*Y BERPOL(4) = YTO2 - (1.0D0/30.0D0) IF (IPM1.EQ.4) GO TO 20 BERPOL(5) = (Y- (1.0D0/3.0D0))*Z IF (IPM1.EQ.5) GO TO 20 BERPOL(6) = Y*YTO2 - 0.5D0*YTO2 + (1.0D0/42.0D0) IF (IPM1.EQ.6) GO TO 20 Y2 = (ARGNOR-0.5D0)* (ARGNOR-0.5D0) Y4 = Y2*Y2 BERPOL(7) = (Y4-1.5D0*Y2+ (31.0D0/48.0D0))*Z IF (IPM1.EQ.7) GO TO 20 X2 = ARGNOR*ARGNOR X4 = X2*X2 X6 = X4*X2 X7 = X6*ARGNOR X8 = X4*X4 TERM = (2.0D0*X2-7.0D0*X4+14.0D0*X6)/3.0D0 BERPOL(8) = X8 - 4.0D0*X7 + TERM - (1.0D0/30.0D0) IF (IPM1.EQ.8) GO TO 20 TERM = X8 - 4.5D0*X7 + 6.0D0*X6 - 4.2D0*X4 + 2.0D0*X2 - 0.3D0 BERPOL(9) = ARGNOR*TERM IF (IPM1.EQ.9) GO TO 20 X9 = ARGNOR*X8 X10 = X4*X6 TERM = X10 - 5.0D0*X9 + 7.5D0*X8 - 7.0D0*X6 + 5.0D0*X4 BERPOL(10) = TERM - 1.5D0*X2 + (5.0D0/66.0D0) IF (IPM1.EQ.10) GO TO 20 TERM = X10 - 5.5D0*X9 - 11.0D0* (X6-X4) - 5.5D0*X2 BERPOL(11) = ARGNOR* (TERM+ ((5.0D0+55.0D0*X8)/6.0D0)) IF (IPM1.EQ.11) GO TO 20 X12 = X4*X8 TERM = X12 - 6.0D0*ARGNOR*X10 + 11.0D0*X10 - 16.5D0* (X8+X4) BERPOL(12) = TERM + 22.0D0*X6 + 5.0D0*X2 - (691.0D0/2730.0D0) IF (IPM1.EQ.12) GO TO 20 TERM = X12 + 13.0D0*X10 - 143.0D0/6.0D0*X8 + 286.0D0/7.0D0*X6 - + 42.9D0*X4 BERPOL(13) = (TERM+65.0D0/3.0D0*X2-691.0D0/210.0D0)*ARGNOR - + 6.5D0*X12 IF (IPM1.EQ.13) GO TO 20 X14 = X12*X2 TERM = X14 + (-7.0D0*ARGNOR+91.0D0/6.0D0)*X12 - 100.1/3.0D0*X10 BERPOL(14) = TERM + 71.5*X8 - 100.1D0*X6 + 455.0D0/6.0D0*X4 - + 69.1D0/3.0D0*X2 + 7.0D0/6.0D0 IF (IPM1.EQ.14) GO TO 20 TERM = X14 + 17.5D0*X12 - 45.5D0*X10 + 715.0D0/6.0D0*X8 - + 214.5D0*X6 BERPOL(15) = (TERM+227.5D0*X4-691.0D0/6.0D0*X2+17.5D0)*ARGNOR - + 7.5D0*X14 C 20 CONTINUE HP = XLAMB(1)*BERPOL(1) HPS = -HP FATT = 1.0D0 ODD = .FALSE. IF (IPM1.EQ.1) RETURN C C DO 40 I = 2,IPM1 FATT = FATT*DBLE(I) TEMP = XLAMB(I)*BERPOL(I)/FATT HP = HP + TEMP IF (ODD) GO TO 30 HPS = HPS + TEMP ODD = .TRUE. GO TO 40 C 30 CONTINUE ODD = .FALSE. HPS = HPS - TEMP 40 CONTINUE C C RETURN C********************************************************************** C C END OF SUBROUTINE ****** HPMIN1 ****** C C********************************************************************** END SUBROUTINE TRASF(VETT,N,NCODE) C C********************************************************************** C C BEGINNING OF SUBROUTINE ****** TRASF ****** C C********************************************************************** C C GIVEN THE INPUT ARRAY VETT OF N+2 COMPONENTS TRASF COMPUTES IN C PLACE THE VALUES (I=1,2,..,N/2+1) C C A(I)=VETT(I)*COS(PI*(I-1)/N)+VETT(I+1)*SIN(PI*(I-1)/N) C B(I)=VETT(I)*COS(PI*(I-1)/N)-VETT(I+1)*SIN(PI*(I-1)/N) IF NCODE=1 C B(I)=-VETT(I)*SIN(PI*(I-1)/N)+VETT(I+1)*COS(PI*(I-1)/N) IF NCODE= C C IT RETURNS A(I) IN VETT(2*I-1) AND B(I) IN VETT(2*I) C C C C INPUT - VETT,N,NCODE C C VETT -DOUBLE PRECISION(N+2)- ARRAY TO BE TRANSFORMED C C N -INTEGER- N+2 IS THE SIZE OF VETT C C NCODE -INTEGER- ALLOWS THE CHOICE OF B(I) C C C OUTPUT - VETT C C VETT -DOUBLE PRECISION(N+2)- TRANSFORMED ARRAY C C*********************************************************************** C*********************************************************************** C FIRST DECLARATION STATEMENT OF TRASF C*********************************************************************** C IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER N,NCODE C INTEGER N2 DIMENSION VETT(*) C * * C * THESE ARE THE VALUES OF PI AND COS(PI / 4.) * C * WITH 18 DECIMAL FIGURES. CHANGE THIS CARD IF A * C * MORE PRECISE ARITHMETIC IS USED. * C * * DATA PI/0.D0/,COSPIQ/0.D0/ C C*********************************************************************** C FIRST EXECUTABLE STATEMENT OF TRASF C*********************************************************************** C **** INITIALIZATION **** C IF (PI.EQ.0.D0) PI = 4.D0*ATAN(1.D0) IF (COSPIQ.EQ.0.D0) COSPIQ = SQRT(0.5D0) PIN = PI/N COSPIN = COS(PIN) SINPIN = SIN(PIN) N2 = N/2 C C *** COMPUTATION OF COS((PI/2N)(K-1)) , SIN((PI/2N)(K-1)) *** C DO 60 K = 1,N2,2 IF (K.NE.1) GO TO 10 C = 1 S = 0.0 GO TO 20 C 10 CONTINUE APP1 = C*COSPIN - S*SINPIN S = S*COSPIN + C*SINPIN C = APP1 20 CONTINUE C **** COMPUTE THE SUM **** C APP1 = VETT(K+1)*S APP2 = VETT(K)*S APP3 = VETT(K)*C VETT(K) = APP3 + APP1 IF(NCODE.EQ.0) GO TO 30 VETT(K+1)=APP3-APP1 GO TO 40 30 CONTINUE VETT(K+1) = VETT(K+1)*C - APP2 40 CONTINUE I = N + 3 - K APP1 = VETT(I)*C APP2 = VETT(I-1)*C APP3 = VETT(I-1)*S VETT(I-1) = APP3 + APP1 IF(NCODE.EQ.0) GO TO 50 VETT(I)= APP3 - APP1 GO TO 60 50 CONTINUE VETT(I) = VETT(I)*S - APP2 60 CONTINUE C I = N2 + 1 APP1 = VETT(I+1)*COSPIQ APP2 = VETT(I)*COSPIQ VETT(I) = APP2 + APP1 IF(NCODE.EQ.0) GO TO 70 VETT(I+1)=APP2-APP1 GO TO 80 70 CONTINUE VETT(I+1) = APP1 - APP2 80 CONTINUE C RETURN C C*********************************************************************** C C END OF SUBROUTINE ****** TRASF ****** C C*********************************************************************** C********************************************************************** END C C COLLECTION OF NUMERICAL DIFFERENTIATION ROUTINES FOR C THE PACKAGE FOURCO. C DOUBLE PRECISION VERSION C SUBROUTINE ENDACE(XVAL,NDER,HBASE,DER,EREST,FUN) C C C *****EVALUATE NUMERICAL DERIVATIVES AND CORRESPONDING ERRORS***** C C C *****PURPOSE***** C C A SUBROUTINE FOR NUMERICAL DIFFERENTIATION AT A POINT. C IT RETURNS A SET OF APPROXIMATIONS TO THE J-TH ORDER DERIVATIVE C (J=1,2,...,14) OF FUN(X) EVALUATED AT X=XVAL AND, FOR EACH C DERIVATIVE, AN ERROR ESTIMATE ( WHICH INCLUDES THE EFFECT OF C AMPLIFICATION OF ROUND-OFF ERRORS ). C C C *** INPUT PARAMETERS *** C C C (1) XVAL DOUBLE PRECISION. THE ABSCISSA AT WHICH THE SET OF C DERIVATIVES IS REQUIRED. C C (2) NDER INTEGER. THE HIGHEST ORDER DERIVATIVE REQUIRED. C IF(NDER.GT.0) ALL DERIVATIVES UP TO MIN(NDER,14) C ARE CALCULATED. C IF(NDER.LT.0 AND NDER EVEN) EVEN ORDER DERIVATIVES C UP TO MIN(-NDER,14) ARE CALCULATED. C IF(NDER.LT.0 AND NDER ODD) ODD ORDER DERIVATIVES C UP TO MIN(-NDER,13) ARE CALCULATED. C C (3) HBASE DOUBLE PRECISION. A STEP LENGTH. C C (6) FUN THE NAME OF A DOUBLE PRECISION FUNCTION SUBPROGRAM, C WHICH IS REQUIRED BY THE ROUTINE AS A SUBPROGRAM C AND WHICH REPRESENTS THE FUNCTION BEING DIFFERENTIATED. C THE ROUTINE REQUIRES 21 FUNCTION EVALUATIONS FUN(X) C LOCATED AT X=XVAL AND AT C X=XVAL+(2*J-1)*HBASE, J=-9,-8,...,9,10. C THE FUNCTION VALUE AT X=XVAL IS DISREGARDED WHEN C ONLY ODD ORDER DERIVATIVES ARE REQUIRED. C C C *** OUTPUT PARAMETERS *** C C C (4) DER(J) J=1,2,...,14. DOUBLE PRECISION. A VECTOR OF C APPROXIMATIONS TO THE J-TH DERIVATIVE OF FUN(X) C EVALUATED AT X=XVAL. C C (5) EREST(J) J=1,2,...,14. DOUBLE PRECISION. A VECTOR OF C ESTIMATES OF THE ABSOLUTE ACCURACY OF DER(J). C THESE ARE NEGATIVE WHEN EREST(J).GT.ABS(DER(J)), OR C WHEN, FOR SOME OTHER REASON THE ROUTINE IS DOUBTFUL C ABOUT THE VALIDITY OF THE RESULT. C C C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION DER(14),EREST(14) DIMENSION RANGE(7),RTAB(10,7),TZEROM(10),ACOF(10) DATA ZERO,ONE,THRON2,TWO / 0.0D0,1.0D0,1.5D0,2.0D0/ C C C QUIT ON ZERO HBASE AND ZERO NDER. SET CONTROL PARAMETER NDPAR. C NDPAR = +1 ODD-ORDER DERIVATIVES ONLY. C NDPAR = 0 ALL DERIVATIVES. C NDPAR = -1 EVEN-ORDER DERIVATIVES ONLY. C C IF(HBASE.EQ.ZERO) GOTO 900 NDERP=NDER NDPAR=0 IF(NDER) 10,900,20 10 NDPAR=4*(NDER/2)-2*NDER-1 NDERP=-NDER 20 CONTINUE IF(NDERP.GT.14) NDERP=14 C C C NEXT, EVALUATE 21 FUNCTION VALUES, SET ACOF(N), AND SET FIRST C COLUMN OF RTAB FOR ODD ORDER DERIVATIVES AND STORE IN TZEROM C THE FIRST COLUMN OF RTAB FOR EVEN ORDER DERIVATIVES. C C FZ=FUN(XVAL) DO 30 N=1,10 NN1=2*N-1 NSQ=NN1*NN1 ACOF(N)=NSQ XNUM=NN1 H=HBASE*XNUM TEMP=FUN(XVAL+H) TERM=FUN(XVAL-H) RTAB(N,1)=(TEMP-TERM)/(TWO*H) TZEROM(N)=(TEMP-FZ-FZ+TERM)/(TWO*H**2) 30 CONTINUE C C C SET UP FOR ODD-ORDER DERIVATIVES. C C IF(NDPAR.EQ.-1) GOTO 40 NPAR=1 NSMAX=(NDERP+1)/2 IF(NSMAX.LE.0) GOTO 701 GOTO 60 C C C SET UP FOR EVEN-ORDER DERIVATIVES C C 40 CONTINUE NPAR=-1 IF(NDPAR.EQ.+1) GOTO 900 NSMAX=NDERP/2 IF(NSMAX.LE.0) GOTO 900 C C C SET THE FIRST COLUMN OF THE T-TABLE FOR EVEN ORDER DERIVATIVES C C DO 50 J=1,10 50 RTAB(J,1)=TZEROM(J) 60 CONTINUE C C ODD-ORDER AND EVEN ORDER PATHS JOIN UP HERE C NEGER=0 HSQ=HBASE*HBASE FACT=ONE HFACT=ONE DO 700 NS= 1,NSMAX C C C FIRST CALCULATE ELEMENTS OF T-TABLE FOR CURRENT NS VALUE C C NPMIN=NS IF(NS.EQ.1) NPMIN=2 DO 250 NP= NPMIN,7 NKTOP=10-NP+1 NPR=NP-NS+1 DO 200 NK= 1,NKTOP J=NK+NP-1 A=ACOF(J) B=ACOF(NK) TERM=ZERO TEMP=ZERO IF(NP.NE.NS) TEMP=A*RTAB(NK,NPR-1)-B*RTAB(NK+1,NPR-1) IF(NS.NE.1) TERM=-RTAB(NK,NPR)+RTAB(NK+1,NPR) RTAB(NK,NPR)= (TERM+TEMP)/(A-B) 200 CONTINUE 250 CONTINUE C C NOW CALCULATE NUP,NLO AND RANMIN. C DO 370 NP=NS,7 NTOP=11-NP NPR=NP-NS+1 XMAX=RTAB(1,NPR) XMIN=RTAB(1,NPR) NMAX=1 NMIN=1 C DO 340 NK=2,NTOP TEMP=RTAB(NK,NPR) IF(TEMP.GT.XMAX) NMAX=NK IF(TEMP.GT.XMAX) XMAX=TEMP IF(TEMP.LT.XMIN) NMIN=NK IF(TEMP.LT.XMIN) XMIN=TEMP 340 CONTINUE C RANGE(NP)=XMAX-XMIN IF(NP.NE.NS) GOTO 350 RANMIN=RANGE(NP) GOTO 355 350 CONTINUE IF(RANGE(NP).GE.RANMIN) GOTO 360 RANMIN=RANGE(NP) 355 CONTINUE NPMIN=NP NUP=NMAX NLO=NMIN 360 CONTINUE 370 CONTINUE C C TERM=ZERO NTOP=11-NPMIN NTOPM2=NTOP-2 XNUM=NTOPM2 J=NPMIN-NS+1 DO 440 NK= 1,NTOP IF(NK.EQ.NUP) GOTO 439 IF(NK.EQ.NLO) GOTO 439 TERM=TERM+RTAB(NK,J) 439 CONTINUE 440 CONTINUE TERM=TERM/XNUM C C JNS=2*NS-1 IF(NPAR.LT.0) JNS=2*NS XJNS=JNS FACT=FACT*XJNS DER(JNS)=TERM*FACT/HFACT IF(JNS.EQ.10) RANMIN=THRON2*RANMIN IF(JNS.EQ.11) RANMIN=THRON2*RANMIN IF(JNS.GE.12) RANMIN=TWO*RANMIN EREST(JNS)=RANMIN*FACT/HFACT C C SET SIGN OF EREST. EREST NEGATIVE EITHER IF IT SWAMPS DER,OR IF C TWO PREVIOUS CONSECUTIVE ERESTS OF THIS PARITY ARE NEGATIVE. C IT MAY ALSO BE SET NEGATIVE AT END (750). C IF(NEGER.GE.2) EREST(JNS)=-EREST(JNS) IF(NEGER.GE.2)GOTO 690 IF(TERM.LT.RANMIN) EREST(JNS)=-EREST(JNS) IF(-TERM.GE.RANMIN) EREST(JNS)=-EREST(JNS) IF(EREST(JNS).LT.ZERO) NEGER=NEGER+1 IF(EREST(JNS).GE.ZERO) NEGER=0 690 CONTINUE C HFACT=HSQ*HFACT FACT=FACT*(XJNS+ONE) 700 CONTINUE 701 CONTINUE IF(NPAR.GT.0)GOTO 40 C C SET SIGN OF EREST NEGATIVE IF TWO PREVIOUS CONSECUTIVE ERESTS C ARE NEGATIVE. C IF(NDPAR.NE.0)GOTO 900 IF(NDERP.LE.2) GOTO 900 NEGER=0 ERPREV=EREST(1) DO 760 J= 2,NDERP ERNOW=EREST(J) IF(NEGER.EQ.2) GOTO 740 IF(ERPREV.GE.ZERO) GOTO 750 IF(ERNOW.GE.ZERO) GOTO 750 740 NEGER=2 IF(ERNOW.GT.ZERO) EREST(J)=-ERNOW 750 ERPREV=ERNOW 760 CONTINUE C 900 CONTINUE RETURN C*** *** *** END OF ENDACE *** *** *** C*********************************************************** END SUBROUTINE CHPEND(XVAL,HBASE,DER,FUN) C C ********************************************************************** C C ***** SUBROUTINE CHPEND ***** C C THIS IS A MODIFIED VERSION OF ENDACE. IT COMPUTES A SET OF C APPROXIMATIONS TO THE J-TH DERIVATIVE (J=1,2,..,5) OF FUN(X) C AT X=XVAL. C NO ERROR ESTIMATE IS RETURNED. C C INPUT - XVAL,HBASE,FUN C C XVAL -DOUBLE PRECISION- THE POINT AT WHICH THE DERIVATIVES C ARE REQUIRED. C C HBASE -DOUBLE PRECISION- THE STEP LENGHT C C FUN -DOUBLE PRECISION FUNCTION- THE FUNCTION BEING DIFFEREN C TIATED. THE ROUTINE REQUIRES 9 EVALUATUIONS OF ..FUN.. C IN THE INTERVAL (XVAL-7*HBASE,XVAL+7*HBASE). C C C OUTPUT - DER C C DER -DOUBLE PRECISION(5)- ARRAY OF FUNCTION DERIVATIVES.D C CONTAINS THE J-TH DERIVATIVE OF FUN. C C C C REFERENCE - J.N.LYNESS - SOFTWARE FOR NUMERICAL DIFFERENTIATION C*********************************************************************** C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION DER(5) DIMENSION RTAB(4,3),ACOF(4),TZEROM(4),RANGE(3) C EXTERNAL FUN * DATA ZERO,ONE,TWO/0.0D0,1.0D0,2.0D0/ C C FZ = FUN(XVAL) DO 10 N = 1,4 NN1 = 2*N - 1 NSQ = NN1*NN1 ACOF(N) = DBLE(NSQ) XNUM = DBLE(NN1) H = HBASE*XNUM TEMP = FUN(XVAL+H) TERM = FUN(XVAL-H) RTAB(N,1) = (TEMP-TERM)/ (TWO*H) TZEROM(N) = (TEMP-FZ-FZ+TERM)/ (TWO*H**2) 10 CONTINUE NPAR = 1 NSMAX = 3 GO TO 40 C 20 CONTINUE NSMAX = 2 NPAR = -1 DO 30 J = 1,4 RTAB(J,1) = TZEROM(J) 30 CONTINUE 40 CONTINUE HSQ = HBASE*HBASE FACT = ONE HFACT = ONE DO 150 NS = 1,NSMAX NPMIN = NS IF (NS.EQ.1) NPMIN = 2 DO 60 NP = NPMIN,3 NKTOP = 5 - NP NPR = NP - NS + 1 DO 50 NK = 1,NKTOP J = NP + NK - 1 A = ACOF(J) B = ACOF(NK) TERM = ZERO TEMP = ZERO IF (NP.NE.NS) TEMP = A*RTAB(NK,NPR-1) - + B*RTAB(NK+1,NPR-1) IF (NS.NE.1) TERM = -RTAB(NK,NPR) + RTAB(NK+1,NPR) RTAB(NK,NPR) = (TERM+TEMP)/ (A-B) 50 CONTINUE 60 CONTINUE C DO 110 NP = NS,3 NTOP = 5 - NP NPR = NP - NS + 1 XMAX = RTAB(1,NPR) XMIN = RTAB(1,NPR) NMAX = 1 NMIN = 1 DO 70 NK = 2,NTOP TEMP = RTAB(NK,NPR) IF (TEMP.GT.XMAX) NMAX = NK IF (TEMP.GT.XMAX) XMAX = TEMP IF (TEMP.LT.XMIN) NMIN = NK IF (TEMP.LT.XMIN) XMIN = TEMP 70 CONTINUE RANGE(NP) = XMAX - XMIN IF (NP.NE.NS) GO TO 80 RANMIN = RANGE(NP) GO TO 90 C 80 CONTINUE IF (RANGE(NP).GE.RANMIN) GO TO 100 RANMIN = RANGE(NP) 90 CONTINUE NPMIN = NP NLO = NMIN 100 CONTINUE 110 CONTINUE C TERM = ZERO NTOP = 5 - NPMIN J = NPMIN - NS + 1 DO 120 NK = 1,NTOP IF (NK.EQ.NLO) GO TO 120 TERM = TERM + RTAB(NK,J) 120 CONTINUE TERM = TERM/ (NTOP-1) JNS = 2*NS - 1 IF (NPAR.LT.0) JNS = 2*NS XJNS = JNS FACT = FACT*XJNS DER(JNS) = TERM*FACT/HFACT 140 CONTINUE HFACT = HSQ*HFACT FACT = FACT* (XJNS+ONE) 150 CONTINUE IF (NPAR.EQ.1) GO TO 20 C RETURN C ******************************************************************* C * * C * END OF SUBROUTINE ****** CHPEND ****** * C * * C ******************************************************************* END SUBROUTINE ONESID(NCODE,XVAL,HBASE,DER,FUN) C C ********************************************************************** C C ***** SUBROUTINE ONESID ***** C C THIS IS A MODIFIED VERSION OF ENDACE. IT COMPUTES A SET OF C APPROXIMATIONS TO THE J-TH DERIVATIVE (J=1,2,..,5) OF FUN(X) C AT X=XVAL. C ONE SIDED DIFFERENCE APPROXIMATIONS ARE USED. C NO ERROR ESTIMATE IS RETURNED. C C INPUT - NCODE,XVAL,HBASE,FUN C C NCODE -INTEGER- IF=1 SELECT RIGHT ONE SIDED DIFFERENCE C IF=-1 SELECT LEFT ONE SIDED DIFFERENCE C C XVAL -DOUBLE PRECISION- THE POINT AT WHICH THE DERIVATIVES C ARE REQUIRED. C C HBASE -DOUBLE PRECISION- THE STEP LENGHT C C FUN -DOUBLE PRECISION FUNCTION- THE FUNCTION BEING DIFFEREN C TIATED. THE ROUTINE REQUIRES 7 EVALUATIONS OF ..FUN.. C IN THE INTERVAL (XVAL,XVAL+6*HBASE) IF NCODE=1 C (XVAL-6*HBASE,XVAL) IF NCODE=-1. C C C OUTPUT - DER C C DER -DOUBLE PRECISION(5)- ARRAY OF FUNCTION DERIVATIVES.DE C CONTAINS THE J-TH DERIVATIVE OF FUN. C C C C REFERENCE - J.N.LYNESS - SOFTWARE FOR NUMERICAL DIFFERENTIATION C*********************************************************************** C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION DER(5) DIMENSION RTAB(7,6),ACOF(7),RANGE(7) C EXTERNAL FUN * DATA ZERO,ONE/0.0D0,1.0D0/ C C RTAB(1,1)=FUN(XVAL) ACOF(1)=ZERO DO 10 N = 2,7 ACOF(N) = DBLE(N-1) H = HBASE*ACOF(N) IF(NCODE.EQ.-1)TEMP=FUN(XVAL-H) IF(NCODE.EQ.1)TEMP=FUN(XVAL+H) RTAB(N,1)=TEMP 10 CONTINUE NSMAX=6 FACT = ONE DO 150 NS = 1,NSMAX NPMIN = NS IF (NS.EQ.1) NPMIN = 2 DO 60 NP = NPMIN,6 NKTOP = 8 - NP NPR = NP - NS + 1 DO 50 NK = 1,NKTOP J = NP + NK - 1 A = ACOF(J)*HBASE B = ACOF(NK)*HBASE IF(NCODE.EQ.-1)A=-A IF(NCODE.EQ.-1)B=-B TERM = ZERO TEMP = ZERO IF (NP.NE.NS) TEMP = A*RTAB(NK,NPR-1) - + B*RTAB(NK+1,NPR-1) IF (NS.NE.1) TERM = -RTAB(NK,NPR) + RTAB(NK+1,NPR) RTAB(NK,NPR) = (TERM+TEMP)/ (A-B) 50 CONTINUE 60 CONTINUE DO 110 NP = NS,6 NTOP = 8 - NP NPR = NP - NS + 1 XMAX = RTAB(1,NPR) XMIN = RTAB(1,NPR) NMAX = 1 NMIN = 1 DO 70 NK = 2,NTOP TEMP = RTAB(NK,NPR) IF (TEMP.GT.XMAX) NMAX = NK IF (TEMP.GT.XMAX) XMAX = TEMP IF (TEMP.LT.XMIN) NMIN = NK IF (TEMP.LT.XMIN) XMIN = TEMP 70 CONTINUE RANGE(NP) = XMAX - XMIN IF (NP.NE.NS) GO TO 80 RANMIN = RANGE(NP) GO TO 90 C 80 CONTINUE IF (RANGE(NP).GE.RANMIN) GO TO 100 RANMIN = RANGE(NP) 90 CONTINUE NPMIN = NP NLO = NMIN 100 CONTINUE 110 CONTINUE C TERM = ZERO NTOP = 8 - NPMIN J = NPMIN - NS + 1 DO 120 NK = 1,NTOP IF (NK.EQ.NLO) GO TO 120 TERM = TERM + RTAB(NK,J) 120 CONTINUE TERM = TERM/ (NTOP-1) JNS=NS XJNS = MAX(JNS-1,1) FACT = FACT*XJNS IF(JNS.EQ.1) GOTO 150 DER(JNS-1)=TERM*FACT 150 CONTINUE C RETURN C ******************************************************************* C * * C * END OF SUBROUTINE ****** ONESIDE ****** * C * * C ******************************************************************* END SUBROUTINE ONESMP(NCODE,FUNVAL,IXVAL,HBASE,DER) C C ********************************************************************** C C ***** SUBROUTINE ONESMP ***** C C THIS IS A MODIFIED VERSION OF ENDACE. IT COMPUTES A SET OF C APPROXIMATIONS TO THE J-TH DERIVATIVE (J=1,2,..,5) OF FUN(X) C AT X=XVAL. C IT REQUIRES A SET OF FUNCTION VALUES STORED IN AN INPUT ARRAY. C ONE SIDED DIFFERENCE APPROXIMATIONS ARE USED. C NO ERROR ESTIMATE IS RETURNED. C C INPUT - NCODE,FUNVAL,IXVAL,HBASE C C NCODE -INTEGER- IF=1 SELECT RIGHT ONE SIDED DIFFERENCE C IF=-1 SELECT LEFT ONE SIDED DIFFERENCE C C FUNVAL -DOUBLE PRECISION(*)- ARRAY OF THE VALUES OF THE C FUNCTION TO BE DIFFERENTIATED. C FUNCTION VALUES MUST BE AT EQUALLY SPACED POINTS C WITH STEP ..HBASE.. C C IXVAL -INTEGER- THE POINT WHERE THE DERIVATIVES ARE C REQUIRED IS STORED IN FUNVAL(IXVAL). C C HBASE -DOUBLE PRECISION- THE STEP LENGHT C C C OUTPUT - DER C C DER -DOUBLE PRECISION(5)- ARRAY OF FUNCTION DERIVATIVES.DE C CONTAINS THE J-TH DERIVATIVE OF FUN. C C C C REFERENCE - J.N.LYNESS - SOFTWARE FOR NUMERICAL DIFFERENTIATION C*********************************************************************** C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION DER(5),FUNVAL(*) DIMENSION RTAB(7,6),ACOF(7),RANGE(7) C EXTERNAL FUN * DATA ZERO,ONE/0.0D0,1.0D0/ C C RTAB(1,1)=FUNVAL(IXVAL) ACOF(1)=ZERO DO 10 N = 2,7 ACOF(N) = DBLE(N-1) H = HBASE*ACOF(N) IF(NCODE.EQ.-1)TEMP=FUNVAL(IXVAL-N+1) IF(NCODE.EQ.1)TEMP=FUNVAL(IXVAL+N-1) RTAB(N,1)=TEMP 10 CONTINUE NSMAX=6 FACT = ONE DO 150 NS = 1,NSMAX NPMIN = NS IF (NS.EQ.1) NPMIN = 2 DO 60 NP = NPMIN,6 NKTOP = 8 - NP NPR = NP - NS + 1 DO 50 NK = 1,NKTOP J = NP + NK - 1 A = ACOF(J)*HBASE B = ACOF(NK)*HBASE IF(NCODE.EQ.-1)A=-A IF(NCODE.EQ.-1)B=-B TERM = ZERO TEMP = ZERO IF (NP.NE.NS) TEMP = A*RTAB(NK,NPR-1) - + B*RTAB(NK+1,NPR-1) IF (NS.NE.1) TERM = -RTAB(NK,NPR) + RTAB(NK+1,NPR) RTAB(NK,NPR) = (TERM+TEMP)/ (A-B) 50 CONTINUE 60 CONTINUE DO 110 NP = NS,6 NTOP = 8 - NP NPR = NP - NS + 1 XMAX = RTAB(1,NPR) XMIN = RTAB(1,NPR) NMAX = 1 NMIN = 1 DO 70 NK = 2,NTOP TEMP = RTAB(NK,NPR) IF (TEMP.GT.XMAX) NMAX = NK IF (TEMP.GT.XMAX) XMAX = TEMP IF (TEMP.LT.XMIN) NMIN = NK IF (TEMP.LT.XMIN) XMIN = TEMP 70 CONTINUE RANGE(NP) = XMAX - XMIN IF (NP.NE.NS) GO TO 80 RANMIN = RANGE(NP) GO TO 90 C 80 CONTINUE IF (RANGE(NP).GE.RANMIN) GO TO 100 RANMIN = RANGE(NP) 90 CONTINUE NPMIN = NP NLO = NMIN 100 CONTINUE 110 CONTINUE C TERM = ZERO NTOP = 8 - NPMIN J = NPMIN - NS + 1 DO 120 NK = 1,NTOP IF (NK.EQ.NLO) GO TO 120 TERM = TERM + RTAB(NK,J) 120 CONTINUE TERM = TERM/ (NTOP-1) JNS=NS XJNS = MAX(JNS-1,1) FACT = FACT*XJNS IF(JNS.EQ.1) GOTO 150 DER(JNS-1)=TERM*FACT 150 CONTINUE C RETURN C ******************************************************************* C * * C * END OF SUBROUTINE ****** ONESMP ****** * C * * C ******************************************************************* END C**************************************************************** C C COLLECTION OF FISHPACK ROUTINES FOR REAL FORWARD FFT C DOUBLE PRECISION VERSION C**************************************************************** C THESE ARE A COPY OF THE FOLLOWING 9 SUBROUTINES C AS THEY APPEAR IN C P.N.SWARZTRAUBER - FISHPACK, A PACKAGE OF FORTRAN SUBPROGRAMS C FOR THE FAST FOURIER TRANSFORM OF PERIODIC AND OTHER C SYMMETRIC SEQUENCES. IN PORTLIB,(1979) C ONLY CHANGES TO IMPROVE PORTABILITY OF THE CODE HAVE BEEN MADE. C C RFFTI,RFFTI1,RFFTF,RFFTF1,RADF2,RADF3,RADF4,RADF5,RADFG. C**************************************************************** SUBROUTINE RFFTI (N,WSAVE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION WSAVE(*) IF (N .EQ. 1) RETURN CALL RFFTI1 (N,WSAVE(N+1),WSAVE(2*N+1)) RETURN END SUBROUTINE RFFTI1 (N,WA,IFAC) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION IFAC DIMENSION WA(*) ,IFAC(*) ,NTRYH(4) DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/ DATA TPI/0.D0/ IF(TPI.EQ.0.D0)TPI=8.D0*ATAN(1.D0) NL = N NF = 0 J = 0 101 J = J+1 IF (J-4) 102,102,103 102 NTRY = NTRYH(J) GO TO 104 103 NTRY = NTRY+2 104 NQ = NL/NTRY NR = NL-NTRY*NQ IF (NR) 101,105,101 105 NF = NF+1 IFAC(NF+2) = NTRY NL = NQ IF (NTRY .NE. 2) GO TO 107 IF (NF .EQ. 1) GO TO 107 DO 106 I=2,NF IB = NF-I+2 IFAC(IB+2) = IFAC(IB+1) 106 CONTINUE IFAC(3) = 2 107 IF (NL .NE. 1) GO TO 104 IFAC(1) = N IFAC(2) = NF ARGH = TPI/DBLE(N) IS = 0 NFM1 = NF-1 L1 = 1 IF (NFM1 .EQ. 0) RETURN DO 110 K1=1,NFM1 IP = IFAC(K1+2) LD = 0 L2 = L1*IP IDO = N/L2 IPM = IP-1 DO 109 J=1,IPM LD = LD+L1 I = IS ARGLD = DBLE(LD)*ARGH FI = 0. DO 108 II=3,IDO,2 I = I+2 FI = FI+1. ARG = FI*ARGLD WA(I-1) = COS(ARG) WA(I) = SIN(ARG) 108 CONTINUE IS = IS+IDO 109 CONTINUE L1 = L2 110 CONTINUE RETURN END SUBROUTINE RFFTF (N,R,WSAVE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION R(*) ,WSAVE(*) IF (N .EQ. 1) RETURN CALL RFFTF1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1)) RETURN END SUBROUTINE RFFTF1 (N,C,CH,WA,IFAC) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION IFAC DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*) NF = IFAC(2) NA = 1 L2 = N IW = N DO 111 K1=1,NF KH = NF-K1 IP = IFAC(KH+3) L1 = L2/IP IDO = N/L2 IDL1 = IDO*L1 IW = IW-(IP-1)*IDO NA = 1-NA IF (IP .NE. 4) GO TO 102 IX2 = IW+IDO IX3 = IX2+IDO IF (NA .NE. 0) GO TO 101 CALL RADF4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3)) GO TO 110 101 CALL RADF4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3)) GO TO 110 102 IF (IP .NE. 2) GO TO 104 IF (NA .NE. 0) GO TO 103 CALL RADF2 (IDO,L1,C,CH,WA(IW)) GO TO 110 103 CALL RADF2 (IDO,L1,CH,C,WA(IW)) GO TO 110 104 IF (IP .NE. 3) GO TO 106 IX2 = IW+IDO IF (NA .NE. 0) GO TO 105 CALL RADF3 (IDO,L1,C,CH,WA(IW),WA(IX2)) GO TO 110 105 CALL RADF3 (IDO,L1,CH,C,WA(IW),WA(IX2)) GO TO 110 106 IF (IP .NE. 5) GO TO 108 IX2 = IW+IDO IX3 = IX2+IDO IX4 = IX3+IDO IF (NA .NE. 0) GO TO 107 CALL RADF5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4)) GO TO 110 107 CALL RADF5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4)) GO TO 110 108 IF (IDO .EQ. 1) NA = 1-NA IF (NA .NE. 0) GO TO 109 CALL RADFG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW)) NA = 1 GO TO 110 109 CALL RADFG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW)) NA = 0 110 L2 = L1 111 CONTINUE IF (NA .EQ. 1) RETURN DO 112 I=1,N C(I) = CH(I) 112 CONTINUE RETURN END SUBROUTINE RADF2 (IDO,L1,CC,CH,WA1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) , 1 WA1(*) DO 101 K=1,L1 CH(1,1,K) = CC(1,K,1)+CC(1,K,2) CH(IDO,2,K) = CC(1,K,1)-CC(1,K,2) 101 CONTINUE IF (IDO-2) 107,105,102 102 IDP2 = IDO+2 DO 104 K=1,L1 DO 103 I=3,IDO,2 IC = IDP2-I TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) CH(I,1,K) = CC(I,K,1)+TI2 CH(IC,2,K) = TI2-CC(I,K,1) CH(I-1,1,K) = CC(I-1,K,1)+TR2 CH(IC-1,2,K) = CC(I-1,K,1)-TR2 103 CONTINUE 104 CONTINUE IF (MOD(IDO,2) .EQ. 1) RETURN 105 DO 106 K=1,L1 CH(1,2,K) = -CC(IDO,K,2) CH(IDO,1,K) = CC(IDO,K,1) 106 CONTINUE 107 RETURN END SUBROUTINE RADF3 (IDO,L1,CC,CH,WA1,WA2) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) , 1 WA1(*) ,WA2(*) DATA TAUR,TAUI /-.5D0,0.D0/ IF(TAUI.EQ.0.D0)TAUI=SQRT(0.75D0) DO 101 K=1,L1 CR2 = CC(1,K,2)+CC(1,K,3) CH(1,1,K) = CC(1,K,1)+CR2 CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2)) CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2 101 CONTINUE IF (IDO .EQ. 1) RETURN IDP2 = IDO+2 DO 103 K=1,L1 DO 102 I=3,IDO,2 IC = IDP2-I DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) CR2 = DR2+DR3 CI2 = DI2+DI3 CH(I-1,1,K) = CC(I-1,K,1)+CR2 CH(I,1,K) = CC(I,K,1)+CI2 TR2 = CC(I-1,K,1)+TAUR*CR2 TI2 = CC(I,K,1)+TAUR*CI2 TR3 = TAUI*(DI2-DI3) TI3 = TAUI*(DR3-DR2) CH(I-1,3,K) = TR2+TR3 CH(IC-1,2,K) = TR2-TR3 CH(I,3,K) = TI2+TI3 CH(IC,2,K) = TI3-TI2 102 CONTINUE 103 CONTINUE RETURN END SUBROUTINE RADF4 (IDO,L1,CC,CH,WA1,WA2,WA3) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) , 1 WA1(*) ,WA2(*) ,WA3(*) DATA HSQT2 /0.D0/ IF(HSQT2.EQ.0.D0)HSQT2=SQRT(0.5D0) DO 101 K=1,L1 TR1 = CC(1,K,2)+CC(1,K,4) TR2 = CC(1,K,1)+CC(1,K,3) CH(1,1,K) = TR1+TR2 CH(IDO,4,K) = TR2-TR1 CH(IDO,2,K) = CC(1,K,1)-CC(1,K,3) CH(1,3,K) = CC(1,K,4)-CC(1,K,2) 101 CONTINUE IF (IDO-2) 107,105,102 102 IDP2 = IDO+2 DO 104 K=1,L1 DO 103 I=3,IDO,2 IC = IDP2-I CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4) CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4) TR1 = CR2+CR4 TR4 = CR4-CR2 TI1 = CI2+CI4 TI4 = CI2-CI4 TI2 = CC(I,K,1)+CI3 TI3 = CC(I,K,1)-CI3 TR2 = CC(I-1,K,1)+CR3 TR3 = CC(I-1,K,1)-CR3 CH(I-1,1,K) = TR1+TR2 CH(IC-1,4,K) = TR2-TR1 CH(I,1,K) = TI1+TI2 CH(IC,4,K) = TI1-TI2 CH(I-1,3,K) = TI4+TR3 CH(IC-1,2,K) = TR3-TI4 CH(I,3,K) = TR4+TI3 CH(IC,2,K) = TR4-TI3 103 CONTINUE 104 CONTINUE IF (MOD(IDO,2) .EQ. 1) RETURN 105 CONTINUE DO 106 K=1,L1 TI1 = -HSQT2*(CC(IDO,K,2)+CC(IDO,K,4)) TR1 = HSQT2*(CC(IDO,K,2)-CC(IDO,K,4)) CH(IDO,1,K) = TR1+CC(IDO,K,1) CH(IDO,3,K) = CC(IDO,K,1)-TR1 CH(1,2,K) = TI1-CC(IDO,K,3) CH(1,4,K) = TI1+CC(IDO,K,3) 106 CONTINUE 107 RETURN END SUBROUTINE RADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,L1,5) ,CH(IDO,5,L1) , 1 WA1(*) ,WA2(*) ,WA3(*) ,WA4(*) DATA TR11,TI11,TR12,TI12/0.D0,0.D0,0.D0,0.D0/ DATA PI/0.D0/ IF(PI.EQ.0.D0)PI=4.D0*ATAN(1.D0) IF(TR11.EQ.0.D0)TR11=COS(2.D0/5.D0*PI) IF(TI11.EQ.0.D0)TI11=SIN(2.D0/5.D0*PI) IF(TR12.EQ.0.D0)TR12=-SIN(3.D0/10.D0*PI) IF(TI12.EQ.0.D0)TI12=COS(3.D0/10.D0*PI) DO 101 K=1,L1 CR2 = CC(1,K,5)+CC(1,K,2) CI5 = CC(1,K,5)-CC(1,K,2) CR3 = CC(1,K,4)+CC(1,K,3) CI4 = CC(1,K,4)-CC(1,K,3) CH(1,1,K) = CC(1,K,1)+CR2+CR3 CH(IDO,2,K) = CC(1,K,1)+TR11*CR2+TR12*CR3 CH(1,3,K) = TI11*CI5+TI12*CI4 CH(IDO,4,K) = CC(1,K,1)+TR12*CR2+TR11*CR3 CH(1,5,K) = TI12*CI5-TI11*CI4 101 CONTINUE IF (IDO .EQ. 1) RETURN IDP2 = IDO+2 DO 103 K=1,L1 DO 102 I=3,IDO,2 IC = IDP2-I DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) DR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4) DI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4) DR5 = WA4(I-2)*CC(I-1,K,5)+WA4(I-1)*CC(I,K,5) DI5 = WA4(I-2)*CC(I,K,5)-WA4(I-1)*CC(I-1,K,5) CR2 = DR2+DR5 CI5 = DR5-DR2 CR5 = DI2-DI5 CI2 = DI2+DI5 CR3 = DR3+DR4 CI4 = DR4-DR3 CR4 = DI3-DI4 CI3 = DI3+DI4 CH(I-1,1,K) = CC(I-1,K,1)+CR2+CR3 CH(I,1,K) = CC(I,K,1)+CI2+CI3 TR2 = CC(I-1,K,1)+TR11*CR2+TR12*CR3 TI2 = CC(I,K,1)+TR11*CI2+TR12*CI3 TR3 = CC(I-1,K,1)+TR12*CR2+TR11*CR3 TI3 = CC(I,K,1)+TR12*CI2+TR11*CI3 TR5 = TI11*CR5+TI12*CR4 TI5 = TI11*CI5+TI12*CI4 TR4 = TI12*CR5-TI11*CR4 TI4 = TI12*CI5-TI11*CI4 CH(I-1,3,K) = TR2+TR5 CH(IC-1,2,K) = TR2-TR5 CH(I,3,K) = TI2+TI5 CH(IC,2,K) = TI5-TI2 CH(I-1,5,K) = TR3+TR4 CH(IC-1,4,K) = TR3-TR4 CH(I,5,K) = TI3+TI4 CH(IC,4,K) = TI4-TI3 102 CONTINUE 103 CONTINUE RETURN END SUBROUTINE RADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 1 C1(IDO,L1,IP) ,C2(IDL1,IP), 2 CH2(IDL1,IP) ,WA(*) DATA TPI/0.D0/ IF(TPI.EQ.0.D0)TPI=8.D0*ATAN(1.D0) ARG = TPI/DBLE(IP) DCP = COS(ARG) DSP = SIN(ARG) IPPH = (IP+1)/2 IPP2 = IP+2 IDP2 = IDO+2 NBD = (IDO-1)/2 IF (IDO .EQ. 1) GO TO 119 DO 101 IK=1,IDL1 CH2(IK,1) = C2(IK,1) 101 CONTINUE DO 103 J=2,IP DO 102 K=1,L1 CH(1,K,J) = C1(1,K,J) 102 CONTINUE 103 CONTINUE IF (NBD .GT. L1) GO TO 107 IS = -IDO DO 106 J=2,IP IS = IS+IDO IDIJ = IS DO 105 I=3,IDO,2 IDIJ = IDIJ+2 DO 104 K=1,L1 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J) CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J) 104 CONTINUE 105 CONTINUE 106 CONTINUE GO TO 111 107 IS = -IDO DO 110 J=2,IP IS = IS+IDO DO 109 K=1,L1 IDIJ = IS DO 108 I=3,IDO,2 IDIJ = IDIJ+2 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J) CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J) 108 CONTINUE 109 CONTINUE 110 CONTINUE 111 IF (NBD .LT. L1) GO TO 115 DO 114 J=2,IPPH JC = IPP2-J DO 113 K=1,L1 DO 112 I=3,IDO,2 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC) C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC) C1(I,K,J) = CH(I,K,J)+CH(I,K,JC) C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J) 112 CONTINUE 113 CONTINUE 114 CONTINUE GO TO 121 115 DO 118 J=2,IPPH JC = IPP2-J DO 117 I=3,IDO,2 DO 116 K=1,L1 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC) C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC) C1(I,K,J) = CH(I,K,J)+CH(I,K,JC) C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J) 116 CONTINUE 117 CONTINUE 118 CONTINUE GO TO 121 119 DO 120 IK=1,IDL1 C2(IK,1) = CH2(IK,1) 120 CONTINUE 121 DO 123 J=2,IPPH JC = IPP2-J DO 122 K=1,L1 C1(1,K,J) = CH(1,K,J)+CH(1,K,JC) C1(1,K,JC) = CH(1,K,JC)-CH(1,K,J) 122 CONTINUE 123 CONTINUE C AR1 = 1. AI1 = 0. DO 127 L=2,IPPH LC = IPP2-L AR1H = DCP*AR1-DSP*AI1 AI1 = DCP*AI1+DSP*AR1 AR1 = AR1H DO 124 IK=1,IDL1 CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2) CH2(IK,LC) = AI1*C2(IK,IP) 124 CONTINUE DC2 = AR1 DS2 = AI1 AR2 = AR1 AI2 = AI1 DO 126 J=3,IPPH JC = IPP2-J AR2H = DC2*AR2-DS2*AI2 AI2 = DC2*AI2+DS2*AR2 AR2 = AR2H DO 125 IK=1,IDL1 CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J) CH2(IK,LC) = CH2(IK,LC)+AI2*C2(IK,JC) 125 CONTINUE 126 CONTINUE 127 CONTINUE DO 129 J=2,IPPH DO 128 IK=1,IDL1 CH2(IK,1) = CH2(IK,1)+C2(IK,J) 128 CONTINUE 129 CONTINUE C IF (IDO .LT. L1) GO TO 132 DO 131 K=1,L1 DO 130 I=1,IDO CC(I,1,K) = CH(I,K,1) 130 CONTINUE 131 CONTINUE GO TO 135 132 DO 134 I=1,IDO DO 133 K=1,L1 CC(I,1,K) = CH(I,K,1) 133 CONTINUE 134 CONTINUE 135 DO 137 J=2,IPPH JC = IPP2-J J2 = J+J DO 136 K=1,L1 CC(IDO,J2-2,K) = CH(1,K,J) CC(1,J2-1,K) = CH(1,K,JC) 136 CONTINUE 137 CONTINUE IF (IDO .EQ. 1) RETURN IF (NBD .LT. L1) GO TO 141 DO 140 J=2,IPPH JC = IPP2-J J2 = J+J DO 139 K=1,L1 DO 138 I=3,IDO,2 IC = IDP2-I CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC) CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC) CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC) CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J) 138 CONTINUE 139 CONTINUE 140 CONTINUE RETURN 141 DO 144 J=2,IPPH JC = IPP2-J J2 = J+J DO 143 I=3,IDO,2 IC = IDP2-I DO 142 K=1,L1 CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC) CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC) CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC) CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J) 142 CONTINUE 143 CONTINUE 144 CONTINUE RETURN END C*********************************************************************** C C ******* END OF THE COLLECTION OF FISHPACK ROUTINES ****** C C*********************************************************************** PROGRAM DRIVER C********************************************************************* C C SIMPLE DRIVER FOR THE PACKAGE FOURCO C DOUBLE PRECISION VERSION C C THIS DRIVER TESTS ALL THE FEATURES OF THE PACKAGE FOURCO C THE OUTPUT ARE THE VALUES OF THE FIRST 100 COSINE/SINE C FOURIER COEFFICIENTS ON THE INTERVAL (0,1) OF THE FUNCTION C F(X) = EXP(3.*X) C COMPUTED FIRST USING THE TOP END FCOEFS WITH 3 DIFFERENT C NUMERICAL DIFFERENTIATION ROUTINES,THE RESTART FEATURE AND C PASSING PRECALCULATED EXACT DERIVATIVES, AND THEN USING THE C TOP END FCOETD. C C********************************************************************* C C THE DRIVER STARTS HERE C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION C(100),S(100),CEXACT(20),SEXACT(20), $ WORK(5000),XLAMB(20),FUNVAL(65) EXTERNAL FUN C C THIS IS THE VALUE OF PI C PI=4.0D0*ATAN(1.0D0) C C SET THE ENDPOINTS OF THE INTERVAL C A=0.D0 B=1.D0 C C SET THE NUMBER OF REQUESTED COEFFICIENTS C NCOEF=100 C C SET THE REQUESTED ABSOLUTE ACCURACY C TOL=1.D-10 C C SET THE PHYSICAL LIMIT PARAMETER C MTOP=32 PRINT 1 1 FORMAT(8X,'--------------------------------------------'/ $ 4X,'THIS RUN TESTS ALL THE FEATURES OF THE PACKAGE FOURCO :'/ $ 4X,'.TOP-END FCOEFS WITH 3 NUM.DIFFERENTIATION ROUTINES,'/ $ 5X,'THE RESTART FEATURE AND PRECALCULATED DERIVATIVES'/ $ 4X,'.TOP-END FCOETD FOR SAMPLED FUNCTION VALUES'/ $ 8X,'--------------------------------------------'// $ 5X,' TEST FUNCTION IS EXP(3.D0*X)'/) C C NOW CALL THE TOP-END FCOEFS USING ENDACE(ACCURATE NUMERICAL C DIFFERENTIATION ROUTINE) C IDER=1 IFIRST=0 C CALL FCOEFS(A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,C,S,IERR,ESTER,MACT) C PRINT 5,IDER,IFIRST,MTOP,TOL 5 FORMAT(1X//4X,'...........................................'/ $ 4X,'THIS IS OUTPUT FROM PACKAGE FOURCO USING TOP-END'/ $ 4X,'FCOEFS WITH ENDACE(ACCURATE NUM.DIFFERENT. ROUTINE)'/ $ 4X,'.........................................'/ $ 14X,'IDER=',I3,' IFIRST=',I3/ $ 14X,'MTOP=',I4,' TOL=',D12.2/) C PRINT THE ESTIMATED ABSOLUTE ERROR, THE ERROR PARAMETER AND C THE VALUE OF MACT C PRINT 10,ESTER,IERR,MACT 10 FORMAT(20X, $ ' ESTIMATED ERROR',5X,D10.3/ $ ' ERROR PARAMETER',10X,I2/ $ ' MACT ',8X,I4) C C PRINT THE VALUE OF THE APPROXIMATED COSINE COEFFICIENTS COMPUTED C BY FCOEFS , THE VALUES OF THE EXACT COEFFICIENTS C AND THE ABSOLUTE ERRORS C PRINT 20 20 FORMAT(1X///8X,' COSINE FOURIER COEFFICIENTS'// $ ' I-1',6X,'COMPUTED',5X,'EXACT',4X, $ 'COMPUTED-EXACT'/) DO 100 I=1,NCOEF,5 J=I/5 +1 CEXACT(J)=(EXP(3.0D0)-1.0D0)*3.0D0/(9.0D0+4.0D0*PI**2*(I-1)**2) C C REMEMBER OUR DEFINITION OF FOURIER COEFFICIENTS C IF (J.NE.1)CEXACT(J)=2.0D0*CEXACT(J) DIFFS=C(I)-CEXACT(J) PRINT 30,I,C(I),CEXACT(J),DIFFS 30 FORMAT(2X,I4,3(2X,D10.3)) 100 CONTINUE C C AND NOW THE SINE COEFFICIENTS ... C PRINT 40 40 FORMAT(1X///8X,' SINE FOURIER COEFFICIENTS'// $ ' I-1',6X,'COMPUTED',5X,'EXACT',4X, $ 'COMPUTED-EXACT'/) DO 200 I=1,NCOEF,5 J=I/5 +1 SEXACT(J)=-(EXP(3.0D0)-1.0D0)*2.0D0*PI*(I-1) $ /(9.0D0+4.0D0*PI**2*(I-1)**2) SEXACT(J)=2.0D0*SEXACT(J) DIFFS=S(I)-SEXACT(J) PRINT 50,I,S(I),SEXACT(J),DIFFS 50 FORMAT(2X,I4,3(2X,D10.3)) 200 CONTINUE C C NOW CALL THE TOP-END FCOEFS USING CHPEND(CHEAP NUMERICAL C DIFFERENTIATION ROUTINE) C IDER=0 IFIRST=0 C CALL FCOEFS(A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,C,S,IERR,ESTER,MACT) C PRINT 45,IDER,IFIRST,MTOP,TOL 45 FORMAT(1X//4X,'...........................................'/ $ 4X,'THIS IS OUTPUT FROM PACKAGE FOURCO USING TOP-END'/ $ 4X,'FCOEFS WITH CHPEND(CHEAP NUM.DIFFERENT. ROUTINE)'/ $ 4X,'.........................................'/ $ 14X,'IDER=',I3,' IFIRST=',I3/ $ 14X,'MTOP=',I4,' TOL=',D12.2/) C PRINT THE ESTIMATED ABSOLUTE ERROR, THE ERROR PARAMETER AND C THE VALUE OF MACT C PRINT 10,ESTER,IERR,MACT C C PRINT THE VALUE OF THE APPROXIMATED COSINE COEFFICIENTS COMPUTED C BY FCOEFS , THE VALUES OF THE EXACT COEFFICIENTS C AND THE ABSOLUTE ERRORS C PRINT 20 DO 101 I=1,NCOEF,5 J=I/5 +1 DIFFS=C(I)-CEXACT(J) PRINT 30,I,C(I),CEXACT(J),DIFFS 101 CONTINUE C C AND NOW THE SINE COEFFICIENTS ... C PRINT 40 DO 201 I=1,NCOEF,5 J=I/5 +1 DIFFS=S(I)-SEXACT(J) PRINT 50,I,S(I),SEXACT(J),DIFFS 201 CONTINUE C C NOW CALL THE TOP-END FCOEFS USING THE RESTART FEATURE. C PREVIOUS RESULTS WILL BE USED. C IFIRST=1 MTOP=512 C CALL FCOEFS(A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,C,S,IERR,ESTER,MACT) C PRINT 65,IDER,IFIRST,MTOP,TOL 65 FORMAT(1X//4X,'...........................................'/ $ 4X,'THIS IS OUTPUT FROM PACKAGE FOURCO USING TOP-END'/ $ 4X,'FCOEFS WITH RESTART FEATURE'/ $ 4X,'.........................................'/ $ 14X,'IDER=',I3,' IFIRST=',I3/ $ 14X,'MTOP=',I4,' TOL=',D12.2/) C PRINT THE ESTIMATED ABSOLUTE ERROR, THE ERROR PARAMETER AND C THE VALUE OF MACT C PRINT 10,ESTER,IERR,MACT C C PRINT THE VALUE OF THE APPROXIMATED COSINE COEFFICIENTS COMPUTED C BY FCOEFS , THE VALUES OF THE EXACT COEFFICIENTS C AND THE ABSOLUTE ERRORS C PRINT 20 DO 102 I=1,NCOEF,5 J=I/5 +1 DIFFS=C(I)-CEXACT(J) PRINT 30,I,C(I),CEXACT(J),DIFFS 102 CONTINUE C C AND NOW THE SINE COEFFICIENTS ... C PRINT 40 DO 202 I=1,NCOEF,5 J=I/5 +1 DIFFS=S(I)-SEXACT(J) PRINT 50,I,S(I),SEXACT(J),DIFFS 202 CONTINUE C C C NOW CALL THE TOP-END FCOEFS USING ONESID(ONE-SIDED NUMERICAL C DIFFERENTIATION ROUTINE) C IDER=2 IFIRST=0 MTOP=32 C CALL FCOEFS(A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,C,S,IERR,ESTER,MACT) C PRINT 75,IDER,IFIRST,MTOP,TOL 75 FORMAT(1X//4X,'...........................................'/ $ 4X,'THIS IS OUTPUT FROM PACKAGE FOURCO USING TOP-END'/ $ 4X,'FCOEFS WITH ONESID(ONE-SIDED NUM.DIFFERENT.ROUTINE)'/ $ 4X,'.........................................'/ $ 14X,'IDER=',I3,' IFIRST=',I3/ $ 14X,'MTOP=',I4,' TOL=',D12.2/) C PRINT THE ESTIMATED ABSOLUTE ERROR, THE ERROR PARAMETER AND C THE VALUE OF MACT C PRINT 10,ESTER,IERR,MACT C C PRINT THE VALUE OF THE APPROXIMATED COSINE COEFFICIENTS COMPUTED C BY FCOEFS , THE VALUES OF THE EXACT COEFFICIENTS C AND THE ABSOLUTE ERRORS C PRINT 20 DO 103 I=1,NCOEF,5 J=I/5 +1 DIFFS=C(I)-CEXACT(J) PRINT 30,I,C(I),CEXACT(J),DIFFS 103 CONTINUE C C AND NOW THE SINE COEFFICIENTS ... C PRINT 40 DO 203 I=1,NCOEF,5 J=I/5 +1 DIFFS=S(I)-SEXACT(J) PRINT 50,I,S(I),SEXACT(J),DIFFS 203 CONTINUE C C C NOW CALL THE TOP-END FCOEFS PASSING PRECALCULATED DERIVATIVES. C XLAMB(I)=I-TH DERIVATIVE OF FUN AT B MINUS I-TH DERIVATIVE OF C FUN AT A. C IDER=-5 DO 121 I=1,ABS(IDER) XLAMB(I) = 3**I * (EXP(3.D0*B) - EXP(3.D0*A)) 121 CONTINUE IFIRST=0 C CALL FCOEFS(A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,C,S,IERR,ESTER,MACT) C PRINT 55,ABS(IDER),IDER,IFIRST,MTOP,TOL 55 FORMAT(1X//4X,'...........................................'/ $ 4X,'THIS IS OUTPUT FROM PACKAGE FOURCO USING TOP-END'/ $ 4X,'FCOEFS PASSING ',I2,' PRECALCULATED EXACT DERIVATIVES'/ $ 4X,'.........................................'/ $ 14X,'IDER=',I3,' IFIRST=',I3/ $ 14X,'MTOP=',I4,' TOL=',D12.2/) C PRINT THE ESTIMATED ABSOLUTE ERROR, THE ERROR PARAMETER AND C THE VALUE OF MACT C PRINT 10,ESTER,IERR,MACT C C PRINT THE VALUE OF THE APPROXIMATED COSINE COEFFICIENTS COMPUTED C BY FCOEFS , THE VALUES OF THE EXACT COEFFICIENTS C AND THE ABSOLUTE ERRORS C PRINT 20 DO 104 I=1,NCOEF,5 J=I/5 +1 DIFFS=C(I)-CEXACT(J) PRINT 30,I,C(I),CEXACT(J),DIFFS 104 CONTINUE C C AND NOW THE SINE COEFFICIENTS ... C PRINT 40 DO 204 I=1,NCOEF,5 J=I/5 +1 DIFFS=S(I)-SEXACT(J) PRINT 50,I,S(I),SEXACT(J),DIFFS 204 CONTINUE C C C NOW CALL THE TOP-END FCOETD THAT REQUIRES SAMPLED FUNCTION C VALUES. C C MTOP=33 STEP=(B-A)/(MTOP-1) DO 105 I=1,MTOP FUNVAL(I)=FUN(A+(I-1)*STEP) 105 CONTINUE CALL FCOETD(A,B,NCOEF,MTOP,FUNVAL,WORK,C,S,IERR,ESTER) C PRINT 115,IDER,IFIRST,MTOP,TOL 115 FORMAT(1X//4X,'...........................................'/ $ 4X,'THIS IS OUTPUT FROM PACKAGE FOURCO USING TOP-END'/ $ 4X,'FCOETD THAT REQUIRES SAMPLED FUNCTION VALUES'/ $ 4X,'.........................................'/ $ 14X,'IDER=',I3,' IFIRST=',I3/ $ 14X,'MTOP=',I4,' TOL=',D12.2/) C PRINT THE ESTIMATED ABSOLUTE ERROR, THE ERROR PARAMETER AND C THE VALUE OF MACT C PRINT 10,ESTER,IERR,MACT C C PRINT THE VALUE OF THE APPROXIMATED COSINE COEFFICIENTS COMPUTED C BY FCOEFS , THE VALUES OF THE EXACT COEFFICIENTS C AND THE ABSOLUTE ERRORS C PRINT 20 DO 106 I=1,NCOEF,5 J=I/5 +1 DIFFS=C(I)-CEXACT(J) PRINT 30,I,C(I),CEXACT(J),DIFFS 106 CONTINUE C C AND NOW THE SINE COEFFICIENTS ... C PRINT 40 DO 206 I=1,NCOEF,5 J=I/5 +1 DIFFS=S(I)-SEXACT(J) PRINT 50,I,S(I),SEXACT(J),DIFFS 206 CONTINUE C C STOP C******************************************************************** END DOUBLE PRECISION FUNCTION FUN(X) C********************************************************************* C C THIS IS THE TEST FUNCTION C DOUBLE PRECISION X FUN=EXP(3.0D0*X) RETURN C C END OF THE TEST FUNCTION C******************************************************************** C END OF THE SIMPLE DRIVER C******************************************************************** END PROGRAM CDRIVE C*********************************************************************** C C DRIVER FOR THE PACKAGE FOURCO C DOUBLE PRECISION VERSION C C*********************************************************************** C C INCLUDES MAIN, TRYFUN, INPUT, EXACT, BERF . C MAIN CALLS INPUT AND THEN FOURCO. C C THIS DRIVER HAS BEEN DESIGNED FOR A STRINGENT TEST OF THE C PACKAGE FOURCO. C C THE VARIOUS OPTIONS ARE DESCRIBED BELOW AND IN THE SUBROUTINE C INPUT. C C*********************************************************************** C C THE LIST OF TEST FUNCTIONS IS IN THE SUBROUTINE INPUT. C THE VALUE OF NFUN IN THE SUBROUTINE INPUT SELECTS A TESTFUNCTION. C THE ARRAY AVEC CONTAINS THE PARAMATERS AMP,ALPHA AND BETA C WHICH MODIFY THE BEHAVIOUR OF THE TEST FUNCTION (SEE SUBROUTINE C INPUT). C C THE VALUE OF THE REQUESTED INPUT ACCURACY IS ALSO SET IN C THE SUBROUTINE INPUT (RELTOL). C C THE VALUE OF THE ARRAY NIX IN THE SUBROUTINE INPUT ALLOWS TO C SELECT OUTPUT TEST PRINTS. C SET NIX(I)=0 I=1,9 IF ONLY DEFAULT PRINTS ARE REQUESTED. C C IF NIX(10)=0 THEN THE FUNCTION SELECTED IN THE SUBROUTINE INPUT C IS USED C IF NIX(10)=1 THEN FOURCO IS TESTED USING 8 TEST FUNCTIONS C IN THE SUBROUTINE INPUT C C THE VALUE OF NIX(11) SELECTS A RUN OF FCOEFS (NIX(11)=0) C OR A RUN OF FCOETD (NIX(11)=1). C THE VALUE OF NIX(12) IS THE VALUE OF THE INPUT ARGUMENT OF C FCOEFS IDER, WHICH SELECTS THE DIFFERENTIATION ROUTINE. C*********************************************************************** C*********************************************************************** C C MAIN C C*********************************************************************** C IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /M360/AVEC(10),AAAA,BBBB,A,B,RELTOL,TOLFAC,SFAC,AVINC, + NIX(12),NVEC(20),NFUN,N2TOP,MODSEQ,MTOP C EXTERNAL TRYFUN * DIMENSION FCOS(240),FSIN(240) DIMENSION WORK(5000),XLAMB(20),FUNVAL(1025) DIMENSION CCOS(240),CSIN(240) C C NFUN = 0 CALL INPUT C 1 CONTINUE STEP = 1.0D0 - AVEC(4) TOL = RELTOL NCOEF = 240 C N2 = 1 10 CONTINUE C C ********** LOOP STARTS HERE. RE-ENTRY FROM 900 BELOW *********** C C SECTION TO SET TOLERANCE C IF (MODSEQ.EQ.1) GO TO 30 IF (MODSEQ.EQ.-1) GO TO 30 BIG = 0.0D0 AA = A DIFF = (B-A)/4.0D0 C DO 20 J = 1,5 TEMP = TRYFUN(AA) BIG = MAX(ABS(TEMP),BIG) AA = AA + DIFF 20 CONTINUE TOL = BIG*RELTOL 30 CONTINUE C 40 CONTINUE C NVEC(2) = N2 C PRINT 9000,AVEC(1),AVEC(2),AVEC(3),AVEC(4),NFUN,NVEC(2) * 9000 FORMAT (/,' FUNCTION SPECIFICATION. AVEC, NVEC ',//,4D16.8,2I4,//) C NVEC(10) = 0 C NVEC(10) IS USED TO COUNT NUMBER OF CALLS TO TRYFUN C IF (NIX(11).EQ.1) GO TO 50 PRINT 9010 * 9010 FORMAT (' WELCOME TO FCOEFS. STANDARD TOP END'// + ' INPUT ARGUMENTS ARE NCOEF,TOL,A,B,MTOP,IDER,IFIRST') PRINT 9015,NCOEF,TOL,A,B,MTOP,NIX(12),IFIRST 9015 FORMAT(1X,I4,3D16.8,1X,2I4,I2) GO TO 60 * 50 CONTINUE PRINT 9020 MTOP=257 STEP=(B-A)/(MTOP-1) DO 85 IJ=1,MTOP FUNVAL(IJ)=TRYFUN(A+(IJ-1)*STEP) 85 CONTINUE * 9020 FORMAT (' WELCOME TO FCOETD. TOP END FOR SAMPLED FUNCTION'// + ' INPUT ARGUMENTS ARE NCOEF,A,B,MTOP') PRINT 9030,NCOEF,A,B,MTOP * 9030 FORMAT (1X,I4,2D16.8,I4,/) C 60 CONTINUE IF (NIX(11).EQ.1) GO TO 80 70 CONTINUE IDER=NIX(12) IFIRST=0 CALL FCOEFS(A,B,TRYFUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,FCOS,FSIN,IERR,ESTER,MACT) GO TO 90 * 80 CONTINUE C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CALL FCOETD(A,B,NCOEF,MTOP,FUNVAL,WORK,C,S,IERR,ESTER) MACT=MTOP 90 CONTINUE C PRINT 9040,IERR,NVEC(10),MACT * 9040 FORMAT (/' IERR=',I2,2X,'GLOBAL N.OF FUN.EVAL=',I5,' MACT= ',I5/) C IF (NFUN.NE.1) GO TO 100 PRINT 9050 * 9050 FORMAT (' THE FUNCTION IS K*EXP(ALPHA(X-BETA)) ') PRINT 9060,AVEC(3),AVEC(1),AVEC(2) * 9060 FORMAT (/,'K=',D16.8,2X,'ALPHA=',D16.8,2X,'BETA=',D16.8,/) GO TO 120 C 100 CONTINUE IF (NFUN.NE.2) GO TO 110 PRINT 9070,N2 * 9070 FORMAT (' THE FUNCTION IS THE BERNOULLI POLYNOMIAL OF DEGREE',I3, + /) GO TO 120 C 110 CONTINUE PRINT 9080,NFUN,AVEC(4) * 9080 FORMAT (' THE FUNCTION IS NFUN',I2,' RHO =',D16.8/) C 120 CONTINUE C PRINT 9090,NFUN * 9090 FORMAT (' NFUN = ',I1,/) PRINT 9100,ESTER * 9100 FORMAT ('********** ESTER =',D16.8,' **********',/) C IF (NIX(1).GT.0) PRINT 9110 * 9110 FORMAT (1X//' //// ACCURACY CONTROL SUSPENDED ////'//) IF (NIX(2).EQ.1) PRINT 9120 * 9120 FORMAT (1X,' **** SINGLE PRECISION FUNCTION VALUES ****') C C CALL EXACT(CCOS,CSIN) C PRINT 9130,NFUN,N2 * 9130 FORMAT (//,' **** THE COSINE COEFFICIENTS **** NFUN,N2 =',2I3,/) PRINT 9140 * 9140 FORMAT (16X,'********** DOUBLE PRECISION **********',/) PRINT 9150 * 9150 FORMAT (2X,'J-1',4X,'J',5X,'CALCULATED',8X,'EXACT',9X, + 'DIFFERENCE') C J = 1 DO 140 JJJ = 1,30 JJ = J - 1 DIFF = FCOS(J) - CCOS(J) NOP = 1 IF (ABS(DIFF).GT.ABS(ESTER)) NOP = -1 PRINT 9160,JJ,J,FCOS(J),CCOS(J),DIFF,NOP GO TO 130 * 9160 FORMAT (2I5,3D16.8,2X,I2) C 130 CONTINUE J = J + 1 IF (J.GE.12) J = J + 9 IF (J.GE.55) J = J + 15 IF (J.GT.240) GO TO 150 140 CONTINUE 150 CONTINUE C NAB = 10 DO 160 I = 1,240 IF (NAB.EQ.0) GO TO 170 DIF = FCOS(I) - CCOS(I) IF (ABS(DIF).LE.ABS(ESTER)) GO TO 160 NAB = NAB - 1 PRINT 9170,I * 9170 FORMAT ('ERROR EXCESSIVE IN COEFFICIENT ',I3) 160 CONTINUE 170 CONTINUE C PRINT 9180,NFUN,N2 * 9180 FORMAT (//,' ***** THE SINE COEFFICIENTS **** NFUN,N2 =',2I3,/) PRINT 9140 PRINT 9150 J = 1 DO 190 JJJ = 1,30 JJ = J - 1 DIFF = FSIN(J) - CSIN(J) NOP = 1 IF (ABS(DIFF).GT.ABS(ESTER)) NOP = -1 PRINT 9160,JJ,J,FSIN(J),CSIN(J),DIFF,NOP GO TO 180 * 180 CONTINUE J = J + 1 IF (J.GE.12) J = J + 9 IF (J.GE.55) J = J + 15 IF (J.GT.240) GO TO 200 190 CONTINUE 200 CONTINUE C NAB = 10 DO 210 I = 1,240 IF (NAB.EQ.0) GO TO 220 DIF = FSIN(I) - CSIN(I) IF (ABS(DIF).LE.ABS(ESTER)) GO TO 210 NAB = NAB - 1 PRINT 9170,I 210 CONTINUE 220 CONTINUE C C THE FUNCTION UPDATE. C ALTER EITHER FUNCTION OR TOLERANCE C N2 = N2 + 1 IF (N2.GT.N2TOP) GO TO 240 IF (MODSEQ.GT.0) GO TO 230 TOL = TOL*TOLFAC GO TO 40 C 230 CONTINUE IF (NFUN.LE.3) AVEC(1) = AVEC(1) + AVINC STEP = SFAC*STEP AVEC(4) = 1.0D0 - STEP GO TO 10 C C 240 CONTINUE C IF(NIX(10).EQ.0) GO TO 250 IF(NFUN.GE.9) GO TO 250 CALL INPUT GOTO 1 C 250 CONTINUE C STOP C C END OF MAIN C*********************************************************************** END DOUBLE PRECISION FUNCTION TRYFUN(X) C*********************************************************************** C C FUNCTION TRYFUN C C*********************************************************************** C*********************************************************************** C C THIS FUNCTION CONTAINS ALL THE TEST FUNCTIONS. C THE VALUE OF NFUN (IN COMMON),WHICH IS SET IN THE C SUBROUTINE INPUT,SELECTS THE DESIRED FUNCTION C C*********************************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Z) REAL SINFUN C TEMPCOM COMMON LIST COMMON /M360/AVEC(10),AAAA,BBBB,A,B,RELTOL,TOLFAC,SFAC,AVINC, + NIX(12),NVEC(20),NFUN,N2TOP,MODSEQ,MTOP C DIMENSION BERFUN(12),BFFACT(12) * TWOPI = 8.0D0*ATAN(1.0D0) ALPHA = AVEC(1) ALSQ = ALPHA*ALPHA BETA = AVEC(2) ARG = ALPHA* (X-BETA) TWOPIX = TWOPI*ARG RHO = AVEC(4) RHOSQ = RHO*RHO NVEC(10) = NVEC(10) + 1 N2 = NVEC(2) C IF (NFUN.LE.0) GO TO 100 10 CONTINUE IF (NFUN.NE.1) GO TO 20 ARG = ALPHA* (X-BETA) TRYFUN = EXP(ARG) GO TO 90 * 20 CONTINUE IF (NFUN.NE.2) GO TO 30 IF (N2.LT.0) GO TO 100 IF (N2.GT.12) GO TO 100 IF (N2.EQ.0) TRYFUN = 1.0D0 IF (N2.EQ.0) GO TO 90 XX = ARG NPER = 0 CALL BERF(XX,NPER,BERZ,BERFUN,BFFACT) TRYFUN = BERFUN(N2) GO TO 90 * 30 CONTINUE IF (NFUN.NE.4) GO TO 40 XNUM = 1.0D0 - RHOSQ DEN = 1.0D0 - 2.0D0*COS(TWOPIX)*RHO + RHOSQ TRYFUN = XNUM/DEN GO TO 90 * 40 CONTINUE IF (NFUN.NE.5) GO TO 50 XNUM = RHO* (SIN(TWOPIX)) DEN = 1.0D0 - RHO* (COS(TWOPIX)) TRYFUN = 2.0D0*ATAN2(XNUM,DEN) GO TO 90 * 50 CONTINUE IF (NFUN.NE.6) GO TO 60 ARG = 1.0D0 - 2.0D0*RHO*COS(TWOPIX) + RHOSQ TERM = 2.0D0*LOG(ARG) TRYFUN = TERM/2.0D0 GO TO 90 * 60 CONTINUE IF (NFUN.NE.7) GO TO 70 XNUM = 2.0D0*RHO*SIN(TWOPIX) DEN = 1.0D0 - RHOSQ TRYFUN = ATAN2(XNUM,DEN) GO TO 90 * 70 CONTINUE IF (NFUN.NE.8) GO TO 80 ARG = RHO*COS(TWOPIX) TERM1 = EXP(ARG) ARG2 = RHO*SIN(TWOPIX) TERM2 = COS(ARG2) TRYFUN = TERM1*TERM2 GO TO 90 * 80 CONTINUE TERM = (1.0D0+RHO*SIN(TWOPIX)) TRYFUN = 1.0D0/TERM 90 CONTINUE TRYFUN = AVEC(3)*TRYFUN IF (NIX(4).EQ.0) GO TO 110 PRINT 9000,NIX(4) * 9000 FORMAT (/,'OUTPUT FROM TRYFUN. NIX(4) =',I2) PRINT 9010,NVEC(10),NFUN,N2,X,TRYFUN * 9010 FORMAT (/,'TRYFUN',3I5,3X,'X =',D12.4,3X,'FUN =',D12.4) NIX(4) = NIX(4) - 1 GO TO 110 C 100 CONTINUE PRINT 9020,NFUN,N2 * 9020 FORMAT (//,' FUNCTION NOT DEFINED. NFUN AND(2) ARE',2I3,//) C 110 CONTINUE IF (NIX(2).NE.1) GO TO 120 SINFUN = TRYFUN TRYFUN = DBLE(SINFUN) 120 CONTINUE C C RETURN C C END OF TRYFUN C*********************************************************************** END SUBROUTINE INPUT C*********************************************************************** C C SUBROUTINE INPUT C C*********************************************************************** C C THIS IS A FREE FORM DATA INPUT SUBROUTINE. IT TAKES THE PLACE C OF A READ STATEMENT. THE DATA IS TRANSMITTED THROUGH COMMON. C ALL THE INPUT ARGUMENTS FOR FCOEFS AND FCOETD ARE SELECTED C IN THIS ROUTINE. C A TEST SESSION SHOULD CONSIST IN CHANGING THOSE INPUT VALUES C (SEE COMMENTS BELOW) AND RUN THE PACKAGE C C*********************************************************************** C IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /M360/AVEC(10),AAAA,BBBB,A,B,RELTOL,TOLFAC,SFAC,AVINC, + NIX(12),NVEC(20),NFUN,N2TOP,MODSEQ,MTOP C C*********************************************************************** C SET THE PARAMETERS NIX C C EACH NIX(J) J LT 10 CONTROLS THE PRINTING OF CODE CHECK OUTPUT IN C SOME PART OF THE PROGRAM. IN GENERAL THE EFFECT IS THIS C IF NIX(J) .LT.0 THERE IS ALWAYS PRINT OUT C IF NIX(J) .EQ.0 THERE IS NEVER PRINT OUT C IF NIX(J) = N .GT.0 THERE IS PRINT OUT THE FIRST N TIMES THE C FIRST N TIMES THE PRINT STATEMENT IS ENCOUNTERED. C C NIX(12) IS THE VALUE OF IDER THE INPUT ARGUMENT OF FCOEFS C WHICH SELECTS THE DIFFERENTIATION ROUTINE USED BY C THE PROCEDURE. C NIX(11) IF EQ 0 CHOOSE STANDARD TOP END. IF EQ 1 NON STANDARD. C NIX(10) IF EQ 0 NFUN SELECTED IN INPUT IS USED C IF EQ 1 A LOOP ON 8 TEST FUNCTIONS IS EXECUTED(ONLY C NFUN=3 IS NOT USED) C NIX(9) NOT USED C NIX(8) NOT USED C NIX(7) NOT USED C NIX(6) NOT USED C NIX(5) IN SPEC. PRINTS PARAMETERS C NIX(4) IN TRYFUN. PRINTS XVALUE, FUNCTION VALUE C NIX(3) NOT USED C NIX(2) IF EQ 1 TRYFUN RETURNS SINGLE PRECISION VALUES C NIX(1) NOT USED C C*********************************************************************** C NIX(12)= 1 NIX(11)= 0 NIX(10)= 1 NIX(9) = 0 NIX(8) = 0 NIX(7) = 0 NIX(6) = 0 NIX(5) = 0 NIX(4) = 0 NIX(3) = 0 NIX(2) = 0 NIX(1) = 0 C C SPECIFY THE INTERVAL C A = 0.0D0 B = 1.0D0 IF (B.LE.A) PRINT 9000 * 9000 FORMAT ('***** WARNING B MUST BE GREATER THAN A *****') C C C*********************************************************************** C THE USER MUST SPECIFY A VALUE OF NFUN FROM THE C FOLLOWING POSSIBLE CASES. ENTER YOUR VALUE BELOW. C C NFUN = 1, A SET OF MORE PRONOUNCED EXPONENTIALS C NFUN = 2, A SET OF BERNOULLI POLYNOMIALS C NFUN = 3, A SET OF STRICTER TOLERANCES USING ONE EXP.FUN C NFUN = 4, F(X)=(1-P**2)/(1-2*P*COS(2*PI*X)+P**2) C NFUN = 5, F(X)= 2*ARCTAN((PSIN(2*PI*X))/(1-PCOS(2*PI*X))) C NFUN = 6, F(X)= LN(1-2PCOS(2*PI*X)+P**2) C NFUN = 7, F(X)=ARCTAN((2PSIN(2*PI*X))/(1-P**2) C NFUN = 8, F(X)=EXP(RHOCOS(2*PI*X))*COS(RHOSIN(2*PI*X)) C NFUN = 9, F(X) = 1 / (1+PSIN(2*PI*X)) C C*********************************************************************** IF(NIX(10).EQ.0)THEN C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C CHANGE NEXT STATEMENT IF YOU WANT TO CHANGE THE TEST FUNCTION C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NFUN = 1 ELSE NFUN = NFUN+1 IF(NFUN.EQ.3) NFUN = 4 ENDIF C*********************************************************************** C NOW SET THE PARAMETERS FOR THE FIRST FUNCTION OF THIS SET. C THESE ARE AVEC(J), J=1,2,3,4. C C THE FUNCTION IS AMP * F(ALPHA X - BETA') C C AVEC(1) REPRESENTS ALPHA C AVEC(2) REPRESENTS BETA C AVEC(3) REPRESENTS THE OVERALL MULTIPLIER AMP. C AVEC(4) REPRESENTS RHO. WHEN SFAC = 0.5 AND AVEC(4) = 0.2 WE GET C RHO =.2,.6,.8,.9,..... C C*********************************************************************** AVEC(1) = 1.0D0 AVEC(2) = 0.125D0 AVEC(2) = 1.0D0 AVEC(2) = 0.0D0 AVEC(3) = 1.0D0 AVEC(4) = 0.2D0 IF (NFUN.GE.4) AVEC(1) = 1.0D0/ (B-A) C C THESE NEXT PARAMETERS GEAR THE NEXT FUNCTIONS C N2TOP = 2 AVINC = 3.0D0 IF (NFUN.NE.1) AVINC = 0.0D0 SFAC = 0.5D0 C C SPECIFY MODSEQ. POSITIVE FOR VARYING FUNCTION. 1 FOR ABS TOLERANCE. C NEGATIVE FOR VARYING TOLERANCE. 2 FOR REL TOLERANCE C MODSEQ = 2 IF (NFUN.EQ.3) MODSEQ = -1 IF (NFUN.EQ.2) MODSEQ = IABS(MODSEQ) C C RELTOL IS THE INITIAL TOLERANCE. (ABSOLUTE WHEN MODSEQ = +1 OR -1.) C (RELATIVE TO FUNAVE WHEN MODSEQ = +2 OR -2.) FUNAVE IS THE C AVERAGE OF FIVE EQUALLY SPACED FUNCTION VALUES. C RELTOL = 1.0D-16 TOLFAC = SQRT(0.1D0) IF (MODSEQ.LT.0) RELTOL = 1.0D-3 C MTOP = 512 C C THERE IS NO NFUN=3 AND NFUN = 2 ARE BERNOULLI POLYNOMIALS.HENCE C IF (NFUN.EQ.2) A = 0.0D0 IF (NFUN.EQ.2) B = 1.0D0 C IF (NFUN.EQ.3) THEN NFUN = 1 RELTOL = 1.0D-3 AVEC(1) = 16.0D0 AVINC = 0.0D0 END IF C PRINT 9010,NFUN * 9010 FORMAT (' THIS IS OUTPUT FROM A RUN OF DRIVER WITH NFUN = ',I3, + //) PRINT 9020,A,B * 9020 FORMAT (' INTERVAL A = ',D8.2,9X,'B = ',D8.2,/) IF (NIX(1).GT.0) PRINT 9030 * 9030 FORMAT (1X//' //// ACCURACY CONTROL SUSPENDED ////'//) IF (NIX(2).EQ.1) PRINT 9040 * 9040 FORMAT (1X,'**** SINGLE PRECISION FUNCTION VALUES ****') PRINT 9050 * 9050 FORMAT (//' THE FIRST FUNCTION TREATED HAS PARAMETERS AS FOLLOWS' + ) PRINT 9060 * 9060 FORMAT (/,' NFUN N2 AVEC(1)=ALPHA',7X,'AVEC(2)=BETA',9X, + 'AVEC(3)=AMP',9X,'AVEC(4)=RHO') PRINT 9070,NFUN, (AVEC(KKK),KKK=1,4) * 9070 FORMAT (I4,' 1 ',D15.7,3D20.7,/) C PRINT 9080,MTOP,SFAC,RELTOL * 9080 FORMAT (' MTOP = ',I3,9X,'SFAC = ',D8.2,5X,'RELTOL = ',D8.2,/) PRINT 9090,N2TOP,AVINC,TOLFAC * 9090 FORMAT (' N2TOP = ',I3,10X,'AVINC = ',D8.2,5X,'TOLFAC = ',D8.2,/) PRINT 9100,MODSEQ * 9100 FORMAT (' MODSEQ =',I3,4X, + 'POSITIVE MEANS FUNCTION VARIES. NEGATIVE',1X, + 'MEANS TOLERANCE VARIES',//) C IF (NIX(5).NE.0) THEN PRINT 9110 * 9110 FORMAT (' CODE-CHECK INPUT PARAMETERS',/,/) PRINT 9120,NIX(10),NIX(9),NIX(8) * 9120 FORMAT (' NIX(10) = ',I3,8X,'NIX(9) = ',I3,8X,'NIX(8) = ',I3,/) PRINT 9130,NIX(6),NIX(4),NIX(3) * 9130 FORMAT (' NIX(6) = ',I3,9X,'NIX(4) = ',I3,8X,'NIX(3) = ',I3,/) END IF C RETURN C C END OF INPUT C************************************************************** END SUBROUTINE EXACT(CCOS,CSIN) C************************************************************** C C SUBROUTINE EXACT C C************************************************************** C C EXACT RETURNS THE EXACT VALUES OF THE FOURIER COEFFICIENTS C OF THE TEST FUNCTION SELECTED IN THE SUBROUTINE INPUT. C THIS SUBROUTINE CONTAINS THE ANALYTIC FORM OF THE FOURIER C COEFFICIENTS OF ALL THE TEST FUNCTIONS PROVIDED BY THE FUNCTION C TRYFUN AND MODIFIED BY MEANS OF PARAMETERS IN THE ARRAY AVEC C C**************************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /M360/AVEC(10),AAAA,BBBB,A,B,RELTOL,TOLFAC,SFAC,AVINC, + NIX(12),NVEC(20),NFUN,N2TOP,MODSEQ,MTOP C DIMENSION CCOS(240),CSIN(240) C CALCULATES THE FIRST 240 FOURIER COEFFICIENTS OF TRYFUN USING C AN ANALYTIC EXPRESSION. C TWOPI = 8.0D0*ATAN(1.0D0) ALPHA = AVEC(1) ALSQ = ALPHA*ALPHA BETA = AVEC(2) CON = 1.0D0 RHO = AVEC(4) RHOSQ = RHO*RHO FACTOR = ALPHA* (A-BETA) C DO 10 J = 1,240 CCOS(J) = 0.0D0 CSIN(J) = 0.0D0 10 CONTINUE C N2 = NVEC(2) IF (NFUN.LE.0) GO TO 160 C IF (NFUN.NE.1) GO TO 30 C BMA = B - A CONS = EXP(-ALPHA*BETA) TERM = EXP(ALPHA*B) - EXP(ALPHA*A) CCOS(1) = (TERM*CONS)/ (ALPHA*BMA) FACT = 2.0D0*TERM*CONS/BMA TPR = TWOPI DO 20 JJ = 1,239 JJJ = JJ + 1 RAT = FACT/ ((TPR*TPR)/ (BMA*BMA)+ALSQ) CCOS(JJJ) = RAT*ALPHA CSIN(JJJ) = -RAT*TPR/BMA TPR = TPR + TWOPI 20 CONTINUE GO TO 170 C 30 CONTINUE IF (NFUN.NE.2) GO TO 60 C C NFUN = 2 FUNCTION IS BERNOULLI POLYNOMIAL OF DEGREE N2. C IF (N2.LT.0) GO TO 160 IF (N2.EQ.0) CCOS(1) = 1.0D0 IF (N2.EQ.0) GO TO 170 C FACTN2 = 1.0D0 XJ = 1.0D0 DO 40 J = 1,N2 FACTN2 = FACTN2*XJ XJ = XJ + 1.0D0 40 CONTINUE C NREM = N2 - 4* (N2/4) FACT = 2.0D0*FACTN2/ (TWOPI**N2) XJJ = 1.0D0 DO 50 J = 2,240 TEMP = XJJ** (-N2) RES = FACT*TEMP IF (NREM.EQ.0) CCOS(J) = -RES IF (NREM.EQ.2) CCOS(J) = RES IF (NREM.EQ.1) CSIN(J) = -RES IF (NREM.EQ.3) CSIN(J) = RES XJJ = XJJ + 1.0D0 50 CONTINUE C GO TO 170 C 60 CONTINUE IF (NFUN.EQ.9) GO TO 140 TERM = 1.0D0 DENOM = 1.0D0 IF (NFUN.EQ.4) CCOS(1) = 1.0D0 IF (NFUN.EQ.8) CCOS(1) = 1.0D0 DO 130 J = 1,239 K = J + 1 XJ = J ARG = TWOPI*XJ*FACTOR CTERM = COS(ARG) STERM = SIN(ARG) TERM = TERM*RHO COF = CTERM*TERM SOF = STERM*TERM 70 CONTINUE IF (NFUN.EQ.4) THEN CCOS(K) = 2.0D0*COF CSIN(K) = -2.0D0*SOF END IF * 80 CONTINUE IF (NFUN.EQ.5) THEN CCOS(K) = 2.0D0*SOF/XJ CSIN(K) = 2.0D0*COF/XJ END IF * 90 CONTINUE IF (NFUN.EQ.6) THEN CCOS(K) = -2.0D0*COF/XJ CSIN(K) = 2.0D0*SOF/XJ END IF * 100 CONTINUE IF (NFUN.EQ.7) THEN NREM = J - 2* (J/2) IF (NREM.EQ.0) GO TO 130 CCOS(K) = 2.0D0*SOF/XJ CSIN(K) = 2.0D0*COF/XJ END IF * 110 CONTINUE IF (NFUN.EQ.8) THEN DENOM = DENOM*XJ CCOS(K) = COF/DENOM CSIN(K) = -SOF/DENOM IF (DENOM.GT.1.0D30) THEN DO 120 L = K,240 CCOS(L) = 0.0D0 120 CONTINUE GO TO 170 * END IF * END IF * 130 CONTINUE C GO TO 170 C 140 CONTINUE IF (NFUN.NE.9) GO TO 160 TERM = 1.0D0 - RHOSQ ADB = SQRT(TERM)/TERM CCOS(1) = ADB CPAR = (SQRT(TERM)-1.0D0)/RHO CSQR = CPAR*CPAR CP = -ADB*CSQR SP = ADB*CPAR DO 150 J = 2,238,2 XJ = J XJM1 = J - 1 ARGE = TWOPI*XJM1*FACTOR ARGO = TWOPI*XJ*FACTOR CCOS(J) = 2.0D0*SP*SIN(ARGE) CSIN(J) = 2.0D0*SP*COS(ARGE) CCOS(J+1) = 2.0D0*CP*COS(ARGO) CSIN(J+1) = -2.0D0*CP*SIN(ARGO) CP = -CP*CSQR SP = -SP*CSQR 150 CONTINUE C GO TO 170 C 160 CONTINUE PRINT 9000,NFUN,N2 * 9000 FORMAT (//,' FUNCTION NOT DEFINED. NFUN AND(2) ARE',2I3,//) C 170 CONTINUE C AMPLIT = AVEC(3) DO 180 J = 1,240 CCOS(J) = AMPLIT*CCOS(J) CSIN(J) = AMPLIT*CSIN(J) 180 CONTINUE C RETURN C ********** END OF EXACT ********** C*********************************************************************** END SUBROUTINE BERF(XX,NPER,BERZ,BERFUN,BFFACT) C*********************************************************************** C C SUBROUTINE BERF C C*********************************************************************** C C BERF COMPUTES BERNOULLI FUNCTIONS OR POLYNOMIAL UP TO THE DEGREE C 12 C C XX ARGUMENT OF BERNOULLI FUNCTION OR POLYNOMIAL C NPER EQ.0 GIVES THE POLYNOMIAL (NON-PERIODIC). C NE.0 GIVES THE (PERIODIC) FUNCTION, WHICH COINCIDES WITH T C POLYNOMIAL IN THE OPEN INTERVAL (0,1). C BERZ = 1.0D0 C BERFUN(12) THE VALUES OF THE BERNOULLI POLYNOMIALS (OR FUNCTIONS) C OF DEGREES 1 TO 12 INCLUSIVE. C BFFACT(12) BFFACT(J) = BERFUN(J) / FACTORIAL(J) C*********************************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION BERFUN(12),BFFACT(12) * X = XX BERZ = 1.0D0 BERFUN(1) = X - 0.5D0 IF (NPER.EQ.0) GO TO 70 C C WE WANT PERIODIC FUNCTIONS 10 CONTINUE IF (X) 20,50,30 20 CONTINUE X = X + 1.0D0 GO TO 10 * 30 CONTINUE IF (X-1.0D0) 60,50,40 40 CONTINUE X = X - 1.0D0 GO TO 30 * 50 CONTINUE C HERE WE HAVE X = 0.0 OR 1.0 BERFUN(1) = 0.0D0 60 CONTINUE C 70 CONTINUE Y = X* (X-1.0D0) Z = Y* (X-0.5D0) BERFUN(2) = Y + (1.0D0/6.0D0) BERFUN(3) = Z BERFUN(4) = Y*Y - (1.0D0/30.0D0) BERFUN(5) = (Y- (1.0D0/3.0D0))*Z BERFUN(6) = Y*Y*Y - 0.5D0*Y*Y + (1.0D0/42.0D0) X2 = X*X X4 = X2*X2 X6 = X4*X2 X8 = X4*X4 X10 = X4*X6 X12 = X4*X8 Y2 = (X-0.5D0)* (X-0.5D0) Y4 = Y2*Y2 BERFUN(7) = (Y4-1.5D0*Y2+ (31.0D0/48.0D0))*Z TERM = (2.0D0*X2-7.0D0*X4+14.0D0*X6)/3.0D0 BERFUN(8) = X8 - 4.0D0*X*X6 + TERM - (1.0D0/30.0D0) TERM = X8 - 4.5D0*X*X6 + 6.0D0*X6 - 4.2D0*X4 + 2.0D0*X2 - 0.3D0 BERFUN(9) = X*TERM TERM = X10 - 5.0D0*X*X8 + 7.5D0*X8 - 7.0D0*X6 + 5.0D0*X4 BERFUN(10) = TERM - 1.5D0*X2 + (5.0D0/66.0D0) TERM = X10 - 5.5D0*X*X8 - 11.0D0* (X6-X4) - 5.5D0*X2 BERFUN(11) = (TERM+ ((5.0D0+55.0D0*X8)/6.0D0))*X TERM = X12 - 6.0D0*X*X10 + 11.0D0*X10 - 16.5D0* (X8+X4) BERFUN(12) = TERM + 22.0D0*X6 + 5.0D0*X2 - (691.0D0/2730.0D0) FACT = 1.0D0 FICT = 1.0D0 DO 80 J = 1,12 BFFACT(J) = BERFUN(J)/FACT FICT = FICT + 1.0D0 FACT = FACT*FICT 80 CONTINUE RETURN C C END OF BERF C*********************************************************************** C C ********* END OF DRIVER ********* C C*********************************************************************** END C C*********************************************************************** C* * C* ******** PACKAGE FOURCO ******** * C* * C* * C* ..... PACKAGE FOR COMPUTING SINE AND COSINE FOURIER ..... * C* ..... COEFFICIENTS OF A FUNCTION OVER AN INTERVAL ..... * C* * C* SINGLE PRECISION VERSION - SEPTEMBER 1986 * C* * C* * C***** AUTHORS : GIULIO GIUNTA AND ALMERICO MURLI * C* DIPARTIMENTO DI MATEMATICA E APPLICAZIONI * C* UNIVERSITA' DI NAPOLI * C* NAPLES - ITALY * C* * C*********************************************************************** C C***** GENERAL PURPOSE - THIS PACKAGE EVALUATES THE FIRST ..NCOEF.. C TRIGONOMETRICAL FOURIER COEFFICIENTS C C C(J)=2/(B-A)*INTEGRAL BETWEEN A,B(FUN(X)* C COS(2*PI*(J-1)*(X-A)/(B-A))) C C S(J)=2/(B-A)*INTEGRAL BETWEEN A,B(FUN(X)* C SIN(2*PI*(J-1)*(X-A)/(B-A))) C C J=1,2,....,NCOEF C C (HOPEFULLY) TO A UNIFORM ABSOLUTE ACCURACY ..TOL.. C C THE PACKAGE HAS TWO TOP-ENDS (FCOEFS,FCOETD) THAT C THE USER CAN SELECT DEPENDING ON HIS INPUT PROBLEM. C C --FCOEFS. TO USE THIS TOP-END THE USER MUST PROVIDE A C FUNCTION SUBPROGRAM FOR EVALUATING ..FUN.. AT ANY C POINT OF THE INTEGRATION INTERVAL (A,B). C C --FCOETD. TO USE THIS TOP-END THE USER MUST PROVIDE C SAMPLED VALUES OF FUN AT EQUALLY SPACED POINTS OF C THE INTEGRATION INTERVAL (A,B). C C ALL THAT THE USER MUST KNOW IN ORDER TO USE THE C PACKAGE IS THE CALLING SEQUENCE OF THE TOP-END C SUITABLE FOR HIS PROBLEM C C C***** STRUCTURE AND PORTABILITY - THE PACKAGE COMPRISES A BASIC C ROUTINE FCOEBA,TWO TOP-ENDS,A SET OF FFT ROUTINES FOR C REAL DATA,A SET OF ROUTINES FOR NUMERICAL DIFFEREN- C TIATION AND TWO AUXILIARY ROUTINES. C ALL IS WRITTEN IN FORTRAN 77. ALL MACHINE DEPENDENT C CONSTANTS ARE SPECIFIED BY R1MACH FORTRAN FUNCTION C USED IN PORT LIBRARY. ALL THE CALLS TO INTRINSIC C FORTRAN FUNCTIONS USE GENERIC NAMES. MATHEMATICAL C CONSTANTS ARE COMPUTED TO FULL MACHINE ACCURACY IN C THE CODE. C C***** FORTRAN SUPPLIED FUNCTIONS - REAL,SIN,COS,ABS,MIN,MAX,SQRT,ATAN C MOD. C***** OTHER REQUIRED FUNCTIONS - R1MACH C C C********************************************************************** C********************************************************************** C STANDARD TOP END FCOEFS * C********************************************************************** SUBROUTINE FCOEFS(A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,C,S,IERR,ESTER,MACT) C C********************************************************************** C C TOP-END FCOEFS ALLOWS TO COMPUTE TRIGONOMETRICAL FOURIER COEFFI- C CIENTS WHEN THE INTEGRAND FUNCTION ..FUN.. IS PROVIDED AS A C FUNCTION SUBPROGRAM C C C INPUT PARAMETERS - A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK,XLAMB. C C A,B -REAL- LOWER AND UPPER LIMIT OF C INTEGRATION. C C FUN -REAL FUNCTION OF ONE REAL C ARGUMENT - THE INTEGRAND FUNCTION. C (EXTERNAL USER SUPPLIED). FUN MUST BE A C SMOOTH REAL VALUED FUNCTION ON THE CLOSED C INTERVAL (A,B). THE ACTUAL NAME FOR C FUN NEEDS TO BE DECLARED EXTERNAL IN THE C CALLING PROGRAM. C C NCOEF -INTEGER- THE REQUIRED NUMBER OF FOURIER COEF C C TOL -REAL - THE REQUIRED ABSOLUTE C ACCURACY OF THE COEFFICIENTS. NOTE THAT C IF A VALUE OF TOL .LT. THE MACHINE EPSILON C IS PROVIDED,THE ROUTINE TRIES TO DO THE C BEST IT CAN (SEE NOTE 1). C HOWEVER ..ESTER.. GREATER THAN THE C MACHINE EPSILON AND ..IERR.. NE 0 SHOULD C BE EXPECTED IN OUTPUT. C C MTOP -INTEGER- THE MAXIMUM NUMBER OF FUNCTION C EVALUATIONS IN FCOEBA IS MTOP+3. MTOP SHOULD C BE .GT. 16 AND AN INTEGER POWER OF 2. C IF MTOP IS ERRONEOUSLY SET .LT. 16 IT IS C SET TO 16 BY DEFAULT (TAKE CARE TO THE SIZE C THE WORKSPACE ARRAY). IF IT IS NOT SET TO C A POWER OF 2 THE ACTUAL VALUE USED IS THE C LARGEST POWER OF 2 LESS THAN MTOP. C IF MTOP IS SET NEGATIVE THEN ONE C OF THE FEATURES OF THE ABNORMAL TERMINATION C CRITERION IS SWITCHED OFF(SEE NOTE 1). C C IDER -INTEGER- ALLOWS THE CHOICE OF DIFFERENTIATIO C ROUTINE. C =0 USE THE STANDARD DIFFERENT.ROUTINE(CHPEN C 18 VALUES OF ..FUN.. IN THE INTERVAL C (A-(B-A)*7/64.,B+(B-A)*7/64.) ARE USED. C =1 USE A MORE ACCURATE AND EXPANSIVE DIFF. C ROUTINE(ENDACE).42 VALUES OF ..FUN.. IN C INTERVAL (A-(B-A)*19/64.,B+(B-A)*19/64.) C ARE USED. CHOOSE THIS OPTION IF THE C ACCURACY REQUIREMENTS ARE STRINGENT. C =2 USE A ONE-SIDED DIFFERENT.ROUTINE. C 14 VALUES OF ..FUN.. IN THE INTERVAL C (A,B) ARE USED. C =-P THE USER PROVIDES THE VALUE OF THE C FIRST P DERIVATIVES OF ..FUN.. AT ..A.. C AND ..B.. (SEE ..XLAMB.. BELOW). C P MUST BE LESS OR EQUAL TO 14,OTHERWISE C ONLY THE FIRST 14 DERIVATIVES ARE USED. C C IFIRST -INTEGER- NORMALLY MUST BE SET TO ZERO. C IF SET NONZERO IT CAUSES THE RESTART OF THE C PROCEDURE USING RESULTS FROM THE PREVIOUS C CALL TO THE TOP-END. IN THIS CASE ..IDER, C WORK,XLAMB,MACT,IERR.. FROM THE PREVIOUS C RUN MUST BE UNCHANGED (SEE ..IERR.. BELOW). C C WORK -REAL ARRAY(4*MTOP+20)- WORK SPACE C ARRAY. C C XLAMB -REAL(15)- C NORMALLY NOT USED. IF ..IDER.. IS SET C NEGATIVE THEN XLAMB MUST CONTAIN C (I) (I) C XLAMB(I)=FUN (B) - FUN (A) C I=1,2,...,ABS(IDER) C C I.E. THE DIFFERENCES OF THE FIRST ABS(IDER) C DERIVATIVES OF ..FUN.. AT ..B.. AND ..A.. C C C OUTPUT PARAMETERS - C,S,IERR,ESTER,MACT C C C -REAL(NCOEF)- ARRAY CONTAINING THE C FIRST NCOEF COSINE FOURIER COEFFICIENTS C C S -REAL(NCOEF)- ARRAY CONTAINING THE C FIRST NCOEF SINE FOURIER COEFFICIENTS C C /// WARNING : C(J),S(J) ARE THE (J-1)-TH /// C /// TRIGONOMETRIC FOURIER COEFFICIENTS /// C C IERR -INTEGER- ERROR INDICATOR C C IERR=-1 INPUT ERROR. ..NCOEF.. LE 0 OR REST C INVOKED WITH IERR FROM THE PREVIOUS C RUN NOT POSITIVE. C IERR=0 NORMAL TERMINATION.TOLERANCE SATISF C IERR.GT.0 PHYSICAL LIMIT TERMINATION. C (POSSIBILITY OF RESTART) C IERR=-2 ABNORMAL TERMINATION C C /// SEE NOTE 1. /// C C ESTER -REAL- THE ESTIMATED ABSOLUTE ACCURACY C NEGATIVE WHEN ..IERR.. NE 0. C C MACT -INTEGER- THE TOTAL NUMBER OF FUNCTION C EVALUATIONS ACTUALLY USED. C C*********************************************************************** C C NOTE 1. C C AFTER STAGE S (AFTER EXPENDING 2**S +3 FUNCTION VALUES) THE C ROUTINE FORMS AN ACCURACY ESTIMATE E(S) C C (1) IF E(S) .LT. TOL THE CALCULATION TERMINATES WITH IERR=0 C (THIS IS NORMAL TERMINATION) C C (2) IF (1) IS NOT SATISFIED, BUT 2**S = MTOP C THE CALCULATION TERMINATES WITH IERR=1 C (THIS IS PHYSICAL LIMIT TERMINATION) C C (3) IF (1) AND (2) ARE NOT SATISFIED BUT BOTH C C (A) E(S) .LT. ABNC1 * MAX(E(J):J=1,2,..,S-1) C C AND C C (B) E(S) .GT. ABNC2 * E(S-1) C C THE CALCULATION TERMINATES WITH IERR=-2 C (THIS IS AN ABNORMAL TERMINATION) C ABNC1=0.01 AND ABNC2=SQRT(0.1) HAVE BEEN CHOSEN AS A RESULT C OF EXPERIMENTATION. C C THE USER CAN ELIMINATE THE ABNORMAL FEATURE (3) BY USING C -ABS(MTOP) IN PLACE OF MTOP AS AN INPUT ARGUMENT.(THIS C REPLACES ABNC1 BY 0 INTERNALLY) C THE USER MAY ENSURE A RESULT BASED ON MTOP FUNCTION VALUES C BY BOTH USING -ABS(MTOP) IN PLACE OF MTOP AND SETTING C TOL=0.0E0. C C COMMENT. UNDER NORMAL CONDITIONS IT IS USUALLY THE CASE THAT C IF IERR.GE.1 ONE MIGHT OBTAIN BETTER ACCURACY BY RAISING THE C VALUE OF MTOP AND RE-RUNNING WITH IFIRST=1 C IF IERR=-2 ONLY MARGINALLY BETTER ACCURACY CAN BE OBTAINED AT C SIGNIFICANT EXPENSE SINCE THE PRESENT ACCURACY IS PROBABLY C GOVERNED BY NOISE LEVEL IN THE FUNCTION VALUES.NOTE THAT C IN THE LATTER CASE RESTART IS NOT ALLOWED. C********************************************************************** C FIRST DECLARATION STATEMENT OF FCOEFS C********************************************************************** IMPLICIT REAL (A-H,O-Z) INTEGER NCOEF,IERR,MTOP,MACT,IDER,IFIRST INTEGER IPM1,IPM2 DIMENSION C(NCOEF),S(NCOEF),WORK(*) DIMENSION XLAMB(*),EREST1(14),EREST2(14) EXTERNAL FUN DATA EPS/0.E0/ C C********************************************************************** C FIRST EXECUTABLE STATEMENT OF FCOEFS C********************************************************************** C IF (EPS.EQ.0.E0) EPS=R1MACH(4) IF(IFIRST.EQ.0)GOTO 20 IF(IERR.GT.0) GO TO 10 IERR = -1 RETURN C 10 CONTINUE IPM1=IERR IF(IDER.LT.0) INITM=MACT-3 IF(IDER.EQ.0) INITM=MACT-21 IF(IDER.EQ.1) INITM=MACT-45 IF(IDER.EQ.2) INITM=MACT-17 DO 15 I=1,INITM+2 WORK(INITM+ABS(MTOP)+5-I) = WORK(2*INITM+5-I) 15 CONTINUE GO TO 130 C 20 CONTINUE INITM=8 IF (NCOEF.GT.0) GO TO 30 IERR = -1 RETURN C 30 CONTINUE IF (IDER.GE.0) GOTO 40 IPM1 = MIN(1-IDER,15) IPM2 = IPM1-1 GO TO 110 C 40 CONTINUE C * * C * HBASE IS THE STEPSIZE FOR THE DIFFERENTIATION ROUTINES * C * * HBASE = (B-A)/64.E0 C C CHOOSE THE DIFFERENTIATION ROUTINE C IF (IDER.NE.1) GO TO 70 CALL ENDACE(B,14,HBASE,WORK,EREST1,FUN) CALL ENDACE(A,14,HBASE,XLAMB,EREST2,FUN) DO 50 K=1,14 IF (EREST1(K).LT.0.E0.OR.EREST2(K).LT.0.E0) GO TO 60 50 CONTINUE C IPM1 = 15 GO TO 90 C 60 CONTINUE IPM1 = K GO TO 90 C 70 CONTINUE IF (IDER.NE.2) GO TO 80 CALL ONESID(-1,B,HBASE,WORK,FUN) CALL ONESID(1,A,HBASE,XLAMB,FUN) IPM1 = 6 GO TO 90 C 80 CONTINUE IPM1 = 6 CALL CHPEND(B,HBASE,WORK,FUN) CALL CHPEND(A,HBASE,XLAMB,FUN) 90 CONTINUE C IPM2 = IPM1 - 1 DO 100 K = 1,IPM2 XLAMB(K) = WORK(K) - XLAMB(K) 100 CONTINUE C C CHECK FOR PERIODIC DERIVATIVES C 110 CONTINUE DO 120 K = 1,IPM2 I = IPM2 - K + 1 IF(ABS(XLAMB(I)).GT.EPS * 1.0D4)GO TO 130 IPM1 = IPM1 - 1 120 CONTINUE C 130 CONTINUE C CALL FCOEBA(A,B,FUN,NCOEF,TOL,MTOP,IPM1,INITM,IFIRST,WORK(1), + WORK(ABS(MTOP)+3),WORK(2*ABS(MTOP)+5),XLAMB,C,S, + IERR,ESTER,MACT) C IF(IFIRST.NE.0) MACT = MACT-1 IF(IDER.EQ.0) MACT = MACT+18 IF(IDER.EQ.1) MACT = MACT+42 IF(IDER.EQ.2) MACT = MACT+14 RETURN C C C*********************************************************************** C C ******* END OF SUBROUTINE FCOEFS ******* C C*********************************************************************** END SUBROUTINE FCOETD(A,B,NCOEF,MTOP,FUNVAL,WORK,C,S,IERR + ,ESTER) C**************************************************************** C********************************************************************** C TOP END FCOETD * C********************************************************************** C C TOP-END FCOETD ALLOWS TO COMPUTE TRIGONOMETRICAL FOURIER C COEFFICIENTS WHEN THE INTEGRAND FUNCTION FUN IS PROVIDED C BY MEANS OF A SET OF FUNCTION VALUES AT EQUALLY SPACED POINTS C OF THE INTEGRATION INTERVAL (A,B). C C C INPUT PARAMETERS - A,B,NCOEF,MTOP,FUNVAL,WORK. C C A,B -REAL- THE BOUNDS OF INTEGRATION C INTERVAL C C NCOEF -INTEGER- THE REQUIRED NUMBER OF FOURIER COE C C MTOP -INTEGER- THE NUMBER OF FUNCTION VALUES C STORED IN THE ARRAY FUNVAL. IT MUST BE C MTOP = 2**K + 1 , K .GE.4 AND LE 16 C IF MTOP IS ERRONEOUSLY SET IERR=-1 IS C RETURNED. C C FUNVAL -REAL(MTOP)- ARRAY CONTAINING THE C FUNCTION VALUES. C THE VALUE OF FUN(A+(I-1)*(B-A)/(MTOP-1)) IS C STORED IN FUNVAL(I) , I=1,2,..,MTOP. C C WORK -REAL(4*MTOP+20)- WORK SPACE ARRAY C C C OUTPUT PARAMETERS - C,S,IERR,ESTER C C C -REAL(NCOEF)- ARRAY CONTAINING THE C FIRST NCOEF COSINE FOURIER COEFFICIENTS C C S -REAL(NCOEF)- ARRAY CONTAINING THE C FIRST NCOEF SINE FOURIER COEFFICIENTS C C /// WARNING : C(J),S(J) ARE THE (J-1)-TH /// C /// TRIGONOMETRIC FOURIER COEFFICIENTS /// C C IERR -INTEGER- ERROR INDICATOR C C IERR=-1 INPUT ERROR.NCOEF.LT.0 OR MTOP C NOT EQUAL TO 2**K +1 , K.GE.4 C IERR.GE.0 NORMAL TERMINATION. C C ESTER - REAL- THE ESTIMATED ABSOLUTE ACCURACY. C C C*********************************************************************** C C*********************************************************************** C FIRST DECLARATION STATEMENT OF FCOETD C*********************************************************************** IMPLICIT REAL (A-H,O-Z) INTEGER NCOEF,MTOP,IERR INTEGER IPOW2,IXVAL,IPM1,IPM2 DIMENSION FUNVAL(MTOP),WORK(*),C(NCOEF),S(NCOEF) DIMENSION XLAMB(6) C EXTERNAL DUMFUN C DATA EPS/0.E0/ C*********************************************************************** C FIRST EXECUTABLE STATEMENT OF FCOETD C*********************************************************************** IF (EPS.EQ.0.E0) EPS=R1MACH(4) C C **** TEST ON INPUT ARGUMENTS **** C IF (NCOEF.GT.0.AND.MTOP.GT.16) GO TO 10 IERR = -1 RETURN C 10 CONTINUE IPOW2 = 16 DO 20 K=1,12 IF (MTOP.EQ.IPOW2+1) GO TO 30 IPOW2 = IPOW2*2 20 CONTINUE C IERR = -1 RETURN C 30 CONTINUE C C **** COMPUTE APPROXIMATED DERIVATIVES **** C NMTOP = MTOP -1 HBASE = (B-A) / NMTOP IXVAL = MTOP CALL ONESMP(-1,FUNVAL,IXVAL,HBASE,WORK) IXVAL = 1 CALL ONESMP(1,FUNVAL,IXVAL,HBASE,XLAMB) IPM1 = 6 IPM2 = 5 DO 40 K=1,IPM2 XLAMB(K) = WORK(K)-XLAMB(K) 40 CONTINUE C C C CHECK FOR PERIODIC DERIVATIVES C DO 50 K=1,IPM2 I = IPM2-K+1 IF (ABS(XLAMB(I)).GT.EPS*1.0D4) GO TO 60 IPM1 = IPM1-1 50 CONTINUE C 60 CONTINUE IFIRST = 0 INITM = -NMTOP DO 70 K=1,MTOP WORK(K) = FUNVAL(K) 70 CONTINUE C TOL=0.E0 CALL FCOEBA(A,B,DUMFUN,NCOEF,TOL,NMTOP,IPM1,INITM,IFIRST + ,WORK(1),WORK(ABS(MTOP)+3),WORK(2*ABS(MTOP)+5) + ,XLAMB,C,S,IERR,ESTER,MACT) C ESTER=-ESTER RETURN C*********************************************************************** C C ******* END OF SUBROUTINE FCOETD ******* C C*********************************************************************** END REAL FUNCTION DUMFUN(X) C THIS IS A DUMMY FUNCTION, JUST TO FILL THE CALLING SEQUENCE C OF FCOEBA IN FCOETD. IN FACT FCOETD USES SAMPLED FUNCTION C VALUES STORED IN FUNVAL. REAL X DUMFUN = X RETURN END SUBROUTINE FCOEBA(A,B,FUN,NCOEF,TOL,MTOP,IPM1,INITM,IFIRST, + WORK1,WORK2,WORK3,XLAMB,C,S,IERR,ESTER,MACT) C*********************************************************************** C C C ******** BEGINNING OF SUBROUTINE **** FCOEBA ************ C C C*********************************************************************** C C C*********************************************************************** C C INPUT PARAMETERS - A,B,FUN,NCOEF,TOL,MTOP,IPM1,INITM,IFIRST,WORK1, C WORK2,WORK3,XLAMB. C C A,B,FUN,NCOEF,TOL,MTOP,IFIRST AS IN FCOEFS C C IPM1 -INTEGER- THE NUMBER OF INPUT DERIVATIVES +1 C C INITM -INTEGER- THE NUMBER OF FUNCTION EVALUATIONS C USED FOR FIRST APPROXIMATION. C IT IS SET NEGATIVE IF FUN IS SAMPLED (TOP-END C FCOETD) C C WORK1,WORK2,WORK3 -REAL ARRAYS- WORKSPACE C C XLAMB -REAL(IPM1)- THE ARRAY OF LENGTH C IPM1 CONTAINING THE FOLLOWING C QUANTITIES C C (I) (I) C XLAMB(I)=FUN (B) - FUN (A) C C I=1, 2, ..., IPM1 -1 C C OUTPUT PARAMETERS - C,S,IERR,ESTER,MACT AS IN FCOEFS C C C*********************************************************************** C FIRST DECLARATION STATEMENT OF FCOEBA C*********************************************************************** C IMPLICIT REAL (A-H,O-Z) LOGICAL RNDLEV INTEGER NCOEF,IPM1,IERR,MTOP,MACT,INITM,IFIRST INTEGER MAXM,IPM2,NHALF,N,NQUART DIMENSION XLAMB(*),C(NCOEF),S(NCOEF),WORK1(*),WORK2(*),WORK3(*) DATA EPS/0.E0/,TWOPI/0.E0/ C * * C******************************************************************* C * FIRST EXECUTABLE STATEMENT OF FCOEBA * C******************************************************************* C **** INITIALIZATION **** C -------------- C IF (EPS.EQ.0.E0) EPS=R1MACH(4) IF (TWOPI.EQ.0.E0) TWOPI=8.E0*ATAN(1.E0) C XL = B - A MACT = 0 C * * C * THE ABNORMAL TERMINATION CRITERION IS BASED ON THE VALUES * C * OF COSTANTS ABNC1,ABNC2 AND IT IS SWITCHED OFF WHEN MTOP * C * IS NEGATIVE * C * * ABNC1 = 0.01E0 ABNC2 = SQRT(0.1E0) IF (MTOP.LT.0) ABNC1 = 0.E0 C RNDLEV = .FALSE. MAXM = MAX(16,ABS(MTOP)) IERR = 0 ESTMAX = 0.E0 N = ABS(INITM) NHALF = N/2 NQUART = N/4 IPM2 = IPM1 -1 IF (INITM.LT.0) GO TO 20 IF (IFIRST.NE.0) GO TO 130 DO 10 I = 1,MAXM+2 WORK1(I) = 0.0 WORK2(I) = 0.0 10 CONTINUE C FA = FUN(A) FB = FUN(B) C C * * C * MODIFICATION OF ARRAY XLAMB * C * * 20 CONTINUE IF (IPM2.EQ.0) GO TO 40 C XLK = XL** (IPM2) C DO 30 I = 1,IPM2 K = IPM1 - I XLAMB(K+1) = XLAMB(K)*XLK XLK = XLK/XL 30 CONTINUE C 40 CONTINUE IF (INITM.GT.0) XLAMB(1) = FB - FA IF (INITM.LT.0) XLAMB(1) = WORK1(N+1) - WORK1(1) IF (ABS(XLAMB(1)).LE.EPS*1.0D4 .AND. IPM1.EQ.1) IPM1 = 0 C C********************************************************************** C C **** FIRST APPROXIMATION **** C ------------------- C DO 60 K = 1,NHALF - 1 ARGNOR = REAL(K)/REAL(N) CALL HPMIN1(ARGNOR,XLAMB,IPM1,HP,HPS) IF (INITM.GT.0) GO TO 50 WORK1(K+1) = WORK1(K+1) - HP WORK1(N-K+1) = WORK1(N-K+1) - HPS GO TO 60 50 CONTINUE APP = ARGNOR*XL WORK1(K+1) = FUN(APP+A) - HP WORK1(N-K+1) = FUN(B-APP) - HPS 60 CONTINUE C CALL HPMIN1(0.0E0,XLAMB,IPM1,HP,HPS) CALL HPMIN1(0.5E0,XLAMB,IPM1,HPHAL,DUMMY) IF (INITM.GT.0) GO TO 70 WORK1(1) = (WORK1(1)+WORK1(N+1)-(HP+HPS)) * 0.5E0 WORK1(NHALF+1) = WORK1(NHALF+1) - HPHAL WORK1(N+1) = 0.E0 GO TO 80 70 CONTINUE WORK1(1) = 0.5E0* (FA+FB-(HP+HPS)) WORK1(NHALF+1) = FUN(A+XL/2.0E0) - HPHAL C 80 CONTINUE C C DFT -N POINTS ON WORK1 - IN PLACE. C CALL RFFTI(N,WORK3) CALL RFFTF(N,WORK1,WORK3) C OVN = 1./REAL(N) WORK1(N+2) = 0.E0 DO 90 K=N+1,3,-1 WORK1(K) = OVN*WORK1(K-1) IF(MOD(K,2).EQ.0) WORK1(K) = -WORK1(K) 90 CONTINUE WORK1(2) = 0.E0 WORK1(1) = OVN*WORK1(1) C********************************************************************** C **** MAIN LOOP **** C ----------- C C (REPEAT UNTILL THE STOPPING CRITERION IS SATISFIED) C 100 CONTINUE C WORK2(1) = WORK1(1) WORK2(2) = WORK1(2) WORK2(N+1) = WORK1(N+1) WORK2(N+2) = WORK1(N+2) DO 110 K=3,N WORK2(K) = 2.0E0*WORK1(K) 110 CONTINUE C CALL TRASF(WORK2,N,1) C TRIGG1=0.E0 TRIGG2=0.E0 DO 120 K=1,N+2,2 TRIGG1=TRIGG1+WORK2(K) TRIGG2=TRIGG2+WORK2(K+1) 120 CONTINUE C 130 CONTINUE NQUART = NHALF NHALF = N N = 2*N C * * C * ESTIMATE THE ACCURACY OF THE APPROXIMATION * C * * IF (INITM.GT.0) GO TO 140 GIBB1 = 0.E0 GIBB2 = 0.E0 GO TO 150 C 140 CONTINUE ARGNOR = 1./REAL(N) AN = ARGNOR*XL CALL HPMIN1(ARGNOR,XLAMB,IPM1,HP,HPS) APP1 = FUN(AN+A) - HP APP2 = FUN(B-AN) - HPS GIBB1 = ABS(APP1-TRIGG1) GIBB2 = ABS(APP2-TRIGG2) 150 CONTINUE EST = SQRT((WORK1(NHALF-1)*2.E0)**2+ (WORK1(NHALF)*2.E0)**2+ + WORK1(NHALF+1)**2) ESTER = MAX(2.E0*GIBB1,2.E0*GIBB2,EST,EPS) C * * C * THIS IS THE STOPPING CRITERION. IF SATISFIED JUMP OUT OF C * THE MAIN LOOP * C * * C * IF THE ACCURACY REQUIREMENT IS ACHIEVED STOP THE UPDATING * C * * IF (ESTER.LE.TOL) GO TO 240 C * * C * * C * IF IT IS IMPOSSIBLE TO ACHIEVE THE USER REQUIRED ACCURACY * C * WITH THE COST MAXM THEN FOURIER COEFFICIENTS OF BEST ACCURACY * C * COMPUTABLE BY THE PROCEDURE ARE RETURNED. * C * * IF (N.GT.MAXM) GO TO 230 C IF (RNDLEV) GO TO 160 IF (NHALF.EQ.4) ESTMAX = ESTER IF (ESTMAX.LT.ESTER) ESTMAX = ESTER IF (NHALF.GE.8 .AND. ESTER.LT.ABNC1*ESTMAX) RNDLEV = .TRUE. GO TO 170 C * * C * IF THE ROUNDOFF LEVEL HAS BEEN ACHIEVED AND THE PRESENT * C * ESTIMATED ERROR IS GT THE PREVIOUS ONE TIMES THE CONSTANT * C * ABCN2 THEN STOP THE UPDATING STAGE * C * * 160 CONTINUE C IF (ESTER.GE.ABNC2*PREVES) GO TO 230 C 170 CONTINUE C * * C * COMPUTE NEW APPROXIMATION * C * * PREVES = ESTER C WORK2(1) = APP1 WORK2(NHALF) = APP2 C C DO 180 K = 2,NQUART ARGNOR = REAL(2*K-1)/REAL(N) AN = ARGNOR*XL CALL HPMIN1(ARGNOR,XLAMB,IPM1,HP,HPS) WORK2(K) = FUN(AN+A) - HP WORK2(NHALF+1-K) = FUN(B-AN) - HPS 180 CONTINUE C * * C * CALLING FFT + TRASF A MIDPOINT TRAPEZOIDAL RULE IS COMPUTED * C * * C DFT -NHALF POINTS ON WORK2 - IN PLACE C TEMP=WORK2(NHALF) DO 190 K=NHALF-1,1,-1 WORK2(K+1)=WORK2(K) 190 CONTINUE C WORK2(1)=TEMP WORK2(NHALF+1)=0.E0 WORK2(NHALF+2)=0.E0 CALL RFFTI(NHALF,WORK3) CALL RFFTF(NHALF,WORK2,WORK3) C OVNH = 1./REAL(NHALF) WORK2(NHALF+2) = 0.E0 DO 200 K=NHALF+1,3,-1 WORK2(K) = OVNH*WORK2(K-1) IF(MOD(K,2).EQ.0) WORK2(K) = -WORK2(K) 200 CONTINUE WORK2(2) = 0.E0 WORK2(1) = OVNH*WORK2(1) C CALL TRASF(WORK2,NHALF,0) C DO 210 K = 1,NHALF,2 IK=N+3-K WORK1(IK) = 0.5* (WORK2(K+1)-WORK1(K+1)) WORK1(IK-1) = 0.5* (WORK1(K)-WORK2(K)) 210 CONTINUE C DO 220 K = 1,NHALF+2 WORK1(K) = 0.5* (WORK1(K)+WORK2(K)) 220 CONTINUE C GO TO 100 C C **** END OF MAIN LOOP **** C*********************************************************************** C **** COMPUTE FINAL RESULTS **** C --------------------- C 230 CONTINUE IERR = IPM1 IF (RNDLEV .AND. ESTER.GE.ABNC2*PREVES) IERR = -2 ESTER = -ESTER C 240 CONTINUE C MACT = NHALF+3 C(1) = WORK1(1) S(1) = 0.0 IF (NCOEF.EQ.1) RETURN C * * C * PERIODIC FUNCTIONS AND NON PERIODIC FUNCTIONS WITH * C * PERIODIC DERIVATIVES ARE TREATED AS A SPECIAL CASE * C * * IF (IPM1.GT.0) GO TO 270 C C DO 260 K = 2,NCOEF IF (K.GT.NQUART) GO TO 250 C(K) = 2.0E0*WORK1(2*K-1) S(K) = 2.0E0*WORK1(2*K) GO TO 260 C 250 CONTINUE C(K) = 0.0E0 IF (K.EQ.NQUART+1) C(K) = WORK1(NHALF+1) S(K) = 0.0E0 260 CONTINUE C IF (IERR.EQ.0.AND.ESTER.LT.0.E0) IERR = 1 C RETURN C 270 CONTINUE IF (IPM1.NE.1) GO TO 300 C C DO 290 K = 2,NCOEF IF (K.GT.NQUART) GO TO 280 C(K) = 2.0E0*WORK1(2*K-1) S(K) = -2.0E0*XLAMB(1)/ (TWOPI*REAL(K-1)) + 2.0E0*WORK1(2*K) GO TO 290 C 280 CONTINUE C(K) = 0.0E0 IF (K.EQ.NQUART+1) C(K) = WORK1(NHALF+1) S(K) = -2.0E0*XLAMB(1)/ (TWOPI*REAL(K-1)) 290 CONTINUE C C RETURN C 300 CONTINUE C * * C * THIS IS THE GENERAL CASE * C * * C DO 340 K = 2,NCOEF C OV2PK = 1.E0 / (TWOPI*REAL(K-1)) OV2PK2 = OV2PK*OV2PK CFACT = OV2PK2 SFACT = -OV2PK CSUM = XLAMB(2) * CFACT SSUM = XLAMB(1) * SFACT IF (IPM1.EQ.2) GO TO 320 DO 310 KK=3,IPM1,2 CFACT = -CFACT * OV2PK2 SFACT = -SFACT * OV2PK2 IF (KK+1.LE.IPM1) CSUM = CSUM + XLAMB(KK+1) * CFACT SSUM = SSUM + XLAMB(KK) * SFACT 310 CONTINUE C 320 CONTINUE IF (K.GT.NQUART) GO TO 330 C(K) = 2.0E0*(CSUM + WORK1(2*K-1)) S(K) = 2.0E0*(SSUM + WORK1(2*K)) GO TO 340 C 330 CONTINUE C(K) = 2.0E0*CSUM S(K) = 2.0E0*SSUM IF (K.EQ.NQUART+1) C(K) = C(K) + WORK1(NHALF+1) C 340 CONTINUE C C C C RETURN C*********************************************************************** C C ************** END OF SUBROUTINE FCOEBA *************** C C*********************************************************************** C END SUBROUTINE HPMIN1(ARGNOR,XLAMB,IPM1,HP,HPS) C*********************************************************************** C C BEGINNING OF SUBROUTINE ****** HPMIN1 ****** C C C*********************************************************************** C C C HPMIN1 EVALUATES THE POLYNOMIAL C C HPMIN1(X)=SUM(I=1,IPM1)(XLAMB(I)*BERNOULLI(I,X)/FACTORIAL(I)) C C WHERE BERNOULLI(I,X) IS THE I-TH BERNOULLI POLYNOMIAL C C AT THE POINTS ..ARGNOR.. AND 1-ARGNOR. C IT RETURNS SUCH THE VALUES IN HP AND HPS RESPECTIVELY C C THE DEGREE OF POLYNOMIAL HPMIN1 IS IPM1. MAXIMUM VALUE FOR C IPM1 IS 15. C C C INPUT - ARGNOR,XLAMB,IPM1 C ARGNOR - REAL - (0,1)-NORMALIZED ARGUMENT C XLAMB -REAL(IPM1)-ARRAY OF FUNCTION DERIVATIVES C SEE INITIAL COMMENTS OF FCOEBA. C IPM1 - INTEGER- THE DEGREE OF POLYNOMIAL HPMIN1 C C OUTPUT - HP,HPS C HP - REAL - THE VALUE OF HPMIN1 AT POINT ARGNOR C HPS - REAL - THE VALUE OF HPMIN1 AT POINT 1-ARGNOR C C*********************************************************************** C*********************************************************************** C FIRST DECLARATION STATEMENT OF HPMIN1 C*********************************************************************** IMPLICIT REAL (A-H,O-Z) INTEGER IPM1 C DIMENSION XLAMB(*) C DIMENSION BERPOL(15) LOGICAL ODD C*********************************************************************** C FIRST EXECUTABLE STATEMENT OF HPMIN1 C*********************************************************************** IF (IPM1.GT.0) GO TO 10 HP = 0.0E0 HPS = 0.0E0 RETURN C 10 CONTINUE BERPOL(1) = ARGNOR - 0.5E0 IF (IPM1.EQ.1) GO TO 20 Y = ARGNOR* (ARGNOR-1.0E0) Z = Y* (ARGNOR-0.5E0) BERPOL(2) = Y + (1.0E0/6.0E0) IF (IPM1.EQ.2) GO TO 20 BERPOL(3) = Z IF (IPM1.EQ.3) GO TO 20 YTO2 = Y*Y BERPOL(4) = YTO2 - (1.0E0/30.0E0) IF (IPM1.EQ.4) GO TO 20 BERPOL(5) = (Y- (1.0E0/3.0E0))*Z IF (IPM1.EQ.5) GO TO 20 BERPOL(6) = Y*YTO2 - 0.5E0*YTO2 + (1.0E0/42.0E0) IF (IPM1.EQ.6) GO TO 20 Y2 = (ARGNOR-0.5E0)* (ARGNOR-0.5E0) Y4 = Y2*Y2 BERPOL(7) = (Y4-1.5E0*Y2+ (31.0E0/48.0E0))*Z IF (IPM1.EQ.7) GO TO 20 X2 = ARGNOR*ARGNOR X4 = X2*X2 X6 = X4*X2 X7 = X6*ARGNOR X8 = X4*X4 TERM = (2.0E0*X2-7.0E0*X4+14.0E0*X6)/3.0E0 BERPOL(8) = X8 - 4.0E0*X7 + TERM - (1.0E0/30.0E0) IF (IPM1.EQ.8) GO TO 20 TERM = X8 - 4.5E0*X7 + 6.0E0*X6 - 4.2E0*X4 + 2.0E0*X2 - 0.3E0 BERPOL(9) = ARGNOR*TERM IF (IPM1.EQ.9) GO TO 20 X9 = ARGNOR*X8 X10 = X4*X6 TERM = X10 - 5.0E0*X9 + 7.5E0*X8 - 7.0E0*X6 + 5.0E0*X4 BERPOL(10) = TERM - 1.5E0*X2 + (5.0E0/66.0E0) IF (IPM1.EQ.10) GO TO 20 TERM = X10 - 5.5E0*X9 - 11.0E0* (X6-X4) - 5.5E0*X2 BERPOL(11) = ARGNOR* (TERM+ ((5.0E0+55.0E0*X8)/6.0E0)) IF (IPM1.EQ.11) GO TO 20 X12 = X4*X8 TERM = X12 - 6.0E0*ARGNOR*X10 + 11.0E0*X10 - 16.5E0* (X8+X4) BERPOL(12) = TERM + 22.0E0*X6 + 5.0E0*X2 - (691.0E0/2730.0E0) IF (IPM1.EQ.12) GO TO 20 TERM = X12 + 13.0E0*X10 - 143.0E0/6.0E0*X8 + 286.0E0/7.0E0*X6 - + 42.9E0*X4 BERPOL(13) = (TERM+65.0E0/3.0E0*X2-691.0E0/210.0E0)*ARGNOR - + 6.5E0*X12 IF (IPM1.EQ.13) GO TO 20 X14 = X12*X2 TERM = X14 + (-7.0E0*ARGNOR+91.0E0/6.0E0)*X12 - 100.1/3.0E0*X10 BERPOL(14) = TERM + 71.5*X8 - 100.1E0*X6 + 455.0E0/6.0E0*X4 - + 69.1E0/3.0E0*X2 + 7.0E0/6.0E0 IF (IPM1.EQ.14) GO TO 20 TERM = X14 + 17.5E0*X12 - 45.5E0*X10 + 715.0E0/6.0E0*X8 - + 214.5E0*X6 BERPOL(15) = (TERM+227.5E0*X4-691.0E0/6.0E0*X2+17.5E0)*ARGNOR - + 7.5E0*X14 C 20 CONTINUE HP = XLAMB(1)*BERPOL(1) HPS = -HP FATT = 1.0E0 ODD = .FALSE. IF (IPM1.EQ.1) RETURN C C DO 40 I = 2,IPM1 FATT = FATT*REAL(I) TEMP = XLAMB(I)*BERPOL(I)/FATT HP = HP + TEMP IF (ODD) GO TO 30 HPS = HPS + TEMP ODD = .TRUE. GO TO 40 C 30 CONTINUE ODD = .FALSE. HPS = HPS - TEMP 40 CONTINUE C C RETURN C********************************************************************** C C END OF SUBROUTINE ****** HPMIN1 ****** C C********************************************************************** END SUBROUTINE TRASF(VETT,N,NCODE) C C********************************************************************** C C BEGINNING OF SUBROUTINE ****** TRASF ****** C C********************************************************************** C C GIVEN THE INPUT ARRAY VETT OF N+2 COMPONENTS TRASF COMPUTES IN C PLACE THE VALUES (I=1,2,..,N/2+1) C C A(I)=VETT(I)*COS(PI*(I-1)/N)+VETT(I+1)*SIN(PI*(I-1)/N) C B(I)=VETT(I)*COS(PI*(I-1)/N)-VETT(I+1)*SIN(PI*(I-1)/N) IF NCODE=1 C B(I)=-VETT(I)*SIN(PI*(I-1)/N)+VETT(I+1)*COS(PI*(I-1)/N) IF NCODE= C C IT RETURNS A(I) IN VETT(2*I-1) AND B(I) IN VETT(2*I) C C C C INPUT - VETT,N,NCODE C C VETT -REAL(N+2)- ARRAY TO BE TRANSFORMED C C N -INTEGER- N+2 IS THE SIZE OF VETT C C NCODE -INTEGER- ALLOWS THE CHOICE OF B(I) C C C OUTPUT - VETT C C VETT -REAL(N+2)- TRANSFORMED ARRAY C C*********************************************************************** C*********************************************************************** C FIRST DECLARATION STATEMENT OF TRASF C*********************************************************************** C IMPLICIT REAL (A-H,O-Z) INTEGER N,NCODE C INTEGER N2 DIMENSION VETT(*) C * * C * THESE ARE THE VALUES OF PI AND COS(PI / 4.) * C * WITH 18 DECIMAL FIGURES. CHANGE THIS CARD IF A * C * MORE PRECISE ARITHMETIC IS USED. * C * * DATA PI/0.E0/,COSPIQ/0.E0/ C C*********************************************************************** C FIRST EXECUTABLE STATEMENT OF TRASF C*********************************************************************** C **** INITIALIZATION **** C IF (PI.EQ.0.E0) PI = 4.E0*ATAN(1.E0) IF (COSPIQ.EQ.0.E0) COSPIQ = SQRT(0.5E0) PIN = PI/N COSPIN = COS(PIN) SINPIN = SIN(PIN) N2 = N/2 C C *** COMPUTATION OF COS((PI/2N)(K-1)) , SIN((PI/2N)(K-1)) *** C DO 60 K = 1,N2,2 IF (K.NE.1) GO TO 10 C = 1 S = 0.0 GO TO 20 C 10 CONTINUE APP1 = C*COSPIN - S*SINPIN S = S*COSPIN + C*SINPIN C = APP1 20 CONTINUE C **** COMPUTE THE SUM **** C APP1 = VETT(K+1)*S APP2 = VETT(K)*S APP3 = VETT(K)*C VETT(K) = APP3 + APP1 IF(NCODE.EQ.0) GO TO 30 VETT(K+1)=APP3-APP1 GO TO 40 30 CONTINUE VETT(K+1) = VETT(K+1)*C - APP2 40 CONTINUE I = N + 3 - K APP1 = VETT(I)*C APP2 = VETT(I-1)*C APP3 = VETT(I-1)*S VETT(I-1) = APP3 + APP1 IF(NCODE.EQ.0) GO TO 50 VETT(I)= APP3 - APP1 GO TO 60 50 CONTINUE VETT(I) = VETT(I)*S - APP2 60 CONTINUE C I = N2 + 1 APP1 = VETT(I+1)*COSPIQ APP2 = VETT(I)*COSPIQ VETT(I) = APP2 + APP1 IF(NCODE.EQ.0) GO TO 70 VETT(I+1)=APP2-APP1 GO TO 80 70 CONTINUE VETT(I+1) = APP1 - APP2 80 CONTINUE C RETURN C C*********************************************************************** C C END OF SUBROUTINE ****** TRASF ****** C C*********************************************************************** C********************************************************************** END C C COLLECTION OF NUMERICAL DIFFERENTIATION ROUTINES FOR C THE PACKAGE FOURCO. C SINGLE PRECISION VERSION C SUBROUTINE ENDACE(XVAL,NDER,HBASE,DER,EREST,FUN) C C C *****EVALUATE NUMERICAL DERIVATIVES AND CORRESPONDING ERRORS***** C C C *****PURPOSE***** C C A SUBROUTINE FOR NUMERICAL DIFFERENTIATION AT A POINT. C IT RETURNS A SET OF APPROXIMATIONS TO THE J-TH ORDER DERIVATIVE C (J=1,2,...,14) OF FUN(X) EVALUATED AT X=XVAL AND, FOR EACH C DERIVATIVE, AN ERROR ESTIMATE ( WHICH INCLUDES THE EFFECT OF C AMPLIFICATION OF ROUND-OFF ERRORS ). C C C *** INPUT PARAMETERS *** C C C (1) XVAL REAL. THE ABSCISSA AT WHICH THE SET OF C DERIVATIVES IS REQUIRED. C C (2) NDER INTEGER. THE HIGHEST ORDER DERIVATIVE REQUIRED. C IF(NDER.GT.0) ALL DERIVATIVES UP TO MIN(NDER,14) C ARE CALCULATED. C IF(NDER.LT.0 AND NDER EVEN) EVEN ORDER DERIVATIVES C UP TO MIN(-NDER,14) ARE CALCULATED. C IF(NDER.LT.0 AND NDER ODD) ODD ORDER DERIVATIVES C UP TO MIN(-NDER,13) ARE CALCULATED. C C (3) HBASE REAL. A STEP LENGTH. C C (6) FUN THE NAME OF A REAL FUNCTION SUBPROGRAM, C WHICH IS REQUIRED BY THE ROUTINE AS A SUBPROGRAM C AND WHICH REPRESENTS THE FUNCTION BEING DIFFERENTIATED. C THE ROUTINE REQUIRES 21 FUNCTION EVALUATIONS FUN(X) C LOCATED AT X=XVAL AND AT C X=XVAL+(2*J-1)*HBASE, J=-9,-8,...,9,10. C THE FUNCTION VALUE AT X=XVAL IS DISREGARDED WHEN C ONLY ODD ORDER DERIVATIVES ARE REQUIRED. C C C *** OUTPUT PARAMETERS *** C C C (4) DER(J) J=1,2,...,14. REAL. A VECTOR OF C APPROXIMATIONS TO THE J-TH DERIVATIVE OF FUN(X) C EVALUATED AT X=XVAL. C C (5) EREST(J) J=1,2,...,14. REAL. A VECTOR OF C ESTIMATES OF THE ABSOLUTE ACCURACY OF DER(J). C THESE ARE NEGATIVE WHEN EREST(J).GT.ABS(DER(J)), OR C WHEN, FOR SOME OTHER REASON THE ROUTINE IS DOUBTFUL C ABOUT THE VALIDITY OF THE RESULT. C C C IMPLICIT REAL (A-H,O-Z) DIMENSION DER(14),EREST(14) DIMENSION RANGE(7),RTAB(10,7),TZEROM(10),ACOF(10) DATA ZERO,ONE,THRON2,TWO / 0.0E0,1.0E0,1.5E0,2.0E0/ C C C QUIT ON ZERO HBASE AND ZERO NDER. SET CONTROL PARAMETER NDPAR. C NDPAR = +1 ODD-ORDER DERIVATIVES ONLY. C NDPAR = 0 ALL DERIVATIVES. C NDPAR = -1 EVEN-ORDER DERIVATIVES ONLY. C C IF(HBASE.EQ.ZERO) GOTO 900 NDERP=NDER NDPAR=0 IF(NDER) 10,900,20 10 NDPAR=4*(NDER/2)-2*NDER-1 NDERP=-NDER 20 CONTINUE IF(NDERP.GT.14) NDERP=14 C C C NEXT, EVALUATE 21 FUNCTION VALUES, SET ACOF(N), AND SET FIRST C COLUMN OF RTAB FOR ODD ORDER DERIVATIVES AND STORE IN TZEROM C THE FIRST COLUMN OF RTAB FOR EVEN ORDER DERIVATIVES. C C FZ=FUN(XVAL) DO 30 N=1,10 NN1=2*N-1 NSQ=NN1*NN1 ACOF(N)=NSQ XNUM=NN1 H=HBASE*XNUM TEMP=FUN(XVAL+H) TERM=FUN(XVAL-H) RTAB(N,1)=(TEMP-TERM)/(TWO*H) TZEROM(N)=(TEMP-FZ-FZ+TERM)/(TWO*H**2) 30 CONTINUE C C C SET UP FOR ODD-ORDER DERIVATIVES. C C IF(NDPAR.EQ.-1) GOTO 40 NPAR=1 NSMAX=(NDERP+1)/2 IF(NSMAX.LE.0) GOTO 701 GOTO 60 C C C SET UP FOR EVEN-ORDER DERIVATIVES C C 40 CONTINUE NPAR=-1 IF(NDPAR.EQ.+1) GOTO 900 NSMAX=NDERP/2 IF(NSMAX.LE.0) GOTO 900 C C C SET THE FIRST COLUMN OF THE T-TABLE FOR EVEN ORDER DERIVATIVES C C DO 50 J=1,10 50 RTAB(J,1)=TZEROM(J) 60 CONTINUE C C ODD-ORDER AND EVEN ORDER PATHS JOIN UP HERE C NEGER=0 HSQ=HBASE*HBASE FACT=ONE HFACT=ONE DO 700 NS= 1,NSMAX C C C FIRST CALCULATE ELEMENTS OF T-TABLE FOR CURRENT NS VALUE C C NPMIN=NS IF(NS.EQ.1) NPMIN=2 DO 250 NP= NPMIN,7 NKTOP=10-NP+1 NPR=NP-NS+1 DO 200 NK= 1,NKTOP J=NK+NP-1 A=ACOF(J) B=ACOF(NK) TERM=ZERO TEMP=ZERO IF(NP.NE.NS) TEMP=A*RTAB(NK,NPR-1)-B*RTAB(NK+1,NPR-1) IF(NS.NE.1) TERM=-RTAB(NK,NPR)+RTAB(NK+1,NPR) RTAB(NK,NPR)= (TERM+TEMP)/(A-B) 200 CONTINUE 250 CONTINUE C C NOW CALCULATE NUP,NLO AND RANMIN. C DO 370 NP=NS,7 NTOP=11-NP NPR=NP-NS+1 XMAX=RTAB(1,NPR) XMIN=RTAB(1,NPR) NMAX=1 NMIN=1 C DO 340 NK=2,NTOP TEMP=RTAB(NK,NPR) IF(TEMP.GT.XMAX) NMAX=NK IF(TEMP.GT.XMAX) XMAX=TEMP IF(TEMP.LT.XMIN) NMIN=NK IF(TEMP.LT.XMIN) XMIN=TEMP 340 CONTINUE C RANGE(NP)=XMAX-XMIN IF(NP.NE.NS) GOTO 350 RANMIN=RANGE(NP) GOTO 355 350 CONTINUE IF(RANGE(NP).GE.RANMIN) GOTO 360 RANMIN=RANGE(NP) 355 CONTINUE NPMIN=NP NUP=NMAX NLO=NMIN 360 CONTINUE 370 CONTINUE C C TERM=ZERO NTOP=11-NPMIN NTOPM2=NTOP-2 XNUM=NTOPM2 J=NPMIN-NS+1 DO 440 NK= 1,NTOP IF(NK.EQ.NUP) GOTO 439 IF(NK.EQ.NLO) GOTO 439 TERM=TERM+RTAB(NK,J) 439 CONTINUE 440 CONTINUE TERM=TERM/XNUM C C JNS=2*NS-1 IF(NPAR.LT.0) JNS=2*NS XJNS=JNS FACT=FACT*XJNS DER(JNS)=TERM*FACT/HFACT IF(JNS.EQ.10) RANMIN=THRON2*RANMIN IF(JNS.EQ.11) RANMIN=THRON2*RANMIN IF(JNS.GE.12) RANMIN=TWO*RANMIN EREST(JNS)=RANMIN*FACT/HFACT C C SET SIGN OF EREST. EREST NEGATIVE EITHER IF IT SWAMPS DER,OR IF C TWO PREVIOUS CONSECUTIVE ERESTS OF THIS PARITY ARE NEGATIVE. C IT MAY ALSO BE SET NEGATIVE AT END (750). C IF(NEGER.GE.2) EREST(JNS)=-EREST(JNS) IF(NEGER.GE.2)GOTO 690 IF(TERM.LT.RANMIN) EREST(JNS)=-EREST(JNS) IF(-TERM.GE.RANMIN) EREST(JNS)=-EREST(JNS) IF(EREST(JNS).LT.ZERO) NEGER=NEGER+1 IF(EREST(JNS).GE.ZERO) NEGER=0 690 CONTINUE C HFACT=HSQ*HFACT FACT=FACT*(XJNS+ONE) 700 CONTINUE 701 CONTINUE IF(NPAR.GT.0)GOTO 40 C C SET SIGN OF EREST NEGATIVE IF TWO PREVIOUS CONSECUTIVE ERESTS C ARE NEGATIVE. C IF(NDPAR.NE.0)GOTO 900 IF(NDERP.LE.2) GOTO 900 NEGER=0 ERPREV=EREST(1) DO 760 J= 2,NDERP ERNOW=EREST(J) IF(NEGER.EQ.2) GOTO 740 IF(ERPREV.GE.ZERO) GOTO 750 IF(ERNOW.GE.ZERO) GOTO 750 740 NEGER=2 IF(ERNOW.GT.ZERO) EREST(J)=-ERNOW 750 ERPREV=ERNOW 760 CONTINUE C 900 CONTINUE RETURN C*** *** *** END OF ENDACE *** *** *** C*********************************************************** END SUBROUTINE CHPEND(XVAL,HBASE,DER,FUN) C C ********************************************************************** C C ***** SUBROUTINE CHPEND ***** C C THIS IS A MODIFIED VERSION OF ENDACE. IT COMPUTES A SET OF C APPROXIMATIONS TO THE J-TH DERIVATIVE (J=1,2,..,5) OF FUN(X) C AT X=XVAL. C NO ERROR ESTIMATE IS RETURNED. C C INPUT - XVAL,HBASE,FUN C C XVAL -REAL- THE POINT AT WHICH THE DERIVATIVES C ARE REQUIRED. C C HBASE -REAL- THE STEP LENGHT C C FUN -REAL FUNCTION- THE FUNCTION BEING DIFFEREN C TIATED. THE ROUTINE REQUIRES 9 EVALUATUIONS OF ..FUN.. C IN THE INTERVAL (XVAL-7*HBASE,XVAL+7*HBASE). C C C OUTPUT - DER C C DER -REAL(5)- ARRAY OF FUNCTION DERIVATIVES.DER(J) C CONTAINS THE J-TH DERIVATIVE OF FUN. C C C C REFERENCE - J.N.LYNESS - SOFTWARE FOR NUMERICAL DIFFERENTIATION C*********************************************************************** C IMPLICIT REAL (A-H,O-Z) DIMENSION DER(5) DIMENSION RTAB(4,3),ACOF(4),TZEROM(4),RANGE(3) C EXTERNAL FUN * DATA ZERO,ONE,TWO/0.0E0,1.0E0,2.0E0/ C C FZ = FUN(XVAL) DO 10 N = 1,4 NN1 = 2*N - 1 NSQ = NN1*NN1 ACOF(N) = REAL(NSQ) XNUM = REAL(NN1) H = HBASE*XNUM TEMP = FUN(XVAL+H) TERM = FUN(XVAL-H) RTAB(N,1) = (TEMP-TERM)/ (TWO*H) TZEROM(N) = (TEMP-FZ-FZ+TERM)/ (TWO*H**2) 10 CONTINUE NPAR = 1 NSMAX = 3 GO TO 40 C 20 CONTINUE NSMAX = 2 NPAR = -1 DO 30 J = 1,4 RTAB(J,1) = TZEROM(J) 30 CONTINUE 40 CONTINUE HSQ = HBASE*HBASE FACT = ONE HFACT = ONE DO 150 NS = 1,NSMAX NPMIN = NS IF (NS.EQ.1) NPMIN = 2 DO 60 NP = NPMIN,3 NKTOP = 5 - NP NPR = NP - NS + 1 DO 50 NK = 1,NKTOP J = NP + NK - 1 A = ACOF(J) B = ACOF(NK) TERM = ZERO TEMP = ZERO IF (NP.NE.NS) TEMP = A*RTAB(NK,NPR-1) - + B*RTAB(NK+1,NPR-1) IF (NS.NE.1) TERM = -RTAB(NK,NPR) + RTAB(NK+1,NPR) RTAB(NK,NPR) = (TERM+TEMP)/ (A-B) 50 CONTINUE 60 CONTINUE C DO 110 NP = NS,3 NTOP = 5 - NP NPR = NP - NS + 1 XMAX = RTAB(1,NPR) XMIN = RTAB(1,NPR) NMAX = 1 NMIN = 1 DO 70 NK = 2,NTOP TEMP = RTAB(NK,NPR) IF (TEMP.GT.XMAX) NMAX = NK IF (TEMP.GT.XMAX) XMAX = TEMP IF (TEMP.LT.XMIN) NMIN = NK IF (TEMP.LT.XMIN) XMIN = TEMP 70 CONTINUE RANGE(NP) = XMAX - XMIN IF (NP.NE.NS) GO TO 80 RANMIN = RANGE(NP) GO TO 90 C 80 CONTINUE IF (RANGE(NP).GE.RANMIN) GO TO 100 RANMIN = RANGE(NP) 90 CONTINUE NPMIN = NP NLO = NMIN 100 CONTINUE 110 CONTINUE C TERM = ZERO NTOP = 5 - NPMIN J = NPMIN - NS + 1 DO 120 NK = 1,NTOP IF (NK.EQ.NLO) GO TO 120 TERM = TERM + RTAB(NK,J) 120 CONTINUE TERM = TERM/ (NTOP-1) JNS = 2*NS - 1 IF (NPAR.LT.0) JNS = 2*NS XJNS = JNS FACT = FACT*XJNS DER(JNS) = TERM*FACT/HFACT 140 CONTINUE HFACT = HSQ*HFACT FACT = FACT* (XJNS+ONE) 150 CONTINUE IF (NPAR.EQ.1) GO TO 20 C RETURN C ******************************************************************* C * * C * END OF SUBROUTINE ****** CHPEND ****** * C * * C ******************************************************************* END SUBROUTINE ONESID(NCODE,XVAL,HBASE,DER,FUN) C C ********************************************************************** C C ***** SUBROUTINE ONESID ***** C C THIS IS A MODIFIED VERSION OF ENDACE. IT COMPUTES A SET OF C APPROXIMATIONS TO THE J-TH DERIVATIVE (J=1,2,..,5) OF FUN(X) C AT X=XVAL. C ONE SIDED DIFFERENCE APPROXIMATIONS ARE USED. C NO ERROR ESTIMATE IS RETURNED. C C INPUT - NCODE,XVAL,HBASE,FUN C C NCODE -INTEGER- IF=1 SELECT RIGHT ONE SIDED DIFFERENCE C IF=-1 SELECT LEFT ONE SIDED DIFFERENCE C C XVAL -REAL- THE POINT AT WHICH THE DERIVATIVES C ARE REQUIRED. C C HBASE -REAL- THE STEP LENGHT C C FUN -REAL FUNCTION- THE FUNCTION BEING DIFFEREN C TIATED. THE ROUTINE REQUIRES 7 EVALUATIONS OF ..FUN.. C IN THE INTERVAL (XVAL,XVAL+6*HBASE) IF NCODE=1 C (XVAL-6*HBASE,XVAL) IF NCODE=-1. C C C OUTPUT - DER C C DER -REAL(5)- ARRAY OF FUNCTION DERIVATIVES.DER(J) C CONTAINS THE J-TH DERIVATIVE OF FUN. C C C C REFERENCE - J.N.LYNESS - SOFTWARE FOR NUMERICAL DIFFERENTIATION C*********************************************************************** C IMPLICIT REAL (A-H,O-Z) DIMENSION DER(5) DIMENSION RTAB(7,6),ACOF(7),RANGE(7) C EXTERNAL FUN * DATA ZERO,ONE/0.0E0,1.0E0/ C C RTAB(1,1)=FUN(XVAL) ACOF(1)=ZERO DO 10 N = 2,7 ACOF(N) = REAL(N-1) H = HBASE*ACOF(N) IF(NCODE.EQ.-1)TEMP=FUN(XVAL-H) IF(NCODE.EQ.1)TEMP=FUN(XVAL+H) RTAB(N,1)=TEMP 10 CONTINUE NSMAX=6 FACT = ONE DO 150 NS = 1,NSMAX NPMIN = NS IF (NS.EQ.1) NPMIN = 2 DO 60 NP = NPMIN,6 NKTOP = 8 - NP NPR = NP - NS + 1 DO 50 NK = 1,NKTOP J = NP + NK - 1 A = ACOF(J)*HBASE B = ACOF(NK)*HBASE IF(NCODE.EQ.-1)A=-A IF(NCODE.EQ.-1)B=-B TERM = ZERO TEMP = ZERO IF (NP.NE.NS) TEMP = A*RTAB(NK,NPR-1) - + B*RTAB(NK+1,NPR-1) IF (NS.NE.1) TERM = -RTAB(NK,NPR) + RTAB(NK+1,NPR) RTAB(NK,NPR) = (TERM+TEMP)/ (A-B) 50 CONTINUE 60 CONTINUE DO 110 NP = NS,6 NTOP = 8 - NP NPR = NP - NS + 1 XMAX = RTAB(1,NPR) XMIN = RTAB(1,NPR) NMAX = 1 NMIN = 1 DO 70 NK = 2,NTOP TEMP = RTAB(NK,NPR) IF (TEMP.GT.XMAX) NMAX = NK IF (TEMP.GT.XMAX) XMAX = TEMP IF (TEMP.LT.XMIN) NMIN = NK IF (TEMP.LT.XMIN) XMIN = TEMP 70 CONTINUE RANGE(NP) = XMAX - XMIN IF (NP.NE.NS) GO TO 80 RANMIN = RANGE(NP) GO TO 90 C 80 CONTINUE IF (RANGE(NP).GE.RANMIN) GO TO 100 RANMIN = RANGE(NP) 90 CONTINUE NPMIN = NP NLO = NMIN 100 CONTINUE 110 CONTINUE C TERM = ZERO NTOP = 8 - NPMIN J = NPMIN - NS + 1 DO 120 NK = 1,NTOP IF (NK.EQ.NLO) GO TO 120 TERM = TERM + RTAB(NK,J) 120 CONTINUE TERM = TERM/ (NTOP-1) JNS=NS XJNS = MAX(JNS-1,1) FACT = FACT*XJNS IF(JNS.EQ.1) GOTO 150 DER(JNS-1)=TERM*FACT 150 CONTINUE C RETURN C ******************************************************************* C * * C * END OF SUBROUTINE ****** ONESIDE ****** * C * * C ******************************************************************* END SUBROUTINE ONESMP(NCODE,FUNVAL,IXVAL,HBASE,DER) C C ********************************************************************** C C ***** SUBROUTINE ONESMP ***** C C THIS IS A MODIFIED VERSION OF ENDACE. IT COMPUTES A SET OF C APPROXIMATIONS TO THE J-TH DERIVATIVE (J=1,2,..,5) OF FUN(X) C AT X=XVAL. C IT REQUIRES A SET OF FUNCTION VALUES STORED IN AN INPUT ARRAY. C ONE SIDED DIFFERENCE APPROXIMATIONS ARE USED. C NO ERROR ESTIMATE IS RETURNED. C C INPUT - NCODE,FUNVAL,IXVAL,HBASE C C NCODE -INTEGER- IF=1 SELECT RIGHT ONE SIDED DIFFERENCE C IF=-1 SELECT LEFT ONE SIDED DIFFERENCE C C FUNVAL -REAL(*)- ARRAY OF THE VALUES OF THE C FUNCTION TO BE DIFFERENTIATED. C FUNCTION VALUES MUST BE AT EQUALLY SPACED POINTS C WITH STEP ..HBASE.. C C IXVAL -INTEGER- THE POINT WHERE THE DERIVATIVES ARE C REQUIRED IS STORED IN FUNVAL(IXVAL). C C HBASE -REAL- THE STEP LENGHT C C C OUTPUT - DER C C DER -REAL(5)- ARRAY OF FUNCTION DERIVATIVES.DER(J) C CONTAINS THE J-TH DERIVATIVE OF FUN. C C C C REFERENCE - J.N.LYNESS - SOFTWARE FOR NUMERICAL DIFFERENTIATION C*********************************************************************** C IMPLICIT REAL (A-H,O-Z) DIMENSION DER(5),FUNVAL(*) DIMENSION RTAB(7,6),ACOF(7),RANGE(7) C EXTERNAL FUN * DATA ZERO,ONE/0.0E0,1.0E0/ C C RTAB(1,1)=FUNVAL(IXVAL) ACOF(1)=ZERO DO 10 N = 2,7 ACOF(N) = REAL(N-1) H = HBASE*ACOF(N) IF(NCODE.EQ.-1)TEMP=FUNVAL(IXVAL-N+1) IF(NCODE.EQ.1)TEMP=FUNVAL(IXVAL+N-1) RTAB(N,1)=TEMP 10 CONTINUE NSMAX=6 FACT = ONE DO 150 NS = 1,NSMAX NPMIN = NS IF (NS.EQ.1) NPMIN = 2 DO 60 NP = NPMIN,6 NKTOP = 8 - NP NPR = NP - NS + 1 DO 50 NK = 1,NKTOP J = NP + NK - 1 A = ACOF(J)*HBASE B = ACOF(NK)*HBASE IF(NCODE.EQ.-1)A=-A IF(NCODE.EQ.-1)B=-B TERM = ZERO TEMP = ZERO IF (NP.NE.NS) TEMP = A*RTAB(NK,NPR-1) - + B*RTAB(NK+1,NPR-1) IF (NS.NE.1) TERM = -RTAB(NK,NPR) + RTAB(NK+1,NPR) RTAB(NK,NPR) = (TERM+TEMP)/ (A-B) 50 CONTINUE 60 CONTINUE DO 110 NP = NS,6 NTOP = 8 - NP NPR = NP - NS + 1 XMAX = RTAB(1,NPR) XMIN = RTAB(1,NPR) NMAX = 1 NMIN = 1 DO 70 NK = 2,NTOP TEMP = RTAB(NK,NPR) IF (TEMP.GT.XMAX) NMAX = NK IF (TEMP.GT.XMAX) XMAX = TEMP IF (TEMP.LT.XMIN) NMIN = NK IF (TEMP.LT.XMIN) XMIN = TEMP 70 CONTINUE RANGE(NP) = XMAX - XMIN IF (NP.NE.NS) GO TO 80 RANMIN = RANGE(NP) GO TO 90 C 80 CONTINUE IF (RANGE(NP).GE.RANMIN) GO TO 100 RANMIN = RANGE(NP) 90 CONTINUE NPMIN = NP NLO = NMIN 100 CONTINUE 110 CONTINUE C TERM = ZERO NTOP = 8 - NPMIN J = NPMIN - NS + 1 DO 120 NK = 1,NTOP IF (NK.EQ.NLO) GO TO 120 TERM = TERM + RTAB(NK,J) 120 CONTINUE TERM = TERM/ (NTOP-1) JNS=NS XJNS = MAX(JNS-1,1) FACT = FACT*XJNS IF(JNS.EQ.1) GOTO 150 DER(JNS-1)=TERM*FACT 150 CONTINUE C RETURN C ******************************************************************* C * * C * END OF SUBROUTINE ****** ONESMP ****** * C * * C ******************************************************************* END C**************************************************************** C C COLLECTION OF FISHPACK ROUTINES FOR REAL FORWARD FFT C SINGLE PRECISION VERSION C**************************************************************** C THESE ARE A COPY OF THE FOLLOWING 9 SUBROUTINES C AS THEY APPEAR IN C P.N.SWARZTRAUBER - FISHPACK, A PACKAGE OF FORTRAN SUBPROGRAMS C FOR THE FAST FOURIER TRANSFORM OF PERIODIC AND OTHER C SYMMETRIC SEQUENCES. IN PORTLIB,(1979) C ONLY CHANGES TO IMPROVE PORTABILITY OF THE CODE HAVE BEEN MADE. C C RFFTI,RFFTI1,RFFTF,RFFTF1,RADF2,RADF3,RADF4,RADF5,RADFG. C**************************************************************** SUBROUTINE RFFTI (N,WSAVE) IMPLICIT REAL (A-H,O-Z) DIMENSION WSAVE(*) IF (N .EQ. 1) RETURN CALL RFFTI1 (N,WSAVE(N+1),WSAVE(2*N+1)) RETURN END SUBROUTINE RFFTI1 (N,WA,IFAC) IMPLICIT REAL (A-H,O-Z) REAL IFAC DIMENSION WA(*) ,IFAC(*) ,NTRYH(4) DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/ DATA TPI/0.E0/ IF(TPI.EQ.0.E0)TPI=8.E0*ATAN(1.E0) NL = N NF = 0 J = 0 101 J = J+1 IF (J-4) 102,102,103 102 NTRY = NTRYH(J) GO TO 104 103 NTRY = NTRY+2 104 NQ = NL/NTRY NR = NL-NTRY*NQ IF (NR) 101,105,101 105 NF = NF+1 IFAC(NF+2) = NTRY NL = NQ IF (NTRY .NE. 2) GO TO 107 IF (NF .EQ. 1) GO TO 107 DO 106 I=2,NF IB = NF-I+2 IFAC(IB+2) = IFAC(IB+1) 106 CONTINUE IFAC(3) = 2 107 IF (NL .NE. 1) GO TO 104 IFAC(1) = N IFAC(2) = NF ARGH = TPI/REAL(N) IS = 0 NFM1 = NF-1 L1 = 1 IF (NFM1 .EQ. 0) RETURN DO 110 K1=1,NFM1 IP = IFAC(K1+2) LD = 0 L2 = L1*IP IDO = N/L2 IPM = IP-1 DO 109 J=1,IPM LD = LD+L1 I = IS ARGLD = REAL(LD)*ARGH FI = 0. DO 108 II=3,IDO,2 I = I+2 FI = FI+1. ARG = FI*ARGLD WA(I-1) = COS(ARG) WA(I) = SIN(ARG) 108 CONTINUE IS = IS+IDO 109 CONTINUE L1 = L2 110 CONTINUE RETURN END SUBROUTINE RFFTF (N,R,WSAVE) IMPLICIT REAL (A-H,O-Z) DIMENSION R(*) ,WSAVE(*) IF (N .EQ. 1) RETURN CALL RFFTF1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1)) RETURN END SUBROUTINE RFFTF1 (N,C,CH,WA,IFAC) IMPLICIT REAL (A-H,O-Z) REAL IFAC DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*) NF = IFAC(2) NA = 1 L2 = N IW = N DO 111 K1=1,NF KH = NF-K1 IP = IFAC(KH+3) L1 = L2/IP IDO = N/L2 IDL1 = IDO*L1 IW = IW-(IP-1)*IDO NA = 1-NA IF (IP .NE. 4) GO TO 102 IX2 = IW+IDO IX3 = IX2+IDO IF (NA .NE. 0) GO TO 101 CALL RADF4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3)) GO TO 110 101 CALL RADF4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3)) GO TO 110 102 IF (IP .NE. 2) GO TO 104 IF (NA .NE. 0) GO TO 103 CALL RADF2 (IDO,L1,C,CH,WA(IW)) GO TO 110 103 CALL RADF2 (IDO,L1,CH,C,WA(IW)) GO TO 110 104 IF (IP .NE. 3) GO TO 106 IX2 = IW+IDO IF (NA .NE. 0) GO TO 105 CALL RADF3 (IDO,L1,C,CH,WA(IW),WA(IX2)) GO TO 110 105 CALL RADF3 (IDO,L1,CH,C,WA(IW),WA(IX2)) GO TO 110 106 IF (IP .NE. 5) GO TO 108 IX2 = IW+IDO IX3 = IX2+IDO IX4 = IX3+IDO IF (NA .NE. 0) GO TO 107 CALL RADF5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4)) GO TO 110 107 CALL RADF5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4)) GO TO 110 108 IF (IDO .EQ. 1) NA = 1-NA IF (NA .NE. 0) GO TO 109 CALL RADFG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW)) NA = 1 GO TO 110 109 CALL RADFG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW)) NA = 0 110 L2 = L1 111 CONTINUE IF (NA .EQ. 1) RETURN DO 112 I=1,N C(I) = CH(I) 112 CONTINUE RETURN END SUBROUTINE RADF2 (IDO,L1,CC,CH,WA1) IMPLICIT REAL (A-H,O-Z) DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) , 1 WA1(*) DO 101 K=1,L1 CH(1,1,K) = CC(1,K,1)+CC(1,K,2) CH(IDO,2,K) = CC(1,K,1)-CC(1,K,2) 101 CONTINUE IF (IDO-2) 107,105,102 102 IDP2 = IDO+2 DO 104 K=1,L1 DO 103 I=3,IDO,2 IC = IDP2-I TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) CH(I,1,K) = CC(I,K,1)+TI2 CH(IC,2,K) = TI2-CC(I,K,1) CH(I-1,1,K) = CC(I-1,K,1)+TR2 CH(IC-1,2,K) = CC(I-1,K,1)-TR2 103 CONTINUE 104 CONTINUE IF (MOD(IDO,2) .EQ. 1) RETURN 105 DO 106 K=1,L1 CH(1,2,K) = -CC(IDO,K,2) CH(IDO,1,K) = CC(IDO,K,1) 106 CONTINUE 107 RETURN END SUBROUTINE RADF3 (IDO,L1,CC,CH,WA1,WA2) IMPLICIT REAL (A-H,O-Z) DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) , 1 WA1(*) ,WA2(*) DATA TAUR,TAUI /-.5E0,0.E0/ IF(TAUI.EQ.0.E0)TAUI=SQRT(0.75E0) DO 101 K=1,L1 CR2 = CC(1,K,2)+CC(1,K,3) CH(1,1,K) = CC(1,K,1)+CR2 CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2)) CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2 101 CONTINUE IF (IDO .EQ. 1) RETURN IDP2 = IDO+2 DO 103 K=1,L1 DO 102 I=3,IDO,2 IC = IDP2-I DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) CR2 = DR2+DR3 CI2 = DI2+DI3 CH(I-1,1,K) = CC(I-1,K,1)+CR2 CH(I,1,K) = CC(I,K,1)+CI2 TR2 = CC(I-1,K,1)+TAUR*CR2 TI2 = CC(I,K,1)+TAUR*CI2 TR3 = TAUI*(DI2-DI3) TI3 = TAUI*(DR3-DR2) CH(I-1,3,K) = TR2+TR3 CH(IC-1,2,K) = TR2-TR3 CH(I,3,K) = TI2+TI3 CH(IC,2,K) = TI3-TI2 102 CONTINUE 103 CONTINUE RETURN END SUBROUTINE RADF4 (IDO,L1,CC,CH,WA1,WA2,WA3) IMPLICIT REAL (A-H,O-Z) DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) , 1 WA1(*) ,WA2(*) ,WA3(*) DATA HSQT2 /0.E0/ IF(HSQT2.EQ.0.E0)HSQT2=SQRT(0.5E0) DO 101 K=1,L1 TR1 = CC(1,K,2)+CC(1,K,4) TR2 = CC(1,K,1)+CC(1,K,3) CH(1,1,K) = TR1+TR2 CH(IDO,4,K) = TR2-TR1 CH(IDO,2,K) = CC(1,K,1)-CC(1,K,3) CH(1,3,K) = CC(1,K,4)-CC(1,K,2) 101 CONTINUE IF (IDO-2) 107,105,102 102 IDP2 = IDO+2 DO 104 K=1,L1 DO 103 I=3,IDO,2 IC = IDP2-I CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4) CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4) TR1 = CR2+CR4 TR4 = CR4-CR2 TI1 = CI2+CI4 TI4 = CI2-CI4 TI2 = CC(I,K,1)+CI3 TI3 = CC(I,K,1)-CI3 TR2 = CC(I-1,K,1)+CR3 TR3 = CC(I-1,K,1)-CR3 CH(I-1,1,K) = TR1+TR2 CH(IC-1,4,K) = TR2-TR1 CH(I,1,K) = TI1+TI2 CH(IC,4,K) = TI1-TI2 CH(I-1,3,K) = TI4+TR3 CH(IC-1,2,K) = TR3-TI4 CH(I,3,K) = TR4+TI3 CH(IC,2,K) = TR4-TI3 103 CONTINUE 104 CONTINUE IF (MOD(IDO,2) .EQ. 1) RETURN 105 CONTINUE DO 106 K=1,L1 TI1 = -HSQT2*(CC(IDO,K,2)+CC(IDO,K,4)) TR1 = HSQT2*(CC(IDO,K,2)-CC(IDO,K,4)) CH(IDO,1,K) = TR1+CC(IDO,K,1) CH(IDO,3,K) = CC(IDO,K,1)-TR1 CH(1,2,K) = TI1-CC(IDO,K,3) CH(1,4,K) = TI1+CC(IDO,K,3) 106 CONTINUE 107 RETURN END SUBROUTINE RADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) IMPLICIT REAL (A-H,O-Z) DIMENSION CC(IDO,L1,5) ,CH(IDO,5,L1) , 1 WA1(*) ,WA2(*) ,WA3(*) ,WA4(*) DATA TR11,TI11,TR12,TI12/0.E0,0.E0,0.E0,0.E0/ DATA PI/0.E0/ IF(PI.EQ.0.E0)PI=4.E0*ATAN(1.E0) IF(TR11.EQ.0.E0)TR11=COS(2.E0/5.E0*PI) IF(TI11.EQ.0.E0)TI11=SIN(2.E0/5.E0*PI) IF(TR12.EQ.0.E0)TR12=-SIN(3.E0/10.E0*PI) IF(TI12.EQ.0.E0)TI12=COS(3.E0/10.E0*PI) DO 101 K=1,L1 CR2 = CC(1,K,5)+CC(1,K,2) CI5 = CC(1,K,5)-CC(1,K,2) CR3 = CC(1,K,4)+CC(1,K,3) CI4 = CC(1,K,4)-CC(1,K,3) CH(1,1,K) = CC(1,K,1)+CR2+CR3 CH(IDO,2,K) = CC(1,K,1)+TR11*CR2+TR12*CR3 CH(1,3,K) = TI11*CI5+TI12*CI4 CH(IDO,4,K) = CC(1,K,1)+TR12*CR2+TR11*CR3 CH(1,5,K) = TI12*CI5-TI11*CI4 101 CONTINUE IF (IDO .EQ. 1) RETURN IDP2 = IDO+2 DO 103 K=1,L1 DO 102 I=3,IDO,2 IC = IDP2-I DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) DR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4) DI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4) DR5 = WA4(I-2)*CC(I-1,K,5)+WA4(I-1)*CC(I,K,5) DI5 = WA4(I-2)*CC(I,K,5)-WA4(I-1)*CC(I-1,K,5) CR2 = DR2+DR5 CI5 = DR5-DR2 CR5 = DI2-DI5 CI2 = DI2+DI5 CR3 = DR3+DR4 CI4 = DR4-DR3 CR4 = DI3-DI4 CI3 = DI3+DI4 CH(I-1,1,K) = CC(I-1,K,1)+CR2+CR3 CH(I,1,K) = CC(I,K,1)+CI2+CI3 TR2 = CC(I-1,K,1)+TR11*CR2+TR12*CR3 TI2 = CC(I,K,1)+TR11*CI2+TR12*CI3 TR3 = CC(I-1,K,1)+TR12*CR2+TR11*CR3 TI3 = CC(I,K,1)+TR12*CI2+TR11*CI3 TR5 = TI11*CR5+TI12*CR4 TI5 = TI11*CI5+TI12*CI4 TR4 = TI12*CR5-TI11*CR4 TI4 = TI12*CI5-TI11*CI4 CH(I-1,3,K) = TR2+TR5 CH(IC-1,2,K) = TR2-TR5 CH(I,3,K) = TI2+TI5 CH(IC,2,K) = TI5-TI2 CH(I-1,5,K) = TR3+TR4 CH(IC-1,4,K) = TR3-TR4 CH(I,5,K) = TI3+TI4 CH(IC,4,K) = TI4-TI3 102 CONTINUE 103 CONTINUE RETURN END SUBROUTINE RADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) IMPLICIT REAL (A-H,O-Z) DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 1 C1(IDO,L1,IP) ,C2(IDL1,IP), 2 CH2(IDL1,IP) ,WA(*) DATA TPI/0.E0/ IF(TPI.EQ.0.E0)TPI=8.E0*ATAN(1.E0) ARG = TPI/REAL(IP) DCP = COS(ARG) DSP = SIN(ARG) IPPH = (IP+1)/2 IPP2 = IP+2 IDP2 = IDO+2 NBD = (IDO-1)/2 IF (IDO .EQ. 1) GO TO 119 DO 101 IK=1,IDL1 CH2(IK,1) = C2(IK,1) 101 CONTINUE DO 103 J=2,IP DO 102 K=1,L1 CH(1,K,J) = C1(1,K,J) 102 CONTINUE 103 CONTINUE IF (NBD .GT. L1) GO TO 107 IS = -IDO DO 106 J=2,IP IS = IS+IDO IDIJ = IS DO 105 I=3,IDO,2 IDIJ = IDIJ+2 DO 104 K=1,L1 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J) CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J) 104 CONTINUE 105 CONTINUE 106 CONTINUE GO TO 111 107 IS = -IDO DO 110 J=2,IP IS = IS+IDO DO 109 K=1,L1 IDIJ = IS DO 108 I=3,IDO,2 IDIJ = IDIJ+2 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J) CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J) 108 CONTINUE 109 CONTINUE 110 CONTINUE 111 IF (NBD .LT. L1) GO TO 115 DO 114 J=2,IPPH JC = IPP2-J DO 113 K=1,L1 DO 112 I=3,IDO,2 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC) C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC) C1(I,K,J) = CH(I,K,J)+CH(I,K,JC) C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J) 112 CONTINUE 113 CONTINUE 114 CONTINUE GO TO 121 115 DO 118 J=2,IPPH JC = IPP2-J DO 117 I=3,IDO,2 DO 116 K=1,L1 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC) C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC) C1(I,K,J) = CH(I,K,J)+CH(I,K,JC) C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J) 116 CONTINUE 117 CONTINUE 118 CONTINUE GO TO 121 119 DO 120 IK=1,IDL1 C2(IK,1) = CH2(IK,1) 120 CONTINUE 121 DO 123 J=2,IPPH JC = IPP2-J DO 122 K=1,L1 C1(1,K,J) = CH(1,K,J)+CH(1,K,JC) C1(1,K,JC) = CH(1,K,JC)-CH(1,K,J) 122 CONTINUE 123 CONTINUE C AR1 = 1. AI1 = 0. DO 127 L=2,IPPH LC = IPP2-L AR1H = DCP*AR1-DSP*AI1 AI1 = DCP*AI1+DSP*AR1 AR1 = AR1H DO 124 IK=1,IDL1 CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2) CH2(IK,LC) = AI1*C2(IK,IP) 124 CONTINUE DC2 = AR1 DS2 = AI1 AR2 = AR1 AI2 = AI1 DO 126 J=3,IPPH JC = IPP2-J AR2H = DC2*AR2-DS2*AI2 AI2 = DC2*AI2+DS2*AR2 AR2 = AR2H DO 125 IK=1,IDL1 CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J) CH2(IK,LC) = CH2(IK,LC)+AI2*C2(IK,JC) 125 CONTINUE 126 CONTINUE 127 CONTINUE DO 129 J=2,IPPH DO 128 IK=1,IDL1 CH2(IK,1) = CH2(IK,1)+C2(IK,J) 128 CONTINUE 129 CONTINUE C IF (IDO .LT. L1) GO TO 132 DO 131 K=1,L1 DO 130 I=1,IDO CC(I,1,K) = CH(I,K,1) 130 CONTINUE 131 CONTINUE GO TO 135 132 DO 134 I=1,IDO DO 133 K=1,L1 CC(I,1,K) = CH(I,K,1) 133 CONTINUE 134 CONTINUE 135 DO 137 J=2,IPPH JC = IPP2-J J2 = J+J DO 136 K=1,L1 CC(IDO,J2-2,K) = CH(1,K,J) CC(1,J2-1,K) = CH(1,K,JC) 136 CONTINUE 137 CONTINUE IF (IDO .EQ. 1) RETURN IF (NBD .LT. L1) GO TO 141 DO 140 J=2,IPPH JC = IPP2-J J2 = J+J DO 139 K=1,L1 DO 138 I=3,IDO,2 IC = IDP2-I CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC) CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC) CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC) CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J) 138 CONTINUE 139 CONTINUE 140 CONTINUE RETURN 141 DO 144 J=2,IPPH JC = IPP2-J J2 = J+J DO 143 I=3,IDO,2 IC = IDP2-I DO 142 K=1,L1 CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC) CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC) CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC) CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J) 142 CONTINUE 143 CONTINUE 144 CONTINUE RETURN END C*********************************************************************** C C ******* END OF THE COLLECTION OF FISHPACK ROUTINES ****** C C*********************************************************************** PROGRAM DRIVER C********************************************************************* C C SIMPLE DRIVER FOR THE PACKAGE FOURCO C SINGLE PRECISION VERSION C THIS DRIVER TESTS ALL THE FEATURES OF THE PACKAGE FOURCO C THE OUTPUT ARE THE VALUES OF THE FIRST 100 COSINE/SINE C FOURIER COEFFICIENTS ON THE INTERVAL (0,1) OF THE FUNCTION C F(X) = EXP(3.*X) C COMPUTED FIRST USING THE TOP END FCOEFS WITH 3 DIFFERENT C NUMERICAL DIFFERENTIATION ROUTINES,THE RESTART FEATURE AND C PASSING PRECALCULATED EXACT DERIVATIVES, AND THEN USING THE C TOP END FCOETD. C C********************************************************************* C C THE DRIVER STARTS HERE C IMPLICIT REAL (A-H,O-Z) DIMENSION C(100),S(100),CEXACT(20),SEXACT(20), $ WORK(5000),XLAMB(20),FUNVAL(65) EXTERNAL FUN C C THIS IS THE VALUE OF PI C PI=4.0E0*ATAN(1.0E0) C C SET THE ENDPOINTS OF THE INTERVAL C A=0.E0 B=1.E0 C C SET THE NUMBER OF REQUESTED COEFFICIENTS C NCOEF=100 C C SET THE REQUESTED ABSOLUTE ACCURACY C TOL=1.E-6 C C SET THE PHYSICAL LIMIT PARAMETER C MTOP=32 PRINT 1 1 FORMAT(8X,'--------------------------------------------'/ $ 4X,'THIS RUN TESTS ALL THE FEATURES OF THE PACKAGE FOURCO :'/ $ 4X,'.TOP-END FCOEFS WITH 3 NUM.DIFFERENTIATION ROUTINES,'/ $ 5X,'THE RESTART FEATURE AND PRECALCULATED DERIVATIVES'/ $ 4X,'.TOP-END FCOETD FOR SAMPLED FUNCTION VALUES'/ $ 8X,'--------------------------------------------'// $ 5X,' TEST FUNCTION IS EXP(3.E0*X)'/) C C NOW CALL THE TOP-END FCOEFS USING ENDACE(ACCURATE NUMERICAL C DIFFERENTIATION ROUTINE) C IDER=1 IFIRST=0 C CALL FCOEFS(A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,C,S,IERR,ESTER,MACT) C PRINT 5,IDER,IFIRST,MTOP,TOL 5 FORMAT(1X//4X,'...........................................'/ $ 4X,'THIS IS OUTPUT FROM PACKAGE FOURCO USING TOP-END'/ $ 4X,'FCOEFS WITH ENDACE(ACCURATE NUM.DIFFERENT. ROUTINE)'/ $ 4X,'.........................................'/ $ 14X,'IDER=',I3,' IFIRST=',I3/ $ 14X,'MTOP=',I4,' TOL=',E12.2/) C PRINT THE ESTIMATED ABSOLUTE ERROR, THE ERROR PARAMETER AND C THE VALUE OF MACT C PRINT 10,ESTER,IERR,MACT 10 FORMAT(20X, $ ' ESTIMATED ERROR',5X,E10.3/ $ ' ERROR PARAMETER',10X,I2/ $ ' MACT ',8X,I4) C C PRINT THE VALUE OF THE APPROXIMATED COSINE COEFFICIENTS COMPUTED C BY FCOEFS , THE VALUES OF THE EXACT COEFFICIENTS C AND THE ABSOLUTE ERRORS C PRINT 20 20 FORMAT(1X///8X,' COSINE FOURIER COEFFICIENTS'// $ ' I-1',6X,'COMPUTED',5X,'EXACT',4X, $ 'COMPUTED-EXACT'/) DO 100 I=1,NCOEF,5 J=I/5 +1 CEXACT(J)=(EXP(3.0E0)-1.0E0)*3.0E0/(9.0E0+4.0E0*PI**2*(I-1)**2) C C REMEMBER OUR DEFINITION OF FOURIER COEFFICIENTS C IF (J.NE.1)CEXACT(J)=2.0E0*CEXACT(J) DIFFS=C(I)-CEXACT(J) PRINT 30,I,C(I),CEXACT(J),DIFFS 30 FORMAT(2X,I4,3(2X,E10.3)) 100 CONTINUE C C AND NOW THE SINE COEFFICIENTS ... C PRINT 40 40 FORMAT(1X///8X,' SINE FOURIER COEFFICIENTS'// $ ' I-1',6X,'COMPUTED',5X,'EXACT',4X, $ 'COMPUTED-EXACT'/) DO 200 I=1,NCOEF,5 J=I/5 +1 SEXACT(J)=-(EXP(3.0E0)-1.0E0)*2.0E0*PI*(I-1) $ /(9.0E0+4.0E0*PI**2*(I-1)**2) SEXACT(J)=2.0E0*SEXACT(J) DIFFS=S(I)-SEXACT(J) PRINT 50,I,S(I),SEXACT(J),DIFFS 50 FORMAT(2X,I4,3(2X,E10.3)) 200 CONTINUE C C NOW CALL THE TOP-END FCOEFS USING CHPEND(CHEAP NUMERICAL C DIFFERENTIATION ROUTINE) C IDER=0 IFIRST=0 C CALL FCOEFS(A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,C,S,IERR,ESTER,MACT) C PRINT 45,IDER,IFIRST,MTOP,TOL 45 FORMAT(1X//4X,'...........................................'/ $ 4X,'THIS IS OUTPUT FROM PACKAGE FOURCO USING TOP-END'/ $ 4X,'FCOEFS WITH CHPEND(CHEAP NUM.DIFFERENT. ROUTINE)'/ $ 4X,'.........................................'/ $ 14X,'IDER=',I3,' IFIRST=',I3/ $ 14X,'MTOP=',I4,' TOL=',E12.2/) C PRINT THE ESTIMATED ABSOLUTE ERROR, THE ERROR PARAMETER AND C THE VALUE OF MACT C PRINT 10,ESTER,IERR,MACT C C PRINT THE VALUE OF THE APPROXIMATED COSINE COEFFICIENTS COMPUTED C BY FCOEFS , THE VALUES OF THE EXACT COEFFICIENTS C AND THE ABSOLUTE ERRORS C PRINT 20 DO 101 I=1,NCOEF,5 J=I/5 +1 DIFFS=C(I)-CEXACT(J) PRINT 30,I,C(I),CEXACT(J),DIFFS 101 CONTINUE C C AND NOW THE SINE COEFFICIENTS ... C PRINT 40 DO 201 I=1,NCOEF,5 J=I/5 +1 DIFFS=S(I)-SEXACT(J) PRINT 50,I,S(I),SEXACT(J),DIFFS 201 CONTINUE C C NOW CALL THE TOP-END FCOEFS USING THE RESTART FEATURE. C PREVIOUS RESULTS WILL BE USED. C IFIRST=1 MTOP=512 C CALL FCOEFS(A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,C,S,IERR,ESTER,MACT) C PRINT 65,IDER,IFIRST,MTOP,TOL 65 FORMAT(1X//4X,'...........................................'/ $ 4X,'THIS IS OUTPUT FROM PACKAGE FOURCO USING TOP-END'/ $ 4X,'FCOEFS WITH RESTART FEATURE'/ $ 4X,'.........................................'/ $ 14X,'IDER=',I3,' IFIRST=',I3/ $ 14X,'MTOP=',I4,' TOL=',E12.2/) C PRINT THE ESTIMATED ABSOLUTE ERROR, THE ERROR PARAMETER AND C THE VALUE OF MACT C PRINT 10,ESTER,IERR,MACT C C PRINT THE VALUE OF THE APPROXIMATED COSINE COEFFICIENTS COMPUTED C BY FCOEFS , THE VALUES OF THE EXACT COEFFICIENTS C AND THE ABSOLUTE ERRORS C PRINT 20 DO 102 I=1,NCOEF,5 J=I/5 +1 DIFFS=C(I)-CEXACT(J) PRINT 30,I,C(I),CEXACT(J),DIFFS 102 CONTINUE C C AND NOW THE SINE COEFFICIENTS ... C PRINT 40 DO 202 I=1,NCOEF,5 J=I/5 +1 DIFFS=S(I)-SEXACT(J) PRINT 50,I,S(I),SEXACT(J),DIFFS 202 CONTINUE C C C NOW CALL THE TOP-END FCOEFS USING ONESID(ONE-SIDED NUMERICAL C DIFFERENTIATION ROUTINE) C IDER=2 IFIRST=0 MTOP=32 C CALL FCOEFS(A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,C,S,IERR,ESTER,MACT) C PRINT 75,IDER,IFIRST,MTOP,TOL 75 FORMAT(1X//4X,'...........................................'/ $ 4X,'THIS IS OUTPUT FROM PACKAGE FOURCO USING TOP-END'/ $ 4X,'FCOEFS WITH ONESID(ONE-SIDED NUM.DIFFERENT.ROUTINE)'/ $ 4X,'.........................................'/ $ 14X,'IDER=',I3,' IFIRST=',I3/ $ 14X,'MTOP=',I4,' TOL=',E12.2/) C PRINT THE ESTIMATED ABSOLUTE ERROR, THE ERROR PARAMETER AND C THE VALUE OF MACT C PRINT 10,ESTER,IERR,MACT C C PRINT THE VALUE OF THE APPROXIMATED COSINE COEFFICIENTS COMPUTED C BY FCOEFS , THE VALUES OF THE EXACT COEFFICIENTS C AND THE ABSOLUTE ERRORS C PRINT 20 DO 103 I=1,NCOEF,5 J=I/5 +1 DIFFS=C(I)-CEXACT(J) PRINT 30,I,C(I),CEXACT(J),DIFFS 103 CONTINUE C C AND NOW THE SINE COEFFICIENTS ... C PRINT 40 DO 203 I=1,NCOEF,5 J=I/5 +1 DIFFS=S(I)-SEXACT(J) PRINT 50,I,S(I),SEXACT(J),DIFFS 203 CONTINUE C C C NOW CALL THE TOP-END FCOEFS PASSING PRECALCULATED DERIVATIVES. C XLAMB(I)=I-TH DERIVATIVE OF FUN AT B MINUS I-TH DERIVATIVE OF C FUN AT A. C IDER=-5 DO 121 I=1,ABS(IDER) XLAMB(I) = 3**I * (EXP(3.E0*B) - EXP(3.E0*A)) 121 CONTINUE IFIRST=0 C CALL FCOEFS(A,B,FUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,C,S,IERR,ESTER,MACT) C PRINT 55,ABS(IDER),IDER,IFIRST,MTOP,TOL 55 FORMAT(1X//4X,'...........................................'/ $ 4X,'THIS IS OUTPUT FROM PACKAGE FOURCO USING TOP-END'/ $ 4X,'FCOEFS PASSING ',I2,' PRECALCULATED EXACT DERIVATIVES'/ $ 4X,'.........................................'/ $ 14X,'IDER=',I3,' IFIRST=',I3/ $ 14X,'MTOP=',I4,' TOL=',E12.2/) C PRINT THE ESTIMATED ABSOLUTE ERROR, THE ERROR PARAMETER AND C THE VALUE OF MACT C PRINT 10,ESTER,IERR,MACT C C PRINT THE VALUE OF THE APPROXIMATED COSINE COEFFICIENTS COMPUTED C BY FCOEFS , THE VALUES OF THE EXACT COEFFICIENTS C AND THE ABSOLUTE ERRORS C PRINT 20 DO 104 I=1,NCOEF,5 J=I/5 +1 DIFFS=C(I)-CEXACT(J) PRINT 30,I,C(I),CEXACT(J),DIFFS 104 CONTINUE C C AND NOW THE SINE COEFFICIENTS ... C PRINT 40 DO 204 I=1,NCOEF,5 J=I/5 +1 DIFFS=S(I)-SEXACT(J) PRINT 50,I,S(I),SEXACT(J),DIFFS 204 CONTINUE C C C NOW CALL THE TOP-END FCOETD THAT REQUIRES SAMPLED FUNCTION C VALUES. C C MTOP=33 STEP=(B-A)/(MTOP-1) DO 105 I=1,MTOP FUNVAL(I)=FUN(A+(I-1)*STEP) 105 CONTINUE CALL FCOETD(A,B,NCOEF,MTOP,FUNVAL,WORK,C,S,IERR,ESTER) C PRINT 115,IDER,IFIRST,MTOP,TOL 115 FORMAT(1X//4X,'...........................................'/ $ 4X,'THIS IS OUTPUT FROM PACKAGE FOURCO USING TOP-END'/ $ 4X,'FCOETD THAT REQUIRES SAMPLED FUNCTION VALUES'/ $ 4X,'.........................................'/ $ 14X,'IDER=',I3,' IFIRST=',I3/ $ 14X,'MTOP=',I4,' TOL=',E12.2/) C PRINT THE ESTIMATED ABSOLUTE ERROR, THE ERROR PARAMETER AND C THE VALUE OF MACT C PRINT 10,ESTER,IERR,MACT C C PRINT THE VALUE OF THE APPROXIMATED COSINE COEFFICIENTS COMPUTED C BY FCOEFS , THE VALUES OF THE EXACT COEFFICIENTS C AND THE ABSOLUTE ERRORS C PRINT 20 DO 106 I=1,NCOEF,5 J=I/5 +1 DIFFS=C(I)-CEXACT(J) PRINT 30,I,C(I),CEXACT(J),DIFFS 106 CONTINUE C C AND NOW THE SINE COEFFICIENTS ... C PRINT 40 DO 206 I=1,NCOEF,5 J=I/5 +1 DIFFS=S(I)-SEXACT(J) PRINT 50,I,S(I),SEXACT(J),DIFFS 206 CONTINUE C C STOP C******************************************************************** END REAL FUNCTION FUN(X) C********************************************************************* C C THIS IS THE TEST FUNCTION C REAL X FUN=EXP(3.0E0*X) RETURN C C END OF THE TEST FUNCTION C******************************************************************** C END OF THE SIMPLE DRIVER C******************************************************************** END PROGRAM CDRIVE C*********************************************************************** C C DRIVER FOR THE PACKAGE FOURCO C SINGLE PRECISION VERSION C*********************************************************************** C C INCLUDES MAIN, TRYFUN, INPUT, EXACT, BERF . C MAIN CALLS INPUT AND THEN FOURCO. C C THIS DRIVER HAS BEEN DESIGNED FOR A STRINGENT TEST OF THE C PACKAGE FOURCO. C C THE VARIOUS OPTIONS ARE DESCRIBED BELOW AND IN THE SUBROUTINE C INPUT. C C*********************************************************************** C C THE LIST OF TEST FUNCTIONS IS IN THE SUBROUTINE INPUT. C THE VALUE OF NFUN IN THE SUBROUTINE INPUT SELECTS A TESTFUNCTION. C THE ARRAY AVEC CONTAINS THE PARAMATERS AMP,ALPHA AND BETA C WHICH MODIFY THE BEHAVIOUR OF THE TEST FUNCTION (SEE SUBROUTINE C INPUT). C C THE VALUE OF THE REQUESTED INPUT ACCURACY IS ALSO SET IN C THE SUBROUTINE INPUT (RELTOL). C C THE VALUE OF THE ARRAY NIX IN THE SUBROUTINE INPUT ALLOWS TO C SELECT OUTPUT TEST PRINTS. C SET NIX(I)=0 I=1,9 IF ONLY DEFAULT PRINTS ARE REQUESTED. C C IF NIX(10)=0 THEN THE FUNCTION SELECTED IN THE SUBROUTINE INPUT C IS USED C IF NIX(10)=1 THEN FOURCO IS TESTED USING 8 TEST FUNCTIONS C IN THE SUBROUTINE INPUT C C THE VALUE OF NIX(11) SELECTS A RUN OF FCOEFS (NIX(11)=0) C OR A RUN OF FCOETD (NIX(11)=1). C THE VALUE OF NIX(12) IS THE VALUE OF THE INPUT ARGUMENT OF C FCOEFS IDER, WHICH SELECTS THE DIFFERENTIATION ROUTINE. C*********************************************************************** C*********************************************************************** C C MAIN C C*********************************************************************** C IMPLICIT REAL (A-H,O-Z) COMMON /M360/AVEC(10),AAAA,BBBB,A,B,RELTOL,TOLFAC,SFAC,AVINC, + NIX(12),NVEC(20),NFUN,N2TOP,MODSEQ,MTOP C EXTERNAL TRYFUN * DIMENSION FCOS(240),FSIN(240) DIMENSION WORK(5000),XLAMB(20),FUNVAL(1025) DIMENSION CCOS(240),CSIN(240) C C NFUN = 0 CALL INPUT C 1 CONTINUE STEP = 1.0E0 - AVEC(4) TOL = RELTOL NCOEF = 240 C N2 = 1 10 CONTINUE C C ********** LOOP STARTS HERE. RE-ENTRY FROM 900 BELOW *********** C C SECTION TO SET TOLERANCE C IF (MODSEQ.EQ.1) GO TO 30 IF (MODSEQ.EQ.-1) GO TO 30 BIG = 0.0E0 AA = A DIFF = (B-A)/4.0E0 C DO 20 J = 1,5 TEMP = TRYFUN(AA) BIG = MAX(ABS(TEMP),BIG) AA = AA + DIFF 20 CONTINUE TOL = BIG*RELTOL 30 CONTINUE C 40 CONTINUE C NVEC(2) = N2 C PRINT 9000,AVEC(1),AVEC(2),AVEC(3),AVEC(4),NFUN,NVEC(2) * 9000 FORMAT (/,' FUNCTION SPECIFICATION. AVEC, NVEC ',//,4E16.8,2I4,//) C NVEC(10) = 0 C NVEC(10) IS USED TO COUNT NUMBER OF CALLS TO TRYFUN C IF (NIX(11).EQ.1) GO TO 50 PRINT 9010 * 9010 FORMAT (' WELCOME TO FCOEFS. STANDARD TOP END'// + ' INPUT ARGUMENTS ARE NCOEF,TOL,A,B,MTOP,IDER,IFIRST') PRINT 9015,NCOEF,TOL,A,B,MTOP,NIX(12),IFIRST 9015 FORMAT(1X,I4,3E16.8,1X,2I4,I2) GO TO 60 * 50 CONTINUE PRINT 9020 MTOP=257 STEP=(B-A)/(MTOP-1) DO 85 IJ=1,MTOP FUNVAL(IJ)=TRYFUN(A+(IJ-1)*STEP) 85 CONTINUE * 9020 FORMAT (' WELCOME TO FCOETD. TOP END FOR SAMPLED FUNCTION'// + ' INPUT ARGUMENTS ARE NCOEF,A,B,MTOP') PRINT 9030,NCOEF,A,B,MTOP * 9030 FORMAT (1X,I4,2E16.8,I4,/) C 60 CONTINUE IF (NIX(11).EQ.1) GO TO 80 70 CONTINUE IDER=NIX(12) IFIRST=0 CALL FCOEFS(A,B,TRYFUN,NCOEF,TOL,MTOP,IDER,IFIRST,WORK, + XLAMB,FCOS,FSIN,IERR,ESTER,MACT) GO TO 90 * 80 CONTINUE C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CALL FCOETD(A,B,NCOEF,MTOP,FUNVAL,WORK,C,S,IERR,ESTER) MACT=MTOP 90 CONTINUE C PRINT 9040,IERR,NVEC(10),MACT * 9040 FORMAT (/' IERR=',I2,2X,'GLOBAL N.OF FUN.EVAL=',I5,' MACT= ',I5/) C IF (NFUN.NE.1) GO TO 100 PRINT 9050 * 9050 FORMAT (' THE FUNCTION IS K*EXP(ALPHA(X-BETA)) ') PRINT 9060,AVEC(3),AVEC(1),AVEC(2) * 9060 FORMAT (/,'K=',E16.8,2X,'ALPHA=',E16.8,2X,'BETA=',E16.8,/) GO TO 120 C 100 CONTINUE IF (NFUN.NE.2) GO TO 110 PRINT 9070,N2 * 9070 FORMAT (' THE FUNCTION IS THE BERNOULLI POLYNOMIAL OF DEGREE',I3, + /) GO TO 120 C 110 CONTINUE PRINT 9080,NFUN,AVEC(4) * 9080 FORMAT (' THE FUNCTION IS NFUN',I2,' RHO =',E16.8/) C 120 CONTINUE C PRINT 9090,NFUN * 9090 FORMAT (' NFUN = ',I1,/) PRINT 9100,ESTER * 9100 FORMAT ('********** ESTER =',E16.8,' **********',/) C IF (NIX(1).GT.0) PRINT 9110 * 9110 FORMAT (1X//' //// ACCURACY CONTROL SUSPENDED ////'//) IF (NIX(2).EQ.1) PRINT 9120 * 9120 FORMAT (1X,' **** SINGLE PRECISION FUNCTION VALUES ****') C C CALL EXACT(CCOS,CSIN) C PRINT 9130,NFUN,N2 * 9130 FORMAT (//,' **** THE COSINE COEFFICIENTS **** NFUN,N2 =',2I3,/) PRINT 9140 * 9140 FORMAT (16X,'********** REAL **********',/) PRINT 9150 * 9150 FORMAT (2X,'J-1',4X,'J',5X,'CALCULATED',8X,'EXACT',9X, + 'DIFFERENCE') C J = 1 DO 140 JJJ = 1,30 JJ = J - 1 DIFF = FCOS(J) - CCOS(J) NOP = 1 IF (ABS(DIFF).GT.ABS(ESTER)) NOP = -1 PRINT 9160,JJ,J,FCOS(J),CCOS(J),DIFF,NOP GO TO 130 * 9160 FORMAT (2I5,3E16.8,2X,I2) C 130 CONTINUE J = J + 1 IF (J.GE.12) J = J + 9 IF (J.GE.55) J = J + 15 IF (J.GT.240) GO TO 150 140 CONTINUE 150 CONTINUE C NAB = 10 DO 160 I = 1,240 IF (NAB.EQ.0) GO TO 170 DIF = FCOS(I) - CCOS(I) IF (ABS(DIF).LE.ABS(ESTER)) GO TO 160 NAB = NAB - 1 PRINT 9170,I * 9170 FORMAT ('ERROR EXCESSIVE IN COEFFICIENT ',I3) 160 CONTINUE 170 CONTINUE C PRINT 9180,NFUN,N2 * 9180 FORMAT (//,' ***** THE SINE COEFFICIENTS **** NFUN,N2 =',2I3,/) PRINT 9140 PRINT 9150 J = 1 DO 190 JJJ = 1,30 JJ = J - 1 DIFF = FSIN(J) - CSIN(J) NOP = 1 IF (ABS(DIFF).GT.ABS(ESTER)) NOP = -1 PRINT 9160,JJ,J,FSIN(J),CSIN(J),DIFF,NOP GO TO 180 * 180 CONTINUE J = J + 1 IF (J.GE.12) J = J + 9 IF (J.GE.55) J = J + 15 IF (J.GT.240) GO TO 200 190 CONTINUE 200 CONTINUE C NAB = 10 DO 210 I = 1,240 IF (NAB.EQ.0) GO TO 220 DIF = FSIN(I) - CSIN(I) IF (ABS(DIF).LE.ABS(ESTER)) GO TO 210 NAB = NAB - 1 PRINT 9170,I 210 CONTINUE 220 CONTINUE C C THE FUNCTION UPDATE. C ALTER EITHER FUNCTION OR TOLERANCE C N2 = N2 + 1 IF (N2.GT.N2TOP) GO TO 240 IF (MODSEQ.GT.0) GO TO 230 TOL = TOL*TOLFAC GO TO 40 C 230 CONTINUE IF (NFUN.LE.3) AVEC(1) = AVEC(1) + AVINC STEP = SFAC*STEP AVEC(4) = 1.0E0 - STEP GO TO 10 C C 240 CONTINUE C IF(NIX(10).EQ.0) GO TO 250 IF(NFUN.GE.9) GO TO 250 CALL INPUT GOTO 1 C 250 CONTINUE C STOP C C END OF MAIN C*********************************************************************** END REAL FUNCTION TRYFUN(X) C*********************************************************************** C C FUNCTION TRYFUN C C*********************************************************************** C*********************************************************************** C C THIS FUNCTION CONTAINS ALL THE TEST FUNCTIONS. C THE VALUE OF NFUN (IN COMMON),WHICH IS SET IN THE C SUBROUTINE INPUT,SELECTS THE DESIRED FUNCTION C C*********************************************************************** IMPLICIT REAL (A-H,O-Z) REAL SINFUN C TEMPCOM COMMON LIST COMMON /M360/AVEC(10),AAAA,BBBB,A,B,RELTOL,TOLFAC,SFAC,AVINC, + NIX(12),NVEC(20),NFUN,N2TOP,MODSEQ,MTOP C DIMENSION BERFUN(12),BFFACT(12) * TWOPI = 8.0E0*ATAN(1.0E0) ALPHA = AVEC(1) ALSQ = ALPHA*ALPHA BETA = AVEC(2) ARG = ALPHA* (X-BETA) TWOPIX = TWOPI*ARG RHO = AVEC(4) RHOSQ = RHO*RHO NVEC(10) = NVEC(10) + 1 N2 = NVEC(2) C IF (NFUN.LE.0) GO TO 100 10 CONTINUE IF (NFUN.NE.1) GO TO 20 ARG = ALPHA* (X-BETA) TRYFUN = EXP(ARG) GO TO 90 * 20 CONTINUE IF (NFUN.NE.2) GO TO 30 IF (N2.LT.0) GO TO 100 IF (N2.GT.12) GO TO 100 IF (N2.EQ.0) TRYFUN = 1.0E0 IF (N2.EQ.0) GO TO 90 XX = ARG NPER = 0 CALL BERF(XX,NPER,BERZ,BERFUN,BFFACT) TRYFUN = BERFUN(N2) GO TO 90 * 30 CONTINUE IF (NFUN.NE.4) GO TO 40 XNUM = 1.0E0 - RHOSQ DEN = 1.0E0 - 2.0E0*COS(TWOPIX)*RHO + RHOSQ TRYFUN = XNUM/DEN GO TO 90 * 40 CONTINUE IF (NFUN.NE.5) GO TO 50 XNUM = RHO* (SIN(TWOPIX)) DEN = 1.0E0 - RHO* (COS(TWOPIX)) TRYFUN = 2.0E0*ATAN2(XNUM,DEN) GO TO 90 * 50 CONTINUE IF (NFUN.NE.6) GO TO 60 ARG = 1.0E0 - 2.0E0*RHO*COS(TWOPIX) + RHOSQ TERM = 2.0E0*LOG(ARG) TRYFUN = TERM/2.0E0 GO TO 90 * 60 CONTINUE IF (NFUN.NE.7) GO TO 70 XNUM = 2.0E0*RHO*SIN(TWOPIX) DEN = 1.0E0 - RHOSQ TRYFUN = ATAN2(XNUM,DEN) GO TO 90 * 70 CONTINUE IF (NFUN.NE.8) GO TO 80 ARG = RHO*COS(TWOPIX) TERM1 = EXP(ARG) ARG2 = RHO*SIN(TWOPIX) TERM2 = COS(ARG2) TRYFUN = TERM1*TERM2 GO TO 90 * 80 CONTINUE TERM = (1.0E0+RHO*SIN(TWOPIX)) TRYFUN = 1.0E0/TERM 90 CONTINUE TRYFUN = AVEC(3)*TRYFUN IF (NIX(4).EQ.0) GO TO 110 PRINT 9000,NIX(4) * 9000 FORMAT (/,'OUTPUT FROM TRYFUN. NIX(4) =',I2) PRINT 9010,NVEC(10),NFUN,N2,X,TRYFUN * 9010 FORMAT (/,'TRYFUN',3I5,3X,'X =',E12.4,3X,'FUN =',E12.4) NIX(4) = NIX(4) - 1 GO TO 110 C 100 CONTINUE PRINT 9020,NFUN,N2 * 9020 FORMAT (//,' FUNCTION NOT DEFINED. NFUN AND(2) ARE',2I3,//) C 110 CONTINUE IF (NIX(2).NE.1) GO TO 120 SINFUN = TRYFUN TRYFUN = REAL(SINFUN) 120 CONTINUE C C RETURN C C END OF TRYFUN C*********************************************************************** END SUBROUTINE INPUT C*********************************************************************** C C SUBROUTINE INPUT C C*********************************************************************** C C THIS IS A FREE FORM DATA INPUT SUBROUTINE. IT TAKES THE PLACE C OF A READ STATEMENT. THE DATA IS TRANSMITTED THROUGH COMMON. C ALL THE INPUT ARGUMENTS FOR FCOEFS AND FCOETD ARE SELECTED C IN THIS ROUTINE. C A TEST SESSION SHOULD CONSIST IN CHANGING THOSE INPUT VALUES C (SEE COMMENTS BELOW) AND RUN THE PACKAGE C C*********************************************************************** C IMPLICIT REAL (A-H,O-Z) COMMON /M360/AVEC(10),AAAA,BBBB,A,B,RELTOL,TOLFAC,SFAC,AVINC, + NIX(12),NVEC(20),NFUN,N2TOP,MODSEQ,MTOP C C*********************************************************************** C SET THE PARAMETERS NIX C C EACH NIX(J) J LT 10 CONTROLS THE PRINTING OF CODE CHECK OUTPUT IN C SOME PART OF THE PROGRAM. IN GENERAL THE EFFECT IS THIS C IF NIX(J) .LT.0 THERE IS ALWAYS PRINT OUT C IF NIX(J) .EQ.0 THERE IS NEVER PRINT OUT C IF NIX(J) = N .GT.0 THERE IS PRINT OUT THE FIRST N TIMES THE C FIRST N TIMES THE PRINT STATEMENT IS ENCOUNTERED. C C NIX(12) IS THE VALUE OF IDER THE INPUT ARGUMENT OF FCOEFS C WHICH SELECTS THE DIFFERENTIATION ROUTINE USED BY C THE PROCEDURE. C NIX(11) IF EQ 0 CHOOSE STANDARD TOP END. IF EQ 1 NON STANDARD. C NIX(10) IF EQ 0 NFUN SELECTED IN INPUT IS USED C IF EQ 1 A LOOP ON 8 TEST FUNCTIONS IS EXECUTED(ONLY C NFUN=3 IS NOT USED) C NIX(9) NOT USED C NIX(8) NOT USED C NIX(7) NOT USED C NIX(6) NOT USED C NIX(5) IN SPEC. PRINTS PARAMETERS C NIX(4) IN TRYFUN. PRINTS XVALUE, FUNCTION VALUE C NIX(3) NOT USED C NIX(2) IF EQ 1 TRYFUN RETURNS SINGLE PRECISION VALUES C NIX(1) NOT USED C C*********************************************************************** C NIX(12)= 1 NIX(11)= 0 NIX(10)= 1 NIX(9) = 0 NIX(8) = 0 NIX(7) = 0 NIX(6) = 0 NIX(5) = 0 NIX(4) = 0 NIX(3) = 0 NIX(2) = 0 NIX(1) = 0 C C SPECIFY THE INTERVAL C A = 0.0E0 B = 1.0E0 IF (B.LE.A) PRINT 9000 * 9000 FORMAT ('***** WARNING B MUST BE GREATER THAN A *****') C C C*********************************************************************** C THE USER MUST SPECIFY A VALUE OF NFUN FROM THE C FOLLOWING POSSIBLE CASES. ENTER YOUR VALUE BELOW. C C NFUN = 1, A SET OF MORE PRONOUNCED EXPONENTIALS C NFUN = 2, A SET OF BERNOULLI POLYNOMIALS C NFUN = 3, A SET OF STRICTER TOLERANCES USING ONE EXP.FUN C NFUN = 4, F(X)=(1-P**2)/(1-2*P*COS(2*PI*X)+P**2) C NFUN = 5, F(X)= 2*ARCTAN((PSIN(2*PI*X))/(1-PCOS(2*PI*X))) C NFUN = 6, F(X)= LN(1-2PCOS(2*PI*X)+P**2) C NFUN = 7, F(X)=ARCTAN((2PSIN(2*PI*X))/(1-P**2) C NFUN = 8, F(X)=EXP(RHOCOS(2*PI*X))*COS(RHOSIN(2*PI*X)) C NFUN = 9, F(X) = 1 / (1+PSIN(2*PI*X)) C C*********************************************************************** IF(NIX(10).EQ.0)THEN C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C CHANGE NEXT STATEMENT IF YOU WANT TO CHANGE THE TEST FUNCTION C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NFUN = 1 ELSE NFUN = NFUN+1 IF(NFUN.EQ.3) NFUN = 4 ENDIF C*********************************************************************** C NOW SET THE PARAMETERS FOR THE FIRST FUNCTION OF THIS SET. C THESE ARE AVEC(J), J=1,2,3,4. C C THE FUNCTION IS AMP * F(ALPHA X - BETA') C C AVEC(1) REPRESENTS ALPHA C AVEC(2) REPRESENTS BETA C AVEC(3) REPRESENTS THE OVERALL MULTIPLIER AMP. C AVEC(4) REPRESENTS RHO. WHEN SFAC = 0.5 AND AVEC(4) = 0.2 WE GET C RHO =.2,.6,.8,.9,..... C C*********************************************************************** AVEC(1) = 1.0E0 AVEC(2) = 0.125E0 AVEC(2) = 1.0E0 AVEC(2) = 0.0E0 AVEC(3) = 1.0E0 AVEC(4) = 0.2E0 IF (NFUN.GE.4) AVEC(1) = 1.0E0/ (B-A) C C THESE NEXT PARAMETERS GEAR THE NEXT FUNCTIONS C N2TOP = 2 AVINC = 3.0E0 IF (NFUN.NE.1) AVINC = 0.0E0 SFAC = 0.5E0 C C SPECIFY MODSEQ. POSITIVE FOR VARYING FUNCTION. 1 FOR ABS TOLERANCE. C NEGATIVE FOR VARYING TOLERANCE. 2 FOR REL TOLERANCE C MODSEQ = 2 IF (NFUN.EQ.3) MODSEQ = -1 IF (NFUN.EQ.2) MODSEQ = IABS(MODSEQ) C C RELTOL IS THE INITIAL TOLERANCE. (ABSOLUTE WHEN MODSEQ = +1 OR -1.) C (RELATIVE TO FUNAVE WHEN MODSEQ = +2 OR -2.) FUNAVE IS THE C AVERAGE OF FIVE EQUALLY SPACED FUNCTION VALUES. C RELTOL = 1.0E-8 TOLFAC = SQRT(0.1E0) IF (MODSEQ.LT.0) RELTOL = 1.0E-3 C MTOP = 512 C C THERE IS NO NFUN=3 AND NFUN = 2 ARE BERNOULLI POLYNOMIALS.HENCE C IF (NFUN.EQ.2) A = 0.0E0 IF (NFUN.EQ.2) B = 1.0E0 C IF (NFUN.EQ.3) THEN NFUN = 1 RELTOL = 1.0E-3 AVEC(1) = 16.0E0 AVINC = 0.0E0 END IF C PRINT 9010,NFUN * 9010 FORMAT (' THIS IS OUTPUT FROM A RUN OF DRIVER WITH NFUN = ',I3, + //) PRINT 9020,A,B * 9020 FORMAT (' INTERVAL A = ',D8.2,9X,'B = ',D8.2,/) IF (NIX(1).GT.0) PRINT 9030 * 9030 FORMAT (1X//' //// ACCURACY CONTROL SUSPENDED ////'//) IF (NIX(2).EQ.1) PRINT 9040 * 9040 FORMAT (1X,'**** SINGLE PRECISION FUNCTION VALUES ****') PRINT 9050 * 9050 FORMAT (//' THE FIRST FUNCTION TREATED HAS PARAMETERS AS FOLLOWS' + ) PRINT 9060 * 9060 FORMAT (/,' NFUN N2 AVEC(1)=ALPHA',7X,'AVEC(2)=BETA',9X, + 'AVEC(3)=AMP',9X,'AVEC(4)=RHO') PRINT 9070,NFUN, (AVEC(KKK),KKK=1,4) * 9070 FORMAT (I4,' 1 ',E15.7,3D20.7,/) C PRINT 9080,MTOP,SFAC,RELTOL * 9080 FORMAT (' MTOP = ',I3,9X,'SFAC = ',D8.2,5X,'RELTOL = ',D8.2,/) PRINT 9090,N2TOP,AVINC,TOLFAC * 9090 FORMAT (' N2TOP = ',I3,10X,'AVINC = ',D8.2,5X,'TOLFAC = ',D8.2,/) PRINT 9100,MODSEQ * 9100 FORMAT (' MODSEQ =',I3,4X, + 'POSITIVE MEANS FUNCTION VARIES. NEGATIVE',1X, + 'MEANS TOLERANCE VARIES',//) C IF (NIX(5).NE.0) THEN PRINT 9110 * 9110 FORMAT (' CODE-CHECK INPUT PARAMETERS',/,/) PRINT 9120,NIX(10),NIX(9),NIX(8) * 9120 FORMAT (' NIX(10) = ',I3,8X,'NIX(9) = ',I3,8X,'NIX(8) = ',I3,/) PRINT 9130,NIX(6),NIX(4),NIX(3) * 9130 FORMAT (' NIX(6) = ',I3,9X,'NIX(4) = ',I3,8X,'NIX(3) = ',I3,/) END IF C RETURN C C END OF INPUT C************************************************************** END SUBROUTINE EXACT(CCOS,CSIN) C************************************************************** C C SUBROUTINE EXACT C C************************************************************** C C EXACT RETURNS THE EXACT VALUES OF THE FOURIER COEFFICIENTS C OF THE TEST FUNCTION SELECTED IN THE SUBROUTINE INPUT. C THIS SUBROUTINE CONTAINS THE ANALYTIC FORM OF THE FOURIER C COEFFICIENTS OF ALL THE TEST FUNCTIONS PROVIDED BY THE FUNCTION C TRYFUN AND MODIFIED BY MEANS OF PARAMETERS IN THE ARRAY AVEC C C**************************************************************** IMPLICIT REAL (A-H,O-Z) COMMON /M360/AVEC(10),AAAA,BBBB,A,B,RELTOL,TOLFAC,SFAC,AVINC, + NIX(12),NVEC(20),NFUN,N2TOP,MODSEQ,MTOP C DIMENSION CCOS(240),CSIN(240) C CALCULATES THE FIRST 240 FOURIER COEFFICIENTS OF TRYFUN USING C AN ANALYTIC EXPRESSION. C TWOPI = 8.0E0*ATAN(1.0E0) ALPHA = AVEC(1) ALSQ = ALPHA*ALPHA BETA = AVEC(2) CON = 1.0E0 RHO = AVEC(4) RHOSQ = RHO*RHO FACTOR = ALPHA* (A-BETA) C DO 10 J = 1,240 CCOS(J) = 0.0E0 CSIN(J) = 0.0E0 10 CONTINUE C N2 = NVEC(2) IF (NFUN.LE.0) GO TO 160 C IF (NFUN.NE.1) GO TO 30 C BMA = B - A CONS = EXP(-ALPHA*BETA) TERM = EXP(ALPHA*B) - EXP(ALPHA*A) CCOS(1) = (TERM*CONS)/ (ALPHA*BMA) FACT = 2.0E0*TERM*CONS/BMA TPR = TWOPI DO 20 JJ = 1,239 JJJ = JJ + 1 RAT = FACT/ ((TPR*TPR)/ (BMA*BMA)+ALSQ) CCOS(JJJ) = RAT*ALPHA CSIN(JJJ) = -RAT*TPR/BMA TPR = TPR + TWOPI 20 CONTINUE GO TO 170 C 30 CONTINUE IF (NFUN.NE.2) GO TO 60 C C NFUN = 2 FUNCTION IS BERNOULLI POLYNOMIAL OF DEGREE N2. C IF (N2.LT.0) GO TO 160 IF (N2.EQ.0) CCOS(1) = 1.0E0 IF (N2.EQ.0) GO TO 170 C FACTN2 = 1.0E0 XJ = 1.0E0 DO 40 J = 1,N2 FACTN2 = FACTN2*XJ XJ = XJ + 1.0E0 40 CONTINUE C NREM = N2 - 4* (N2/4) FACT = 2.0E0*FACTN2/ (TWOPI**N2) XJJ = 1.0E0 DO 50 J = 2,240 TEMP = XJJ** (-N2) RES = FACT*TEMP IF (NREM.EQ.0) CCOS(J) = -RES IF (NREM.EQ.2) CCOS(J) = RES IF (NREM.EQ.1) CSIN(J) = -RES IF (NREM.EQ.3) CSIN(J) = RES XJJ = XJJ + 1.0E0 50 CONTINUE C GO TO 170 C 60 CONTINUE IF (NFUN.EQ.9) GO TO 140 TERM = 1.0E0 DENOM = 1.0E0 IF (NFUN.EQ.4) CCOS(1) = 1.0E0 IF (NFUN.EQ.8) CCOS(1) = 1.0E0 DO 130 J = 1,239 K = J + 1 XJ = J ARG = TWOPI*XJ*FACTOR CTERM = COS(ARG) STERM = SIN(ARG) TERM = TERM*RHO COF = CTERM*TERM SOF = STERM*TERM 70 CONTINUE IF (NFUN.EQ.4) THEN CCOS(K) = 2.0E0*COF CSIN(K) = -2.0E0*SOF END IF * 80 CONTINUE IF (NFUN.EQ.5) THEN CCOS(K) = 2.0E0*SOF/XJ CSIN(K) = 2.0E0*COF/XJ END IF * 90 CONTINUE IF (NFUN.EQ.6) THEN CCOS(K) = -2.0E0*COF/XJ CSIN(K) = 2.0E0*SOF/XJ END IF * 100 CONTINUE IF (NFUN.EQ.7) THEN NREM = J - 2* (J/2) IF (NREM.EQ.0) GO TO 130 CCOS(K) = 2.0E0*SOF/XJ CSIN(K) = 2.0E0*COF/XJ END IF * 110 CONTINUE IF (NFUN.EQ.8) THEN DENOM = DENOM*XJ CCOS(K) = COF/DENOM CSIN(K) = -SOF/DENOM IF (DENOM.GT.1.0D30) THEN DO 120 L = K,240 CCOS(L) = 0.0E0 120 CONTINUE GO TO 170 * END IF * END IF * 130 CONTINUE C GO TO 170 C 140 CONTINUE IF (NFUN.NE.9) GO TO 160 TERM = 1.0E0 - RHOSQ ADB = SQRT(TERM)/TERM CCOS(1) = ADB CPAR = (SQRT(TERM)-1.0E0)/RHO CSQR = CPAR*CPAR CP = -ADB*CSQR SP = ADB*CPAR DO 150 J = 2,238,2 XJ = J XJM1 = J - 1 ARGE = TWOPI*XJM1*FACTOR ARGO = TWOPI*XJ*FACTOR CCOS(J) = 2.0E0*SP*SIN(ARGE) CSIN(J) = 2.0E0*SP*COS(ARGE) CCOS(J+1) = 2.0E0*CP*COS(ARGO) CSIN(J+1) = -2.0E0*CP*SIN(ARGO) CP = -CP*CSQR SP = -SP*CSQR 150 CONTINUE C GO TO 170 C 160 CONTINUE PRINT 9000,NFUN,N2 * 9000 FORMAT (//,' FUNCTION NOT DEFINED. NFUN AND(2) ARE',2I3,//) C 170 CONTINUE C AMPLIT = AVEC(3) DO 180 J = 1,240 CCOS(J) = AMPLIT*CCOS(J) CSIN(J) = AMPLIT*CSIN(J) 180 CONTINUE C RETURN C ********** END OF EXACT ********** C*********************************************************************** END SUBROUTINE BERF(XX,NPER,BERZ,BERFUN,BFFACT) C*********************************************************************** C C SUBROUTINE BERF C C*********************************************************************** C C BERF COMPUTES BERNOULLI FUNCTIONS OR POLYNOMIAL UP TO THE DEGREE C 12 C C XX ARGUMENT OF BERNOULLI FUNCTION OR POLYNOMIAL C NPER EQ.0 GIVES THE POLYNOMIAL (NON-PERIODIC). C NE.0 GIVES THE (PERIODIC) FUNCTION, WHICH COINCIDES WITH T C POLYNOMIAL IN THE OPEN INTERVAL (0,1). C BERZ = 1.0E0 C BERFUN(12) THE VALUES OF THE BERNOULLI POLYNOMIALS (OR FUNCTIONS) C OF DEGREES 1 TO 12 INCLUSIVE. C BFFACT(12) BFFACT(J) = BERFUN(J) / FACTORIAL(J) C*********************************************************************** IMPLICIT REAL (A-H,O-Z) DIMENSION BERFUN(12),BFFACT(12) * X = XX BERZ = 1.0E0 BERFUN(1) = X - 0.5E0 IF (NPER.EQ.0) GO TO 70 C C WE WANT PERIODIC FUNCTIONS 10 CONTINUE IF (X) 20,50,30 20 CONTINUE X = X + 1.0E0 GO TO 10 * 30 CONTINUE IF (X-1.0E0) 60,50,40 40 CONTINUE X = X - 1.0E0 GO TO 30 * 50 CONTINUE C HERE WE HAVE X = 0.0 OR 1.0 BERFUN(1) = 0.0E0 60 CONTINUE C 70 CONTINUE Y = X* (X-1.0E0) Z = Y* (X-0.5E0) BERFUN(2) = Y + (1.0E0/6.0E0) BERFUN(3) = Z BERFUN(4) = Y*Y - (1.0E0/30.0E0) BERFUN(5) = (Y- (1.0E0/3.0E0))*Z BERFUN(6) = Y*Y*Y - 0.5E0*Y*Y + (1.0E0/42.0E0) X2 = X*X X4 = X2*X2 X6 = X4*X2 X8 = X4*X4 X10 = X4*X6 X12 = X4*X8 Y2 = (X-0.5E0)* (X-0.5E0) Y4 = Y2*Y2 BERFUN(7) = (Y4-1.5E0*Y2+ (31.0E0/48.0E0))*Z TERM = (2.0E0*X2-7.0E0*X4+14.0E0*X6)/3.0E0 BERFUN(8) = X8 - 4.0E0*X*X6 + TERM - (1.0E0/30.0E0) TERM = X8 - 4.5E0*X*X6 + 6.0E0*X6 - 4.2E0*X4 + 2.0E0*X2 - 0.3E0 BERFUN(9) = X*TERM TERM = X10 - 5.0E0*X*X8 + 7.5E0*X8 - 7.0E0*X6 + 5.0E0*X4 BERFUN(10) = TERM - 1.5E0*X2 + (5.0E0/66.0E0) TERM = X10 - 5.5E0*X*X8 - 11.0E0* (X6-X4) - 5.5E0*X2 BERFUN(11) = (TERM+ ((5.0E0+55.0E0*X8)/6.0E0))*X TERM = X12 - 6.0E0*X*X10 + 11.0E0*X10 - 16.5E0* (X8+X4) BERFUN(12) = TERM + 22.0E0*X6 + 5.0E0*X2 - (691.0E0/2730.0E0) FACT = 1.0E0 FICT = 1.0E0 DO 80 J = 1,12 BFFACT(J) = BERFUN(J)/FACT FICT = FICT + 1.0E0 FACT = FACT*FICT 80 CONTINUE RETURN C C END OF BERF C*********************************************************************** C C ********* END OF DRIVER ********* C C*********************************************************************** END