diff --git a/pFUnittests/CMakeLists.txt b/pFUnittests/CMakeLists.txt
index 121e05a..39df6b6 100644
--- a/pFUnittests/CMakeLists.txt
+++ b/pFUnittests/CMakeLists.txt
@@ -40,3 +40,7 @@ add_pfunit_ctest (slam_error_handling_test
add_pfunit_ctest (slam_io_test
TEST_SOURCES test_slam_io.pf
LINK_LIBRARIES slam-Fortran)
+
+add_pfunit_ctest (slam_moon_reduction_test
+ TEST_SOURCES test_moon_reduction.pf
+ LINK_LIBRARIES slam-Fortran)
diff --git a/pFUnittests/test_moon_reduction.pf b/pFUnittests/test_moon_reduction.pf
new file mode 100644
index 0000000..4e1820c
--- /dev/null
+++ b/pFUnittests/test_moon_reduction.pf
@@ -0,0 +1,209 @@
+!==============================================================================
+!
+!> @anchor test_moon_reduction
+!!
+!> @brief Unit tests for Moon reference frame transformations
+!!
+!> @author Christopher Kebschull (CHK)
+!!
+!> @date
+!! - CHK: 2026-05-18 (initial design)
+!!
+!!
+!> @details Tests cover:
+!! 1. Moon physical constants (slam_moon_astro)
+!! 2. Rotation matrix orthogonality at J2000.0
+!! 3. Round-trip position: GCRF -> MOON_FIXED -> GCRF
+!! 4. Round-trip position+velocity: GCRF -> MOON_FIXED -> GCRF
+!! 5. MCRF frame ID and name lookups (slam_rframes)
+!
+!------------------------------------------------------------------------------
+module test_moon_reduction
+
+ use funit
+ use slam_types, only: dp
+ use slam_moon_astro, only: getMoonRadius, getMoonGravity, getMoonRotation
+ use slam_moon_reduction, only: Reduction_moon_type
+ use slam_rframes, only: REF_FRAME_MCRF, REF_FRAME_MOON_FIXED, &
+ FRAME_CENTER_MOON, &
+ C_REF_FRAME_MCRF, C_REF_FRAME_MOON_FIXED, &
+ C_FRAME_CENTER_MOON, &
+ getFrameId, getFrameName, &
+ getFrameCenterId, getFrameCenterName
+
+ implicit none
+
+contains
+
+ !--------------------------------------------------------------------------
+ ! Constants tests
+ !--------------------------------------------------------------------------
+
+ @test
+ subroutine test_moonRadius()
+ @assertEqual(1737.4d0, getMoonRadius(), tolerance=1.d-10)
+ end subroutine test_moonRadius
+
+ @test
+ subroutine test_moonGravity()
+ ! DE-430 value, check to 3 decimal places
+ @assertEqual(4902.800066d0, getMoonGravity(), tolerance=1.d-6)
+ end subroutine test_moonGravity
+
+ @test
+ subroutine test_moonRotationRatePositive()
+ @assertTrue(getMoonRotation() > 0.d0)
+ end subroutine test_moonRotationRatePositive
+
+ !--------------------------------------------------------------------------
+ ! Frame ID / name lookup tests
+ !--------------------------------------------------------------------------
+
+ @test
+ subroutine test_frameIds()
+ @assertEqual(11, REF_FRAME_MCRF)
+ @assertEqual(12, REF_FRAME_MOON_FIXED)
+ @assertEqual(3, FRAME_CENTER_MOON)
+ end subroutine test_frameIds
+
+ @test
+ subroutine test_frameNameLookup()
+ @assertEqual(REF_FRAME_MCRF, getFrameId(C_REF_FRAME_MCRF))
+ @assertEqual(REF_FRAME_MOON_FIXED, getFrameId(C_REF_FRAME_MOON_FIXED))
+ @assertEqual(FRAME_CENTER_MOON, getFrameCenterId(C_FRAME_CENTER_MOON))
+ end subroutine test_frameNameLookup
+
+ !--------------------------------------------------------------------------
+ ! Rotation matrix orthogonality at J2000.0
+ !--------------------------------------------------------------------------
+
+ @test
+ subroutine test_rotationMatrixOrthogonal()
+
+ type(Reduction_moon_type) :: moon_red
+ real(dp), dimension(3,3) :: M, should_be_identity
+ real(dp), dimension(3,3) :: eye
+ real(dp) :: time_mjd
+ integer :: i, j
+
+ time_mjd = 51544.5d0 ! J2000.0
+
+ call moon_red%getMoonFixedRotationMatrix(time_mjd)
+ M = moon_red%R_gcrf2moonFixed
+
+ !** M * M^T must equal identity
+ should_be_identity = matmul(M, transpose(M))
+
+ eye = 0.d0
+ do i = 1, 3
+ eye(i,i) = 1.d0
+ end do
+
+ do i = 1, 3
+ do j = 1, 3
+ @assertEqual(eye(i,j), should_be_identity(i,j), tolerance=1.d-12)
+ end do
+ end do
+
+ end subroutine test_rotationMatrixOrthogonal
+
+ !--------------------------------------------------------------------------
+ ! Round-trip position (GCRF -> MOON_FIXED -> GCRF)
+ !--------------------------------------------------------------------------
+
+ @test
+ subroutine test_roundTripPosition()
+
+ type(Reduction_moon_type) :: moon_red
+ real(dp), dimension(3) :: r_gcrf, r_moon_gcrf, r_moon_fixed, r_gcrf_recovered
+ real(dp) :: time_mjd
+ integer :: i
+
+ time_mjd = 58600.d0 ! 2019 Apr 27
+
+ !** Moon position from DE-421 ephemeris (reference: test_solarsystem)
+ r_moon_gcrf = (/242037.66047749575d0, -291029.18822841981d0, -137196.91059160809d0/)
+
+ !** Arbitrary satellite position near Moon
+ r_gcrf = (/243000.d0, -290500.d0, -136800.d0/)
+
+ call moon_red%gcrf2moonFixed(r_gcrf, r_moon_gcrf, time_mjd, r_moon_fixed)
+ call moon_red%moonFixed2gcrf(r_moon_fixed, r_moon_gcrf, time_mjd, r_gcrf_recovered)
+
+ do i = 1, 3
+ @assertEqual(r_gcrf(i), r_gcrf_recovered(i), tolerance=1.d-8)
+ end do
+
+ end subroutine test_roundTripPosition
+
+ !--------------------------------------------------------------------------
+ ! Round-trip position and velocity (GCRF -> MOON_FIXED -> GCRF)
+ !--------------------------------------------------------------------------
+
+ @test
+ subroutine test_roundTripPositionVelocity()
+
+ type(Reduction_moon_type) :: moon_red
+ real(dp), dimension(3) :: r_gcrf, v_gcrf
+ real(dp), dimension(3) :: r_moon_gcrf, v_moon_gcrf
+ real(dp), dimension(3) :: r_moon_fixed, v_moon_fixed
+ real(dp), dimension(3) :: r_gcrf_rec, v_gcrf_rec
+ real(dp) :: time_mjd
+ integer :: i
+
+ time_mjd = 58600.d0
+
+ r_moon_gcrf = (/242037.66047749575d0, -291029.18822841981d0, -137196.91059160809d0/)
+ v_moon_gcrf = (/-0.9d0, 0.7d0, 0.3d0/) ! approximate (km/s)
+
+ r_gcrf = (/243000.d0, -290500.d0, -136800.d0/)
+ v_gcrf = (/-0.85d0, 0.65d0, 0.28d0 /)
+
+ call moon_red%gcrf2moonFixed(r_gcrf, v_gcrf, r_moon_gcrf, v_moon_gcrf, &
+ time_mjd, r_moon_fixed, v_moon_fixed)
+ call moon_red%moonFixed2gcrf(r_moon_fixed, v_moon_fixed, r_moon_gcrf, v_moon_gcrf, &
+ time_mjd, r_gcrf_rec, v_gcrf_rec)
+
+ do i = 1, 3
+ @assertEqual(r_gcrf(i), r_gcrf_rec(i), tolerance=1.d-8)
+ @assertEqual(v_gcrf(i), v_gcrf_rec(i), tolerance=1.d-12)
+ end do
+
+ end subroutine test_roundTripPositionVelocity
+
+ !--------------------------------------------------------------------------
+ ! Verify north pole direction in Moon body-fixed frame is +z
+ !--------------------------------------------------------------------------
+
+ @test
+ subroutine test_northPoleAlongZaxis()
+ use slam_math, only: deg2rad
+
+ type(Reduction_moon_type) :: moon_red
+ real(dp), dimension(3) :: pole_gcrf, pole_fixed
+ real(dp) :: ra0, dec0, time_mjd
+ real(dp) :: dummy_W
+
+ !** IAU 2015 values at J2000.0 (approximate)
+ ra0 = 269.9949d0 * deg2rad
+ dec0 = 66.5392d0 * deg2rad
+ time_mjd = 51544.5d0
+
+ !** North pole direction unit vector in GCRF
+ pole_gcrf(1) = cos(dec0) * cos(ra0)
+ pole_gcrf(2) = cos(dec0) * sin(ra0)
+ pole_gcrf(3) = sin(dec0)
+
+ call moon_red%getMoonFixedRotationMatrix(time_mjd)
+ pole_fixed = matmul(moon_red%R_gcrf2moonFixed, pole_gcrf)
+
+ !** In body-fixed frame the pole should point nearly along +z
+ !** (small deviation because the libration terms shift the exact pole slightly
+ !** from the secular value, so we use a loose tolerance here)
+ @assertEqual(0.d0, pole_fixed(1), tolerance=0.1d0)
+ @assertEqual(0.d0, pole_fixed(2), tolerance=0.1d0)
+ @assertEqual(1.d0, pole_fixed(3), tolerance=0.1d0)
+
+ end subroutine test_northPoleAlongZaxis
+
+end module test_moon_reduction
diff --git a/src/astro/CMakeLists.txt b/src/astro/CMakeLists.txt
index 4840953..7b698d1 100644
--- a/src/astro/CMakeLists.txt
+++ b/src/astro/CMakeLists.txt
@@ -6,10 +6,12 @@ add_library(astro-Fortran OBJECT
slam_orbit_types.f90
slam_astro.f90
slam_astro_conversions.f90
+ slam_moon_astro.f90
slam_txys.f90
slam_iaupn.f
slam_iaupn76.f
slam_reduction.f90
+ slam_moon_reduction.f90
slam_units.f90
slam_parsetle.f90
)
diff --git a/src/astro/slam_astro.f90 b/src/astro/slam_astro.f90
index e90c35b..035f846 100644
--- a/src/astro/slam_astro.f90
+++ b/src/astro/slam_astro.f90
@@ -37,6 +37,8 @@ module slam_astro
real(dp) :: massEarth ! earth's mass in kg
real(dp) :: mu ! gravity constant km^3/s^2
real(dp) :: muorekm2 ! mu/rekm2
+ real(dp) :: mu_moon = 4902.8d0 ! Moon's gravity constant km³/s²
+ real(dp) :: geo_rekm_moon = 1738.0d0 ! Moon's geopotential reference radius km
real(dp) :: omega_earth ! earth's rotational rate rad/s
real(dp) :: flat ! flattening
real(dp) :: sun_radius ! radius of the sun
@@ -79,6 +81,12 @@ module slam_astro
public :: setEarthGeopotentialRadius
public :: setEarthGravity
public :: setEarthRadius
+ public :: setLunarGravity
+ public :: setLunarGeopotentialRadius
+
+ !** lunar getter
+ public :: getLunarGravity
+ public :: getLunarGeopotentialRadius
contains
@@ -578,6 +586,78 @@ subroutine setEarthGeopotentialRadius(rekm_in)
end subroutine setEarthGeopotentialRadius
+ !=================================================================
+ !> @anchor setLunarGravity
+ !> @brief Set the Moon's gravity constant
+ !> @param[in] mu_in Moon's gravity constant in km³/s²
+ !--------------------------------------------------------
+ subroutine setLunarGravity(mu_in)
+
+ real(dp), intent(in) :: mu_in
+
+ character(len=*), parameter :: csubid = 'setLunarGravity'
+
+ if(isControlled()) then
+ if(hasToReturn()) return
+ call checkIn(csubid)
+ end if
+
+ if(mu_in < 4800.d0 .or. mu_in > 5100.d0) then
+ call setError(E_MOON_GRAVITY, WARNING)
+ return
+ end if
+
+ mu_moon = mu_in
+
+ if(isControlled()) call checkOut(csubid)
+
+ end subroutine setLunarGravity
+
+ !=================================================================
+ !> @anchor getLunarGravity
+ !> @brief Get the Moon's gravity constant
+ !> @returns Moon's gravity constant in km³/s²
+ !--------------------------------------------------------
+ real(dp) function getLunarGravity()
+ getLunarGravity = mu_moon
+ end function getLunarGravity
+
+ !=================================================================
+ !> @anchor setLunarGeopotentialRadius
+ !> @brief Set the Moon's geopotential reference radius
+ !> @param[in] rekm_in Moon's reference radius in km
+ !--------------------------------------------------------
+ subroutine setLunarGeopotentialRadius(rekm_in)
+
+ real(dp), intent(in) :: rekm_in
+
+ character(len=*), parameter :: csubid = 'setLunarGeopotentialRadius'
+
+ if(isControlled()) then
+ if(hasToReturn()) return
+ call checkIn(csubid)
+ end if
+
+ if(rekm_in < 1700.d0 .or. rekm_in > 1800.d0) then
+ call setError(E_MOON_RADIUS, WARNING)
+ return
+ end if
+
+ geo_rekm_moon = rekm_in
+
+ if(isControlled()) call checkOut(csubid)
+
+ end subroutine setLunarGeopotentialRadius
+
+ !=================================================================
+ !> @anchor getLunarGeopotentialRadius
+ !> @brief Get the Moon's geopotential reference radius
+ !> @returns Moon's reference radius in km
+ !--------------------------------------------------------
+ real(dp) function getLunarGeopotentialRadius()
+ getLunarGeopotentialRadius = geo_rekm_moon
+ end function getLunarGeopotentialRadius
+
!=========================================================================
!
!> @brief Get Speed of light for this universe
diff --git a/src/astro/slam_moon_astro.f90 b/src/astro/slam_moon_astro.f90
new file mode 100644
index 0000000..79f8719
--- /dev/null
+++ b/src/astro/slam_moon_astro.f90
@@ -0,0 +1,312 @@
+!>------------------------------------------------------------------------------------
+!!
+!> @brief Physical constants for the Moon
+!!
+!! @anchor slam_moon_astro
+!!
+!> @author Christopher Kebschull (CHK)
+!!
+!> @date
+!! - CHK: 2026-05-18 (initial implementation)
+!!
+!!
+!> @details This module provides lunar physical constants analogous to the Earth
+!! constants in slam_astro. Values are taken from the IAU 2015 report
+!! (Archinal et al. 2018) and the DE-430 ephemeris.
+!!
+!> @copyright OKAPI:Orbits
+!!
+!!------------------------------------------------------------------------------------
+module slam_moon_astro
+
+ use slam_types, only: dp
+ use slam_orbit_types, only: kepler_t
+ use slam_math, only: angle, cross, deg2rad, eps6, eps9, halfpi, infinite, mag, pi, redang, twopi, undefined
+ use slam_astro_conversions, only: ELLIPTICAL_INCLINED, CIRCULAR_INCLINED, ELLIPTICAL_EQUATORIAL, CIRCULAR_EQUATORIAL
+ use slam_units, only: UNIT_KM, UNIT_RAD
+
+ implicit none
+
+ private
+
+ !** Moon physical constants
+ !-----------------------------------------------------------
+ real(dp), parameter :: moon_radius_km = 1737.4d0 ! mean radius (km), IAU 2015
+ real(dp), parameter :: moon_mu_km3s2 = 4902.800066d0 ! GM (km^3/s^2), DE-430
+ real(dp), parameter :: moon_rot_rate = 2.6617d-6 ! mean rotation rate (rad/s)
+
+ public :: getMoonRadius
+ public :: getMoonGravity
+ public :: getMoonRotation
+ public :: moon_coe2rv
+ public :: moon_rv2coe
+
+contains
+
+ !=========================================================================
+ !> @anchor getMoonRadius
+ !> @brief Get Moon's mean radius in km
+ !> @returns real(dp) Moon's mean radius in km
+ !-------------------------------------------------------------
+ real(dp) function getMoonRadius()
+ getMoonRadius = moon_radius_km
+ end function getMoonRadius
+
+ !=========================================================================
+ !> @anchor getMoonGravity
+ !> @brief Get Moon's gravitational parameter GM in km^3/s^2
+ !> @returns real(dp) Moon's GM in km^3/s^2
+ !-------------------------------------------------------------
+ real(dp) function getMoonGravity()
+ getMoonGravity = moon_mu_km3s2
+ end function getMoonGravity
+
+ !=========================================================================
+ !> @anchor getMoonRotation
+ !> @brief Get Moon's mean rotation rate in rad/s
+ !> @returns real(dp) Moon's mean rotation rate in rad/s
+ !-------------------------------------------------------------
+ real(dp) function getMoonRotation()
+ getMoonRotation = moon_rot_rate
+ end function getMoonRotation
+
+ !=========================================================================
+ !> @anchor moon_coe2rv
+ !> @brief Convert selenocentric osculating Keplerian elements to Cartesian
+ !!
+ !> @param[in] kep Keplerian elements (angles in radians, sma in km)
+ !> @param[out] r_mcrf Position in MCRF (km)
+ !> @param[out] v_mcrf Velocity in MCRF (km/s)
+ !!
+ !> @details Uses Moon's GM. The MCRF (Moon-Centred Reference Frame) has
+ !! axes aligned with GCRF, so the output can be translated to
+ !! GCRF by adding the Moon's GCRF position and velocity.
+ !! Conversion via the perifocal (PQW) frame:
+ !! r_pqw = [r*cos(nu), r*sin(nu), 0]
+ !! v_pqw = sqrt(GM/p) * [-sin(nu), e+cos(nu), 0]
+ !! followed by rotation Rz(-Om)*Rx(-i)*Rz(-om) to MCRF.
+ !---------------------------------------------------------------------------
+ subroutine moon_coe2rv(kep, r_mcrf, v_mcrf)
+
+ type(kepler_t), intent(in) :: kep
+ real(dp), dimension(3), intent(out) :: r_mcrf
+ real(dp), dimension(3), intent(out) :: v_mcrf
+
+ real(dp) :: p, r_mag, vfac
+ real(dp) :: cosRaan, sinRaan, cosInc, sinInc, cosAop, sinAop, cosNu, sinNu
+ real(dp) :: r_pqw(3), v_pqw(3)
+ real(dp) :: Q(3,3)
+
+ cosNu = cos(kep%tran); sinNu = sin(kep%tran)
+ p = kep%sma * (1.d0 - kep%ecc**2)
+ r_mag = p / (1.d0 + kep%ecc * cosNu)
+
+ !** position and velocity in perifocal frame
+ r_pqw = r_mag * (/cosNu, sinNu, 0.d0/)
+ vfac = sqrt(moon_mu_km3s2 / p)
+ v_pqw = vfac * (/-sinNu, kep%ecc + cosNu, 0.d0/)
+
+ !** rotation matrix: perifocal -> MCRF (inertial)
+ !** Q = Rz(-Omega) * Rx(-i) * Rz(-omega)
+ cosRaan = cos(kep%raan); sinRaan = sin(kep%raan)
+ cosInc = cos(kep%inc); sinInc = sin(kep%inc)
+ cosAop = cos(kep%aop); sinAop = sin(kep%aop)
+
+ Q(1,1) = cosRaan*cosAop - sinRaan*sinAop*cosInc
+ Q(1,2) = -cosRaan*sinAop - sinRaan*cosAop*cosInc
+ Q(1,3) = sinRaan*sinInc
+ Q(2,1) = sinRaan*cosAop + cosRaan*sinAop*cosInc
+ Q(2,2) = -sinRaan*sinAop + cosRaan*cosAop*cosInc
+ Q(2,3) = -cosRaan*sinInc
+ Q(3,1) = sinInc*sinAop
+ Q(3,2) = sinInc*cosAop
+ Q(3,3) = cosInc
+
+ r_mcrf = matmul(Q, r_pqw)
+ v_mcrf = matmul(Q, v_pqw)
+
+ return
+
+ end subroutine moon_coe2rv
+
+ !=========================================================================
+ !> @anchor moon_rv2coe
+ !> @brief Convert Moon-centered Cartesian state (MCRF) to osculating Keplerian elements
+ !!
+ !> @param[in] r_mcrf Position in MCRF (km)
+ !> @param[in] v_mcrf Velocity in MCRF (km/s)
+ !> @param[out] kep Keplerian elements (angles in radians, sma in km)
+ !> @param[out] otype Orbit type (ELLIPTICAL_INCLINED, CIRCULAR_INCLINED, etc.)
+ !!
+ !> @details Uses Moon's GM. The algorithm mirrors rv2coe from slam_astro_conversions
+ !! but uses moon_mu_km3s2 instead of Earth's GM.
+ !---------------------------------------------------------------------------
+ subroutine moon_rv2coe(r_mcrf, v_mcrf, kep, otype)
+
+ real(dp), dimension(3), intent(in) :: r_mcrf
+ real(dp), dimension(3), intent(in) :: v_mcrf
+ type(kepler_t), intent(out) :: kep
+ integer, intent(out) :: otype
+
+ real(dp) :: c1, rdotv, hk, sme, semipar, temp
+ real(dp) :: hbar(3), ebar(3), nbar(3), ecc_anom
+ real(dp) :: maghbar, magnbar, magr, magv
+ real(dp) :: sine, cose
+ character(len=2) :: typeorbit
+
+ magr = mag(r_mcrf)
+ magv = mag(v_mcrf)
+
+ hbar = cross(r_mcrf, v_mcrf)
+ maghbar = mag(hbar)
+
+ if (maghbar > eps9) then
+
+ nbar(1) = -hbar(2)
+ nbar(2) = hbar(1)
+ nbar(3) = 0.d0
+ magnbar = mag(nbar)
+
+ c1 = magv**2 - moon_mu_km3s2 / magr
+ rdotv = dot_product(r_mcrf, v_mcrf)
+ ebar = (c1 * r_mcrf - rdotv * v_mcrf) / moon_mu_km3s2
+
+ kep%ecc = mag(ebar)
+
+ sme = 0.5d0 * magv**2 - moon_mu_km3s2 / magr
+ if (abs(sme) > eps9) then
+ kep%sma = -moon_mu_km3s2 / (2.d0 * sme)
+ else
+ kep%sma = infinite
+ end if
+
+ semipar = maghbar**2 / moon_mu_km3s2
+
+ hk = hbar(3) / maghbar
+ if (abs(abs(hk) - 1.d0) < eps9) hk = sign(1.d0, hbar(3))
+ kep%inc = acos(hk)
+
+ !** orbit type
+ typeorbit = 'EI'
+ otype = ELLIPTICAL_INCLINED
+ if (kep%ecc < eps9) then
+ if (kep%inc < eps9 .or. abs(kep%inc - pi) < eps9) then
+ typeorbit = 'CE'; otype = CIRCULAR_EQUATORIAL
+ else
+ typeorbit = 'CI'; otype = CIRCULAR_INCLINED
+ end if
+ else if (kep%inc < eps9 .or. abs(kep%inc - pi) < eps9) then
+ typeorbit = 'EE'; otype = ELLIPTICAL_EQUATORIAL
+ end if
+
+ !** RAAN
+ if (magnbar > eps9) then
+ temp = nbar(1) / magnbar
+ if (abs(temp) > 1.d0) temp = sign(1.d0, temp)
+ kep%raan = acos(temp)
+ if (nbar(2) < 0.d0) kep%raan = twopi - kep%raan
+ kep%raan = redang(kep%raan, 2, 1, .false.)
+ else
+ kep%raan = undefined
+ end if
+
+ !** argument of perigee
+ if (typeorbit == 'EI') then
+ call angle(nbar, ebar, kep%aop)
+ if (ebar(3) < 0.d0) kep%aop = twopi - kep%aop
+ kep%aop = redang(kep%aop, 2, 1, .false.)
+ else
+ kep%aop = undefined
+ end if
+
+ !** true anomaly
+ if (typeorbit(1:1) == 'E') then
+ call angle(ebar, r_mcrf, kep%tran)
+ if (rdotv < 0.d0) kep%tran = twopi - kep%tran
+ kep%tran = redang(kep%tran, 2, 1, .false.)
+ else
+ kep%tran = undefined
+ end if
+
+ !** argument of latitude (circular inclined)
+ if (typeorbit == 'CI') then
+ call angle(nbar, r_mcrf, kep%arglat)
+ if (r_mcrf(3) < 0.d0) kep%arglat = twopi - kep%arglat
+ kep%arglat = redang(kep%arglat, 2, 1, .false.)
+ else
+ kep%arglat = undefined
+ end if
+
+ !** longitude of periapsis (elliptical equatorial)
+ if (kep%ecc > eps9 .and. typeorbit == 'EE') then
+ temp = ebar(1) / kep%ecc
+ if (abs(temp) > 1.d0) temp = sign(1.d0, temp)
+ kep%lonper = acos(temp)
+ if (ebar(2) < 0.d0) kep%lonper = twopi - kep%lonper
+ if (kep%inc > halfpi) kep%lonper = twopi - kep%lonper
+ kep%lonper = redang(kep%lonper, 2, 1, .false.)
+ else
+ kep%lonper = undefined
+ end if
+
+ !** true longitude (circular equatorial)
+ if (magr > eps9 .and. typeorbit == 'CE') then
+ temp = r_mcrf(1) / magr
+ if (abs(temp) > 1.d0) temp = sign(1.d0, temp)
+ kep%truelon = acos(temp)
+ if (r_mcrf(2) < 0.d0) kep%truelon = twopi - kep%truelon
+ if (kep%inc > halfpi) kep%truelon = twopi - kep%truelon
+ kep%truelon = redang(kep%truelon, 2, 1, .false.)
+ else
+ kep%truelon = undefined
+ end if
+
+ !** mean anomaly (inline true2mean)
+ if (typeorbit(1:1) == 'E') then
+ if (abs(kep%ecc) < eps6) then
+ kep%man = kep%tran
+ ecc_anom = kep%tran
+ else if (kep%ecc < 0.999d0) then
+ sine = sqrt(1.d0 - kep%ecc**2) * sin(kep%tran) / (1.d0 + kep%ecc * cos(kep%tran))
+ cose = (kep%ecc + cos(kep%tran)) / (1.d0 + kep%ecc * cos(kep%tran))
+ ecc_anom = atan2(sine, cose)
+ kep%man = ecc_anom - kep%ecc * sin(ecc_anom)
+ else if (kep%ecc > 1.0001d0) then
+ sine = sqrt(kep%ecc**2 - 1.d0) * sin(kep%tran) / (1.d0 + kep%ecc * cos(kep%tran))
+ ecc_anom = asinh(sine)
+ kep%man = kep%ecc * sinh(ecc_anom) - ecc_anom
+ else
+ if (abs(kep%tran) < 168.d0 * deg2rad) then
+ ecc_anom = tan(kep%tran * 0.5d0)
+ kep%man = ecc_anom + ecc_anom**3 / 3.d0
+ end if
+ end if
+ if (kep%ecc < 1.d0) then
+ kep%man = mod(kep%man, twopi)
+ if (kep%man < 0.d0) kep%man = kep%man + twopi
+ end if
+ else if (typeorbit == 'CI') then
+ kep%man = kep%arglat
+ else if (typeorbit == 'CE') then
+ kep%man = kep%truelon
+ end if
+
+ else
+
+ kep%sma = undefined; kep%ecc = undefined
+ kep%inc = undefined; kep%raan = undefined
+ kep%aop = undefined; kep%tran = undefined
+ kep%man = undefined; kep%arglat = undefined
+ kep%truelon = undefined; kep%lonper = undefined
+ otype = ELLIPTICAL_INCLINED
+
+ end if
+
+ kep%sma_unit = UNIT_KM
+ kep%angles_unit = UNIT_RAD
+
+ return
+
+ end subroutine moon_rv2coe
+
+end module slam_moon_astro
diff --git a/src/astro/slam_moon_reduction.f90 b/src/astro/slam_moon_reduction.f90
new file mode 100644
index 0000000..f074b8d
--- /dev/null
+++ b/src/astro/slam_moon_reduction.f90
@@ -0,0 +1,362 @@
+!>------------------------------------------------------------------------------------
+!!
+!> @brief Coordinate frame transformations for Moon-centered reference frames
+!!
+!! @anchor slam_moon_reduction
+!!
+!> @author Christopher Kebschull (CHK)
+!!
+!> @date
+!! - CHK: 2026-05-18 (initial implementation)
+!!
+!!
+!> @details Implements transformations between GCRF and Moon-centered frames:
+!!
+!! MCRF (Moon-Centered Reference Frame): inertial frame with origin at the
+!! Moon's centre of mass, axes aligned with GCRF. Obtained by translating
+!! the GCRF origin to the Moon.
+!!
+!! MOON_FIXED: Moon body-fixed rotating frame. The orientation follows the
+!! IAU 2015 report (Archinal et al. 2018, CeMDA 130:22). The north pole
+!! direction and prime meridian angle are computed from analytical series
+!! with 13 argument angles derived from the Moon's mean orbital elements.
+!!
+!! Rotation from MCRF to MOON_FIXED:
+!! M = R_z(W) * R_x(90 - delta0) * R_z(90 + alpha0)
+!!
+!! where alpha0 and delta0 are the right ascension and declination of the
+!! lunar north pole in GCRF, and W is the prime meridian angle.
+!!
+!! Velocity transformation accounts for the instantaneous rotation using
+!! the mean rate dW/dt = 13.17635815 deg/day.
+!!
+!> @copyright OKAPI:Orbits
+!!
+!!------------------------------------------------------------------------------------
+module slam_moon_reduction
+
+ use slam_types, only: dp
+ use slam_math, only: pi, halfPi, deg2rad, cross
+ use slam_error_handling, only: isControlled, hasToReturn, checkIn, checkOut, &
+ setError, FATAL, E_SPECIAL
+
+ implicit none
+
+ private
+
+ !** Mean lunar rotation rate: 13.17635815 deg/day -> rad/s
+ real(dp), parameter :: moon_rot_rate_rads = 13.17635815d0 * deg2rad / 86400.d0
+
+ !** MJD of J2000.0
+ real(dp), parameter :: mjd_j2000 = 51544.5d0
+
+ public :: Reduction_moon_type
+
+ !=================================================================
+ !
+ !> @brief Type for Moon frame transformations
+ !!
+ !! Usage:
+ !! type(Reduction_moon_type) :: moon_red
+ !! call moon_red%gcrf2moonFixed(r_gcrf, r_moon_gcrf, time_mjd, r_moon_fixed)
+ !!
+ !-----------------------------------------------------------------
+ type :: Reduction_moon_type
+
+ real(dp), dimension(3,3) :: R_gcrf2moonFixed = 0.d0 ! rotation MCRF -> MOON_FIXED
+ real(dp) :: rotMatrixDate = -1.d30 ! MJD epoch of stored matrix
+
+ contains
+
+ procedure :: getMoonFixedRotationMatrix
+
+ procedure :: gcrf2moonFixed_r
+ procedure :: gcrf2moonFixed_rv
+ generic :: gcrf2moonFixed => gcrf2moonFixed_r, gcrf2moonFixed_rv
+
+ procedure :: moonFixed2gcrf_r
+ procedure :: moonFixed2gcrf_rv
+ generic :: moonFixed2gcrf => moonFixed2gcrf_r, moonFixed2gcrf_rv
+
+ end type Reduction_moon_type
+
+contains
+
+ !=========================================================================
+ !> @anchor getMoonIAUAngles
+ !> @brief Compute IAU 2015 Moon orientation angles at a given epoch
+ !!
+ !> @param[in] time_mjd Epoch in MJD (TDB)
+ !> @param[out] ra0 Right ascension of north pole (rad)
+ !> @param[out] dec0 Declination of north pole (rad)
+ !> @param[out] W Prime meridian angle (rad)
+ !!
+ !> @details Source: Archinal et al. 2018, CeMDA 130:22, Table 2 (Moon).
+ !---------------------------------------------------------------------------
+ subroutine getMoonIAUAngles(time_mjd, ra0, dec0, W)
+
+ real(dp), intent(in) :: time_mjd
+ real(dp), intent(out) :: ra0, dec0, W
+
+ real(dp) :: T, d
+ real(dp) :: E1, E2, E3, E4, E5, E6, E7, E8, E9, E10, E11, E12, E13
+
+ d = time_mjd - mjd_j2000
+ T = d / 36525.d0
+
+ !** IAU 2015 argument angles (degrees -> radians)
+ E1 = (125.045d0 - 0.0529921d0 * d) * deg2rad
+ E2 = (250.089d0 - 0.1059842d0 * d) * deg2rad
+ E3 = (260.008d0 + 13.012009d0 * d) * deg2rad
+ E4 = (176.625d0 + 13.340716d0 * d) * deg2rad
+ E5 = (357.529d0 + 0.985600d0 * d) * deg2rad
+ E6 = (311.589d0 + 26.4057084d0 * d) * deg2rad
+ E7 = (134.963d0 + 13.064993d0 * d) * deg2rad
+ E8 = (276.617d0 + 0.3287146d0 * d) * deg2rad
+ E9 = ( 34.226d0 + 1.7484877d0 * d) * deg2rad
+ E10 = ( 15.134d0 - 0.1589763d0 * d) * deg2rad
+ E11 = (119.743d0 + 0.0036096d0 * d) * deg2rad
+ E12 = (239.961d0 + 0.1643573d0 * d) * deg2rad
+ E13 = ( 25.053d0 + 12.9590088d0 * d) * deg2rad
+
+ !** North pole right ascension (degrees)
+ ra0 = 269.9949d0 + 0.0013d0*T &
+ - 3.8787d0*sin(E1) - 0.1204d0*sin(E2) + 0.0700d0*sin(E3) &
+ - 0.0172d0*sin(E4) + 0.0072d0*sin(E6) - 0.0052d0*sin(E10) &
+ + 0.0043d0*sin(E13)
+
+ !** North pole declination (degrees)
+ dec0 = 66.5392d0 + 0.0130d0*T &
+ + 1.5419d0*cos(E1) + 0.0239d0*cos(E2) - 0.0278d0*cos(E3) &
+ + 0.0068d0*cos(E4) - 0.0029d0*cos(E6) + 0.0009d0*cos(E7) &
+ + 0.0008d0*cos(E10) - 0.0009d0*cos(E13)
+
+ !** Prime meridian angle (degrees)
+ W = 38.3213d0 + 13.17635815d0*d - 1.4d-12*d*d &
+ + 3.5610d0*sin(E1) + 0.1208d0*sin(E2) - 0.0642d0*sin(E3) &
+ + 0.0158d0*sin(E4) + 0.0252d0*sin(E5) - 0.0066d0*sin(E6) &
+ - 0.0047d0*sin(E7) - 0.0046d0*sin(E8) + 0.0028d0*sin(E9) &
+ + 0.0052d0*sin(E10) + 0.0040d0*sin(E11) + 0.0019d0*sin(E12) &
+ - 0.0044d0*sin(E13)
+
+ ra0 = ra0 * deg2rad
+ dec0 = dec0 * deg2rad
+ W = W * deg2rad
+
+ return
+
+ end subroutine getMoonIAUAngles
+
+ !=========================================================================
+ !> @anchor getMoonFixedRotationMatrix
+ !> @brief Build and cache the MCRF -> MOON_FIXED rotation matrix
+ !!
+ !> @param[in] time_mjd Epoch in MJD (TDB)
+ !!
+ !> @details Uses the IAU 2015 series. The matrix is only recomputed when
+ !! the requested epoch differs from the cached epoch.
+ !! Convention: r_moon_fixed = R * r_mcrf
+ !! R = Rz(W) * Rx(90 - dec0) * Rz(90 + ra0)
+ !---------------------------------------------------------------------------
+ subroutine getMoonFixedRotationMatrix(this, time_mjd)
+
+ class(Reduction_moon_type), intent(inout) :: this
+ real(dp), intent(in) :: time_mjd
+
+ real(dp) :: ra0, dec0, W
+ real(dp) :: ang_node, ang_pole
+ real(dp) :: cn, sn, cp, sp, cW, sW
+ real(dp), dimension(3,3) :: Rz_node, Rx_pole, Rz_W
+
+ if (abs(time_mjd - this%rotMatrixDate) < 1.d-10) return
+
+ call getMoonIAUAngles(time_mjd, ra0, dec0, W)
+
+ !** Rz(90 + ra0)
+ ang_node = halfPi + ra0
+ cn = cos(ang_node); sn = sin(ang_node)
+ Rz_node(1,:) = (/ cn, sn, 0.d0/)
+ Rz_node(2,:) = (/-sn, cn, 0.d0/)
+ Rz_node(3,:) = (/0.d0, 0.d0, 1.d0/)
+
+ !** Rx(90 - dec0)
+ ang_pole = halfPi - dec0
+ cp = cos(ang_pole); sp = sin(ang_pole)
+ Rx_pole(1,:) = (/1.d0, 0.d0, 0.d0/)
+ Rx_pole(2,:) = (/0.d0, cp, sp /)
+ Rx_pole(3,:) = (/0.d0, -sp, cp /)
+
+ !** Rz(W)
+ cW = cos(W); sW = sin(W)
+ Rz_W(1,:) = (/ cW, sW, 0.d0/)
+ Rz_W(2,:) = (/-sW, cW, 0.d0/)
+ Rz_W(3,:) = (/0.d0, 0.d0, 1.d0/)
+
+ !** R = Rz(W) * Rx(90-dec0) * Rz(90+ra0)
+ this%R_gcrf2moonFixed = matmul(Rz_W, matmul(Rx_pole, Rz_node))
+ this%rotMatrixDate = time_mjd
+
+ return
+
+ end subroutine getMoonFixedRotationMatrix
+
+ !=========================================================================
+ !> @anchor gcrf2moonFixed_r
+ !> @brief Transform a position from GCRF to Moon body-fixed
+ !!
+ !> @param[in] r_gcrf Position in GCRF (km)
+ !> @param[in] r_moon_gcrf Moon centre position in GCRF (km)
+ !> @param[in] time_mjd Epoch in MJD (TDB)
+ !> @param[out] r_moon_fixed Position in MOON_FIXED (km)
+ !---------------------------------------------------------------------------
+ subroutine gcrf2moonFixed_r(this, r_gcrf, r_moon_gcrf, time_mjd, r_moon_fixed)
+
+ class(Reduction_moon_type), intent(inout) :: this
+ real(dp), dimension(3), intent(in) :: r_gcrf
+ real(dp), dimension(3), intent(in) :: r_moon_gcrf
+ real(dp), intent(in) :: time_mjd
+ real(dp), dimension(3), intent(out) :: r_moon_fixed
+
+ character(len=*), parameter :: csubid = 'gcrf2moonFixed_r'
+
+ if(isControlled()) then
+ if(hasToReturn()) return
+ call checkIn(csubid)
+ end if
+
+ call this%getMoonFixedRotationMatrix(time_mjd)
+ r_moon_fixed = matmul(this%R_gcrf2moonFixed, r_gcrf - r_moon_gcrf)
+
+ if(isControlled()) call checkOut(csubid)
+ return
+
+ end subroutine gcrf2moonFixed_r
+
+ !=========================================================================
+ !> @anchor gcrf2moonFixed_rv
+ !> @brief Transform position and velocity from GCRF to Moon body-fixed
+ !!
+ !> @param[in] r_gcrf Position in GCRF (km)
+ !> @param[in] v_gcrf Velocity in GCRF (km/s)
+ !> @param[in] r_moon_gcrf Moon centre position in GCRF (km)
+ !> @param[in] v_moon_gcrf Moon centre velocity in GCRF (km/s)
+ !> @param[in] time_mjd Epoch in MJD (TDB)
+ !> @param[out] r_moon_fixed Position in MOON_FIXED (km)
+ !> @param[out] v_moon_fixed Velocity in MOON_FIXED (km/s)
+ !!
+ !> @details The velocity accounts for the frame rotation:
+ !! v_fixed = R * v_mcrf - omega x r_fixed
+ !! where omega = (0, 0, dW/dt) in the body-fixed frame.
+ !---------------------------------------------------------------------------
+ subroutine gcrf2moonFixed_rv(this, r_gcrf, v_gcrf, r_moon_gcrf, v_moon_gcrf, &
+ time_mjd, r_moon_fixed, v_moon_fixed)
+
+ class(Reduction_moon_type), intent(inout) :: this
+ real(dp), dimension(3), intent(in) :: r_gcrf, v_gcrf
+ real(dp), dimension(3), intent(in) :: r_moon_gcrf, v_moon_gcrf
+ real(dp), intent(in) :: time_mjd
+ real(dp), dimension(3), intent(out) :: r_moon_fixed, v_moon_fixed
+
+ real(dp), dimension(3) :: omega
+
+ character(len=*), parameter :: csubid = 'gcrf2moonFixed_rv'
+
+ if(isControlled()) then
+ if(hasToReturn()) return
+ call checkIn(csubid)
+ end if
+
+ call this%getMoonFixedRotationMatrix(time_mjd)
+
+ r_moon_fixed = matmul(this%R_gcrf2moonFixed, r_gcrf - r_moon_gcrf)
+
+ !** omega = (0, 0, dW/dt) in Moon body-fixed frame
+ omega = (/0.d0, 0.d0, moon_rot_rate_rads/)
+ v_moon_fixed = matmul(this%R_gcrf2moonFixed, v_gcrf - v_moon_gcrf) &
+ - cross(omega, r_moon_fixed)
+
+ if(isControlled()) call checkOut(csubid)
+ return
+
+ end subroutine gcrf2moonFixed_rv
+
+ !=========================================================================
+ !> @anchor moonFixed2gcrf_r
+ !> @brief Transform a position from Moon body-fixed to GCRF
+ !!
+ !> @param[in] r_moon_fixed Position in MOON_FIXED (km)
+ !> @param[in] r_moon_gcrf Moon centre position in GCRF (km)
+ !> @param[in] time_mjd Epoch in MJD (TDB)
+ !> @param[out] r_gcrf Position in GCRF (km)
+ !---------------------------------------------------------------------------
+ subroutine moonFixed2gcrf_r(this, r_moon_fixed, r_moon_gcrf, time_mjd, r_gcrf)
+
+ class(Reduction_moon_type), intent(inout) :: this
+ real(dp), dimension(3), intent(in) :: r_moon_fixed
+ real(dp), dimension(3), intent(in) :: r_moon_gcrf
+ real(dp), intent(in) :: time_mjd
+ real(dp), dimension(3), intent(out) :: r_gcrf
+
+ character(len=*), parameter :: csubid = 'moonFixed2gcrf_r'
+
+ if(isControlled()) then
+ if(hasToReturn()) return
+ call checkIn(csubid)
+ end if
+
+ call this%getMoonFixedRotationMatrix(time_mjd)
+ r_gcrf = matmul(transpose(this%R_gcrf2moonFixed), r_moon_fixed) + r_moon_gcrf
+
+ if(isControlled()) call checkOut(csubid)
+ return
+
+ end subroutine moonFixed2gcrf_r
+
+ !=========================================================================
+ !> @anchor moonFixed2gcrf_rv
+ !> @brief Transform position and velocity from Moon body-fixed to GCRF
+ !!
+ !> @param[in] r_moon_fixed Position in MOON_FIXED (km)
+ !> @param[in] v_moon_fixed Velocity in MOON_FIXED (km/s)
+ !> @param[in] r_moon_gcrf Moon centre position in GCRF (km)
+ !> @param[in] v_moon_gcrf Moon centre velocity in GCRF (km/s)
+ !> @param[in] time_mjd Epoch in MJD (TDB)
+ !> @param[out] r_gcrf Position in GCRF (km)
+ !> @param[out] v_gcrf Velocity in GCRF (km/s)
+ !!
+ !> @details Inverse of gcrf2moonFixed_rv:
+ !! v_mcrf = R^T * (v_fixed + omega x r_fixed)
+ !---------------------------------------------------------------------------
+ subroutine moonFixed2gcrf_rv(this, r_moon_fixed, v_moon_fixed, r_moon_gcrf, v_moon_gcrf, &
+ time_mjd, r_gcrf, v_gcrf)
+
+ class(Reduction_moon_type), intent(inout) :: this
+ real(dp), dimension(3), intent(in) :: r_moon_fixed, v_moon_fixed
+ real(dp), dimension(3), intent(in) :: r_moon_gcrf, v_moon_gcrf
+ real(dp), intent(in) :: time_mjd
+ real(dp), dimension(3), intent(out) :: r_gcrf, v_gcrf
+
+ real(dp), dimension(3) :: omega
+
+ character(len=*), parameter :: csubid = 'moonFixed2gcrf_rv'
+
+ if(isControlled()) then
+ if(hasToReturn()) return
+ call checkIn(csubid)
+ end if
+
+ call this%getMoonFixedRotationMatrix(time_mjd)
+
+ r_gcrf = matmul(transpose(this%R_gcrf2moonFixed), r_moon_fixed) + r_moon_gcrf
+
+ omega = (/0.d0, 0.d0, moon_rot_rate_rads/)
+ v_gcrf = matmul(transpose(this%R_gcrf2moonFixed), &
+ v_moon_fixed + cross(omega, r_moon_fixed)) &
+ + v_moon_gcrf
+
+ if(isControlled()) call checkOut(csubid)
+ return
+
+ end subroutine moonFixed2gcrf_rv
+
+end module slam_moon_reduction
diff --git a/src/astro/slam_rframes.f90 b/src/astro/slam_rframes.f90
index deb4abe..bd87dd0 100644
--- a/src/astro/slam_rframes.f90
+++ b/src/astro/slam_rframes.f90
@@ -43,6 +43,8 @@ module slam_rframes
integer, parameter :: REF_FRAME_OCRF = 8 ! Orbit Centered Reference Frame
integer, parameter :: REF_FRAME_RSW = 9 ! Radial, Normal and Binormal satellite based system
integer, parameter :: REF_FRAME_TEME = 10 ! True Equator Mean Equinox
+ integer, parameter :: REF_FRAME_MCRF = 11 ! Moon-Centered Reference Frame (inertial, aligned with GCRF)
+ integer, parameter :: REF_FRAME_MOON_FIXED = 12 ! Moon body-fixed rotating frame (IAU 2015)
!================================================
!
@@ -51,6 +53,7 @@ module slam_rframes
!-----------------------------------------
integer, parameter :: FRAME_CENTER_EARTH = 1 ! Earth
integer, parameter :: FRAME_CENTER_EARTH_BARYCENTER = 2 ! Earth Barycenter
+ integer, parameter :: FRAME_CENTER_MOON = 3 ! Moon
!================================================
!
@@ -59,6 +62,7 @@ module slam_rframes
!-----------------------------------------
character(len=*), parameter :: C_FRAME_CENTER_EARTH = "EARTH"
character(len=*), parameter :: C_FRAME_CENTER_EARTH_BARYCENTER = "EARTH BARYCENTER"
+ character(len=*), parameter :: C_FRAME_CENTER_MOON = "MOON"
!================================================
@@ -78,8 +82,10 @@ module slam_rframes
character(len=*), parameter :: C_REF_FRAME_WGS84 = "WGS84" ! WGS84 reference frame
character(len=*), parameter :: C_REF_FRAME_J2000 = "J2000" ! J2000.0 reference frame
character(len=*), parameter :: C_REF_FRAME_OCRF = "OCRF" ! Orbit Centered Reference Frame
- character(len=*), parameter :: C_REF_FRAME_RSW = "RSW" ! Radial, Normal and Binormal satellite based system
- character(len=*), parameter :: C_REF_FRAME_TEME = "TEME" ! True Equator Mean Equinox
+ character(len=*), parameter :: C_REF_FRAME_RSW = "RSW" ! Radial, Normal and Binormal satellite based system
+ character(len=*), parameter :: C_REF_FRAME_TEME = "TEME" ! True Equator Mean Equinox
+ character(len=*), parameter :: C_REF_FRAME_MCRF = "MCRF" ! Moon-Centered Reference Frame (inertial)
+ character(len=*), parameter :: C_REF_FRAME_MOON_FIXED = "MOON_FIXED" ! Moon body-fixed rotating frame
contains
@@ -174,6 +180,9 @@ integer function getFrameCenterId(cname) result(iout)
case (C_FRAME_CENTER_EARTH_BARYCENTER) !** earth barycenter
iout = FRAME_CENTER_EARTH_BARYCENTER
+ case (C_FRAME_CENTER_MOON) !** moon
+ iout = FRAME_CENTER_MOON
+
case default !** unknown
call setError(E_FRAME_CENTER, FATAL, (/cname/))
@@ -228,6 +237,9 @@ character(len=len(C_FRAME_CENTER_EARTH_BARYCENTER)) function getFrameCenterName(
case(FRAME_CENTER_EARTH_BARYCENTER)
ccenter = C_FRAME_CENTER_EARTH_BARYCENTER
+ case(FRAME_CENTER_MOON)
+ ccenter = C_FRAME_CENTER_MOON
+
case default
write(ctemp,'(i3)') icenter
@@ -307,6 +319,12 @@ end function getFrameCenterName
case(REF_FRAME_TEME)
getFrameName = C_REF_FRAME_TEME
+ case(REF_FRAME_MCRF)
+ getFrameName = C_REF_FRAME_MCRF
+
+ case(REF_FRAME_MOON_FIXED)
+ getFrameName = C_REF_FRAME_MOON_FIXED
+
case default
call setError(E_UNKNOWN_PARAMETER, FATAL, (/cid/))
return
@@ -385,6 +403,12 @@ integer function getFrameId(cname)
case(C_REF_FRAME_TEME)
getFrameId = REF_FRAME_TEME
+ case(C_REF_FRAME_MCRF)
+ getFrameId = REF_FRAME_MCRF
+
+ case(C_REF_FRAME_MOON_FIXED)
+ getFrameId = REF_FRAME_MOON_FIXED
+
case default
call setError(E_UNKNOWN_PARAMETER, FATAL, (/cname/))
return
diff --git a/src/inout/slam_error_handling.f90 b/src/inout/slam_error_handling.f90
index e756cbe..f76f166 100644
--- a/src/inout/slam_error_handling.f90
+++ b/src/inout/slam_error_handling.f90
@@ -144,6 +144,8 @@ module slam_error_handling
integer, parameter, public :: E_EARTH_GRAVITY = 205 !< earth gravity constant value not accepted
integer, parameter, public :: E_EOP_INIT = 206 !< Earth orientation parameters data not
integer, parameter, public :: E_EOP_INDEX = 207 !< Earth orientation parameters data index too small
+ integer, parameter, public :: E_MOON_RADIUS = 208 !< lunar reference radius value not accepted
+ integer, parameter, public :: E_MOON_GRAVITY = 209 !< lunar gravity constant value not accepted
!** time/coordinate conversions
integer, parameter, public :: E_UTC = 300 !< UTC is not defined for dates earlier than Jan 1, 1961
integer, parameter, public :: E_LEAP_SECOND = 301 !< Leap seconds for propagations too far into future can not be considered
@@ -1061,6 +1063,18 @@ subroutine getErrorMessage(code, message, par)
write(message(1:len(message)),'(a)') "Given parameter for Earth's gravity constant seems unrealistic."
end select
+ case(E_MOON_RADIUS)
+ select case(errorLanguage)
+ case default
+ write(message(1:len(message)),'(a)') "Given parameter for the Moon's reference radius seems unrealistic."
+ end select
+
+ case(E_MOON_GRAVITY)
+ select case(errorLanguage)
+ case default
+ write(message(1:len(message)),'(a)') "Given parameter for the Moon's gravity constant seems unrealistic."
+ end select
+
case(E_EOP_INIT)
select case(errorLanguage)
case default
diff --git a/src/math/slam_math.f90 b/src/math/slam_math.f90
index ffe1eb5..66d7931 100644
--- a/src/math/slam_math.f90
+++ b/src/math/slam_math.f90
@@ -125,7 +125,7 @@ end subroutine angle
!> @anchor cross
!!
!----------------------------------------------------
- function cross( &
+ pure function cross( &
vec1, & ! <-- DBL() first vector
vec2 & ! <-- DBL() second vector
)