Skip to content

Instantly share code, notes, and snippets.

@certik

certik/a.f90

Created Oct 4, 2012
Embed
What would you like to do?
real(dp) function Inu_formula(k, x) result(r)
integer, intent(in) :: k
real(dp), intent(in) :: x
select case (k)
case (0)
r = esinh(x)
case (1)
if (x > 0.9_dp) then
r = -esinh(x)/x + ecosh(x)
else
r = x**2/3 + x**4/30 + x**6/840 + x**8/45360 + x**10/3991680 + &
x**12/518918400 + x**14/93405312e3_dp
r = r / exp(x)
end if
case (2)
if (x > 0.9_dp) then
r = (3/x**2 + 1)*esinh(x) - 3*ecosh(x)/x
else
r = x**3/15 + x**5/210 + x**7/7560 + x**9/498960 + &
x**11/51891840 + x**13/7783776e3_dp
r = r / exp(x)
end if
case (3)
if (x > 0.9_dp) then
r = -(15/x**3 + 6/x)*esinh(x) + (15/x**2 + 1)*ecosh(x)
else
r = x**4/105 + x**6/1890 + x**8/83160 + x**10/6486480 + &
x**12/778377600 + x**14/132324192e3_dp
r = r / exp(x)
end if
case (4)
if (x > 0.9_dp) then
r = (105/x**4 + 45/x**2 + 1)*esinh(x) - (105/x**3 + 10/x)*ecosh(x)
else
r = x**5/945 + x**7/20790 + x**9/1081080 + x**11/97297200 + &
x**13/132324192e2_dp
r = r / exp(x)
end if
case default
call stop_error("k = " // str(k) // " not implemented.")
end select
r = r * sqrt(2/(pi*x))
end function
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment