safe_sin_polynomial Function

public elemental function safe_sin_polynomial(x) result(y)

Public only for the unit-test suite (no production module imports it); ignore when developing production code in other modules. sin(x) via Cody-Waite reduction mod π/2 and a 9-term Horner Taylor on the reduced argument. The quadrant k mod 4 picks sin vs cos and a sign:

k mod 4 = 0: sin(x) = sin(r) k mod 4 = 1: sin(x) = cos(r) k mod 4 = 2: sin(x) = -sin(r) k mod 4 = 3: sin(x) = -cos(r)

Cody-Waite splitting of π/2 keeps the reduction accurate for |x| ≲ 2¹⁷. Beyond that argument-precision starts to erode; callers reducing huge arguments (e.g. multi-year tidal phases) should pre-reduce.

sin(0) == 0 exact (the polynomial collapses at r=0 to r·SIN_C0 = 0).

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: x

Return Value real(kind=wp)


Calls

proc~~safe_sin_polynomial~~CallsGraph proc~safe_sin_polynomial safe_sin_polynomial proc~cos_reduced cos_reduced proc~safe_sin_polynomial->proc~cos_reduced proc~sin_reduced sin_reduced proc~safe_sin_polynomial->proc~sin_reduced

Variables

Type Visibility Attributes Name Initial
integer, private :: k_int
real(kind=wp), private :: k_real
integer, private :: quadrant
real(kind=wp), private :: r

Source Code

   elemental function safe_sin_polynomial(x) result(y)
      !! Public only for the unit-test suite (no production module imports it);
      !! ignore when developing production code in other modules.
      !! `sin(x)` via Cody-Waite reduction mod π/2 and a 9-term
      !! Horner Taylor on the reduced argument.  The quadrant
      !! `k mod 4` picks sin vs cos and a sign:
      !!
      !!   k mod 4 = 0:  sin(x) =  sin(r)
      !!   k mod 4 = 1:  sin(x) =  cos(r)
      !!   k mod 4 = 2:  sin(x) = -sin(r)
      !!   k mod 4 = 3:  sin(x) = -cos(r)
      !!
      !! Cody-Waite splitting of π/2 keeps the reduction accurate
      !! for |x| ≲ 2¹⁷.  Beyond that argument-precision starts to
      !! erode; callers reducing huge arguments (e.g. multi-year
      !! tidal phases) should pre-reduce.
      !!
      !! `sin(0) == 0` exact (the polynomial collapses at r=0 to
      !! r·SIN_C0 = 0).
      real(wp), intent(in) :: x
      real(wp) :: y
      real(wp) :: k_real, r
      integer :: k_int, quadrant

      k_real = anint(x*INV_PI_OVER_2)
      k_int = nint(k_real)
      r = (x - k_real*PI_OVER_2_HI) - k_real*PI_OVER_2_LO
      quadrant = modulo(k_int, 4)

      select case (quadrant)
      case (0)
         y = sin_reduced(r)
      case (1)
         y = cos_reduced(r)
      case (2)
         y = -sin_reduced(r)
      case default  ! 3
         y = -cos_reduced(r)
      end select
   end function safe_sin_polynomial