icrf_to_iau_moon Function

    pure function icrf_to_iau_moon(et) result(rotmat)

    implicit none

    real(wp),intent(in)     :: et        !! ephemeris time [sec]
    real(wp),dimension(3,3) :: rotmat    !! the rotation matrix

    real(wp) :: w,dec,ra,d,t
    real(wp) :: e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13

    d = et * sec2day     !interval in days from the J2000 epoch
    t = et * sec2century !interval in Julian centuries (36525 days) from the J2000 epoch

    ! Nutation/precession angles:
    e1  = ( 125.045_wp - 0.0529921_wp  * d ) * deg2rad
    e2  = ( 250.089_wp - 0.1059842_wp  * d ) * deg2rad
    e3  = ( 260.008_wp + 13.0120009_wp * d ) * deg2rad
    e4  = ( 176.625_wp + 13.3407154_wp * d ) * deg2rad
    e5  = ( 357.529_wp + 0.9856003_wp  * d ) * deg2rad
    e6  = ( 311.589_wp + 26.4057084_wp * d ) * deg2rad
    e7  = ( 134.963_wp + 13.0649930_wp * d ) * deg2rad
    e8  = ( 276.617_wp + 0.3287146_wp  * d ) * deg2rad
    e9  = ( 34.226_wp  + 1.7484877_wp  * d ) * deg2rad
    e10 = ( 15.134_wp  - 0.1589763_wp  * d ) * deg2rad
    e11 = ( 119.743_wp + 0.0036096_wp  * d ) * deg2rad
    e12 = ( 239.961_wp + 0.1643573_wp  * d ) * deg2rad
    e13 = ( 25.053_wp  + 12.9590088_wp * d ) * deg2rad

    ra  = ( 269.9949_wp + 0.0031_wp * t &
            - 3.8787_wp * sin(E1 ) &
            - 0.1204_wp * sin(E2 ) &
            + 0.0700_wp * sin(E3 ) &
            - 0.0172_wp * sin(E4 ) &
            + 0.0072_wp * sin(E6 ) &
            - 0.0052_wp * sin(E10) &
            + 0.0043_wp * sin(E13) ) * deg2rad

    dec = ( 66.5392_wp + 0.0130_wp * t &
            + 1.5419_wp * cos(E1 ) &
            + 0.0239_wp * cos(E2 ) &
            - 0.0278_wp * cos(E3 ) &
            + 0.0068_wp * cos(E4 ) &
            - 0.0029_wp * cos(E6 ) &
            + 0.0009_wp * cos(E7 ) &
            + 0.0008_wp * cos(E10) &
            - 0.0009_wp * cos(E13) ) * deg2rad

    w = ( 38.3213_wp + 13.17635815_wp * d - 1.4e-12_wp * d**2 &
          + 3.5610_wp * sin(E1 ) &
          + 0.1208_wp * sin(E2 ) &
          - 0.0642_wp * sin(E3 ) &
          + 0.0158_wp * sin(E4 ) &
          + 0.0252_wp * sin(E5 ) &
          - 0.0066_wp * sin(E6 ) &
          - 0.0047_wp * sin(E7 ) &
          - 0.0046_wp * sin(E8 ) &
          + 0.0028_wp * sin(E9 ) &
          + 0.0052_wp * sin(E10) &
          + 0.0040_wp * sin(E11) &
          + 0.0019_wp * sin(E12) &
          - 0.0044_wp * sin(E13) ) * deg2rad

    rotmat = iau_rotation_matrix(w,dec,ra)

    end function icrf_to_iau_moon