FUNCTION dcosint(xvalue) RESULT(fn_val) ! This program calculates the value of the cosine-integral ! defined as ! DCOSINT(x) = Gamma + Ln(x) + Integral (0 to x) [cos(t)-1]/t dt ! where Gamma is Euler's constant. ! The code uses rational approximations with a maximum accuracy of 20sf. ! INPUT PARAMETER: ! XVALUE - DOUBLE PRECISION - The argument to the function ! MACHINE-DEPENDENT PARAMETERS: ! XLOW - DOUBLE PRECISION - The absolute value below which ! DCOSINT( x ) = gamma + LN(x) , ! to machine precision. ! The recommended value is SQRT(2*EPSNEG) ! XHIGH1 - DOUBLE PRECISION - The value above which ! DCOSINT(x) = sin(x)/x - cos(x)/x^2 ! to machine precision. ! The recommended value is SQRT(6/EPSNEG) ! XHIGH2 - DOUBLE PRECISION - The value above which the trig. functions ! cannot be accurately determined. ! The value of the function is ! DCOSINT(x) = 0.0 ! The recommended value is pi/EPS. ! Values of EPS and EPSNEG for certain machine/compiler ! combinations can be found in the paper ! W.J. CODY Algorithm 665: MACHAR: A subroutine to dynamically ! determine machine parameters, ACM Trans. Math. Soft. 14 (1988) 303-311. ! The current code gives numerical values for XLOW,XHIGH1,XHIGH2 ! suitable for machines whose arithmetic conforms to the IEEE ! standard. The codes will probably work on other machines ! but might overflow or underflow for certain arguments. ! AUTHOR: Allan MacLeod ! Dept. of Mathematics and Statistics ! University of Paisley ! Scotland ! (e-mail: macl_ms0@paisley.ac.uk) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) REAL (dp), INTENT(IN) :: xvalue REAL (dp) :: fn_val ! Local variables INTEGER :: i REAL (dp) :: cx, dif, fival, gival, logval, root, sum, sumden, sumnum, & sx, x, xsq ! DATA VALUES REAL (dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp, three = 3.0_dp, & six = 6.0_dp, twelve = 12.0_dp REAL (dp), PARAMETER :: gamma = 0.57721566490153286060_dp REAL (dp), PARAMETER :: logl1 = 0.046875_dp, logl2 = 0.25_dp ! MACHINE-DEPENDENT PARAMETERS (SUITABLE FOR IEEE MACHINES) REAL (dp), PARAMETER :: xlow = 1.48996E-8_dp, xhigh1 = 2.324953E8_dp, & xhigh2 = 1.4148475E16_dp ! VALUES FOR COS-INTEGRAL FOR 0 < X <= 3 REAL (dp), PARAMETER :: ac1n(0:5) = (/ -0.24607411378767540707_dp, & 0.72113492241301534559E-2_dp, -0.11867127836204767056E-3_dp, & 0.90542655466969866243E-6_dp, -0.34322242412444409037E-8_dp, & 0.51950683460656886834E-11_dp /) REAL (dp), PARAMETER :: ac1d(0:5) = (/ 1.0_dp, 0.12670095552700637845E-1_dp, & 0.78168450570724148921E-4_dp, 0.29959200177005821677E-6_dp, & 0.73191677761328838216E-9_dp, 0.94351174530907529061E-12_dp /) ! VALUES FOR COS-INTEGRAL FOR 3 < X <= 6 REAL (dp), PARAMETER :: ac2n(0:7) = (/ -0.15684781827145408780_dp, & 0.66253165609605468916E-2_dp, -0.12822297297864512864E-3_dp, & 0.12360964097729408891E-5_dp, -0.66450975112876224532E-8_dp, & 0.20326936466803159446E-10_dp, -0.33590883135343844613E-13_dp, & 0.23686934961435015119E-16_dp /) REAL (dp), PARAMETER :: ac2d(0:6) = (/ 1.0_dp, 0.96166044388828741188E-2_dp, & 0.45257514591257035006E-4_dp, 0.13544922659627723233E-6_dp, & 0.27715365686570002081E-9_dp, 0.37718676301688932926E-12_dp, & 0.27706844497155995398E-15_dp /) ! VALUES FOR FI(X) FOR 6 <= X <= 12 REAL (dp), PARAMETER :: afn1(0:7) = (/ 0.99999999962173909991E0_dp, & 0.36451060338631902917E3_dp, 0.44218548041288440874E5_dp, & 0.22467569405961151887E7_dp, 0.49315316723035561922E8_dp, & 0.43186795279670283193E9_dp, 0.11847992519956804350E10_dp, & 0.45573267593795103181E9_dp /) REAL (dp), PARAMETER :: afd1(0:7) = (/ 1.0_dp, 0.36651060273229347594E3_dp, & 0.44927569814970692777E5_dp, 0.23285354882204041700E7_dp, & 0.53117852017228262911E8_dp, 0.50335310667241870372E9_dp, & 0.16575285015623175410E10_dp, 0.11746532837038341076E10_dp /) ! VALUES OF GI(X) FOR 6 <= X <=12 REAL (dp), PARAMETER :: agn1(0:8) = (/ 0.99999999920484901956E0_dp, & 0.51385504875307321394E3_dp, 0.92293483452013810811E5_dp, & 0.74071341863359841727E7_dp, 0.28142356162841356551E9_dp, & 0.49280890357734623984E10_dp, 0.35524762685554302472E11_dp, & 0.79194271662085049376E11_dp, 0.17942522624413898907E11_dp /) REAL (dp), PARAMETER :: agd1(0:8) = (/ 1.0_dp, 0.51985504708814870209E3_dp, & 0.95292615508125947321E5_dp, 0.79215459679762667578E7_dp, & 0.31977567790733781460E9_dp, 0.62273134702439012114E10_dp, & 0.54570971054996441467E11_dp, 0.18241750166645704670E12_dp, & 0.15407148148861454434E12_dp /) ! VALUES FOR FI(X) FOR X > 12 REAL (dp), PARAMETER :: afn2(0:7) = (/ 0.19999999999999978257E1_dp, & 0.22206119380434958727E4_dp, 0.84749007623988236808E6_dp, & 0.13959267954823943232E9_dp, 0.10197205463267975592E11_dp, & 0.30229865264524075951E12_dp, 0.27504053804288471142E13_dp, & 0.21818989704686874983E13_dp /) REAL (dp), PARAMETER :: afd2(0:7) = (/ 1.0_dp, 0.11223059690217167788E4_dp, & 0.43685270974851313242E6_dp, 0.74654702140658116258E8_dp, & 0.58580034751805687471E10_dp, 0.20157980379272098841E12_dp, & 0.26229141857684496445E13_dp, 0.87852907334918467516E13_dp /) ! VALUES FOR GI(X) FOR X > 12 REAL (dp), PARAMETER :: agn2(0:8) = (/ 0.59999999999999993089E1_dp, & 0.96527746044997139158E4_dp, 0.56077626996568834185E7_dp, & 0.15022667718927317198E10_dp, 0.19644271064733088465E12_dp, & 0.12191368281163225043E14_dp, 0.31924389898645609533E15_dp, & 0.25876053010027485934E16_dp, 0.12754978896268878403E16_dp /) REAL (dp), PARAMETER :: agd2(0:8) = (/ 1.0_dp, 0.16287957674166143196E4_dp, & 0.96636303195787870963E6_dp, 0.26839734750950667021E9_dp, & 0.37388510548029219241E11_dp, 0.26028585666152144496E13_dp, & 0.85134283716950697226E14_dp, 0.11304079361627952930E16_dp, & 0.42519841479489798424E16_dp /) ! VALUES FOR AN APPROXIMATION TO LN(X/ROOT) REAL (dp), PARAMETER :: p(0:2) = (/ 0.83930008362695945726E1_dp, & -0.65306663899493304675E1_dp, 0.569155722227490223_dp /) REAL (dp), PARAMETER :: q(0:1) = (/ 0.41965004181347972847E1_dp, & -0.46641666676862479585E1_dp /) ! VALUES OF THE FIRST TWO ROOTS OF THE COSINE-INTEGRAL REAL (dp), PARAMETER :: rt1n = 631.0_dp, rt1d = 1024.0_dp, & rt1r = 0.29454812071623379711E-3_dp REAL (dp), PARAMETER :: rt2n = 3465.0_dp, rt2d = 1024.0_dp, & rt2r = 0.39136005118642639785E-3_dp ! START COMPUTATION x = xvalue IF ( x <= zero ) THEN fn_val = zero RETURN END IF IF ( x <= six ) THEN ! CODE FOR 3 < X < = 6 IF ( x > three ) THEN sumnum = zero sumden = zero xsq = x * x DO i = 7 , 0 , -1 sumnum = sumnum * xsq + ac2n( i ) END DO DO i = 6 , 0 , -1 sumden = sumden * xsq + ac2d( i ) END DO root = rt2n / rt2d dif = ( x - root ) - rt2r sum = root + rt2r IF ( ABS(dif) < logl2 ) THEN cx = dif / ( sum + x ) xsq = cx * cx sx = p(0) + xsq * ( p(1) + xsq * p(2) ) sx = sx / ( q(0) + xsq * ( q(1) + xsq ) ) logval = cx * sx ELSE logval = LOG( x / sum ) END IF fn_val = logval + dif * ( x + sum ) * sumnum / sumden ELSE ! CODE FOR 0 < X < = 3 IF ( x > xlow ) THEN sumnum = zero sumden = zero xsq = x * x DO i = 5 , 0 , -1 sumnum = sumnum * xsq + ac1n( i ) sumden = sumden * xsq + ac1d( i ) END DO root = rt1n / rt1d dif = ( x - root ) - rt1r sum = root + rt1r IF ( ABS(dif) < logl1 ) THEN cx = dif / ( sum + x ) xsq = cx * cx sx = p(0) + xsq * ( p(1) + xsq * p(2) ) sx = sx / ( q(0) + xsq * ( q(1) + xsq ) ) logval = cx * sx ELSE logval = LOG( x / sum ) END IF fn_val = logval + dif * ( x + sum ) * sumnum / sumden ELSE fn_val = gamma + LOG( x ) END IF END IF END IF ! CODE FOR 6 < X < = 12 IF ( x > six .AND. x <= twelve ) THEN sumnum = zero sumden = zero xsq = one / ( x * x ) DO i = 7 , 0 , -1 sumnum = sumnum * xsq + afn1( i ) sumden = sumden * xsq + afd1( i ) END DO fival = sumnum / ( x * sumden ) sumnum = zero sumden = zero DO i = 8 , 0 , -1 sumnum = sumnum * xsq + agn1( i ) sumden = sumden * xsq + agd1( i ) END DO gival = xsq * sumnum / sumden fn_val = fival * SIN( x ) - gival * COS( x ) END IF ! CODE FOR X > 12 IF ( x > twelve ) THEN IF ( x > xhigh2 ) THEN fn_val = zero ELSE cx = COS( x ) sx = SIN( x ) xsq = one / ( x * x ) IF ( x > xhigh1 ) THEN fn_val = sx / x - cx * xsq ELSE sumnum = zero sumden = zero DO i = 7 , 0 , -1 sumnum = sumnum * xsq + afn2( i ) sumden = sumden * xsq + afd2( i ) END DO fival = ( one - xsq * sumnum / sumden ) / x sumnum = zero sumden = zero DO i = 8 , 0 , -1 sumnum = sumnum * xsq + agn2( i ) sumden = sumden * xsq + agd2( i ) END DO gival = ( one - xsq * sumnum / sumden ) * xsq fn_val = sx * fival - cx * gival END IF END IF END IF RETURN END FUNCTION dcosint