Accurate computation of ee - e*sin(ee)
when (e, ee) is close to (1, 0)
Note
must not be used for large ee (absolute) as then rounding worse not better
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=wp), | intent(in) | :: | e | |||
real(kind=wp), | intent(in) | :: | ee |
pure function emkepl(e, ee) implicit none real(wp) :: emkepl real(wp),intent(in) :: e real(wp),intent(in) :: ee real(wp) :: x, ee2, term, d, x0 x = (one - e)*sin(ee) ee2 = -ee*ee term = ee d = zero do d = d + two term = term*ee2/(d*(d + one)) x0 = x x = x - term if (x==x0) exit end do emkepl = x end function emkepl