/************************************************************************/ /* */ /* ALGORITHM 353 - FILON QUADRATURE */ /* */ /* Stephen M. Chase and Lloyd D. Fosdick 7 July 1967. */ /* With modifications suggested by Bo Einarsson, */ /* Research Institute of National Defence, Tumba, Sweden, */ /* 8 December 1969. */ /* */ /************************************************************************/ /* Translated to Fortran 77 by R. A. Vowels. 27 June 2006. */ /* Translated to Fortran 90 by R. A. Vowels, 13 May 2007. */ /* Translated to PL/I by R. A. Vowels, 13 May 2007. */ /* Fortran 77 version works 27 June 2006. */ /* PL/I version works 2 June 2007. */ /* FSER1 evaluates the integrals */ /* 1 1 */ /* C = S F(X) cos (M pi X) dX, S = S F(X) sin (M pi X) dX */ /* 0 0 */ /* using the Filon quadrature algorithm. */ /* The user may request an evaluation of C only, S only, or both C and S. */ /* With the function F(x) = x**2*(1-x), and M = 64, MAX = 20, LC = 1, */ /* LS = 1 and EPS = 1e-10: */ /* Exact solutions are : */ /* C = -1/(64*pi)**2 = -2.473661710e-5 */ /* S = -6/(64*pi)**2 = -7.381790413e-7 */ /* Actual results obtained were: */ /* C = -2.473661709349E-05 */ /* S = -7.381790410980E-07 */ /* LC = 9 LS = 9 */ /* With EPS changed to 1E-20, the actual results were: */ /* C = -2.473661710018E-05 */ /* S= -7.381790412808E-07 */ /* LC = 20 LS = 20 */ (SUBSCRIPTRANGE, FIXEDOVERFLOW, SIZE): FSER1: PROCEDURE (F, EPS, MAX, M, C, S, LC, LS) OPTIONS (REORDER); /* F = name of a function to evaluate F(X), supplied by the user. */ /* EPS = error control for algorithm. */ /* MAX = maximum number of halvings of the step size that are allowed. */ /* M = parameter in formula. */ /* C = value of the cosine integral computed by this subroutine. */ /* S = value of the sine integral computed by this subroutine. */ /* LC = 1 if the cosine integral is to be computed (value returned via */ /* C). The value of step size used to compute C is returned in LC. */ /* (The step size h = 2**(-LC). */ /* = 0 if the cosine integral is not to be computed. */ /* LS = 1 if the sine integral is to be computed (value returned via S). */ /* The value of step size used to compute S is returned in LS. */ /* (The step size h = 2**(-LS). */ /* LS = 0 if the sine integral is not to be computed. */ DECLARE F ENTRY (FLOAT (18)) RETURNS (FLOAT (18)); DECLARE ( EPS, C, S ) FLOAT (18); DECLARE ( MAX, M, LC, LS ) FIXED BINARY (31); DECLARE ( PI, XM, F0, F1, SUMSIN, SUMCOS, B1, TMAX, H, T, TP, ODSIN, ODCOS, XI, THA, FRAC, COSIP, TEMP1, TEMP2, TSQ, A, B2, B3, B4, B5, B, G, S1, S1SQ, C1, P, PVT1, PVT2, T1, T2 ) FLOAT (18); DECLARE ( I, N, NSTOP, NST, LN, IN, M_ON_2, IN_ON_2 ) FIXED BINARY (31); DECLARE ( LLC, LLS, MSWITCH ) FIXED BINARY (31); PI = 3.14159265358979324E0; XM = M; M_ON_2 = M/2; F1 = 1 - 2 * (M - M_ON_2 * 2); F0 = F(0.00000000000000000E0); F1 = F(1.00000000000000000E0) * F1; SUMCOS = (F1 * F0)/2; SUMSIN = 0; B1 = 2.00000000000000000E0/3.00000000000000000E0; /* TMAX is the switchover point in the angle T. */ /* Our analysis indicates that TMAX = 1/6 is the best for a machine */ /* having a 44 bit floating-point mantissa. */ TMAX = 0.16600000000000000E0; /* N is the number of the iteration. Note that we start at the */ /* fourth iteration step. */ /* Actually, the first evaluation of an integral is at N = 5, and */ /* therefore, the first comparison of values is at N = 6. */ N = LOG(2*XM) / 0.69300000000000000E0; PUT SKIP LIST ( N); /* Both TMAX and N may be changed if the machine for which this */ /* routine is intended has greater or less accuracy than a 44-bit mantissa. */ /* If N is changed, then corresponding changes must be made */ /* in the assignments of H and NSTOP. */ H = 1/2.00000000000000000E0**N; NSTOP = ISLL(1, N) - 1; /* A simple way to do NSTOP = 2**N-1; */ T = H * XM; TP = T * PI; NST = 1; MSWITCH = 67; /* LLC and LLS are used by the routine in computed-go-to statements. */ /* As soon as LLS and LLC have been defined, we can use LS and LC */ /* as return values (see description of dummy arguments). */ IF LS <= 0 THEN LLS = 2; ELSE DO; LLS = 1; LS = MAX; END; IF LC <= 0 THEN LLC = 2; ELSE DO; LLC = 1; LC = MAX; END; LN = 1; /* All of the above is executed once per call. */ /* Now the iteration begins. */ L10: ODCOS, ODSIN = 0; /* Begin summation for ODCOS and ODSIN. */ DO I = 1 TO NSTOP BY NST; XI = I; THA = XI * T; /* THA*PI is the angle used in this I-th term. */ /* CIR(I*T*PI) is calculated here using the identity */ /* CIR (integer multiple of PI + fractional mult of PI ) */ /* = COS(integer*PI) * CIR(FRAC*PI) */ /* = (+ or -) * CIR(FRAC*PI). */ FRAC = THA; IN = THA; THA = IN; FRAC = (FRAC - THA) * PI; /* THA is a floating-point integer, FRAC is the fractional part * PI. */ IN_ON_2 = IN/2; COSIP = 1 - 2*(IN - 2*IN_ON_2); TEMP1 = COSIP * F(XI*H); /* TEMP1 = COS(integer part) * F(I*H). */ IF LLS = 2 THEN GO TO L55; ODSIN = TEMP1 * SIN (FRAC) + ODSIN; L55: IF LLC = 1 THEN ODCOS = TEMP1 * COS(FRAC) + ODCOS; END; IF MSWITCH = 67 THEN GO TO L67; ELSE GO TO L70; L67: NST = 2; /* Now have made up for the first 4 iteration steps, so reset these */ /* two numbers to look like the general case. */ NSTOP = ISLL(1, N); /* A simple way to get NSTOP = 2**N; */ MSWITCH = 70; GO TO L92; L70: TSQ = TP * TP; IF T <= TMAX THEN DO; /* Power series for small T. */ /* (The closed form is used with larger values of T; see the ELSE */ /* below.) */ /* The power series are (with 'TN' = TP ** N) */ /* A = (2/45)*T3 - (2/315)*T5 + (2/4725)*T7 */ /* B = (2/3) + (2/15)*T2 - (4/105)*T4 + (2/567)*T6 - (4/22275)*T8 */ /* G = (4/3) - (2/15)*T2 + (1/210)*T4 - (1/11340)*T6 */ A = TP * TSQ * (1 - TSQ * (1 - TSQ / 15) / 7.00000000000000000E0) / 22.5000000000000000E0; B2 = B1 * TSQ / 5; B3 = B2 * TSQ * 2 / 7.00000000000000000E0; B4 = B3 * TSQ / 10.8000000000000000E0; B5 = B4 * TSQ * 14.0000000000000000E0 / 275.000000000000000E0; B = B1 + B2 - B3 + B4 - B5; G = 2*B1 - B2 + B3/8 - B4/40 + 5*B5/896.000000000000000E0; GO TO L80; END; ELSE DO; /* Closed form of the coefficients, where again 'TN' means TP**N. */ /* A = 1/TP + COS(TP)*SIN(TP)/T2 - 2*(SIN(TP))**2/TP */ /* B = 2*((1 + (COS(TP))**2/T2 - 2*SIN(TP)*COS(TP)/T3) */ /* G = 4*(SIN(TP)/T3 - COS(TP)/T2) */ IN = T; IN_ON_2 = IN/2; TEMP1 = 1 - 2 * (IN - 2*IN_ON_2); TEMP2 = IN; /* TEMP1 is COS (integer part of TP), TEMP2 is the fractional part of TP. */ TEMP2 = (T - TEMP2) * PI; S1 = TEMP1 * SIN (TEMP2); C1 = TEMP1 * COS (TEMP2); P = S1 * C1; S1SQ = S1 * S1; A = (((-2*S1SQ/TP) + P)/ TP + 1) / TP; B = 2 * ((-2*P/TP) + 2 - S1SQ) / TSQ; G = 4 * (S1 / TP - C1) / TSQ; END; L80: IF LLS = 2 THEN GO TO L85; /* Have calculated the coefficients, now ready for the integration formulas. */ T2 = H * (A * (F0 - F1) + B * SUMSIN + G * ODSIN); /* ENDT1 is a subroutine that checks for the convergence of the */ /* previous value to within EPS2, where */ /* EPS2 = (1.0 + ABSF(present value))*EPS */ /* EPS is supplied by the user. */ CALL ENDT1 (PVT2, T2, EPS, S, LLS, LN); IF LLS = 1 THEN GO TO L85; L84: LS = N; L85: IF LLC = 2 THEN GO TO L90; /* This is the cosine integral. */ T1 = H * (B * SUMCOS + G * ODCOS); CALL ENDT1 (PVT1, T1, EPS, C, LLC, LN); IF LLC = 1 THEN GO TO L90; LC = N; L90: LN = 2; IF LLC + LLS > 3 THEN GO TO L100; L92: N = N + 1; IF N > MAX THEN GO TO L100; H = H / 2; T = T / 2; TP = TP / 2; NSTOP = 2 * NSTOP; SUMSIN = SUMSIN + ODSIN; SUMCOS = SUMCOS + ODCOS; GO TO L10; L100: S = T2; C = T1; RETURN; END FSER1; ENDT1: PROCEDURE (PREVQT, QUANT, EPS, VALUE, L1, L2) OPTIONS (REORDER); DECLARE ( PREVQT, QUANT, EPS, VALUE ) FLOAT (18); DECLARE ( L1, L2 ) FIXED BINARY (31); DECLARE ( REPS ) FLOAT (18); IF L2 = 1 THEN GO TO L29; REPS = EPS * (1.00000000000000000E0 + ABS (QUANT)); IF ABS(PREVQT - QUANT) <= REPS THEN DO; VALUE = QUANT; L1 = 2; RETURN; END; L29: PREVQT = QUANT; END ENDT1;