add oxlog

source/f90/o.f90

+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

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

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