From d28f5ecff33458d96047dec4be835cccf32395b4 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sun, 22 May 2022 21:16:58 -0400 Subject: [PATCH] PY gravsim: passes unit tests The Python version of the GravitySimulator class is now passing all unit tests. This completes the initial coding. I still need to review documentation across all the language implementations. --- demo/python/astronomy.py | 21 ++- generate/template/astronomy.py | 21 ++- generate/test.js | 2 +- generate/test.py | 185 +++++++++++++++++++++++++-- source/python/astronomy/astronomy.py | 21 ++- 5 files changed, 228 insertions(+), 22 deletions(-) diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index 0920370e..e1c4049d 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -3461,6 +3461,13 @@ class _body_state_t: self.r = r self.v = v + def clone(self): + '''Make a copy of this body state.''' + return _body_state_t(self.tt, self.r.clone(), self.v.clone()) + + def __sub__(self, other): + return _body_state_t(self.tt, self.r - other.r, self.v - other.v) + def _CalcVsopPosVel(model, tt): t = tt / _DAYS_PER_MILLENNIUM @@ -3602,6 +3609,10 @@ class _TerseVector: self.y = y self.z = z + def clone(self): + '''Create a copy of this vector.''' + return _TerseVector(self.x, self.y, self.z) + @staticmethod def zero(): '''Return a zero vector.''' @@ -3692,6 +3703,10 @@ class _body_grav_calc_t: self.v = v # velocity [au/day] self.a = a # acceleration [au/day^2] + def clone(self): + '''Creates a copy of this gravity simulation state.''' + return _body_grav_calc_t(self.tt, self.r.clone(), self.v.clone(), self.a.clone()) + class _grav_sim_t: def __init__(self, bary, grav): @@ -10095,7 +10110,7 @@ class GravitySimulator: raise Error('Invalid body: {}'.format(body)) def _CalcBodyAccelerations(self): - for b in self.curr.bodies.values(): + for b in self.curr.bodies: b.a = _TerseVector.zero() _AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Sun ].r, _SUN_GM) _AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Mercury].r, _MERCURY_GM) @@ -10110,8 +10125,8 @@ class GravitySimulator: def _Duplicate(self): # Copy the current stateinto the previous state, so that both become the same moment in time. gravitators = {} - for body in self.curr.gravitators.values(): - gravitators[body] = self.curr.gravitators[body].clone() + for body, grav in self.curr.gravitators.items(): + gravitators[body] = grav.clone() bodies = [] for b in self.curr.bodies: diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 2d006a3a..c8fed3da 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -1610,6 +1610,13 @@ class _body_state_t: self.r = r self.v = v + def clone(self): + '''Make a copy of this body state.''' + return _body_state_t(self.tt, self.r.clone(), self.v.clone()) + + def __sub__(self, other): + return _body_state_t(self.tt, self.r - other.r, self.v - other.v) + def _CalcVsopPosVel(model, tt): t = tt / _DAYS_PER_MILLENNIUM @@ -1693,6 +1700,10 @@ class _TerseVector: self.y = y self.z = z + def clone(self): + '''Create a copy of this vector.''' + return _TerseVector(self.x, self.y, self.z) + @staticmethod def zero(): '''Return a zero vector.''' @@ -1783,6 +1794,10 @@ class _body_grav_calc_t: self.v = v # velocity [au/day] self.a = a # acceleration [au/day^2] + def clone(self): + '''Creates a copy of this gravity simulation state.''' + return _body_grav_calc_t(self.tt, self.r.clone(), self.v.clone(), self.a.clone()) + class _grav_sim_t: def __init__(self, bary, grav): @@ -7602,7 +7617,7 @@ class GravitySimulator: raise Error('Invalid body: {}'.format(body)) def _CalcBodyAccelerations(self): - for b in self.curr.bodies.values(): + for b in self.curr.bodies: b.a = _TerseVector.zero() _AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Sun ].r, _SUN_GM) _AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Mercury].r, _MERCURY_GM) @@ -7617,8 +7632,8 @@ class GravitySimulator: def _Duplicate(self): # Copy the current stateinto the previous state, so that both become the same moment in time. gravitators = {} - for body in self.curr.gravitators.values(): - gravitators[body] = self.curr.gravitators[body].clone() + for body, grav in self.curr.gravitators.items(): + gravitators[body] = grav.clone() bodies = [] for b in self.curr.bodies: diff --git a/generate/test.js b/generate/test.js index 8dd35bbf..f6e96c63 100644 --- a/generate/test.js +++ b/generate/test.js @@ -2908,7 +2908,7 @@ function GravSimFile(filename, originBody, nsteps, rthresh, vthresh) { max_vdiff = vdiff; prev = rec; } - Debug(`JS GravSimFile(${filename}): PASS - max pos error = ${max_rdiff.toFixed(4)} arcmin, max = ${max_vdiff.toFixed(4)} arcmin.`); + Debug(`JS GravSimFile(${filename}): PASS - max pos error = ${max_rdiff.toFixed(4)} arcmin, max vel error = ${max_vdiff.toFixed(4)} arcmin.`); return 0; } diff --git a/generate/test.py b/generate/test.py index 9ec33213..13a1660d 100755 --- a/generate/test.py +++ b/generate/test.py @@ -1948,13 +1948,18 @@ def VerifyState(func, stats, filename, lnum, time, pos, vel, r_thresh, v_thresh) return 0 -def VerifyStateBody(func, filename, r_thresh, v_thresh): +class JplStateRecord: + def __init__(self, lnum, state): + self.lnum = lnum + self.state = state + + +def JplHorizonsStateVectors(filename): + stateList = [] 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 @@ -1972,22 +1977,32 @@ def VerifyStateBody(func, filename, r_thresh, v_thresh): # 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 VerifyStateBody({} line {}): cannot parse position vector.'.format(filename, lnum)) + print('PY JplHorizonsStateVectors({} line {}): cannot parse position vector.'.format(filename, lnum)) return 1 - pos = [ float(match.group(1)), float(match.group(2)), float(match.group(3)) ] + rx, ry, rz = 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 VerifyStateBody({} line {}): cannot parse velocity vector.'.format(filename, lnum)) + print('PY JplHorizonsStateVectors({} line {}): cannot parse velocity vector.'.format(filename, lnum)) return 1 - vel = [ float(match.group(1)), float(match.group(2)), float(match.group(3)) ] - if VerifyState(func, stats, filename, lnum, time, pos, vel, r_thresh, v_thresh): - print('PY VerifyStateBody({} line {}): FAILED VERIFICATION.'.format(filename, lnum)) - return 1 - count += 1 + vx, vy, vz = float(match.group(1)), float(match.group(2)), float(match.group(3)) + yield JplStateRecord(lnum, astronomy.StateVector(rx, ry, rz, vx, vy, vz, time)) part = (part + 1) % 3 - Debug('PY VerifyStateBody({}): PASS - Tested {} cases. max rdiff={:0.3e}, vdiff={:0.3e}'.format(filename, count, stats.max_rdiff, stats.max_vdiff)) + + +def VerifyStateBody(func, filename, r_thresh, v_thresh): + stats = _bary_stats_t() + count = 0 + for rec in JplHorizonsStateVectors(filename): + time = rec.state.t + pos = [rec.state.x, rec.state.y, rec.state.z] + vel = [rec.state.vx, rec.state.vy, rec.state.vz] + if VerifyState(func, stats, filename, rec.lnum, time, pos, vel, r_thresh, v_thresh): + print('PY VerifyStateBody({} line {}): FAILED VERIFICATION.'.format(filename, rec.lnum)) + return 1 + count += 1 + Debug('PY VerifyStateBody({}): PASS - Tested {} cases. max rdiff={:0.3e}, vdiff={:0.3e}'.format(filename, count, stats.max_rdiff, stats.max_vdiff)) return 0 #----------------------------------------------------------------------------------------------------------- @@ -2508,6 +2523,151 @@ def Repr(): #----------------------------------------------------------------------------------------------------------- +def GravSimTest(): + Debug("") + + if 0 != GravSimEmpty("barystate/Sun.txt", astronomy.Body.SSB, astronomy.Body.Sun, 0.0269, 1.9635): return 1 + if 0 != GravSimEmpty("barystate/Mercury.txt", astronomy.Body.SSB, astronomy.Body.Mercury, 0.5725, 0.9332): return 1 + if 0 != GravSimEmpty("barystate/Venus.txt", astronomy.Body.SSB, astronomy.Body.Venus, 0.1433, 0.1458): return 1 + if 0 != GravSimEmpty("barystate/Earth.txt", astronomy.Body.SSB, astronomy.Body.Earth, 0.0651, 0.2098): return 1 + if 0 != GravSimEmpty("barystate/Mars.txt", astronomy.Body.SSB, astronomy.Body.Mars, 0.1150, 0.1896): return 1 + if 0 != GravSimEmpty("barystate/Jupiter.txt", astronomy.Body.SSB, astronomy.Body.Jupiter, 0.2546, 0.8831): return 1 + if 0 != GravSimEmpty("barystate/Saturn.txt", astronomy.Body.SSB, astronomy.Body.Saturn, 0.3660, 1.0818): return 1 + if 0 != GravSimEmpty("barystate/Uranus.txt", astronomy.Body.SSB, astronomy.Body.Uranus, 0.3107, 0.9321): return 1 + if 0 != GravSimEmpty("barystate/Neptune.txt", astronomy.Body.SSB, astronomy.Body.Neptune, 0.3382, 1.5586): return 1 + if 0 != GravSimEmpty("heliostate/Mercury.txt", astronomy.Body.Sun, astronomy.Body.Mercury, 0.5087, 0.9473): return 1 + if 0 != GravSimEmpty("heliostate/Venus.txt", astronomy.Body.Sun, astronomy.Body.Venus, 0.1214, 0.1543): return 1 + if 0 != GravSimEmpty("heliostate/Earth.txt", astronomy.Body.Sun, astronomy.Body.Earth, 0.0508, 0.2099): return 1 + if 0 != GravSimEmpty("heliostate/Mars.txt", astronomy.Body.Sun, astronomy.Body.Mars, 0.1085, 0.1927): return 1 + if 0 != GravSimEmpty("heliostate/Jupiter.txt", astronomy.Body.Sun, astronomy.Body.Jupiter, 0.2564, 0.8805): return 1 + if 0 != GravSimEmpty("heliostate/Saturn.txt", astronomy.Body.Sun, astronomy.Body.Saturn, 0.3664, 1.0826): return 1 + if 0 != GravSimEmpty("heliostate/Uranus.txt", astronomy.Body.Sun, astronomy.Body.Uranus, 0.3106, 0.9322): return 1 + if 0 != GravSimEmpty("heliostate/Neptune.txt", astronomy.Body.Sun, astronomy.Body.Neptune, 0.3381, 1.5584): return 1 + + Debug("") + nsteps = 20 + + if 0 != GravSimFile("barystate/Ceres.txt", astronomy.Body.SSB, nsteps, 0.6640, 0.6226): return 1 + if 0 != GravSimFile("barystate/Pallas.txt", astronomy.Body.SSB, nsteps, 0.4687, 0.3474): return 1 + if 0 != GravSimFile("barystate/Vesta.txt", astronomy.Body.SSB, nsteps, 0.5806, 0.5462): return 1 + if 0 != GravSimFile("barystate/Juno.txt", astronomy.Body.SSB, nsteps, 0.6760, 0.5750): return 1 + if 0 != GravSimFile("barystate/Bennu.txt", astronomy.Body.SSB, nsteps, 3.7444, 2.6581): return 1 + if 0 != GravSimFile("barystate/Halley.txt", astronomy.Body.SSB, nsteps, 0.0539, 0.0825): return 1 + if 0 != GravSimFile("heliostate/Ceres.txt", astronomy.Body.Sun, nsteps, 0.0445, 0.0355): return 1 + if 0 != GravSimFile("heliostate/Pallas.txt", astronomy.Body.Sun, nsteps, 0.1062, 0.0854): return 1 + if 0 != GravSimFile("heliostate/Vesta.txt", astronomy.Body.Sun, nsteps, 0.1432, 0.1308): return 1 + if 0 != GravSimFile("heliostate/Juno.txt", astronomy.Body.Sun, nsteps, 0.1554, 0.1328): return 1 + if 0 != GravSimFile("geostate/Ceres.txt", astronomy.Body.Earth, nsteps, 6.5689, 6.4797): return 1 + if 0 != GravSimFile("geostate/Pallas.txt", astronomy.Body.Earth, nsteps, 9.3288, 7.3533): return 1 + if 0 != GravSimFile("geostate/Vesta.txt", astronomy.Body.Earth, nsteps, 3.2980, 3.8863): return 1 + if 0 != GravSimFile("geostate/Juno.txt", astronomy.Body.Earth, nsteps, 6.0962, 7.7147): return 1 + + Debug("") + print("PY GravSimTest: PASS") + + return 0 + + +def GravSimEmpty(filename, origin, body, rthresh, vthresh): + max_rdiff = 0.0 + max_vdiff = 0.0 + sim = None + for rec in JplHorizonsStateVectors(filename): + if sim is None: + sim = astronomy.GravitySimulator(origin, rec.state.t, []) + sim.Update(rec.state.t) + calc = sim.SolarSystemBodyState(body) + if origin == astronomy.Body.SSB and body == astronomy.Body.Sun: + rdiff = SsbArcminPosError(rec.state, calc) + else: + rdiff = ArcminPosError(rec.state, calc) + if rdiff > rthresh: + print('PY GravSimEmpty({} line {}): excessive position error = {} arcmin.'.format(filename, rec.lnum, rdiff)) + return 1 + if rdiff > max_rdiff: + max_rdiff = rdiff + + vdiff = ArcminVelError(rec.state, calc) + if vdiff > vthresh: + print('PY GravSimEmpty({} line {}): excessive velocity error = {} arcmin.'.format(filename, rec.lnum, vdiff)) + return 1 + if vdiff > max_vdiff: + max_vdiff = vdiff + Debug('PY GravSimEmpty({:22s}): PASS - max pos error = {:0.4f} arcmin, max vel error = {:0.4f} arcmin.'.format(filename, max_rdiff, max_vdiff)) + return 0 + + +def GravSimFile(filename, originBody, nsteps, rthresh, vthresh): + sim = None + max_rdiff = 0.0 + max_vdiff = 0.0 + for rec in JplHorizonsStateVectors(filename): + if sim == None: + sim = astronomy.GravitySimulator(originBody, rec.state.t, [rec.state]) + time = rec.state.t + smallBodyArray = sim.Update(time) + else: + tt1 = prev.state.t.tt + tt2 = rec.state.t.tt + dt = (tt2 - tt1) / nsteps + for k in range(1, nsteps+1): + time = astronomy.Time.FromTerrestrialTime(tt1 + k*dt) + smallBodyArray = sim.Update(time) + if len(smallBodyArray) != 1: + print('PY GravSimFile({} line {}): unexpected smallBodyArray.length = {}'.format(filename, rec.lnum, len(smallBodyArray))) + return 1 + if time.tt != sim.Time().tt: + print('PY GravSimFile({} line {}): expected {} but simulator reports {}'.format(filename, rec.lnum, time, sim.Time())) + return 1 + rdiff = ArcminPosError(rec.state, smallBodyArray[0]) + if rdiff > rthresh: + print('PY GravSimFile({} line {}): excessive position error = {}'.format(filename, rec.lnum, rdiff)) + return 1 + if rdiff > max_rdiff: + max_rdiff = rdiff + vdiff = ArcminVelError(rec.state, smallBodyArray[0]) + if vdiff > vthresh: + print('PY GravSimFile({} line {}): excessive position error = {}'.format(filename, rec.lnum, vdiff)) + return 1 + if vdiff > max_vdiff: + max_vdiff = vdiff + prev = rec + Debug('PY GravSimFile({:22s}): PASS - max pos error = {:0.4f} arcmin, max vel error = {:0.4f} arcmin.'.format(filename, max_rdiff, max_vdiff)) + return 0 + + +def SsbArcminPosError(correct, calc): + # Scale the SSB based on 1 AU, not on its absolute magnitude, which can become very close to zero. + dx = calc.x - correct.x + dy = calc.y - correct.y + dz = calc.z - correct.z + diffSquared = dx*dx + dy*dy + dz*dz + radians = sqrt(diffSquared) + return 60.0 * math.degrees(radians) + + +def ArcminPosError(correct, calc): + dx = calc.x - correct.x + dy = calc.y - correct.y + dz = calc.z - correct.z + diffSquared = dx*dx + dy*dy + dz*dz + magSquared = correct.x*correct.x + correct.y*correct.y + correct.z*correct.z + radians = sqrt(diffSquared / magSquared) + return 60.0 * math.degrees(radians) + + +def ArcminVelError(correct, calc): + dx = calc.vx - correct.vx + dy = calc.vy - correct.vy + dz = calc.vz - correct.vz + diffSquared = dx*dx + dy*dy + dz*dz + magSquared = correct.vx*correct.vx + correct.vy*correct.vy + correct.vz*correct.vz + radians = sqrt(diffSquared / magSquared) + return 60.0 * math.degrees(radians) + + +#----------------------------------------------------------------------------------------------------------- + UnitTests = { 'aberration': Aberration, 'axis': Axis, @@ -2516,6 +2676,7 @@ UnitTests = { 'elongation': Elongation, 'geoid': Geoid, 'global_solar_eclipse': GlobalSolarEclipse, + 'gravsim': GravSimTest, 'heliostate': HelioState, 'issue_103': Issue103, 'jupiter_moons': JupiterMoons, diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py index 0920370e..e1c4049d 100644 --- a/source/python/astronomy/astronomy.py +++ b/source/python/astronomy/astronomy.py @@ -3461,6 +3461,13 @@ class _body_state_t: self.r = r self.v = v + def clone(self): + '''Make a copy of this body state.''' + return _body_state_t(self.tt, self.r.clone(), self.v.clone()) + + def __sub__(self, other): + return _body_state_t(self.tt, self.r - other.r, self.v - other.v) + def _CalcVsopPosVel(model, tt): t = tt / _DAYS_PER_MILLENNIUM @@ -3602,6 +3609,10 @@ class _TerseVector: self.y = y self.z = z + def clone(self): + '''Create a copy of this vector.''' + return _TerseVector(self.x, self.y, self.z) + @staticmethod def zero(): '''Return a zero vector.''' @@ -3692,6 +3703,10 @@ class _body_grav_calc_t: self.v = v # velocity [au/day] self.a = a # acceleration [au/day^2] + def clone(self): + '''Creates a copy of this gravity simulation state.''' + return _body_grav_calc_t(self.tt, self.r.clone(), self.v.clone(), self.a.clone()) + class _grav_sim_t: def __init__(self, bary, grav): @@ -10095,7 +10110,7 @@ class GravitySimulator: raise Error('Invalid body: {}'.format(body)) def _CalcBodyAccelerations(self): - for b in self.curr.bodies.values(): + for b in self.curr.bodies: b.a = _TerseVector.zero() _AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Sun ].r, _SUN_GM) _AddAcceleration(b.a, b.r, self.curr.gravitators[Body.Mercury].r, _MERCURY_GM) @@ -10110,8 +10125,8 @@ class GravitySimulator: def _Duplicate(self): # Copy the current stateinto the previous state, so that both become the same moment in time. gravitators = {} - for body in self.curr.gravitators.values(): - gravitators[body] = self.curr.gravitators[body].clone() + for body, grav in self.curr.gravitators.items(): + gravitators[body] = grav.clone() bodies = [] for b in self.curr.bodies: