From 9869fd8bfc97794ae72ee40e2a5e53fc64da0b73 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sat, 13 Nov 2021 19:34:13 -0500 Subject: [PATCH] C BaryState: Calculate Moon's barycentric position and velocity vectors. The C function Astronomy_BaryState() now supports BODY_MOON. Because of the complexity of the CalcMoon() function, I ended up calculating two positions close together in time, and using dr/dt to estimate the velocity vector. --- generate/ctest.c | 1 + generate/template/astronomy.c | 54 +++++++++++++++++++++++++++++++---- source/c/astronomy.c | 54 +++++++++++++++++++++++++++++++---- 3 files changed, 99 insertions(+), 10 deletions(-) diff --git a/generate/ctest.c b/generate/ctest.c index d01a8495..7a2815a1 100644 --- a/generate/ctest.c +++ b/generate/ctest.c @@ -4559,6 +4559,7 @@ static int BaryStateTest(void) CHECK(BaryStateBody(BODY_URANUS, "barystate/Uranus.txt", 1.71e-3, 1.03e-6)); CHECK(BaryStateBody(BODY_NEPTUNE, "barystate/Neptune.txt", 2.95e-3, 1.39e-6)); CHECK(BaryStateBody(BODY_PLUTO, "barystate/Pluto.txt", 2.05e-3, 1.91e-7)); + CHECK(BaryStateBody(BODY_MOON, "barystate/Moon.txt", 2.35e-5, 1.13e-6)); printf("C BaryStateTest: PASS\n"); fail: diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 977b142a..49daea9f 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -1865,6 +1865,40 @@ astro_vector_t Astronomy_GeoMoon(astro_time_t time) } +body_state_t GeoMoonState(astro_time_t time) +{ + /* + This is a hack, because trying to figure out how to derive a time + derivative for CalcMoon() would be extremely painful! + Calculate just before and just after the given time. + Average to find position, subtract to find velocity. + */ + const double dt = 1.0e-5; /* 0.864 seconds */ + astro_vector_t r1, r2; + astro_time_t t1, t2; + body_state_t s; + + t1 = Astronomy_AddDays(time, -dt); + t2 = Astronomy_AddDays(time, +dt); + + r1 = Astronomy_GeoMoon(t1); + r2 = Astronomy_GeoMoon(t2); + + /* The desired position is the average of the two calculated positions. */ + s.r.x = (r1.x + r2.x) / 2; + s.r.y = (r1.y + r2.y) / 2; + s.r.z = (r1.z + r2.z) / 2; + + /* The difference of the position vectors divided by the time span gives the velocity vector. */ + s.v.x = (r2.x - r1.x) / (2 * dt); + s.v.y = (r2.y - r1.y) / (2 * dt); + s.v.z = (r2.z - r1.z) / (2 * dt); + s.tt = time.tt; + + return s; +} + + /** * @brief Calculates the Moon's libration angles at a given moment in time. * @@ -3036,7 +3070,7 @@ astro_state_vector_t Astronomy_BaryState(astro_body_t body, astro_time_t time) { astro_state_vector_t state; body_state_t bary[5]; - body_state_t planet; + body_state_t planet, earth; if (body == BODY_SSB) { @@ -3076,14 +3110,12 @@ astro_state_vector_t Astronomy_BaryState(astro_body_t body, astro_time_t time) case BODY_NEPTUNE: return ExportState(bary[4], time); /* Handle the remaining VSOP bodies: Mercury, Venus, Earth, Mars. */ - /* Otherwise, we need to calculate the heliocentric state of the given body */ - /* and add the Sun's heliocentric state to obtain the body's barycentric state. */ - /* BarySun + HelioBody = BaryBody */ case BODY_MERCURY: case BODY_VENUS: case BODY_EARTH: case BODY_MARS: planet = CalcVsopPosVel(&vsop[body], time.tt); + /* BarySun + HelioBody = BaryBody */ state.x = bary[0].r.x + planet.r.x; state.y = bary[0].r.y + planet.r.y; state.z = bary[0].r.z + planet.r.z; @@ -3094,7 +3126,19 @@ astro_state_vector_t Astronomy_BaryState(astro_body_t body, astro_time_t time) state.status = ASTRO_SUCCESS; return state; - /* FIXFIXFIX - add support for BODY_MOON, BODY_EMB. */ + case BODY_MOON: + planet = GeoMoonState(time); + earth = CalcVsopPosVel(&vsop[BODY_EARTH], time.tt); + state.x = bary[0].r.x + earth.r.x + planet.r.x; + state.y = bary[0].r.y + earth.r.y + planet.r.y; + state.z = bary[0].r.z + earth.r.z + planet.r.z; + state.vx = bary[0].v.x + earth.v.x + planet.v.x; + state.vy = bary[0].v.y + earth.v.y + planet.v.y; + state.vz = bary[0].v.z + earth.v.z + planet.v.z; + state.t = time; + state.status = ASTRO_SUCCESS; + return state; + default: return StateVecError(ASTRO_INVALID_BODY, time); } diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 284e77e1..7812a1a8 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -2054,6 +2054,40 @@ astro_vector_t Astronomy_GeoMoon(astro_time_t time) } +body_state_t GeoMoonState(astro_time_t time) +{ + /* + This is a hack, because trying to figure out how to derive a time + derivative for CalcMoon() would be extremely painful! + Calculate just before and just after the given time. + Average to find position, subtract to find velocity. + */ + const double dt = 1.0e-5; /* 0.864 seconds */ + astro_vector_t r1, r2; + astro_time_t t1, t2; + body_state_t s; + + t1 = Astronomy_AddDays(time, -dt); + t2 = Astronomy_AddDays(time, +dt); + + r1 = Astronomy_GeoMoon(t1); + r2 = Astronomy_GeoMoon(t2); + + /* The desired position is the average of the two calculated positions. */ + s.r.x = (r1.x + r2.x) / 2; + s.r.y = (r1.y + r2.y) / 2; + s.r.z = (r1.z + r2.z) / 2; + + /* The difference of the position vectors divided by the time span gives the velocity vector. */ + s.v.x = (r2.x - r1.x) / (2 * dt); + s.v.y = (r2.y - r1.y) / (2 * dt); + s.v.z = (r2.z - r1.z) / (2 * dt); + s.tt = time.tt; + + return s; +} + + /** * @brief Calculates the Moon's libration angles at a given moment in time. * @@ -4275,7 +4309,7 @@ astro_state_vector_t Astronomy_BaryState(astro_body_t body, astro_time_t time) { astro_state_vector_t state; body_state_t bary[5]; - body_state_t planet; + body_state_t planet, earth; if (body == BODY_SSB) { @@ -4315,14 +4349,12 @@ astro_state_vector_t Astronomy_BaryState(astro_body_t body, astro_time_t time) case BODY_NEPTUNE: return ExportState(bary[4], time); /* Handle the remaining VSOP bodies: Mercury, Venus, Earth, Mars. */ - /* Otherwise, we need to calculate the heliocentric state of the given body */ - /* and add the Sun's heliocentric state to obtain the body's barycentric state. */ - /* BarySun + HelioBody = BaryBody */ case BODY_MERCURY: case BODY_VENUS: case BODY_EARTH: case BODY_MARS: planet = CalcVsopPosVel(&vsop[body], time.tt); + /* BarySun + HelioBody = BaryBody */ state.x = bary[0].r.x + planet.r.x; state.y = bary[0].r.y + planet.r.y; state.z = bary[0].r.z + planet.r.z; @@ -4333,7 +4365,19 @@ astro_state_vector_t Astronomy_BaryState(astro_body_t body, astro_time_t time) state.status = ASTRO_SUCCESS; return state; - /* FIXFIXFIX - add support for BODY_MOON, BODY_EMB. */ + case BODY_MOON: + planet = GeoMoonState(time); + earth = CalcVsopPosVel(&vsop[BODY_EARTH], time.tt); + state.x = bary[0].r.x + earth.r.x + planet.r.x; + state.y = bary[0].r.y + earth.r.y + planet.r.y; + state.z = bary[0].r.z + earth.r.z + planet.r.z; + state.vx = bary[0].v.x + earth.v.x + planet.v.x; + state.vy = bary[0].v.y + earth.v.y + planet.v.y; + state.vz = bary[0].v.z + earth.v.z + planet.v.z; + state.t = time; + state.status = ASTRO_SUCCESS; + return state; + default: return StateVecError(ASTRO_INVALID_BODY, time); }