Skip to content

Commit

Permalink
Removed skyfield dependency and replaced time eval with spiceypy
Browse files Browse the repository at this point in the history
  • Loading branch information
@itsmoosh
itsmoosh committed Nov 13, 2022
1 parent 7b870d1 commit ff1ef35
Show file tree
Hide file tree
Showing 2 changed files with 159 additions and 10 deletions.
17 changes: 7 additions & 10 deletions MoonMag/eval_induced_field.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,10 +14,10 @@
import matplotlib.pyplot as plt
import matplotlib as mpl

from skyfield import api as skyapi
import spiceypy as spice
from datetime import datetime as dtime
from datetime import timezone

from MoonMag import _excitation
from MoonMag.config import *
import MoonMag.symmetry_funcs as sym
import MoonMag.asymmetry_funcs as asym
Expand Down Expand Up @@ -86,16 +86,13 @@ def run_calcs(bname, comp, recalc, plot_field, plot_asym, synodic_only=True,
n_max_main = nprm_max_main + p_max_main

# Set time to evaluate magnetic fields
tscale = skyapi.load.timescale()
J2000 = tscale.J(2000)
spiceTLS = 'naif0012.tls'
spice.furnsh(os.path.join(_excitation, spiceTLS))

tzone = timezone.utc # Timezone code readable by skyfield module to apply to eval_datetime string in config file -- see https://rhodesmill.org/skyfield/time.html#utc-and-your-timezone
t = dtime.fromisoformat(eval_datetime).replace(tzinfo=tzone) # Evaluate at the specified eval_datetime in config file, in the above time zone
#t = tscale.now().utc_datetime() # Evaluate as the bodies are right NOW
#t = J2000.utc_datetime() # Evaluate AT J2000
#tsec = spice.str2et(dtime.now().isoformat()) # Evaluate as the bodies are right NOW
#tsec = 0 # Evaluate AT J2000
tsec = spice.str2et(eval_datetime)

tdiff_jd = ( tscale.from_datetime(t) - J2000 )
tsec = tdiff_jd * 86400

# Linear arrays of values to loop over for nprm,mprm,p,q,n,m
nprmvals = [ nprm for nprm in range(1,nprm_max_main+1) for _ in range(-nprm,nprm+1) ]
Expand Down
152 changes: 152 additions & 0 deletions MoonMag/excitation/naif0012.tls
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,152 @@
KPL/LSK


LEAPSECONDS KERNEL FILE
===========================================================================

Modifications:
--------------

2016, Jul. 14 NJB Modified file to account for the leapsecond that
will occur on December 31, 2016.

2015, Jan. 5 NJB Modified file to account for the leapsecond that
will occur on June 30, 2015.

2012, Jan. 5 NJB Modified file to account for the leapsecond that
will occur on June 30, 2012.

2008, Jul. 7 NJB Modified file to account for the leapsecond that
will occur on December 31, 2008.

2005, Aug. 3 NJB Modified file to account for the leapsecond that
will occur on December 31, 2005.

1998, Jul 17 WLT Modified file to account for the leapsecond that
will occur on December 31, 1998.

1997, Feb 22 WLT Modified file to account for the leapsecond that
will occur on June 30, 1997.

1995, Dec 14 KSZ Corrected date of last leapsecond from 1-1-95
to 1-1-96.

1995, Oct 25 WLT Modified file to account for the leapsecond that
will occur on Dec 31, 1995.

1994, Jun 16 WLT Modified file to account for the leapsecond on
June 30, 1994.

1993, Feb. 22 CHA Modified file to account for the leapsecond on
June 30, 1993.

1992, Mar. 6 HAN Modified file to account for the leapsecond on
June 30, 1992.

1990, Oct. 8 HAN Modified file to account for the leapsecond on
Dec. 31, 1990.


Explanation:
------------

The contents of this file are used by the routine DELTET to compute the
time difference

[1] DELTA_ET = ET - UTC

the increment to be applied to UTC to give ET.

The difference between UTC and TAI,

[2] DELTA_AT = TAI - UTC

is always an integral number of seconds. The value of DELTA_AT was 10
seconds in January 1972, and increases by one each time a leap second
is declared. Combining [1] and [2] gives

[3] DELTA_ET = ET - (TAI - DELTA_AT)

= (ET - TAI) + DELTA_AT

The difference (ET - TAI) is periodic, and is given by

[4] ET - TAI = DELTA_T_A + K sin E

where DELTA_T_A and K are constant, and E is the eccentric anomaly of the
heliocentric orbit of the Earth-Moon barycenter. Equation [4], which ignores
small-period fluctuations, is accurate to about 0.000030 seconds.

The eccentric anomaly E is given by

[5] E = M + EB sin M

where M is the mean anomaly, which in turn is given by

[6] M = M + M t
0 1

where t is the number of ephemeris seconds past J2000.

Thus, in order to compute DELTA_ET, the following items are necessary.

DELTA_TA
K
EB
M0
M1
DELTA_AT after each leap second.

The numbers, and the formulation, are taken from the following sources.

1) Moyer, T.D., Transformation from Proper Time on Earth to
Coordinate Time in Solar System Barycentric Space-Time Frame
of Reference, Parts 1 and 2, Celestial Mechanics 23 (1981),
33-56 and 57-68.

2) Moyer, T.D., Effects of Conversion to the J2000 Astronomical
Reference System on Algorithms for Computing Time Differences
and Clock Rates, JPL IOM 314.5--942, 1 October 1985.

The variable names used above are consistent with those used in the
Astronomical Almanac.

\begindata

DELTET/DELTA_T_A = 32.184
DELTET/K = 1.657D-3
DELTET/EB = 1.671D-2
DELTET/M = ( 6.239996D0 1.99096871D-7 )

DELTET/DELTA_AT = ( 10, @1972-JAN-1
11, @1972-JUL-1
12, @1973-JAN-1
13, @1974-JAN-1
14, @1975-JAN-1
15, @1976-JAN-1
16, @1977-JAN-1
17, @1978-JAN-1
18, @1979-JAN-1
19, @1980-JAN-1
20, @1981-JUL-1
21, @1982-JUL-1
22, @1983-JUL-1
23, @1985-JUL-1
24, @1988-JAN-1
25, @1990-JAN-1
26, @1991-JAN-1
27, @1992-JUL-1
28, @1993-JUL-1
29, @1994-JUL-1
30, @1996-JAN-1
31, @1997-JUL-1
32, @1999-JAN-1
33, @2006-JAN-1
34, @2009-JAN-1
35, @2012-JUL-1
36, @2015-JUL-1
37, @2017-JAN-1 )

\begintext


0 comments on commit ff1ef35

Please sign in to comment.