mirror/astronomy
1
0
Fork 0
mirror of https://github.com/cosinekitty/astronomy.git synced 2026-05-19 06:17:03 -04:00
Code Issues Packages Projects Releases Wiki Activity

C BaryState: Calculate Moon's barycentric position and velocity vectors.

Browse Source 
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.
This commit is contained in:
Don Cross
2021-11-13 19:34:13 -05:00
parent 71cb92df08
commit 9869fd8bfc
3 changed files with 99 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
generate/ctest.c
View File

@@ -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:

54
generate/template/astronomy.c
View File

@@ -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);
}

54
source/c/astronomy.c
View File

@@ -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);
}