mean_obliquity_of_ecliptic_iau2006 Function

    pure function mean_obliquity_of_ecliptic_iau2006(et) result(e)

    implicit none

    real(wp),intent(in) :: et  !! ephemeris time (sec)
    real(wp)            :: e   !! obliquity of ecliptic (deg)

    real(wp) :: t  !! time in centuries from the J2000 epoch

    real(wp),parameter,dimension(6) :: c = [84381.406_wp,&
                                            -46.836769_wp,&
                                            -0.0001831_wp,&
                                            0.00200340_wp,&
                                            -0.000000576_wp,&
                                            -0.0000000434_wp] !! coefficients

    ! convert input time to centuries:
    t = et*sec2day*day2century

    ! use horner's rule:
    e = (c(1)+t*(c(2)+t*(c(3)+t*(c(4)+t*(c(5)+t*c(6))))))*arcsec2deg

    end function mean_obliquity_of_ecliptic_iau2006