SUBROUTINE FSPL2 ( F, A, B, FPA, FPB, FBA, FBB, W, EPS, MAX, & LC, LS, C, S ) c*********************************************************************72 C C THIS SUBROUTINE HAS ADAPTED SEVERAL IDEAS FROM C ALGORITHM 353, FILON QUADRATURE BY CHASE AND FOSDICK, C COMM. ACM, VOL. 12, P.457-458, 1969. C C F(X)=THE FUNCTION TO BE INTEGRATED, SUPPLIED BY THE USER C AND DECLARED 'EXTERNAL' IN THE CALLING PROGRAM. C DATA PI / 3.141592653589793 / C C A=LOWER QUADRATURE LIMIT AND B=UPPER QUADRATURE LIMIT. C IF A.GE.B THE COMPUTATION IS BYPASSED AND THE SIGNS OF C LC, LS, AND EPS ARE CHANGED. C C FPA AND FPB ARE THE VALUES OF THE DERIVATIVE OF F(X). C FBA AND FBB ARE THE CORRESPONDING VALUES OF THE SECOND C DERIVATIVE AT THE POINTS A AND B. C W=THE ANGULAR FREQUENCY. C C EPS = REQUIRED ACCURACY, DEFINED BY C |ERROR| <= EPS*(1.+|C|) C AND C |ERROR| <= EPS*(1.+|S|) C IF CONVERGENCE IS NOT OBTAINED, THE VALUE C OF EPS IS RETURNED WITH NEGATIVE SIGN. C C MAX=THE MAXIMUM NUMBER OF PARTITIONS (INTERVAL HALVINGS). C IN THIS ROUTINE THE INTERNAL VARIABLE MXN DEFINED BELOW C IS USED INSTEAD OF MAX. C LC POSITIVE ON ENTRY INDICATES THAT C IS WANTED. C LS POSITIVE ON ENTRY INDICATES THAT S IS WANTED. C ON EXIT LC AND LS GIVE THE NUMBER OF PARTITIONS USED C FOR THE COMPUTATION OF C AND S. C C THIS ROUTINE CALLS THE SUBROUTINE ENDT2. C DIMENSION PVTC(7),PVTS(7) IF(EPS.LT.0.) GOTO 5 IF(A.LT.B) GOTO 10 EPS=-EPS 5 LC=-LC LS=-LS RETURN 10 N=1 W1=ABS(W) TEMP=2.0*(B-A)*W1/PI IF(TEMP.GT.2.0) N=ALOG(TEMP)/0.693 C 0.693=ALOG(2.) ROUNDED DOWNWARDS. MXN=MAX0(MAX,N+1) FA=F(A) FB=F(B) COSA=COS(W1*A) SINA=SIN(W1*A) COSB=COS(W1*B) SINB=SIN(W1*B) H=(B-A)/FLOAT(2**N) NSTOP=2**N-1 NST=1 C TMAX IS THE SWITCH-OVER POINT FOR TETA. C ANALYSIS SHOWS THAT WITH A 56 BIT FLOATING POINT MANTISSA TMAX=0.2 C IS SUITABLE, WHILE WITH A 24 BIT MANTISSA WE PREFER TMAX=1. C TMAXB IS THE SWITCH-OVER POINT IN BETA, WHERE THE C CANCELLATION IS STRONGEST. TMAXB=5.*TMAX C LLC AND LLS ARE USED BY THE ROUTINE IN COMPUTED-GO-TO C STATEMENTS. AS SOON AS LLS AND LLC HAVE BEEN DEFINED, C WE CAN USE LS AND LC AS RETURN PARAMETERS (SEE ABOVE). IF(LS)11,11,12 11 LLS=2 GOTO 13 12 LLS=1 13 IF(LC)14,14,15 14 LLC=2 GOTO 17 15 LLC=1 17 CONTINUE SUMCOS=0.5*(FA*COSA+FB*COSB) SUMSIN=0.5*(FA*SINA+FB*SINB) C ALL OF THE ABOVE IS EXECUTED ONLY ONCE PER CALL. C NOW THE ITERATION BEGINS. C THE CONSTANT 'M' IS USED IN THE RICHARDSON EXTRAPOLATION. C M-1 IS THE NUMBER OF TIMES THE ORIGINAL STEP LENGTH 'H' C HAS BEEN DIVIDED BY TWO. M=1 20 CONTINUE H2=H*H TETA=W1*H DO 65 I=1,NSTOP,NST X=A+H*FLOAT(I) WX=W1*X GOTO (50,55),LLS 50 SUMSIN=SUMSIN+F(X)*SIN(WX) 55 GOTO (60,65),LLC 60 SUMCOS=SUMCOS+F(X)*COS(WX) 65 CONTINUE T2=TETA*TETA TEMP=1.0-SIN(0.5*TETA)**2/1.5 IF(TETA-TMAX)70,70,75 C 70 IS THE POWER SERIES FOR SMALL TETA, 75 IS THE CLOSED C FORM USED WITH LARGER VALUES OF TETA. C THE COEFFICIENTS OF THE DIFFERENT POWER SERIES BELOW ARE C GIVEN IN EXACT FORM, COMPARE WITH THE REFERENCE ABOVE. 70 ALFA=TETA*(1.0-T2*(2.0/15.0-T2*(19.0/1680.0- & T2*(13.0/25200.0-T2*(293.0/19958400.0- & T2*181.0/619164000.0)))))/12.0 DELTA=-1.0/12.0+T2*(1.0/90.0-T2*(5.0/12096.0- & T2*(1.0/129600.0-T2/11404800.0))) EPSIL=1.0-T2*(1.0/6.0-T2*(0.0125-T2*(17.0/30240.0- & T2*(31.0/1814400.0-T2/2661120.0)))) T3=T2 72 BETA=TETA*H2*(1.0-T2/21.0*(1.0-T2*(1.0/48.0- & T2*(1.0/3960.0-T3/494208.0))))/180.0 GOTO 80 C CLOSED FORM OF THE COEFFICIENTS. 75 TEMP1=(0.5*TETA)**2 TEMP2=SIN(0.5*TETA)**2/TEMP1 TEMP3=SIN(TETA)/TETA ALFA=(TEMP-TEMP2*TEMP3)/TETA DELTA=(TEMP-TEMP2)/T2 EPSIL=TEMP2*TEMP2 IF(TETA-TMAXB) 76,76,78 76 T3=T2*(1.-T2*(1./175.-T2*(1./40800.-T2/12209400.))) GOTO 72 78 BETA=(TEMP-TEMP3)/(TETA*W1*W1) C HAVE CALCULATED THE COEFFICIENTS, NOW READY FOR THE C INTEGRATION FORMULAS. 80 GOTO (81,85),LLS 81 TS=H*((BETA*FBB-ALFA*FB)*COSB+(ALFA*FA-BETA*FBA)*COSA+ & DELTA*H*(FPB*SINB-FPA*SINA)+EPSIL*SUMSIN)/TEMP CALL ENDT2(PVTS,TS,EPS,S,LLS,M) LS=N 85 GOTO (86,90),LLC 86 TC=H*((ALFA*FB-BETA*FBB)*SINB+(BETA*FBA-ALFA*FA)*SINA+ & DELTA*H*(FPB*COSB-FPA*COSA)+EPSIL*SUMCOS)/TEMP CALL ENDT2(PVTC,TC,EPS,C,LLC,M) LC=N 90 CONTINUE C NOW TEST TO SEE IF DONE. IF(LLC+LLS-3) 92,92,100 92 N=N+1 C THIS IS THE BEGINNING OF THE ITERATION. IF(N-MXN) 95,95,99 95 H=0.5*H NST=2 NSTOP=2**N M=M+1 GOTO 20 99 EPS=-EPS 100 CONTINUE IF(LS.GT.0.AND.W.LT.0.0) S=-S RETURN END SUBROUTINE ENDT2(PREVOT,QUANT,EPS,VALUE,L,M) c*********************************************************************72 C C ENDT2 IS A SUBROUTINE THAT PERFORMS RICHARDSON EXTRA- C POLATION OF THE VALUES 'QUANT' WHICH ARE INTRODUCED INTO C THE ROUTINE EACH TIME IT IS CALLED, EACH TIME WITH C INCREASING VALUE OF 'M', STARTING WITH M = 1. THE CURRENT C VALUES ARE STORED IN THE ARRAY 'PREVOT', WHERE 'PREVOT(1)' C AT EXIT IS EQUAL TO 'QUANT'. THE BEST VALUE FOR THE MOMENT C IS GIVEN IN 'VALUE'. ENDT2 REQUIRES THE PRESENT VALUE TO C AGREE WITH THE PREVIOUS VALUE TO WITHIN EPS2, WHERE C EPS2 = EPS*(1.0 + ABS(PRESENT VALUE)). C EPS IS SUPPLIED BY THE USER. C THE ERROR EXPANSION IS OF THE FORM C ERROR = C4*H**4 + C6*H**6 + C8*H**8 + ... + CN*H**N + ... C DIMENSION PREVOT(7),RICH(7) DATA RICH(1) / 0.0/, RICH(2) / 15.0/, & RICH(3) / 63.0/, RICH(4) / 255.0/, & RICH(5) / 1023.0/, RICH(6) /4095.0/, & RICH(7) /16383.0/ C RICH(1) IS NOT USED. C RICH(K) = 2**(2*K) - 1, K = 2,3,4,5,6,7. TEMP2=PREVOT(1) PREVOT(1)=QUANT TEMP1=QUANT IF(M.EQ.1)GOTO 30 20 REPS=EPS*(1.0+ABS(QUANT)) DO 23 K=2,M DIFF=TEMP1-TEMP2 IF(ABS(DIFF)-REPS) 25,25,22 22 IF(K.EQ.8) GOTO 30 TEMP1=TEMP1+DIFF/RICH(K) TEMP2=PREVOT(K) PREVOT(K)=TEMP1 23 CONTINUE GO TO 30 25 L=2 30 VALUE=TEMP1 RETURN END