diff --git a/source/f90/o.f90 b/source/f90/o.f90 new file mode 100644 index 0000000000000000000000000000000000000000..a9f4545dc2766b627f5ff8f208c71440d2a623c5 --- /dev/null +++ b/source/f90/o.f90 @@ -0,0 +1,23 @@ +double precision function o(x) + +! *** rational approximations for ei(x) * exp(-x) + e1(x) * exp(+x) +! *** in the intervals (0, 1.3) and (1.3, infinity) <accuracy : 4.e-04> +! *** used for fourier transforming half-lorentzian functions + + implicit none + + double precision, intent(in) :: x + + if ( x < 0 ) then + error stop '*** Program abort. Function o is undefined for x < 0 ***' + elseif ( x == 0.0d0 ) then + o = 0.0d0 + elseif ( x <= 1.3d0) then + o = -dsinh(x) * dlog(x**2) + x * ((0.03114d0 * x**2 + 0.41666d0) * x**2 + 0.84557d0) + else + o = (((202.91d0 / x**2 + 932.21d0) / x**2 + 41.740d0) / x**2 + 2.0d0) / & + (((540.88d0 / x**2 + 345.67d0) / x**2 + 18.961d0) / x**2 + 1.0d0) / x + endif + + return +end function o diff --git a/tests/oxlog/o1.f90 b/tests/oxlog/o1.f90 index 16c9bf04a05f45027c587d7847e544f00e09c88d..a9c9ea77b32378179160b7f175066d4a37e51d13 100644 --- a/tests/oxlog/o1.f90 +++ b/tests/oxlog/o1.f90 @@ -12,5 +12,7 @@ double precision function o1(u) u2 = u**2 o1 = -dsinh(u) * dlog(u2) + u * ((0.03114d0 * u2 + 0.41666d0) * u2 + 0.84557d0) + write (*,*) u, o1 + return end function o1 diff --git a/tests/oxlog/o2.f90 b/tests/oxlog/o2.f90 index 76c562b2d807f2e00fe70f38ca4045f44023bd2b..463402e4fe43ded715ea410c47d2313f539c7d97 100644 --- a/tests/oxlog/o2.f90 +++ b/tests/oxlog/o2.f90 @@ -12,6 +12,6 @@ double precision function o2(u) u2 = 1 / u**2 o2 = (((202.91d0 * u2 + 932.21d0) * u2 + 41.740d0) * u2 + 2.0d0) / & (((540.88d0 * u2 + 345.67d0) * u2 + 18.961d0) * u2 + 1.0d0) / u - + write (*,*) u, o2 return end function o2