SUBROUTINE AUSINC(F, A, B, LAMBDA, EPSREL, EPSABS, MAXITR, * S, NFCALL, IER ) C*********************** PURPOSE: ******************************* C GIVEN F(X), LAMBDA, THIS SUBROUTINE COMPUTES S, THE SINC C APPROXIMATION TO THE INTEGRAL OF F(X)/(X-LAMBDA) OVER [A,B], C WHERE A < LAMBDA < B C A SEQUENCE OF APPROXIMATIONS IS COMPUTED UNTIL CONVERGENCE IS C ATTAINED. A MIXED ABSOLUTE/RELATIVE ERROR BOUND IS USED, WHERE C THE CURRENT VALUE IS ACCEPTED WHEN THE DIFFERENCE BETWEEN TWO C CONSECUTIVE VALUES IN THE SEQUENCE IS LESS THAN THE LARGER OF C EPSABS AND EPSREL*(ABSOLUTE VALUE OF CURRENT APPROXIMATION). C IF THE USER HAS NO IDEA OF THE MAGNITUDE OF THE INTEGRAL THEN C IT IS RECOMMENDED THAT EPSABS AND EPSREL BE CHOSEN TO BE EQUAL. C IF EPSREL IS SET TO ZERO, THAT IS A PURE ABSOLUTE ERROR TEST C FOR CONVERGENCE IS APPLIED, THEN THE CODE CHECKS TO SEE IF C THE ABSOLUTE VALUE OF THE APPROXIMATION TO THE INTEGRAL IS > 1. C IF SO, THEN EPSREL IS RESET EQUAL TO EPSABS, AND A MIXED C ERROR TEST IS APPLIED. C*********************** PARAMETERS: ******************************* C INPUT: C F(X) DOUBLE PRECISION FUNCTION: OF THE FORM F(X). C A,B DOUBLE PRECISION: THE LIMITS OF INTEGRATION. C LAMBDA DOUBLE PRECISION. THE LOCATION OF THE SINGULARITY. C EPSREL DOUBLE PRECISION: THE REQUESTED RELATIVE ERROR. C EPSREL .GE. 0.0D0. SEE NOTE ABOVE. C EPSABS DOUBLE PRECISION: THE REQUESTED ABSOLUTE ERROR. C EPSABS .GE. 0.0D0 AND C EPSABS + EPSREL .GT. 0.0D0. C MAXITR INTEGER: THE MAXIMUM NUMBER OF ITERATIONS TO BE C CARRIED OUT. AN INITIAL VALUE OF H IS CHOSEN C EMPIRICALLY, AND THEN A SEQUENCE OF AT MOST MAXITR C SINC APPROXIMATIONS WITH H REPLACED BY H/2 IS CARRIED C OUT. IF AGREEMENT BETWEEN TWO OF THESE TO WITHIN A C MIXED RELATIVE AND ABSOLUTE ERROR TEST IS REACHED, C THE CODE RETURNS THE LAST APPROXIMATION, TOGETHER C WITH THE THE LAST ABSOLUTE AND RELATIVE DIFFERENCES C AS ERROR ESTIMATES. C C OUTPUT: C EPSREL DOUBLE PRECISION: THE RELATIVE ERROR ESTIMATE. C EPSABS DOUBLE PRECISION: THE ABSOLUTE ERROR ESTIMATE. C S DOUBLE PRECISION: THE APPROXIMATION. C NFCALL INTEGER: THE NUMBER OF FUNCTION EVALUATIONS. C IER INTEGER: C IER = 0: ERROR REQUEST HAS BEEN SATISFIED. C IER = 1: AFTER MAXITR ATTEMPTS TO FIND TWO C APPROXIMATIONS IN AGREEMENT, THE CODE C FAILED. IT MAY HELP TO INCREASE MAXITR C AND TRY AGAIN. C THE CRITERION USED FOR TERMINATING THE SUM OF THE SERIES IS: C IF ( H*(B-A)*(ABS(TERM)+ABS(TZ)) .LT. C MAX(EPSABS,EPSREL)*HALF ) TERMINATE, C WHERE C C EK = EXP(K*H+SHIFT) C ZK = (B*EK+A)/(EK+1) C FZ = F(ZK) C TZ = EK/((EK+1)*(EK+1)*(ZK-LAMBDA)) C TERM = TZ*FZ C THE QUANTITY "SHIFT" IS CHOSEN TO ENSURE THAT NO SINC POINT CAN C COINCIDE WITH THE SINGULARITY AT LAMBDA. C A SEQUENCE OF VALUES OF H IS USED, H <-- H/2, AND WHEN AGREEMENT C IS REACHED BETWEEN TWO APPROXIMATIONS, SO THAT C C ABS(OLD - NEW) < MAX(EPSABS,EPSREL*ABS(NEW)), C C THE ITERATION STOPS. C*********************************************************************** C********************* DECLARATIONS: ************************** C PARAMETERS: DOUBLE PRECISION F, A, B, H, LAMBDA, EPSREL, EPSABS, S INTEGER MAXITR, NFCALL, IER EXTERNAL F C LOCAL VARIABLES AND INTRINSIC FUNCTIONS USED: INTEGER N, NMAXP, NMAXN DOUBLE PRECISION ZN, PI, ZERO, ONE, HALF, FZ, TZ, TERM, OLD, * Q, E, EP1, MAX, TWO, SHIFT, AL, THREE, TH, * T2H INTEGER COUNT, NBAR, ADD INTRINSIC TAN, LOG, ACOS, ABS, EXP, MAX, MOD C********************* CONSTANTS DEFINED **************************** PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0, * TWO = 2.0D0, THREE = 3.0D0 ) C******************* FIRST EXECUTABLE STATEMENT: ********************** C FIRST, CHOOSE THE INITIAL H. PI = ACOS(-ONE) IF ( EPSREL .GT. ZERO ) THEN H = TWO*PI*PI/(-LOG(EPSREL)) ELSE H = TWO*PI*PI/(-LOG(EPSABS)) ENDIF C NOW CHOOSE THE SHIFT VALUE TO ENSURE THAT EVALUATION POINTS C AVOID THE SINGULARITY: Q = LOG((LAMBDA-A)/(B-LAMBDA)) IF ( MOD(MAXITR,2) .NE. 0 ) THEN AL = (ONE + TWO**(MAXITR+2))/THREE ELSE AL = (ONE + TWO**(MAXITR+1))/THREE ENDIF SHIFT = Q/H - AL/TWO**(MAXITR+1) NBAR = SHIFT SHIFT = (SHIFT - NBAR)*H C INITIALIZATIONS PRIOR TO BEGINNING THE ITERATION: S = ZERO T2H = ZERO COUNT = 0 NFCALL = 0 C THE MAIN LOOP PERFORMING MAXITR REFINEMENTS: 5 CONTINUE COUNT = COUNT + 1 IF ( COUNT .GT. MAXITR ) THEN C TOO MANY ITERATIONS. RETURN WITH IER = 1. IER = 1 GO TO 30 ENDIF OLD = S TH = T2H C FIRST, SUM OVER POSITIVE INDICES: N = -1 IF ( COUNT .EQ. 1 ) THEN ADD = 1 ELSE ADD = 2 ENDIF 10 CONTINUE N = N + ADD E = EXP(N*H+SHIFT) EP1 = E+ONE ZN = (B*E+A)/EP1 FZ = F(ZN) NFCALL = NFCALL + 1 TZ = E/EP1 TZ = TZ/(EP1*(ZN-LAMBDA)) TERM = TZ*FZ C ON THE FIRST PASS, WITH THE INITIAL H, CHECK THE TERMS C TO SEE WHEN TO STOP SUMMING. C STORE THE VALUE OF N WHERE THE TERM IS SMALL ENOUGH C IN NMAXP. C ON SUBSEQUENT PASSES, STOP WHEN N > 2*NMAXP - 1. IF ( COUNT .EQ. 1 ) THEN IF ( H*(B-A)*(ABS(TERM)+ABS(TZ)) * .LT. * MAX(EPSABS,EPSREL)*HALF ) THEN NMAXP = N - 1 GO TO 15 ENDIF ELSE IF ( N .GT. 2*NMAXP-1) THEN TH = TH + TERM NMAXP = 2*NMAXP + 1 GO TO 15 ENDIF ENDIF TH = TH + TERM GO TO 10 15 CONTINUE C NOW SUM OVER NEGATIVE INDICES. IF ( COUNT .EQ. 1 ) THEN N = 0 ELSE N = 1 ENDIF 20 CONTINUE N = N - ADD E = EXP(N*H+SHIFT) EP1 = E+ONE ZN = (B*E+A)/EP1 FZ = F(ZN) NFCALL = NFCALL + 1 TZ = E/EP1 TZ = TZ/(EP1*(ZN-LAMBDA)) TERM = TZ*FZ C ON THE FIRST PASS, WITH THE INITIAL H, CHECK THE TERMS C TO SEE WHEN TO STOP SUMMING. C STORE THE VALUE OF N WHERE THE TERM IS SMALL ENOUGH C IN NMAXN. C ON SUBSEQUENT PASSES, STOP WHEN N < 2*NMAX + 1. IF ( COUNT .EQ. 1 ) THEN IF ( H*(B-A)*(ABS(TERM)+ABS(TZ)) * .LT. * MAX(EPSABS,EPSREL)*HALF ) THEN NMAXN = N + 1 GO TO 25 ENDIF ELSE IF ( N .LT. 2*NMAXN+1) THEN TH = TH + TERM NMAXN = 2*NMAXN - 1 GO TO 25 ENDIF ENDIF TH = TH + TERM GO TO 20 25 CONTINUE S = TH*H*(B-A) S = S + * PI*F(LAMBDA)/TAN((PI/H)*(LOG((LAMBDA-A)/(B-LAMBDA)) * -SHIFT)) NFCALL = NFCALL + 1 IF ( EPSREL .EQ. ZERO .AND. ABS(S) .GT. ONE ) THEN EPSREL = EPSABS ENDIF IF ( ABS(OLD - S ) .GT. MAX(EPSABS,EPSREL*ABS(S)) ) THEN H = H*HALF T2H = TH GO TO 5 ELSE IER = 0 GO TO 30 ENDIF 30 CONTINUE EPSABS = ABS(OLD-S) EPSREL = ABS((OLD-S)/S) RETURN END