diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index 3b7e37b7..428ef6eb 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -4199,6 +4199,81 @@ def GeoVector(body, time, aberration): raise Error('Light-travel time solver did not converge: dt={}'.format(dt)) +def _ExportState(terse, time): + return StateVector( + terse.r.x, terse.r.y, terse.r.z, + terse.v.x, terse.v.y, terse.v.z, + time + ) + + +def BaryState(body, time): + """Calculates barycentric position and velocity vectors for the given body. + + Given a body and a time, calculates the barycentric position and velocity + vectors for the center of that body at that time. + The vectors are expressed in equatorial J2000 coordinates (EQJ). + + Parameters + ---------- + body : Body + The celestial body whose barycentric state vector is to be calculated. + Supported values are `Body.Sun`, `Body.SSB`, and all planets except Pluto: + `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, `Body.Jupiter`, + `Body.Saturn`, `Body.Uranus`, `Body.Neptune`. + time : Time + The date and time for which to calculate position and velocity. + + Returns + ------- + StateVector + An object that contains barycentric position and velocity vectors. + """ + # Trivial case: the solar sytem barycenter itself. + if body == Body.SSB: + return StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time) + + # Find the barycentric positions and velocities for the 5 major bodies. + bary = _major_bodies_t(time.tt) + + # If the caller is asking for one of the major bodies, + # we can immediately return the answer. + if body == Body.Sun: + return _ExportState(bary.Sun, time) + + if body == Body.Jupiter: + return _ExportState(bary.Jupiter, time) + + if body == Body.Saturn: + return _ExportState(bary.Saturn, time) + + if body == Body.Uranus: + return _ExportState(bary.Uranus, time) + + if body == Body.Neptune: + return _ExportState(bary.Neptune, time) + + if 0 <= body.value < len(_vsop): + # Handle the remaining VSOP bodies: Mercury, Venus, Earth, Mars. + # 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 + planet = _CalcVsopPosVel(_vsop[body.value], time.tt) + return StateVector( + bary.Sun.r.x + planet.r.x, + bary.Sun.r.y + planet.r.y, + bary.Sun.r.z + planet.r.z, + bary.Sun.v.x + planet.v.x, + bary.Sun.v.y + planet.v.y, + bary.Sun.v.z + planet.v.z, + time + ) + + # FIXFIXFIX: later, we can add support for Pluto, Moon, EMB, etc. + raise InvalidBodyError() + + def Equator(body, time, observer, ofdate, aberration): """Calculates equatorial coordinates of a celestial body as seen by an observer on the Earth's surface. diff --git a/generate/dotnet/csharp_test/csharp_test.cs b/generate/dotnet/csharp_test/csharp_test.cs index b90c51b9..29c7150e 100644 --- a/generate/dotnet/csharp_test/csharp_test.cs +++ b/generate/dotnet/csharp_test/csharp_test.cs @@ -2839,7 +2839,7 @@ namespace csharp_test { if (0 != BaryStateBody(Body.Sun, "../../barystate/Sun.txt", 1.23e-5, 1.14e-7)) return 1; if (0 != BaryStateBody(Body.Mercury, "../../barystate/Mercury.txt", 5.24e-5, 8.22e-6)) return 1; - if (0 != BaryStateBody(Body.Venus, "../../barystate/Venus.txt", 2.98e-5, 8.22e-6)) return 1; + if (0 != BaryStateBody(Body.Venus, "../../barystate/Venus.txt", 2.98e-5, 8.78e-7)) return 1; if (0 != BaryStateBody(Body.Earth, "../../barystate/Earth.txt", 2.30e-5, 1.09e-6)) return 1; if (0 != BaryStateBody(Body.Mars, "../../barystate/Mars.txt", 4.34e-5, 8.23e-7)) return 1; if (0 != BaryStateBody(Body.Jupiter, "../../barystate/Jupiter.txt", 3.74e-4, 1.78e-6)) return 1; diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 67e33f57..cf5ffd91 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -2928,10 +2928,10 @@ astro_state_vector_t Astronomy_BaryState(astro_body_t body, astro_time_t time) case BODY_URANUS: return ExportState(bary[3], 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 */ - /* Handle the remaining VSOP bodies: Mercury, Venus, Earth, Mars. */ case BODY_MERCURY: case BODY_VENUS: case BODY_EARTH: diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs index 0c44e2d1..f0e27478 100644 --- a/generate/template/astronomy.cs +++ b/generate/template/astronomy.cs @@ -3560,10 +3560,10 @@ $ASTRO_IAU_DATA() case Body.Neptune: return ExportState(bary.Neptune, 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 - // Handle the remaining VSOP bodies: Mercury, Venus, Earth, Mars. int bindex = (int)body; if (bindex >= 0 && bindex < vsop.Length) { @@ -3579,7 +3579,7 @@ $ASTRO_IAU_DATA() ); } - // FIXFIXFIX: later, we can add support for Pluton, Moon, EMB, etc. + // FIXFIXFIX: later, we can add support for Pluto, Moon, EMB, etc. throw new InvalidBodyException(body); } diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 1b832b85..fbb4a34d 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -2169,6 +2169,81 @@ def GeoVector(body, time, aberration): raise Error('Light-travel time solver did not converge: dt={}'.format(dt)) +def _ExportState(terse, time): + return StateVector( + terse.r.x, terse.r.y, terse.r.z, + terse.v.x, terse.v.y, terse.v.z, + time + ) + + +def BaryState(body, time): + """Calculates barycentric position and velocity vectors for the given body. + + Given a body and a time, calculates the barycentric position and velocity + vectors for the center of that body at that time. + The vectors are expressed in equatorial J2000 coordinates (EQJ). + + Parameters + ---------- + body : Body + The celestial body whose barycentric state vector is to be calculated. + Supported values are `Body.Sun`, `Body.SSB`, and all planets except Pluto: + `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, `Body.Jupiter`, + `Body.Saturn`, `Body.Uranus`, `Body.Neptune`. + time : Time + The date and time for which to calculate position and velocity. + + Returns + ------- + StateVector + An object that contains barycentric position and velocity vectors. + """ + # Trivial case: the solar sytem barycenter itself. + if body == Body.SSB: + return StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time) + + # Find the barycentric positions and velocities for the 5 major bodies. + bary = _major_bodies_t(time.tt) + + # If the caller is asking for one of the major bodies, + # we can immediately return the answer. + if body == Body.Sun: + return _ExportState(bary.Sun, time) + + if body == Body.Jupiter: + return _ExportState(bary.Jupiter, time) + + if body == Body.Saturn: + return _ExportState(bary.Saturn, time) + + if body == Body.Uranus: + return _ExportState(bary.Uranus, time) + + if body == Body.Neptune: + return _ExportState(bary.Neptune, time) + + if 0 <= body.value < len(_vsop): + # Handle the remaining VSOP bodies: Mercury, Venus, Earth, Mars. + # 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 + planet = _CalcVsopPosVel(_vsop[body.value], time.tt) + return StateVector( + bary.Sun.r.x + planet.r.x, + bary.Sun.r.y + planet.r.y, + bary.Sun.r.z + planet.r.z, + bary.Sun.v.x + planet.v.x, + bary.Sun.v.y + planet.v.y, + bary.Sun.v.z + planet.v.z, + time + ) + + # FIXFIXFIX: later, we can add support for Pluto, Moon, EMB, etc. + raise InvalidBodyError() + + def Equator(body, time, observer, ofdate, aberration): """Calculates equatorial coordinates of a celestial body as seen by an observer on the Earth's surface. diff --git a/generate/test.py b/generate/test.py index 08e9886a..86eb8962 100755 --- a/generate/test.py +++ b/generate/test.py @@ -1882,7 +1882,99 @@ def Issue103(): #----------------------------------------------------------------------------------------------------------- +class _bary_stats_t: + def __init__(self): + self.max_rdiff = 0.0 + self.max_vdiff = 0.0 + +def StateVectorDiff(vec, x, y, z): + dx = v(vec[0] - x) + dy = v(vec[1] - y) + dz = v(vec[2] - z) + ds = sqrt(dx*dx + dy*dy + dz*dz) + return ds + +def VerifyBaryState(stats, body, filename, lnum, time, pos, vel, r_thresh, v_thresh): + state = astronomy.BaryState(body, time) + rdiff = StateVectorDiff(pos, state.x, state.y, state.z) + if rdiff > stats.max_rdiff: + stats.max_rdiff = rdiff + + vdiff = StateVectorDiff(vel, state.vx, state.vy, state.vz) + if vdiff > stats.max_vdiff: + stats.max_vdiff = vdiff + + if rdiff > r_thresh: + print('PY VerifyBaryState({} line {}): EXCESSIVE position error = {:0.4e}'.format(filename, lnum, rdiff)) + return 1 + + if vdiff > v_thresh: + print('PY VerifyBaryState({} line {}): EXCESSIVE velocity error = {:0.4e}'.format(filename, lnum, vdiff)) + return 1 + + return 0 + + +def BaryStateBody(body, filename, r_thresh, v_thresh): + with open(filename, 'rt') as infile: + lnum = 0 + part = 0 + count = 0 + found_begin = False + stats = _bary_stats_t() + for line in infile: + line = line.rstrip() + lnum += 1 + if not found_begin: + if line == '$$SOE': + found_begin = True + elif line == '$$EOE': + break + else: + if part == 0: + # 2446545.000000000 = A.D. 1986-Apr-24 12:00:00.0000 TDB + tt = float(line.split()[0]) - 2451545.0 # convert JD to J2000 TT + time = astronomy.Time.FromTerrestrialTime(tt) + elif part == 1: + # X = 1.134408131605554E-03 Y =-2.590904586750408E-03 Z =-7.490427225904720E-05 + match = re.match(r'\s*X =\s*(\S+) Y =\s*(\S+) Z =\s*(\S+)', line) + if not match: + print('PY BaryStateBody({} line {}): cannot parse position vector.'.format(filename, lnum)) + return 1 + pos = [ float(match.group(1)), float(match.group(2)), float(match.group(3)) ] + else: # part == 2 + # VX= 9.148038778472862E-03 VY= 3.973823407182510E-03 VZ= 2.765660368640458E-04 + match = re.match(r'\s*VX=\s*(\S+) VY=\s*(\S+) VZ=\s*(\S+)', line) + if not match: + print('PY BaryStateBody({} line {}): cannot parse velocity vector.'.format(filename, lnum)) + return 1 + vel = [ float(match.group(1)), float(match.group(2)), float(match.group(3)) ] + if VerifyBaryState(stats, body, filename, lnum, time, pos, vel, r_thresh, v_thresh): + print('PY BaryStateBody({} line {}): FAILED VERIFICATION.'.format(filename, lnum)) + return 1 + count += 1 + part = (part + 1) % 3 + Debug('PY BaryStateBody({}): PASS - Tested {} cases. max rdiff={:0.3e}, vdiff={:0.3e}'.format(filename, count, stats.max_rdiff, stats.max_vdiff)) + return 0 + + +def BaryState(): + if BaryStateBody(astronomy.Body.Sun, 'barystate/Sun.txt', 1.23e-5, 1.14e-7): return 1 + if BaryStateBody(astronomy.Body.Mercury, 'barystate/Mercury.txt', 5.24e-5, 8.22e-6): return 1 + if BaryStateBody(astronomy.Body.Venus, 'barystate/Venus.txt', 2.98e-5, 8.78e-7): return 1 + if BaryStateBody(astronomy.Body.Earth, 'barystate/Earth.txt', 2.30e-5, 1.09e-6): return 1 + if BaryStateBody(astronomy.Body.Mars, 'barystate/Mars.txt', 4.34e-5, 8.23e-7): return 1 + if BaryStateBody(astronomy.Body.Jupiter, 'barystate/Jupiter.txt', 3.74e-4, 1.78e-6): return 1 + if BaryStateBody(astronomy.Body.Saturn, 'barystate/Saturn.txt', 1.07e-3, 1.71e-6): return 1 + if BaryStateBody(astronomy.Body.Uranus, 'barystate/Uranus.txt', 1.71e-3, 1.03e-6): return 1 + if BaryStateBody(astronomy.Body.Neptune, 'barystate/Neptune.txt', 2.95e-3, 1.39e-6): return 1 + print('PY BaryState: PASS') + return 0 + +#----------------------------------------------------------------------------------------------------------- + UnitTests = { + 'barystate': BaryState, 'constellation': Constellation, 'elongation': Elongation, 'geoid': Geoid, diff --git a/source/c/astronomy.c b/source/c/astronomy.c index b47438f2..bef507c8 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -4154,10 +4154,10 @@ astro_state_vector_t Astronomy_BaryState(astro_body_t body, astro_time_t time) case BODY_URANUS: return ExportState(bary[3], 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 */ - /* Handle the remaining VSOP bodies: Mercury, Venus, Earth, Mars. */ case BODY_MERCURY: case BODY_VENUS: case BODY_EARTH: diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs index 4f9d189a..0cd2a0ac 100644 --- a/source/csharp/astronomy.cs +++ b/source/csharp/astronomy.cs @@ -4760,10 +4760,10 @@ namespace CosineKitty case Body.Neptune: return ExportState(bary.Neptune, 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 - // Handle the remaining VSOP bodies: Mercury, Venus, Earth, Mars. int bindex = (int)body; if (bindex >= 0 && bindex < vsop.Length) { @@ -4779,7 +4779,7 @@ namespace CosineKitty ); } - // FIXFIXFIX: later, we can add support for Pluton, Moon, EMB, etc. + // FIXFIXFIX: later, we can add support for Pluto, Moon, EMB, etc. throw new InvalidBodyException(body); } diff --git a/source/python/README.md b/source/python/README.md index 068a112a..9cdbccb4 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -999,6 +999,25 @@ and the specified body as seen from the center of the Earth. --- + +### BaryState(body, time) + +**Calculates barycentric position and velocity vectors for the given body.** + +Given a body and a time, calculates the barycentric position and velocity +vectors for the center of that body at that time. +The vectors are expressed in equatorial J2000 coordinates (EQJ). + +| Type | Parameter | Description | +| --- | --- | --- | +| [`Body`](#Body) | `body` | The celestial body whose barycentric state vector is to be calculated. Supported values are `Body.Sun`, `Body.SSB`, and all planets except Pluto: `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, `Body.Jupiter`, `Body.Saturn`, `Body.Uranus`, `Body.Neptune`. | +| [`Time`](#Time) | `time` | The date and time for which to calculate position and velocity. | + +### Returns: [`StateVector`](#StateVector) +An object that contains barycentric position and velocity vectors. + +--- + ### BodyCode(name) diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 3b7e37b7..428ef6eb 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -4199,6 +4199,81 @@ def GeoVector(body, time, aberration): raise Error('Light-travel time solver did not converge: dt={}'.format(dt)) +def _ExportState(terse, time): + return StateVector( + terse.r.x, terse.r.y, terse.r.z, + terse.v.x, terse.v.y, terse.v.z, + time + ) + + +def BaryState(body, time): + """Calculates barycentric position and velocity vectors for the given body. + + Given a body and a time, calculates the barycentric position and velocity + vectors for the center of that body at that time. + The vectors are expressed in equatorial J2000 coordinates (EQJ). + + Parameters + ---------- + body : Body + The celestial body whose barycentric state vector is to be calculated. + Supported values are `Body.Sun`, `Body.SSB`, and all planets except Pluto: + `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, `Body.Jupiter`, + `Body.Saturn`, `Body.Uranus`, `Body.Neptune`. + time : Time + The date and time for which to calculate position and velocity. + + Returns + ------- + StateVector + An object that contains barycentric position and velocity vectors. + """ + # Trivial case: the solar sytem barycenter itself. + if body == Body.SSB: + return StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time) + + # Find the barycentric positions and velocities for the 5 major bodies. + bary = _major_bodies_t(time.tt) + + # If the caller is asking for one of the major bodies, + # we can immediately return the answer. + if body == Body.Sun: + return _ExportState(bary.Sun, time) + + if body == Body.Jupiter: + return _ExportState(bary.Jupiter, time) + + if body == Body.Saturn: + return _ExportState(bary.Saturn, time) + + if body == Body.Uranus: + return _ExportState(bary.Uranus, time) + + if body == Body.Neptune: + return _ExportState(bary.Neptune, time) + + if 0 <= body.value < len(_vsop): + # Handle the remaining VSOP bodies: Mercury, Venus, Earth, Mars. + # 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 + planet = _CalcVsopPosVel(_vsop[body.value], time.tt) + return StateVector( + bary.Sun.r.x + planet.r.x, + bary.Sun.r.y + planet.r.y, + bary.Sun.r.z + planet.r.z, + bary.Sun.v.x + planet.v.x, + bary.Sun.v.y + planet.v.y, + bary.Sun.v.z + planet.v.z, + time + ) + + # FIXFIXFIX: later, we can add support for Pluto, Moon, EMB, etc. + raise InvalidBodyError() + + def Equator(body, time, observer, ofdate, aberration): """Calculates equatorial coordinates of a celestial body as seen by an observer on the Earth's surface.