subroutine fser1 ( f, eps, max, m, c, s, lc, ls ) c*********************************************************************72 c cc fser1 c implicit none real a real b real b1 real b2 real b3 real b4 real b5 real c real c1 real cosip real eps real f real f0 real f1 real frac real g real h integer i integer in integer lc integer llc integer lls integer ln integer ls integer m integer max integer mswtch integer n integer nst integer nstop real odcos real odsin real p real pi real pvt1 real pvt2 real s real s1 real s1sq real sumcos real sumsin real t real t1 real t2 real temp1 real temp2 real tha real tmax real tp real tsq real xi real xm pi = 3.1415926535898e+00 xm = m c c f1 = cos ( m * pi ) temporary. c f1 = 1.0e+00 - 2.0e+00 * ( m - ( m / 2 ) * 2 ) f0 = f ( 0.0e+00 ) f1 = f ( 1.0e+00 ) * f1 c c 'cir' will be used throughout these comments to stand for 'sin' or c 'cos' whenever those two symbols may occur. c now define sumcir of the endpoints. c sumcos = ( f1 + f0 ) * 0.5e+00 sumsin = 0.0e+00 b1 = 2.0e+00 / 3.0e+00 c c tmax is the switch-over point in the angle t. c our analysis indicates that tmax = 1/6 is the best for the iliac ii, c which has a 44 bit floating point mantissa. c tmax = 0.166e+00 c c n is the number of iterations. c n = int ( alog ( 2.0e+00 * xm ) / 0.693e+00 ) c c both tmax and n may be changed if the machine for which this c routine is intended has greater or less accuracy than iliac ii. c if n is changed, then the corresponding changes must be made c in the assignments of h and nstop. c h = 1.0e+00 / float ( 2**n ) c c h = 1 / 2**n. c nstop = 2**n - 1 c c nstop = 2**n - 1. c t = h * xm tp = t * pi nst = 1 assign 67 to mswtch c c llc and lls are used by the routine in computed-go-to statements. c as soon as lls and llc have been defined, we can use ls and lc c as return parameters (see above). c if ( ls ) 1, 1, 2 1 lls = 2 go to 3 2 lls = 1 ls = max 3 if ( lc ) 4, 4, 5 4 llc = 2 go to 7 5 llc = 1 lc = max 7 ln = 1 c c all of the above is executed only once per call. c now the iteration begins. c 10 odcos = 0.0e+00 odsin = 0.0e+00 c c begin summation for odcos and odsin. c do 65 i = 1, nstop, nst xi = i tha = xi * t c c tha * pi is the angle used in this i-th term. c c cir(i*t*pi) is calculated here using the identity c cir ( integer multiple of pi + fractional multiple of pi ) c = cos ( integer * pi ) * cir ( frac * pi ) c = ( + 0r - ) * cir ( frac * pi ). c frac = tha in = int ( tha ) tha = in frac = ( frac - tha ) * pi c c tha is a floating point integer, frac is the fractional part * pi. c cosip = 1 - 2 * ( in - 2 * ( in / 2 ) ) temp1 = cosip * f ( xi * h ) c c temp1 = cos ( integer part ) * f ( i * h ). c go to ( 50, 55 ) lls 50 odsin = temp1 * sin ( frac ) + odsin 55 go to ( 60, 65 ) llc 60 odcos = temp1 * cos ( frac ) + odcos 65 continue go to mswtch, ( 67, 70 ) 67 nst = 2 c c now have made up for the first 4 iteration steps, so reset these c two numbers to look like the general case. c nstop = 2**n c c nstop = 2**n, in case you change starting value of n. c assign 70 to mswtch go to 92 70 continue tsq = tp * tp if ( t - tmax ) 74, 74, 75 c c 74 is the power series for small t, 75 is the closed form used with c larger values of t. c the power series (with 'tn' = tp**n) c a = (2/45)*t3 - (2/315)*t5 + (2/4725)*t7 c b = (2/3) + (2/15)*t2 - (4/105)*t4 + (2/567)*t6 - (4/22275)*t8 c g = (4/3) - (2/15)*t2 - (1/210)*t4 + (1/11340)*t6 c the next term in g is too small. it is (1/997920)*t8. c 74 a = tp * tsq * ( 1.0e+00 - tsq * ( 1.0e+00 - tsq / 15.0e+00 ) & / 7.0e+00 ) / 22.5e+00 b2 = b1 * tsq * 0.2e+00 b3 = b2 * tsq * 2.0e+00 / 7.0e+00 b4 = b3 * tsq / 10.8e+00 b5 = b4 * tsq * 14.0e+00 / 275.0e+00 b = b1 + b2 - b3 + b4 - b5 g = 2.0e+00 * b1 - b2 + b3 / 8.0e+00 - b4 / 40.0e+00 c c g = 2*b1 - b2 + b3/8 - b4/40 + 5*b5/896. if you want the t8 c term included in g. c go to 80 c c closed form of the coefficients, where again 'tn' means tp**n. c a = 1/tp + cos(tp)*sin(tp)/t2 - 2*(sin(tp))**2/t3 c b = 2*((1+(cos(tp))**2)/t2 - 2*sin(tp)*cos(tp)/t3). c g = 4*( sin(tp)/t3 - cos(tp)/t2 ). c 75 in = int ( t ) temp1 = 1 - 2 * ( in - 2 * ( in / 2 ) ) temp2 = in c c temp1 is cos ( integer part of tp), temp2 is fractional part of tp. c temp2 = ( t - temp2 ) * pi s1 = temp1 * sin ( temp2 ) c c s1 = sin(tp). c c1 = temp1 * cos ( temp2 ) c c c1 = cos(tp). c p = s1 * c1 s1sq = s1 * s1 a = ((( -2.0e+00 * s1sq / tp ) + p ) / tp + 1.0e+00 ) / tp b = 2.0e+00 * ( ( -2.0e+00 * p / tp ) + 2.0e+00 - s1sq ) / tsq g = 4.0e+00 * ( s1 / tp - c1 ) / tsq 80 go to ( 81, 85 ) lls c c have calculated the coefficients. now ready for the integration c formulas. c 81 t2 = h * ( a * ( f0 - f1 ) + b * sumsin + g * odsin ) c c endt1 is a subroutine which checks for the convergence of the c iterations. endt1 requires the present value to agree with the c previous value to within eps2, where c eps2 = ( 1.0 + abs ( present value ) * eps. c eps is supplied by the user. c call endt1 ( pvt2, t2, eps, s, lls, ln ) go to ( 85, 84 ), lls 84 ls = n 85 go to ( 86, 90 ), llc c c this is the cosine integral. c 86 t1 = h * ( b * sumcos + g * odcos ) call endt1 ( pvt1, t1, eps, c, llc, ln ) go to ( 90, 89 ), llc 89 lc = n 90 ln = 2 c c now test to see if done. c if ( llc + lls - 3 ) 92, 92, 100 92 n = n + 1 c c this is the beginning of the iteration. c if ( n - max ) 95, 95, 100 95 h = 0.5e+00 * h t = 0.5e+00 * t tp = 0.5e+00 * tp nstop = 2 * nstop sumsin = sumsin + odsin sumcos = sumcos + odcos go to 10 100 s = t2 c = t1 return end subroutine endt1 ( prevqt, quant, eps, value, l1, l2 ) c*********************************************************************72 c cc endt1 c implicit none real eps integer l1 integer l2 real prevqt real quant real reps real value go to ( 29, 20 ), l2 20 reps = eps * ( 1.0e+00 + abs ( quant ) ) if ( abs ( prevqt - quant ) - reps ) 25, 25, 29 25 value = quant l1 = 2 go to 30 29 prevqt = quant 30 return end