Accurate computation of s - (1 - g1)*asinh(s)
when (g1, s) is close to (0, 0)
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=wp), | intent(in) | :: | g1 | |||
real(kind=wp), | intent(in) | :: | s |
pure function shmkep (g1, s) implicit none real(wp) :: shmkep real(wp),intent(in) :: g1 real(wp),intent(in) :: s real(wp) :: g,t,tsq,x,term,twoi1,x0 g = one - g1 t = s/(one + sqrt(one + s*s)) tsq = t*t x = s*(g1 + g*tsq) term = two*g*t twoi1 = one do twoi1 = twoi1 + two term = term*tsq x0 = x x = x - term/twoi1 if (x==x0) exit end do shmkep = x end function shmkep