From 23b45dd075e149ea1e0f727519685dfc71fa2dc8 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Tue, 21 Feb 2023 10:56:05 -0500 Subject: [PATCH] Improved Pluto state table in Python code. The generated code for the Pluto state table in Python now uses a class `_pstate` for better type checking. It also makes the code easier to understand. Moved class _TerseVector higher in the source file to reduce the need for quoted forward type declarations. --- demo/python/astronomy.py | 227 ++++++++++++++------------- generate/codegen.c | 4 +- generate/template/astronomy.py | 123 ++++++++------- source/python/astronomy/astronomy.py | 227 ++++++++++++++------------- 4 files changed, 292 insertions(+), 289 deletions(-) diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index 3b999760..77d6810f 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -43,6 +43,50 @@ def _cbrt(x: float) -> float: return -((-x) ** (1.0 / 3.0)) return x ** (1.0 / 3.0) +class _TerseVector: + '''A 3D vector that is not attached to a time. Used privately inside this module for conciseness.''' + + def __init__(self, x: float, y: float, z: float) -> None: + self.x = x + self.y = y + self.z = z + + def clone(self) -> "_TerseVector": + '''Create a copy of this vector.''' + return _TerseVector(self.x, self.y, self.z) + + @staticmethod + def zero() -> "_TerseVector": + '''Return a zero vector.''' + return _TerseVector(0.0, 0.0, 0.0) + + def ToAstroVector(self, time: "Time") -> "Vector": + '''Convert _TerseVector object to Vector object.''' + return Vector(self.x, self.y, self.z, time) + + def quadrature(self) -> float: + '''Return magnitude squared of this vector.''' + return self.x**2 + self.y**2 + self.z**2 + + def mean(self, other: "_TerseVector") -> "_TerseVector": + '''Return the average of this vector and another vector.''' + return _TerseVector((self.x + other.x)/2.0, (self.y + other.y)/2.0, (self.z + other.z)/2.0) + + def __add__(self, other: "_TerseVector") -> "_TerseVector": + return _TerseVector(self.x + other.x, self.y + other.y, self.z + other.z) + + def __sub__(self, other: "_TerseVector") -> "_TerseVector": + return _TerseVector(self.x - other.x, self.y - other.y, self.z - other.z) + + def __mul__(self, scalar: float) -> "_TerseVector": + return _TerseVector(scalar * self.x, scalar * self.y, scalar * self.z) + + def __rmul__(self, scalar: float) -> "_TerseVector": + return _TerseVector(scalar * self.x, scalar * self.y, scalar * self.z) + + def __truediv__(self, scalar: float) -> "_TerseVector": + return _TerseVector(self.x / scalar, self.y / scalar, self.z / scalar) + KM_PER_AU = 1.4959787069098932e+8 # The number of kilometers per astronomical unit. C_AUDAY = 173.1446326846693 # The speed of light expressed in astronomical units per day. AU_PER_LY = 63241.07708807546 # The number of astronomical units in one light-year. @@ -3096,14 +3140,14 @@ def _VsopDeriv(formula: _vsop_formula_t, t: float) -> float: _DAYS_PER_MILLENNIUM = 365250.0 -def _VsopRotate(eclip: "_TerseVector") -> "_TerseVector": +def _VsopRotate(eclip: _TerseVector) -> _TerseVector: # Convert ecliptic cartesian coordinates to equatorial cartesian coordinates. x = eclip.x + 0.000000440360*eclip.y - 0.000000190919*eclip.z y = -0.000000479966*eclip.x + 0.917482137087*eclip.y - 0.397776982902*eclip.z z = 0.397776982902*eclip.y + 0.917482137087*eclip.z return _TerseVector(x, y, z) -def _VsopSphereToRect(lon: float, lat: float, rad: float) -> "_TerseVector": +def _VsopSphereToRect(lon: float, lat: float, rad: float) -> _TerseVector: # Convert spherical coordinates to cartesian coordinates. r_coslat = rad * math.cos(lat) return _TerseVector( @@ -3121,7 +3165,7 @@ def _CalcVsop(model: _vsop_model_t, time: Time) -> Vector: return _VsopRotate(eclip).ToAstroVector(time) class _body_state_t: - def __init__(self, tt: float, r: "_TerseVector", v: "_TerseVector") -> None: + def __init__(self, tt: float, r: _TerseVector, v: _TerseVector) -> None: self.tt = tt self.r = r self.v = v @@ -3210,114 +3254,72 @@ def _CalcEarth(time: Time) -> Vector: #---------------------------------------------------------------------------- # BEGIN Pluto Integrator +class _pstate: + def __init__(self, tt: float, pos: Tuple[float, float, float], vel: Tuple[float, float, float]) -> None: + self.tt = tt + self.pos = _TerseVector(pos[0], pos[1], pos[2]) + self.vel = _TerseVector(vel[0], vel[1], vel[2]) + _PLUTO_NUM_STATES = 51 _PLUTO_TIME_STEP = 29200 _PLUTO_DT = 146 _PLUTO_NSTEPS = 201 -_PlutoStateTable: List[Tuple[float, Tuple[float, float, float], Tuple[float, float, float]]] = [ - ( -730000.0, (-26.118207232108, -14.376168177825, 3.384402515299), ( 1.6339372163656e-03, -2.7861699588508e-03, -1.3585880229445e-03)) -, ( -700800.0, ( 41.974905202127, -0.448502952929, -12.770351505989), ( 7.3458569351457e-04, 2.2785014891658e-03, 4.8619778602049e-04)) -, ( -671600.0, ( 14.706930780744, 44.269110540027, 9.353698474772), (-2.1000147999800e-03, 2.2295915939915e-04, 7.0143443551414e-04)) -, ( -642400.0, (-29.441003929957, -6.430161530570, 6.858481011305), ( 8.4495803960544e-04, -3.0783914758711e-03, -1.2106305981192e-03)) -, ( -613200.0, ( 39.444396946234, -6.557989760571, -13.913760296463), ( 1.1480029005873e-03, 2.2400006880665e-03, 3.5168075922288e-04)) -, ( -584000.0, ( 20.230380950700, 43.266966657189, 7.382966091923), (-1.9754081700585e-03, 5.3457141292226e-04, 7.5929169129793e-04)) -, ( -554800.0, (-30.658325364620, 2.093818874552, 9.880531138071), ( 6.1010603013347e-05, -3.1326500935382e-03, -9.9346125151067e-04)) -, ( -525600.0, ( 35.737703251673, -12.587706024764, -14.677847247563), ( 1.5802939375649e-03, 2.1347678412429e-03, 1.9074436384343e-04)) -, ( -496400.0, ( 25.466295188546, 41.367478338417, 5.216476873382), (-1.8054401046468e-03, 8.3283083599510e-04, 8.0260156912107e-04)) -, ( -467200.0, (-29.847174904071, 10.636426313081, 12.297904180106), (-6.3257063052907e-04, -2.9969577578221e-03, -7.4476074151596e-04)) -, ( -438000.0, ( 30.774692107687, -18.236637015304, -14.945535879896), ( 2.0113162005465e-03, 1.9353827024189e-03, -2.0937793168297e-06)) -, ( -408800.0, ( 30.243153324028, 38.656267888503, 2.938501750218), (-1.6052508674468e-03, 1.1183495337525e-03, 8.3333973416824e-04)) -, ( -379600.0, (-27.288984772533, 18.643162147874, 14.023633623329), (-1.1856388898191e-03, -2.7170609282181e-03, -4.9015526126399e-04)) -, ( -350400.0, ( 24.519605196774, -23.245756064727, -14.626862367368), ( 2.4322321483154e-03, 1.6062008146048e-03, -2.3369181613312e-04)) -, ( -321200.0, ( 34.505274805875, 35.125338586954, 0.557361475637), (-1.3824391637782e-03, 1.3833397561817e-03, 8.4823598806262e-04)) -, ( -292000.0, (-23.275363915119, 25.818514298769, 15.055381588598), (-1.6062295460975e-03, -2.3395961498533e-03, -2.4377362639479e-04)) -, ( -262800.0, ( 17.050384798092, -27.180376290126, -13.608963321694), ( 2.8175521080578e-03, 1.1358749093955e-03, -4.9548725258825e-04)) -, ( -233600.0, ( 38.093671910285, 30.880588383337, -1.843688067413), (-1.1317697153459e-03, 1.6128814698472e-03, 8.4177586176055e-04)) -, ( -204400.0, (-18.197852930878, 31.932869934309, 15.438294826279), (-1.9117272501813e-03, -1.9146495909842e-03, -1.9657304369835e-05)) -, ( -175200.0, ( 8.528924039997, -29.618422200048, -11.805400994258), ( 3.1034370787005e-03, 5.1393633292430e-04, -7.7293066202546e-04)) -, ( -146000.0, ( 40.946857258640, 25.904973592021, -4.256336240499), (-8.3652705194051e-04, 1.8129497136404e-03, 8.1564228273060e-04)) -, ( -116800.0, (-12.326958895325, 36.881883446292, 15.217158258711), (-2.1166103705038e-03, -1.4814420035990e-03, 1.7401209844705e-04)) -, ( -87600.0, ( -0.633258375909, -30.018759794709, -9.171932874950), ( 3.2016994581737e-03, -2.5279858672148e-04, -1.0411088271861e-03)) -, ( -58400.0, ( 42.936048423883, 20.344685584452, -6.588027007912), (-5.0525450073192e-04, 1.9910074335507e-03, 7.7440196540269e-04)) -, ( -29200.0, ( -5.975910552974, 40.611809958460, 14.470131723673), (-2.2184202156107e-03, -1.0562361130164e-03, 3.3652250216211e-04)) -, ( 0.0, ( -9.875369580774, -27.978926224737, -5.753711824704), ( 3.0287533248818e-03, -1.1276087003636e-03, -1.2651326732361e-03)) -, ( 29200.0, ( 43.958831986165, 14.214147973292, -8.808306227163), (-1.4717608981871e-04, 2.1404187242141e-03, 7.1486567806614e-04)) -, ( 58400.0, ( 0.678136763520, 43.094461639362, 13.243238780721), (-2.2358226110718e-03, -6.3233636090933e-04, 4.7664798895648e-04)) -, ( 87600.0, (-18.282602096834, -23.305039586660, -1.766620508028), ( 2.5567245263557e-03, -1.9902940754171e-03, -1.3943491701082e-03)) -, ( 116800.0, ( 43.873338744526, 7.700705617215, -10.814273666425), ( 2.3174803055677e-04, 2.2402163127924e-03, 6.2988756452032e-04)) -, ( 146000.0, ( 7.392949027906, 44.382678951534, 11.629500214854), (-2.1932815453830e-03, -2.1751799585364e-04, 5.9556516201114e-04)) -, ( 175200.0, (-24.981690229261, -16.204012851426, 2.466457544298), ( 1.8193989149580e-03, -2.6765419531201e-03, -1.3848283502247e-03)) -, ( 204400.0, ( 42.530187039511, 0.845935508021, -12.554907527683), ( 6.5059779150669e-04, 2.2725657282262e-03, 5.1133743202822e-04)) -, ( 233600.0, ( 13.999526486822, 44.462363044894, 9.669418486465), (-2.1079296569252e-03, 1.7533423831993e-04, 6.9128485798076e-04)) -, ( 262800.0, (-29.184024803031, -7.371243995762, 6.493275957928), ( 9.3581363109681e-04, -3.0610357109184e-03, -1.2364201089345e-03)) -, ( 292000.0, ( 39.831980671753, -6.078405766765, -13.909815358656), ( 1.1117769689167e-03, 2.2362097830152e-03, 3.6230548231153e-04)) -, ( 321200.0, ( 20.294955108476, 43.417190420251, 7.450091985932), (-1.9742157451535e-03, 5.3102050468554e-04, 7.5938408813008e-04)) -, ( 350400.0, (-30.669992302160, 2.318743558955, 9.973480913858), ( 4.5605107450676e-05, -3.1308219926928e-03, -9.9066533301924e-04)) -, ( 379600.0, ( 35.626122155983, -12.897647509224, -14.777586508444), ( 1.6015684949743e-03, 2.1171931182284e-03, 1.8002516202204e-04)) -, ( 408800.0, ( 26.133186148561, 41.232139187599, 5.006401326220), (-1.7857704419579e-03, 8.6046232702817e-04, 8.0614690298954e-04)) -, ( 438000.0, (-29.576740229230, 11.863535943587, 12.631323039872), (-7.2292830060955e-04, -2.9587820140709e-03, -7.0824296450300e-04)) -, ( 467200.0, ( 29.910805787391, -19.159019294000, -15.013363865194), ( 2.0871080437997e-03, 1.8848372554514e-03, -3.8528655083926e-05)) -, ( 496400.0, ( 31.375957451819, 38.050372720763, 2.433138343754), (-1.5546055556611e-03, 1.1699815465629e-03, 8.3565439266001e-04)) -, ( 525600.0, (-26.360071336928, 20.662505904952, 14.414696258958), (-1.3142373118349e-03, -2.6236647854842e-03, -4.2542017598193e-04)) -, ( 554800.0, ( 22.599441488648, -24.508879898306, -14.484045731468), ( 2.5454108304806e-03, 1.4917058755191e-03, -3.0243665086079e-04)) -, ( 584000.0, ( 35.877864013014, 33.894226366071, -0.224524636277), (-1.2941245730845e-03, 1.4560427668319e-03, 8.4762160640137e-04)) -, ( 613200.0, (-21.538149762417, 28.204068269761, 15.321973799534), (-1.7312117409010e-03, -2.1939631314577e-03, -1.6316913275180e-04)) -, ( 642400.0, ( 13.971521374415, -28.339941764789, -13.083792871886), ( 2.9334630526035e-03, 9.1860931752944e-04, -5.9939422488627e-04)) -, ( 671600.0, ( 39.526942044143, 28.939897360110, -2.872799527539), (-1.0068481658095e-03, 1.7021132888090e-03, 8.3578230511981e-04)) -, ( 700800.0, (-15.576200701394, 34.399412961275, 15.466033737854), (-2.0098814612884e-03, -1.7191109825989e-03, 7.0414782780416e-05)) -, ( 730000.0, ( 4.243252837090, -30.118201690825, -10.707441231349), ( 3.1725847067411e-03, 1.6098461202270e-04, -9.0672150593868e-04)) +_PlutoStateTable: List[_pstate] = [ + _pstate( -730000.0, (-26.118207232108, -14.376168177825, 3.384402515299), ( 1.6339372163656e-03, -2.7861699588508e-03, -1.3585880229445e-03)) +, _pstate( -700800.0, ( 41.974905202127, -0.448502952929, -12.770351505989), ( 7.3458569351457e-04, 2.2785014891658e-03, 4.8619778602049e-04)) +, _pstate( -671600.0, ( 14.706930780744, 44.269110540027, 9.353698474772), (-2.1000147999800e-03, 2.2295915939915e-04, 7.0143443551414e-04)) +, _pstate( -642400.0, (-29.441003929957, -6.430161530570, 6.858481011305), ( 8.4495803960544e-04, -3.0783914758711e-03, -1.2106305981192e-03)) +, _pstate( -613200.0, ( 39.444396946234, -6.557989760571, -13.913760296463), ( 1.1480029005873e-03, 2.2400006880665e-03, 3.5168075922288e-04)) +, _pstate( -584000.0, ( 20.230380950700, 43.266966657189, 7.382966091923), (-1.9754081700585e-03, 5.3457141292226e-04, 7.5929169129793e-04)) +, _pstate( -554800.0, (-30.658325364620, 2.093818874552, 9.880531138071), ( 6.1010603013347e-05, -3.1326500935382e-03, -9.9346125151067e-04)) +, _pstate( -525600.0, ( 35.737703251673, -12.587706024764, -14.677847247563), ( 1.5802939375649e-03, 2.1347678412429e-03, 1.9074436384343e-04)) +, _pstate( -496400.0, ( 25.466295188546, 41.367478338417, 5.216476873382), (-1.8054401046468e-03, 8.3283083599510e-04, 8.0260156912107e-04)) +, _pstate( -467200.0, (-29.847174904071, 10.636426313081, 12.297904180106), (-6.3257063052907e-04, -2.9969577578221e-03, -7.4476074151596e-04)) +, _pstate( -438000.0, ( 30.774692107687, -18.236637015304, -14.945535879896), ( 2.0113162005465e-03, 1.9353827024189e-03, -2.0937793168297e-06)) +, _pstate( -408800.0, ( 30.243153324028, 38.656267888503, 2.938501750218), (-1.6052508674468e-03, 1.1183495337525e-03, 8.3333973416824e-04)) +, _pstate( -379600.0, (-27.288984772533, 18.643162147874, 14.023633623329), (-1.1856388898191e-03, -2.7170609282181e-03, -4.9015526126399e-04)) +, _pstate( -350400.0, ( 24.519605196774, -23.245756064727, -14.626862367368), ( 2.4322321483154e-03, 1.6062008146048e-03, -2.3369181613312e-04)) +, _pstate( -321200.0, ( 34.505274805875, 35.125338586954, 0.557361475637), (-1.3824391637782e-03, 1.3833397561817e-03, 8.4823598806262e-04)) +, _pstate( -292000.0, (-23.275363915119, 25.818514298769, 15.055381588598), (-1.6062295460975e-03, -2.3395961498533e-03, -2.4377362639479e-04)) +, _pstate( -262800.0, ( 17.050384798092, -27.180376290126, -13.608963321694), ( 2.8175521080578e-03, 1.1358749093955e-03, -4.9548725258825e-04)) +, _pstate( -233600.0, ( 38.093671910285, 30.880588383337, -1.843688067413), (-1.1317697153459e-03, 1.6128814698472e-03, 8.4177586176055e-04)) +, _pstate( -204400.0, (-18.197852930878, 31.932869934309, 15.438294826279), (-1.9117272501813e-03, -1.9146495909842e-03, -1.9657304369835e-05)) +, _pstate( -175200.0, ( 8.528924039997, -29.618422200048, -11.805400994258), ( 3.1034370787005e-03, 5.1393633292430e-04, -7.7293066202546e-04)) +, _pstate( -146000.0, ( 40.946857258640, 25.904973592021, -4.256336240499), (-8.3652705194051e-04, 1.8129497136404e-03, 8.1564228273060e-04)) +, _pstate( -116800.0, (-12.326958895325, 36.881883446292, 15.217158258711), (-2.1166103705038e-03, -1.4814420035990e-03, 1.7401209844705e-04)) +, _pstate( -87600.0, ( -0.633258375909, -30.018759794709, -9.171932874950), ( 3.2016994581737e-03, -2.5279858672148e-04, -1.0411088271861e-03)) +, _pstate( -58400.0, ( 42.936048423883, 20.344685584452, -6.588027007912), (-5.0525450073192e-04, 1.9910074335507e-03, 7.7440196540269e-04)) +, _pstate( -29200.0, ( -5.975910552974, 40.611809958460, 14.470131723673), (-2.2184202156107e-03, -1.0562361130164e-03, 3.3652250216211e-04)) +, _pstate( 0.0, ( -9.875369580774, -27.978926224737, -5.753711824704), ( 3.0287533248818e-03, -1.1276087003636e-03, -1.2651326732361e-03)) +, _pstate( 29200.0, ( 43.958831986165, 14.214147973292, -8.808306227163), (-1.4717608981871e-04, 2.1404187242141e-03, 7.1486567806614e-04)) +, _pstate( 58400.0, ( 0.678136763520, 43.094461639362, 13.243238780721), (-2.2358226110718e-03, -6.3233636090933e-04, 4.7664798895648e-04)) +, _pstate( 87600.0, (-18.282602096834, -23.305039586660, -1.766620508028), ( 2.5567245263557e-03, -1.9902940754171e-03, -1.3943491701082e-03)) +, _pstate( 116800.0, ( 43.873338744526, 7.700705617215, -10.814273666425), ( 2.3174803055677e-04, 2.2402163127924e-03, 6.2988756452032e-04)) +, _pstate( 146000.0, ( 7.392949027906, 44.382678951534, 11.629500214854), (-2.1932815453830e-03, -2.1751799585364e-04, 5.9556516201114e-04)) +, _pstate( 175200.0, (-24.981690229261, -16.204012851426, 2.466457544298), ( 1.8193989149580e-03, -2.6765419531201e-03, -1.3848283502247e-03)) +, _pstate( 204400.0, ( 42.530187039511, 0.845935508021, -12.554907527683), ( 6.5059779150669e-04, 2.2725657282262e-03, 5.1133743202822e-04)) +, _pstate( 233600.0, ( 13.999526486822, 44.462363044894, 9.669418486465), (-2.1079296569252e-03, 1.7533423831993e-04, 6.9128485798076e-04)) +, _pstate( 262800.0, (-29.184024803031, -7.371243995762, 6.493275957928), ( 9.3581363109681e-04, -3.0610357109184e-03, -1.2364201089345e-03)) +, _pstate( 292000.0, ( 39.831980671753, -6.078405766765, -13.909815358656), ( 1.1117769689167e-03, 2.2362097830152e-03, 3.6230548231153e-04)) +, _pstate( 321200.0, ( 20.294955108476, 43.417190420251, 7.450091985932), (-1.9742157451535e-03, 5.3102050468554e-04, 7.5938408813008e-04)) +, _pstate( 350400.0, (-30.669992302160, 2.318743558955, 9.973480913858), ( 4.5605107450676e-05, -3.1308219926928e-03, -9.9066533301924e-04)) +, _pstate( 379600.0, ( 35.626122155983, -12.897647509224, -14.777586508444), ( 1.6015684949743e-03, 2.1171931182284e-03, 1.8002516202204e-04)) +, _pstate( 408800.0, ( 26.133186148561, 41.232139187599, 5.006401326220), (-1.7857704419579e-03, 8.6046232702817e-04, 8.0614690298954e-04)) +, _pstate( 438000.0, (-29.576740229230, 11.863535943587, 12.631323039872), (-7.2292830060955e-04, -2.9587820140709e-03, -7.0824296450300e-04)) +, _pstate( 467200.0, ( 29.910805787391, -19.159019294000, -15.013363865194), ( 2.0871080437997e-03, 1.8848372554514e-03, -3.8528655083926e-05)) +, _pstate( 496400.0, ( 31.375957451819, 38.050372720763, 2.433138343754), (-1.5546055556611e-03, 1.1699815465629e-03, 8.3565439266001e-04)) +, _pstate( 525600.0, (-26.360071336928, 20.662505904952, 14.414696258958), (-1.3142373118349e-03, -2.6236647854842e-03, -4.2542017598193e-04)) +, _pstate( 554800.0, ( 22.599441488648, -24.508879898306, -14.484045731468), ( 2.5454108304806e-03, 1.4917058755191e-03, -3.0243665086079e-04)) +, _pstate( 584000.0, ( 35.877864013014, 33.894226366071, -0.224524636277), (-1.2941245730845e-03, 1.4560427668319e-03, 8.4762160640137e-04)) +, _pstate( 613200.0, (-21.538149762417, 28.204068269761, 15.321973799534), (-1.7312117409010e-03, -2.1939631314577e-03, -1.6316913275180e-04)) +, _pstate( 642400.0, ( 13.971521374415, -28.339941764789, -13.083792871886), ( 2.9334630526035e-03, 9.1860931752944e-04, -5.9939422488627e-04)) +, _pstate( 671600.0, ( 39.526942044143, 28.939897360110, -2.872799527539), (-1.0068481658095e-03, 1.7021132888090e-03, 8.3578230511981e-04)) +, _pstate( 700800.0, (-15.576200701394, 34.399412961275, 15.466033737854), (-2.0098814612884e-03, -1.7191109825989e-03, 7.0414782780416e-05)) +, _pstate( 730000.0, ( 4.243252837090, -30.118201690825, -10.707441231349), ( 3.1725847067411e-03, 1.6098461202270e-04, -9.0672150593868e-04)) ] -class _TerseVector: - def __init__(self, x: float, y: float, z: float) -> None: - self.x = x - self.y = y - self.z = z - - def clone(self) -> "_TerseVector": - '''Create a copy of this vector.''' - return _TerseVector(self.x, self.y, self.z) - - @staticmethod - def zero() -> "_TerseVector": - '''Return a zero vector.''' - return _TerseVector(0.0, 0.0, 0.0) - - def ToAstroVector(self, time: Time) -> Vector: - '''Convert _TerseVector object to Vector object.''' - return Vector(self.x, self.y, self.z, time) - - def quadrature(self) -> float: - '''Return magnitude squared of this vector.''' - return self.x**2 + self.y**2 + self.z**2 - - def mean(self, other: "_TerseVector") -> "_TerseVector": - '''Return the average of this vector and another vector.''' - return _TerseVector((self.x + other.x)/2.0, (self.y + other.y)/2.0, (self.z + other.z)/2.0) - - def __add__(self, other: "_TerseVector") -> "_TerseVector": - return _TerseVector(self.x + other.x, self.y + other.y, self.z + other.z) - - def __sub__(self, other: "_TerseVector") -> "_TerseVector": - return _TerseVector(self.x - other.x, self.y - other.y, self.z - other.z) - - def __mul__(self, scalar: float) -> "_TerseVector": - return _TerseVector(scalar * self.x, scalar * self.y, scalar * self.z) - - def __rmul__(self, scalar: float) -> "_TerseVector": - return _TerseVector(scalar * self.x, scalar * self.y, scalar * self.z) - - def __truediv__(self, scalar: float) -> "_TerseVector": - return _TerseVector(self.x / scalar, self.y / scalar, self.z / scalar) - - -def _BodyStateFromTable(entry: Tuple[float, Tuple[float, float, float], Tuple[float, float, float]]) -> _body_state_t: - ( tt, (rx, ry, rz), (vx, vy, vz) ) = entry - return _body_state_t(tt, _TerseVector(rx, ry, rz), _TerseVector(vx, vy, vz)) - - def _AdjustBarycenterPosVel(ssb: _body_state_t, tt: float, body: Body, planet_gm: float) -> _body_state_t: shift = planet_gm / (planet_gm + _SUN_GM) planet = _CalcVsopPosVel(_vsop[body.value], tt) @@ -3429,22 +3431,21 @@ def _ClampIndex(frac: float, nsteps: int) -> int: return index -def _GravFromState(entry: Tuple[float, Tuple[float, float, float], Tuple[float, float, float]]) -> _grav_sim_t: - state = _BodyStateFromTable(entry) - bary = _major_bodies_t(state.tt) - r = state.r + bary.Sun.r - v = state.v + bary.Sun.v +def _GravFromState(entry: _pstate) -> _grav_sim_t: + bary = _major_bodies_t(entry.tt) + r = entry.pos + bary.Sun.r + v = entry.vel + bary.Sun.v a = bary.Acceleration(r) - grav = _body_grav_calc_t(state.tt, r, v, a) + grav = _body_grav_calc_t(entry.tt, r, v, a) return _grav_sim_t(bary, grav) def _GetSegment(cache: List[Optional[List[_body_grav_calc_t]]], tt: float) -> Optional[List[_body_grav_calc_t]]: - if (tt < _PlutoStateTable[0][0]) or (tt > _PlutoStateTable[_PLUTO_NUM_STATES-1][0]): + if (tt < _PlutoStateTable[0].tt) or (tt > _PlutoStateTable[_PLUTO_NUM_STATES-1].tt): # Don't bother calculating a segment. Let the caller crawl backward/forward to this time. return None - seg_index = _ClampIndex((tt - _PlutoStateTable[0][0]) / _PLUTO_TIME_STEP, _PLUTO_NUM_STATES-1) + seg_index = _ClampIndex((tt - _PlutoStateTable[0].tt) / _PLUTO_TIME_STEP, _PLUTO_NUM_STATES-1) if cache[seg_index] is None: seg = cache[seg_index] = [ _GravFromState(_PlutoStateTable[seg_index]).grav ] @@ -3481,7 +3482,7 @@ def _GetSegment(cache: List[Optional[List[_body_grav_calc_t]]], tt: float) -> Op return cache[seg_index] -def _CalcPlutoOneWay(entry: Tuple[float, Tuple[float, float, float], Tuple[float, float, float]], target_tt: float, dt: float) -> _grav_sim_t: +def _CalcPlutoOneWay(entry: _pstate, target_tt: float, dt: float) -> _grav_sim_t: sim = _GravFromState(entry) n = math.ceil((target_tt - sim.grav.tt) / dt) for i in range(n): @@ -3496,7 +3497,7 @@ def _CalcPluto(time: Time, helio: bool) -> StateVector: # The target time is outside the year range 0000..4000. # Calculate it by crawling backward from 0000 or forward from 4000. # FIXFIXFIX - This is super slow. Could optimize this with extra caching if needed. - if time.tt < _PlutoStateTable[0][0]: + if time.tt < _PlutoStateTable[0].tt: sim = _CalcPlutoOneWay(_PlutoStateTable[0], time.tt, -_PLUTO_DT) else: sim = _CalcPlutoOneWay(_PlutoStateTable[_PLUTO_NUM_STATES-1], time.tt, +_PLUTO_DT) diff --git a/generate/codegen.c b/generate/codegen.c index 11aa7ce5..6341922c 100644 --- a/generate/codegen.c +++ b/generate/codegen.c @@ -1572,7 +1572,7 @@ static int PlutoStateTable_Python(cg_context_t *context, const top_model_t *mode fprintf(context->outfile, "_PLUTO_DT = %d\n", PLUTO_DT); fprintf(context->outfile, "_PLUTO_NSTEPS = %d\n", PLUTO_NSTEPS); fprintf(context->outfile, "\n"); - fprintf(context->outfile, "_PlutoStateTable: List[Tuple[float, Tuple[float, float, float], Tuple[float, float, float]]] = [\n"); + fprintf(context->outfile, "_PlutoStateTable: List[_pstate] = [\n"); for (i=0; i < PLUTO_NUM_STATES; ++i) { @@ -1580,7 +1580,7 @@ static int PlutoStateTable_Python(cg_context_t *context, const top_model_t *mode CHECK(TopPosition(model, tt, &equ)); fprintf(context->outfile, - "%c (%10.1lf, (%16.12lf, %16.12lf, %16.12lf), (%20.13le, %20.13le, %20.13le))\n", + "%c _pstate(%10.1lf, (%16.12lf, %16.12lf, %16.12lf), (%20.13le, %20.13le, %20.13le))\n", (i==0 ? ' ' : ','), tt, equ.x, equ.y, equ.z, equ.vx, equ.vy, equ.vz); } diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index fe0b0d55..5fab56b3 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -43,6 +43,50 @@ def _cbrt(x: float) -> float: return -((-x) ** (1.0 / 3.0)) return x ** (1.0 / 3.0) +class _TerseVector: + '''A 3D vector that is not attached to a time. Used privately inside this module for conciseness.''' + + def __init__(self, x: float, y: float, z: float) -> None: + self.x = x + self.y = y + self.z = z + + def clone(self) -> "_TerseVector": + '''Create a copy of this vector.''' + return _TerseVector(self.x, self.y, self.z) + + @staticmethod + def zero() -> "_TerseVector": + '''Return a zero vector.''' + return _TerseVector(0.0, 0.0, 0.0) + + def ToAstroVector(self, time: "Time") -> "Vector": + '''Convert _TerseVector object to Vector object.''' + return Vector(self.x, self.y, self.z, time) + + def quadrature(self) -> float: + '''Return magnitude squared of this vector.''' + return self.x**2 + self.y**2 + self.z**2 + + def mean(self, other: "_TerseVector") -> "_TerseVector": + '''Return the average of this vector and another vector.''' + return _TerseVector((self.x + other.x)/2.0, (self.y + other.y)/2.0, (self.z + other.z)/2.0) + + def __add__(self, other: "_TerseVector") -> "_TerseVector": + return _TerseVector(self.x + other.x, self.y + other.y, self.z + other.z) + + def __sub__(self, other: "_TerseVector") -> "_TerseVector": + return _TerseVector(self.x - other.x, self.y - other.y, self.z - other.z) + + def __mul__(self, scalar: float) -> "_TerseVector": + return _TerseVector(scalar * self.x, scalar * self.y, scalar * self.z) + + def __rmul__(self, scalar: float) -> "_TerseVector": + return _TerseVector(scalar * self.x, scalar * self.y, scalar * self.z) + + def __truediv__(self, scalar: float) -> "_TerseVector": + return _TerseVector(self.x / scalar, self.y / scalar, self.z / scalar) + KM_PER_AU = 1.4959787069098932e+8 # The number of kilometers per astronomical unit. C_AUDAY = 173.1446326846693 # The speed of light expressed in astronomical units per day. AU_PER_LY = 63241.07708807546 # The number of astronomical units in one light-year. @@ -1781,14 +1825,14 @@ def _VsopDeriv(formula: _vsop_formula_t, t: float) -> float: _DAYS_PER_MILLENNIUM = 365250.0 -def _VsopRotate(eclip: "_TerseVector") -> "_TerseVector": +def _VsopRotate(eclip: _TerseVector) -> _TerseVector: # Convert ecliptic cartesian coordinates to equatorial cartesian coordinates. x = eclip.x + 0.000000440360*eclip.y - 0.000000190919*eclip.z y = -0.000000479966*eclip.x + 0.917482137087*eclip.y - 0.397776982902*eclip.z z = 0.397776982902*eclip.y + 0.917482137087*eclip.z return _TerseVector(x, y, z) -def _VsopSphereToRect(lon: float, lat: float, rad: float) -> "_TerseVector": +def _VsopSphereToRect(lon: float, lat: float, rad: float) -> _TerseVector: # Convert spherical coordinates to cartesian coordinates. r_coslat = rad * math.cos(lat) return _TerseVector( @@ -1806,7 +1850,7 @@ def _CalcVsop(model: _vsop_model_t, time: Time) -> Vector: return _VsopRotate(eclip).ToAstroVector(time) class _body_state_t: - def __init__(self, tt: float, r: "_TerseVector", v: "_TerseVector") -> None: + def __init__(self, tt: float, r: _TerseVector, v: _TerseVector) -> None: self.tt = tt self.r = r self.v = v @@ -1895,56 +1939,14 @@ def _CalcEarth(time: Time) -> Vector: #---------------------------------------------------------------------------- # BEGIN Pluto Integrator +class _pstate: + def __init__(self, tt: float, pos: Tuple[float, float, float], vel: Tuple[float, float, float]) -> None: + self.tt = tt + self.pos = _TerseVector(pos[0], pos[1], pos[2]) + self.vel = _TerseVector(vel[0], vel[1], vel[2]) + $ASTRO_PLUTO_TABLE() -class _TerseVector: - def __init__(self, x: float, y: float, z: float) -> None: - self.x = x - self.y = y - self.z = z - - def clone(self) -> "_TerseVector": - '''Create a copy of this vector.''' - return _TerseVector(self.x, self.y, self.z) - - @staticmethod - def zero() -> "_TerseVector": - '''Return a zero vector.''' - return _TerseVector(0.0, 0.0, 0.0) - - def ToAstroVector(self, time: Time) -> Vector: - '''Convert _TerseVector object to Vector object.''' - return Vector(self.x, self.y, self.z, time) - - def quadrature(self) -> float: - '''Return magnitude squared of this vector.''' - return self.x**2 + self.y**2 + self.z**2 - - def mean(self, other: "_TerseVector") -> "_TerseVector": - '''Return the average of this vector and another vector.''' - return _TerseVector((self.x + other.x)/2.0, (self.y + other.y)/2.0, (self.z + other.z)/2.0) - - def __add__(self, other: "_TerseVector") -> "_TerseVector": - return _TerseVector(self.x + other.x, self.y + other.y, self.z + other.z) - - def __sub__(self, other: "_TerseVector") -> "_TerseVector": - return _TerseVector(self.x - other.x, self.y - other.y, self.z - other.z) - - def __mul__(self, scalar: float) -> "_TerseVector": - return _TerseVector(scalar * self.x, scalar * self.y, scalar * self.z) - - def __rmul__(self, scalar: float) -> "_TerseVector": - return _TerseVector(scalar * self.x, scalar * self.y, scalar * self.z) - - def __truediv__(self, scalar: float) -> "_TerseVector": - return _TerseVector(self.x / scalar, self.y / scalar, self.z / scalar) - - -def _BodyStateFromTable(entry: Tuple[float, Tuple[float, float, float], Tuple[float, float, float]]) -> _body_state_t: - ( tt, (rx, ry, rz), (vx, vy, vz) ) = entry - return _body_state_t(tt, _TerseVector(rx, ry, rz), _TerseVector(vx, vy, vz)) - - def _AdjustBarycenterPosVel(ssb: _body_state_t, tt: float, body: Body, planet_gm: float) -> _body_state_t: shift = planet_gm / (planet_gm + _SUN_GM) planet = _CalcVsopPosVel(_vsop[body.value], tt) @@ -2056,22 +2058,21 @@ def _ClampIndex(frac: float, nsteps: int) -> int: return index -def _GravFromState(entry: Tuple[float, Tuple[float, float, float], Tuple[float, float, float]]) -> _grav_sim_t: - state = _BodyStateFromTable(entry) - bary = _major_bodies_t(state.tt) - r = state.r + bary.Sun.r - v = state.v + bary.Sun.v +def _GravFromState(entry: _pstate) -> _grav_sim_t: + bary = _major_bodies_t(entry.tt) + r = entry.pos + bary.Sun.r + v = entry.vel + bary.Sun.v a = bary.Acceleration(r) - grav = _body_grav_calc_t(state.tt, r, v, a) + grav = _body_grav_calc_t(entry.tt, r, v, a) return _grav_sim_t(bary, grav) def _GetSegment(cache: List[Optional[List[_body_grav_calc_t]]], tt: float) -> Optional[List[_body_grav_calc_t]]: - if (tt < _PlutoStateTable[0][0]) or (tt > _PlutoStateTable[_PLUTO_NUM_STATES-1][0]): + if (tt < _PlutoStateTable[0].tt) or (tt > _PlutoStateTable[_PLUTO_NUM_STATES-1].tt): # Don't bother calculating a segment. Let the caller crawl backward/forward to this time. return None - seg_index = _ClampIndex((tt - _PlutoStateTable[0][0]) / _PLUTO_TIME_STEP, _PLUTO_NUM_STATES-1) + seg_index = _ClampIndex((tt - _PlutoStateTable[0].tt) / _PLUTO_TIME_STEP, _PLUTO_NUM_STATES-1) if cache[seg_index] is None: seg = cache[seg_index] = [ _GravFromState(_PlutoStateTable[seg_index]).grav ] @@ -2108,7 +2109,7 @@ def _GetSegment(cache: List[Optional[List[_body_grav_calc_t]]], tt: float) -> Op return cache[seg_index] -def _CalcPlutoOneWay(entry: Tuple[float, Tuple[float, float, float], Tuple[float, float, float]], target_tt: float, dt: float) -> _grav_sim_t: +def _CalcPlutoOneWay(entry: _pstate, target_tt: float, dt: float) -> _grav_sim_t: sim = _GravFromState(entry) n = math.ceil((target_tt - sim.grav.tt) / dt) for i in range(n): @@ -2123,7 +2124,7 @@ def _CalcPluto(time: Time, helio: bool) -> StateVector: # The target time is outside the year range 0000..4000. # Calculate it by crawling backward from 0000 or forward from 4000. # FIXFIXFIX - This is super slow. Could optimize this with extra caching if needed. - if time.tt < _PlutoStateTable[0][0]: + if time.tt < _PlutoStateTable[0].tt: sim = _CalcPlutoOneWay(_PlutoStateTable[0], time.tt, -_PLUTO_DT) else: sim = _CalcPlutoOneWay(_PlutoStateTable[_PLUTO_NUM_STATES-1], time.tt, +_PLUTO_DT) diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py index 3b999760..77d6810f 100644 --- a/source/python/astronomy/astronomy.py +++ b/source/python/astronomy/astronomy.py @@ -43,6 +43,50 @@ def _cbrt(x: float) -> float: return -((-x) ** (1.0 / 3.0)) return x ** (1.0 / 3.0) +class _TerseVector: + '''A 3D vector that is not attached to a time. Used privately inside this module for conciseness.''' + + def __init__(self, x: float, y: float, z: float) -> None: + self.x = x + self.y = y + self.z = z + + def clone(self) -> "_TerseVector": + '''Create a copy of this vector.''' + return _TerseVector(self.x, self.y, self.z) + + @staticmethod + def zero() -> "_TerseVector": + '''Return a zero vector.''' + return _TerseVector(0.0, 0.0, 0.0) + + def ToAstroVector(self, time: "Time") -> "Vector": + '''Convert _TerseVector object to Vector object.''' + return Vector(self.x, self.y, self.z, time) + + def quadrature(self) -> float: + '''Return magnitude squared of this vector.''' + return self.x**2 + self.y**2 + self.z**2 + + def mean(self, other: "_TerseVector") -> "_TerseVector": + '''Return the average of this vector and another vector.''' + return _TerseVector((self.x + other.x)/2.0, (self.y + other.y)/2.0, (self.z + other.z)/2.0) + + def __add__(self, other: "_TerseVector") -> "_TerseVector": + return _TerseVector(self.x + other.x, self.y + other.y, self.z + other.z) + + def __sub__(self, other: "_TerseVector") -> "_TerseVector": + return _TerseVector(self.x - other.x, self.y - other.y, self.z - other.z) + + def __mul__(self, scalar: float) -> "_TerseVector": + return _TerseVector(scalar * self.x, scalar * self.y, scalar * self.z) + + def __rmul__(self, scalar: float) -> "_TerseVector": + return _TerseVector(scalar * self.x, scalar * self.y, scalar * self.z) + + def __truediv__(self, scalar: float) -> "_TerseVector": + return _TerseVector(self.x / scalar, self.y / scalar, self.z / scalar) + KM_PER_AU = 1.4959787069098932e+8 # The number of kilometers per astronomical unit. C_AUDAY = 173.1446326846693 # The speed of light expressed in astronomical units per day. AU_PER_LY = 63241.07708807546 # The number of astronomical units in one light-year. @@ -3096,14 +3140,14 @@ def _VsopDeriv(formula: _vsop_formula_t, t: float) -> float: _DAYS_PER_MILLENNIUM = 365250.0 -def _VsopRotate(eclip: "_TerseVector") -> "_TerseVector": +def _VsopRotate(eclip: _TerseVector) -> _TerseVector: # Convert ecliptic cartesian coordinates to equatorial cartesian coordinates. x = eclip.x + 0.000000440360*eclip.y - 0.000000190919*eclip.z y = -0.000000479966*eclip.x + 0.917482137087*eclip.y - 0.397776982902*eclip.z z = 0.397776982902*eclip.y + 0.917482137087*eclip.z return _TerseVector(x, y, z) -def _VsopSphereToRect(lon: float, lat: float, rad: float) -> "_TerseVector": +def _VsopSphereToRect(lon: float, lat: float, rad: float) -> _TerseVector: # Convert spherical coordinates to cartesian coordinates. r_coslat = rad * math.cos(lat) return _TerseVector( @@ -3121,7 +3165,7 @@ def _CalcVsop(model: _vsop_model_t, time: Time) -> Vector: return _VsopRotate(eclip).ToAstroVector(time) class _body_state_t: - def __init__(self, tt: float, r: "_TerseVector", v: "_TerseVector") -> None: + def __init__(self, tt: float, r: _TerseVector, v: _TerseVector) -> None: self.tt = tt self.r = r self.v = v @@ -3210,114 +3254,72 @@ def _CalcEarth(time: Time) -> Vector: #---------------------------------------------------------------------------- # BEGIN Pluto Integrator +class _pstate: + def __init__(self, tt: float, pos: Tuple[float, float, float], vel: Tuple[float, float, float]) -> None: + self.tt = tt + self.pos = _TerseVector(pos[0], pos[1], pos[2]) + self.vel = _TerseVector(vel[0], vel[1], vel[2]) + _PLUTO_NUM_STATES = 51 _PLUTO_TIME_STEP = 29200 _PLUTO_DT = 146 _PLUTO_NSTEPS = 201 -_PlutoStateTable: List[Tuple[float, Tuple[float, float, float], Tuple[float, float, float]]] = [ - ( -730000.0, (-26.118207232108, -14.376168177825, 3.384402515299), ( 1.6339372163656e-03, -2.7861699588508e-03, -1.3585880229445e-03)) -, ( -700800.0, ( 41.974905202127, -0.448502952929, -12.770351505989), ( 7.3458569351457e-04, 2.2785014891658e-03, 4.8619778602049e-04)) -, ( -671600.0, ( 14.706930780744, 44.269110540027, 9.353698474772), (-2.1000147999800e-03, 2.2295915939915e-04, 7.0143443551414e-04)) -, ( -642400.0, (-29.441003929957, -6.430161530570, 6.858481011305), ( 8.4495803960544e-04, -3.0783914758711e-03, -1.2106305981192e-03)) -, ( -613200.0, ( 39.444396946234, -6.557989760571, -13.913760296463), ( 1.1480029005873e-03, 2.2400006880665e-03, 3.5168075922288e-04)) -, ( -584000.0, ( 20.230380950700, 43.266966657189, 7.382966091923), (-1.9754081700585e-03, 5.3457141292226e-04, 7.5929169129793e-04)) -, ( -554800.0, (-30.658325364620, 2.093818874552, 9.880531138071), ( 6.1010603013347e-05, -3.1326500935382e-03, -9.9346125151067e-04)) -, ( -525600.0, ( 35.737703251673, -12.587706024764, -14.677847247563), ( 1.5802939375649e-03, 2.1347678412429e-03, 1.9074436384343e-04)) -, ( -496400.0, ( 25.466295188546, 41.367478338417, 5.216476873382), (-1.8054401046468e-03, 8.3283083599510e-04, 8.0260156912107e-04)) -, ( -467200.0, (-29.847174904071, 10.636426313081, 12.297904180106), (-6.3257063052907e-04, -2.9969577578221e-03, -7.4476074151596e-04)) -, ( -438000.0, ( 30.774692107687, -18.236637015304, -14.945535879896), ( 2.0113162005465e-03, 1.9353827024189e-03, -2.0937793168297e-06)) -, ( -408800.0, ( 30.243153324028, 38.656267888503, 2.938501750218), (-1.6052508674468e-03, 1.1183495337525e-03, 8.3333973416824e-04)) -, ( -379600.0, (-27.288984772533, 18.643162147874, 14.023633623329), (-1.1856388898191e-03, -2.7170609282181e-03, -4.9015526126399e-04)) -, ( -350400.0, ( 24.519605196774, -23.245756064727, -14.626862367368), ( 2.4322321483154e-03, 1.6062008146048e-03, -2.3369181613312e-04)) -, ( -321200.0, ( 34.505274805875, 35.125338586954, 0.557361475637), (-1.3824391637782e-03, 1.3833397561817e-03, 8.4823598806262e-04)) -, ( -292000.0, (-23.275363915119, 25.818514298769, 15.055381588598), (-1.6062295460975e-03, -2.3395961498533e-03, -2.4377362639479e-04)) -, ( -262800.0, ( 17.050384798092, -27.180376290126, -13.608963321694), ( 2.8175521080578e-03, 1.1358749093955e-03, -4.9548725258825e-04)) -, ( -233600.0, ( 38.093671910285, 30.880588383337, -1.843688067413), (-1.1317697153459e-03, 1.6128814698472e-03, 8.4177586176055e-04)) -, ( -204400.0, (-18.197852930878, 31.932869934309, 15.438294826279), (-1.9117272501813e-03, -1.9146495909842e-03, -1.9657304369835e-05)) -, ( -175200.0, ( 8.528924039997, -29.618422200048, -11.805400994258), ( 3.1034370787005e-03, 5.1393633292430e-04, -7.7293066202546e-04)) -, ( -146000.0, ( 40.946857258640, 25.904973592021, -4.256336240499), (-8.3652705194051e-04, 1.8129497136404e-03, 8.1564228273060e-04)) -, ( -116800.0, (-12.326958895325, 36.881883446292, 15.217158258711), (-2.1166103705038e-03, -1.4814420035990e-03, 1.7401209844705e-04)) -, ( -87600.0, ( -0.633258375909, -30.018759794709, -9.171932874950), ( 3.2016994581737e-03, -2.5279858672148e-04, -1.0411088271861e-03)) -, ( -58400.0, ( 42.936048423883, 20.344685584452, -6.588027007912), (-5.0525450073192e-04, 1.9910074335507e-03, 7.7440196540269e-04)) -, ( -29200.0, ( -5.975910552974, 40.611809958460, 14.470131723673), (-2.2184202156107e-03, -1.0562361130164e-03, 3.3652250216211e-04)) -, ( 0.0, ( -9.875369580774, -27.978926224737, -5.753711824704), ( 3.0287533248818e-03, -1.1276087003636e-03, -1.2651326732361e-03)) -, ( 29200.0, ( 43.958831986165, 14.214147973292, -8.808306227163), (-1.4717608981871e-04, 2.1404187242141e-03, 7.1486567806614e-04)) -, ( 58400.0, ( 0.678136763520, 43.094461639362, 13.243238780721), (-2.2358226110718e-03, -6.3233636090933e-04, 4.7664798895648e-04)) -, ( 87600.0, (-18.282602096834, -23.305039586660, -1.766620508028), ( 2.5567245263557e-03, -1.9902940754171e-03, -1.3943491701082e-03)) -, ( 116800.0, ( 43.873338744526, 7.700705617215, -10.814273666425), ( 2.3174803055677e-04, 2.2402163127924e-03, 6.2988756452032e-04)) -, ( 146000.0, ( 7.392949027906, 44.382678951534, 11.629500214854), (-2.1932815453830e-03, -2.1751799585364e-04, 5.9556516201114e-04)) -, ( 175200.0, (-24.981690229261, -16.204012851426, 2.466457544298), ( 1.8193989149580e-03, -2.6765419531201e-03, -1.3848283502247e-03)) -, ( 204400.0, ( 42.530187039511, 0.845935508021, -12.554907527683), ( 6.5059779150669e-04, 2.2725657282262e-03, 5.1133743202822e-04)) -, ( 233600.0, ( 13.999526486822, 44.462363044894, 9.669418486465), (-2.1079296569252e-03, 1.7533423831993e-04, 6.9128485798076e-04)) -, ( 262800.0, (-29.184024803031, -7.371243995762, 6.493275957928), ( 9.3581363109681e-04, -3.0610357109184e-03, -1.2364201089345e-03)) -, ( 292000.0, ( 39.831980671753, -6.078405766765, -13.909815358656), ( 1.1117769689167e-03, 2.2362097830152e-03, 3.6230548231153e-04)) -, ( 321200.0, ( 20.294955108476, 43.417190420251, 7.450091985932), (-1.9742157451535e-03, 5.3102050468554e-04, 7.5938408813008e-04)) -, ( 350400.0, (-30.669992302160, 2.318743558955, 9.973480913858), ( 4.5605107450676e-05, -3.1308219926928e-03, -9.9066533301924e-04)) -, ( 379600.0, ( 35.626122155983, -12.897647509224, -14.777586508444), ( 1.6015684949743e-03, 2.1171931182284e-03, 1.8002516202204e-04)) -, ( 408800.0, ( 26.133186148561, 41.232139187599, 5.006401326220), (-1.7857704419579e-03, 8.6046232702817e-04, 8.0614690298954e-04)) -, ( 438000.0, (-29.576740229230, 11.863535943587, 12.631323039872), (-7.2292830060955e-04, -2.9587820140709e-03, -7.0824296450300e-04)) -, ( 467200.0, ( 29.910805787391, -19.159019294000, -15.013363865194), ( 2.0871080437997e-03, 1.8848372554514e-03, -3.8528655083926e-05)) -, ( 496400.0, ( 31.375957451819, 38.050372720763, 2.433138343754), (-1.5546055556611e-03, 1.1699815465629e-03, 8.3565439266001e-04)) -, ( 525600.0, (-26.360071336928, 20.662505904952, 14.414696258958), (-1.3142373118349e-03, -2.6236647854842e-03, -4.2542017598193e-04)) -, ( 554800.0, ( 22.599441488648, -24.508879898306, -14.484045731468), ( 2.5454108304806e-03, 1.4917058755191e-03, -3.0243665086079e-04)) -, ( 584000.0, ( 35.877864013014, 33.894226366071, -0.224524636277), (-1.2941245730845e-03, 1.4560427668319e-03, 8.4762160640137e-04)) -, ( 613200.0, (-21.538149762417, 28.204068269761, 15.321973799534), (-1.7312117409010e-03, -2.1939631314577e-03, -1.6316913275180e-04)) -, ( 642400.0, ( 13.971521374415, -28.339941764789, -13.083792871886), ( 2.9334630526035e-03, 9.1860931752944e-04, -5.9939422488627e-04)) -, ( 671600.0, ( 39.526942044143, 28.939897360110, -2.872799527539), (-1.0068481658095e-03, 1.7021132888090e-03, 8.3578230511981e-04)) -, ( 700800.0, (-15.576200701394, 34.399412961275, 15.466033737854), (-2.0098814612884e-03, -1.7191109825989e-03, 7.0414782780416e-05)) -, ( 730000.0, ( 4.243252837090, -30.118201690825, -10.707441231349), ( 3.1725847067411e-03, 1.6098461202270e-04, -9.0672150593868e-04)) +_PlutoStateTable: List[_pstate] = [ + _pstate( -730000.0, (-26.118207232108, -14.376168177825, 3.384402515299), ( 1.6339372163656e-03, -2.7861699588508e-03, -1.3585880229445e-03)) +, _pstate( -700800.0, ( 41.974905202127, -0.448502952929, -12.770351505989), ( 7.3458569351457e-04, 2.2785014891658e-03, 4.8619778602049e-04)) +, _pstate( -671600.0, ( 14.706930780744, 44.269110540027, 9.353698474772), (-2.1000147999800e-03, 2.2295915939915e-04, 7.0143443551414e-04)) +, _pstate( -642400.0, (-29.441003929957, -6.430161530570, 6.858481011305), ( 8.4495803960544e-04, -3.0783914758711e-03, -1.2106305981192e-03)) +, _pstate( -613200.0, ( 39.444396946234, -6.557989760571, -13.913760296463), ( 1.1480029005873e-03, 2.2400006880665e-03, 3.5168075922288e-04)) +, _pstate( -584000.0, ( 20.230380950700, 43.266966657189, 7.382966091923), (-1.9754081700585e-03, 5.3457141292226e-04, 7.5929169129793e-04)) +, _pstate( -554800.0, (-30.658325364620, 2.093818874552, 9.880531138071), ( 6.1010603013347e-05, -3.1326500935382e-03, -9.9346125151067e-04)) +, _pstate( -525600.0, ( 35.737703251673, -12.587706024764, -14.677847247563), ( 1.5802939375649e-03, 2.1347678412429e-03, 1.9074436384343e-04)) +, _pstate( -496400.0, ( 25.466295188546, 41.367478338417, 5.216476873382), (-1.8054401046468e-03, 8.3283083599510e-04, 8.0260156912107e-04)) +, _pstate( -467200.0, (-29.847174904071, 10.636426313081, 12.297904180106), (-6.3257063052907e-04, -2.9969577578221e-03, -7.4476074151596e-04)) +, _pstate( -438000.0, ( 30.774692107687, -18.236637015304, -14.945535879896), ( 2.0113162005465e-03, 1.9353827024189e-03, -2.0937793168297e-06)) +, _pstate( -408800.0, ( 30.243153324028, 38.656267888503, 2.938501750218), (-1.6052508674468e-03, 1.1183495337525e-03, 8.3333973416824e-04)) +, _pstate( -379600.0, (-27.288984772533, 18.643162147874, 14.023633623329), (-1.1856388898191e-03, -2.7170609282181e-03, -4.9015526126399e-04)) +, _pstate( -350400.0, ( 24.519605196774, -23.245756064727, -14.626862367368), ( 2.4322321483154e-03, 1.6062008146048e-03, -2.3369181613312e-04)) +, _pstate( -321200.0, ( 34.505274805875, 35.125338586954, 0.557361475637), (-1.3824391637782e-03, 1.3833397561817e-03, 8.4823598806262e-04)) +, _pstate( -292000.0, (-23.275363915119, 25.818514298769, 15.055381588598), (-1.6062295460975e-03, -2.3395961498533e-03, -2.4377362639479e-04)) +, _pstate( -262800.0, ( 17.050384798092, -27.180376290126, -13.608963321694), ( 2.8175521080578e-03, 1.1358749093955e-03, -4.9548725258825e-04)) +, _pstate( -233600.0, ( 38.093671910285, 30.880588383337, -1.843688067413), (-1.1317697153459e-03, 1.6128814698472e-03, 8.4177586176055e-04)) +, _pstate( -204400.0, (-18.197852930878, 31.932869934309, 15.438294826279), (-1.9117272501813e-03, -1.9146495909842e-03, -1.9657304369835e-05)) +, _pstate( -175200.0, ( 8.528924039997, -29.618422200048, -11.805400994258), ( 3.1034370787005e-03, 5.1393633292430e-04, -7.7293066202546e-04)) +, _pstate( -146000.0, ( 40.946857258640, 25.904973592021, -4.256336240499), (-8.3652705194051e-04, 1.8129497136404e-03, 8.1564228273060e-04)) +, _pstate( -116800.0, (-12.326958895325, 36.881883446292, 15.217158258711), (-2.1166103705038e-03, -1.4814420035990e-03, 1.7401209844705e-04)) +, _pstate( -87600.0, ( -0.633258375909, -30.018759794709, -9.171932874950), ( 3.2016994581737e-03, -2.5279858672148e-04, -1.0411088271861e-03)) +, _pstate( -58400.0, ( 42.936048423883, 20.344685584452, -6.588027007912), (-5.0525450073192e-04, 1.9910074335507e-03, 7.7440196540269e-04)) +, _pstate( -29200.0, ( -5.975910552974, 40.611809958460, 14.470131723673), (-2.2184202156107e-03, -1.0562361130164e-03, 3.3652250216211e-04)) +, _pstate( 0.0, ( -9.875369580774, -27.978926224737, -5.753711824704), ( 3.0287533248818e-03, -1.1276087003636e-03, -1.2651326732361e-03)) +, _pstate( 29200.0, ( 43.958831986165, 14.214147973292, -8.808306227163), (-1.4717608981871e-04, 2.1404187242141e-03, 7.1486567806614e-04)) +, _pstate( 58400.0, ( 0.678136763520, 43.094461639362, 13.243238780721), (-2.2358226110718e-03, -6.3233636090933e-04, 4.7664798895648e-04)) +, _pstate( 87600.0, (-18.282602096834, -23.305039586660, -1.766620508028), ( 2.5567245263557e-03, -1.9902940754171e-03, -1.3943491701082e-03)) +, _pstate( 116800.0, ( 43.873338744526, 7.700705617215, -10.814273666425), ( 2.3174803055677e-04, 2.2402163127924e-03, 6.2988756452032e-04)) +, _pstate( 146000.0, ( 7.392949027906, 44.382678951534, 11.629500214854), (-2.1932815453830e-03, -2.1751799585364e-04, 5.9556516201114e-04)) +, _pstate( 175200.0, (-24.981690229261, -16.204012851426, 2.466457544298), ( 1.8193989149580e-03, -2.6765419531201e-03, -1.3848283502247e-03)) +, _pstate( 204400.0, ( 42.530187039511, 0.845935508021, -12.554907527683), ( 6.5059779150669e-04, 2.2725657282262e-03, 5.1133743202822e-04)) +, _pstate( 233600.0, ( 13.999526486822, 44.462363044894, 9.669418486465), (-2.1079296569252e-03, 1.7533423831993e-04, 6.9128485798076e-04)) +, _pstate( 262800.0, (-29.184024803031, -7.371243995762, 6.493275957928), ( 9.3581363109681e-04, -3.0610357109184e-03, -1.2364201089345e-03)) +, _pstate( 292000.0, ( 39.831980671753, -6.078405766765, -13.909815358656), ( 1.1117769689167e-03, 2.2362097830152e-03, 3.6230548231153e-04)) +, _pstate( 321200.0, ( 20.294955108476, 43.417190420251, 7.450091985932), (-1.9742157451535e-03, 5.3102050468554e-04, 7.5938408813008e-04)) +, _pstate( 350400.0, (-30.669992302160, 2.318743558955, 9.973480913858), ( 4.5605107450676e-05, -3.1308219926928e-03, -9.9066533301924e-04)) +, _pstate( 379600.0, ( 35.626122155983, -12.897647509224, -14.777586508444), ( 1.6015684949743e-03, 2.1171931182284e-03, 1.8002516202204e-04)) +, _pstate( 408800.0, ( 26.133186148561, 41.232139187599, 5.006401326220), (-1.7857704419579e-03, 8.6046232702817e-04, 8.0614690298954e-04)) +, _pstate( 438000.0, (-29.576740229230, 11.863535943587, 12.631323039872), (-7.2292830060955e-04, -2.9587820140709e-03, -7.0824296450300e-04)) +, _pstate( 467200.0, ( 29.910805787391, -19.159019294000, -15.013363865194), ( 2.0871080437997e-03, 1.8848372554514e-03, -3.8528655083926e-05)) +, _pstate( 496400.0, ( 31.375957451819, 38.050372720763, 2.433138343754), (-1.5546055556611e-03, 1.1699815465629e-03, 8.3565439266001e-04)) +, _pstate( 525600.0, (-26.360071336928, 20.662505904952, 14.414696258958), (-1.3142373118349e-03, -2.6236647854842e-03, -4.2542017598193e-04)) +, _pstate( 554800.0, ( 22.599441488648, -24.508879898306, -14.484045731468), ( 2.5454108304806e-03, 1.4917058755191e-03, -3.0243665086079e-04)) +, _pstate( 584000.0, ( 35.877864013014, 33.894226366071, -0.224524636277), (-1.2941245730845e-03, 1.4560427668319e-03, 8.4762160640137e-04)) +, _pstate( 613200.0, (-21.538149762417, 28.204068269761, 15.321973799534), (-1.7312117409010e-03, -2.1939631314577e-03, -1.6316913275180e-04)) +, _pstate( 642400.0, ( 13.971521374415, -28.339941764789, -13.083792871886), ( 2.9334630526035e-03, 9.1860931752944e-04, -5.9939422488627e-04)) +, _pstate( 671600.0, ( 39.526942044143, 28.939897360110, -2.872799527539), (-1.0068481658095e-03, 1.7021132888090e-03, 8.3578230511981e-04)) +, _pstate( 700800.0, (-15.576200701394, 34.399412961275, 15.466033737854), (-2.0098814612884e-03, -1.7191109825989e-03, 7.0414782780416e-05)) +, _pstate( 730000.0, ( 4.243252837090, -30.118201690825, -10.707441231349), ( 3.1725847067411e-03, 1.6098461202270e-04, -9.0672150593868e-04)) ] -class _TerseVector: - def __init__(self, x: float, y: float, z: float) -> None: - self.x = x - self.y = y - self.z = z - - def clone(self) -> "_TerseVector": - '''Create a copy of this vector.''' - return _TerseVector(self.x, self.y, self.z) - - @staticmethod - def zero() -> "_TerseVector": - '''Return a zero vector.''' - return _TerseVector(0.0, 0.0, 0.0) - - def ToAstroVector(self, time: Time) -> Vector: - '''Convert _TerseVector object to Vector object.''' - return Vector(self.x, self.y, self.z, time) - - def quadrature(self) -> float: - '''Return magnitude squared of this vector.''' - return self.x**2 + self.y**2 + self.z**2 - - def mean(self, other: "_TerseVector") -> "_TerseVector": - '''Return the average of this vector and another vector.''' - return _TerseVector((self.x + other.x)/2.0, (self.y + other.y)/2.0, (self.z + other.z)/2.0) - - def __add__(self, other: "_TerseVector") -> "_TerseVector": - return _TerseVector(self.x + other.x, self.y + other.y, self.z + other.z) - - def __sub__(self, other: "_TerseVector") -> "_TerseVector": - return _TerseVector(self.x - other.x, self.y - other.y, self.z - other.z) - - def __mul__(self, scalar: float) -> "_TerseVector": - return _TerseVector(scalar * self.x, scalar * self.y, scalar * self.z) - - def __rmul__(self, scalar: float) -> "_TerseVector": - return _TerseVector(scalar * self.x, scalar * self.y, scalar * self.z) - - def __truediv__(self, scalar: float) -> "_TerseVector": - return _TerseVector(self.x / scalar, self.y / scalar, self.z / scalar) - - -def _BodyStateFromTable(entry: Tuple[float, Tuple[float, float, float], Tuple[float, float, float]]) -> _body_state_t: - ( tt, (rx, ry, rz), (vx, vy, vz) ) = entry - return _body_state_t(tt, _TerseVector(rx, ry, rz), _TerseVector(vx, vy, vz)) - - def _AdjustBarycenterPosVel(ssb: _body_state_t, tt: float, body: Body, planet_gm: float) -> _body_state_t: shift = planet_gm / (planet_gm + _SUN_GM) planet = _CalcVsopPosVel(_vsop[body.value], tt) @@ -3429,22 +3431,21 @@ def _ClampIndex(frac: float, nsteps: int) -> int: return index -def _GravFromState(entry: Tuple[float, Tuple[float, float, float], Tuple[float, float, float]]) -> _grav_sim_t: - state = _BodyStateFromTable(entry) - bary = _major_bodies_t(state.tt) - r = state.r + bary.Sun.r - v = state.v + bary.Sun.v +def _GravFromState(entry: _pstate) -> _grav_sim_t: + bary = _major_bodies_t(entry.tt) + r = entry.pos + bary.Sun.r + v = entry.vel + bary.Sun.v a = bary.Acceleration(r) - grav = _body_grav_calc_t(state.tt, r, v, a) + grav = _body_grav_calc_t(entry.tt, r, v, a) return _grav_sim_t(bary, grav) def _GetSegment(cache: List[Optional[List[_body_grav_calc_t]]], tt: float) -> Optional[List[_body_grav_calc_t]]: - if (tt < _PlutoStateTable[0][0]) or (tt > _PlutoStateTable[_PLUTO_NUM_STATES-1][0]): + if (tt < _PlutoStateTable[0].tt) or (tt > _PlutoStateTable[_PLUTO_NUM_STATES-1].tt): # Don't bother calculating a segment. Let the caller crawl backward/forward to this time. return None - seg_index = _ClampIndex((tt - _PlutoStateTable[0][0]) / _PLUTO_TIME_STEP, _PLUTO_NUM_STATES-1) + seg_index = _ClampIndex((tt - _PlutoStateTable[0].tt) / _PLUTO_TIME_STEP, _PLUTO_NUM_STATES-1) if cache[seg_index] is None: seg = cache[seg_index] = [ _GravFromState(_PlutoStateTable[seg_index]).grav ] @@ -3481,7 +3482,7 @@ def _GetSegment(cache: List[Optional[List[_body_grav_calc_t]]], tt: float) -> Op return cache[seg_index] -def _CalcPlutoOneWay(entry: Tuple[float, Tuple[float, float, float], Tuple[float, float, float]], target_tt: float, dt: float) -> _grav_sim_t: +def _CalcPlutoOneWay(entry: _pstate, target_tt: float, dt: float) -> _grav_sim_t: sim = _GravFromState(entry) n = math.ceil((target_tt - sim.grav.tt) / dt) for i in range(n): @@ -3496,7 +3497,7 @@ def _CalcPluto(time: Time, helio: bool) -> StateVector: # The target time is outside the year range 0000..4000. # Calculate it by crawling backward from 0000 or forward from 4000. # FIXFIXFIX - This is super slow. Could optimize this with extra caching if needed. - if time.tt < _PlutoStateTable[0][0]: + if time.tt < _PlutoStateTable[0].tt: sim = _CalcPlutoOneWay(_PlutoStateTable[0], time.tt, -_PLUTO_DT) else: sim = _CalcPlutoOneWay(_PlutoStateTable[_PLUTO_NUM_STATES-1], time.tt, +_PLUTO_DT)