notebook.community

Edit and run



In [1]:

    
from nustar_pysolar import planning
import astropy.units as u
from astropy.time import Time
from astropy.coordinates import get_body, solar_system_ephemeris, get_body_barycentric, SkyCoord









    



WARNING: AstropyDeprecationWarning: astropy.utils.compat.odict.OrderedDict is now deprecated - import OrderedDict from the collections module instead [astropy.utils.compat.odict]



In [2]:

    
fname = planning.download_occultation_times(outdir='../data/')
print(fname)









    



../data/NUSTAR.2017_129.SHADOW_ANALYSIS.txt



In [3]:

    
tstart = '2017-05-19T01:43:00'
tend = '2017-05-19T11:43:00'
orbits = planning.sunlight_periods(fname, tstart, tend)
solar_system_ephemeris.set('jpl')









    Out[3]:





<ScienceState solar_system_ephemeris: 'jpl'>



In [4]:

    
from skyfield.api import Loader
load = Loader('../data')
ts = load.timescale()
planets = load('jup310.bsp')









    



[#################################] 100% deltat.data
[#################################] 100% deltat.preds
[#################################] 100% Leap_Second.dat
[#################################] 100% jup310.bsp

Use the astropy interface to get the location of Jupiter as the time that you want to use.



In [5]:

    
dt = 0.

# Using JPL Horizons web interface at 2017-05-19T01:34:40
horizon_ephem = SkyCoord(*[193.1535, -4.01689]*u.deg)


for orbit in orbits:
    tstart = orbit[0]
    tend = orbit[1]
    print()
#    print('Orbit duration: ', tstart.isoformat(), tend.isoformat())
    on_time = (tend - tstart).total_seconds()
    
    point_time = tstart + 0.5*(tend - tstart)
    print('Time used for ephemeris: ', point_time.isoformat())
    
    astro_time = Time(point_time)
    solar_system_ephemeris.set('jpl')


    jupiter = get_body('Jupiter', astro_time)
    
    jplephem = SkyCoord(jupiter.ra.deg*u.deg, jupiter.dec.deg*u.deg)
    
    # Switch to the built in ephemris
    solar_system_ephemeris.set('builtin')
    jupiter = get_body('Jupiter', astro_time)
    
    builtin_ephem = SkyCoord(jupiter.ra.deg*u.deg, jupiter.dec.deg*u.deg)
    
    t = ts.from_astropy(astro_time)
    jupiter, earth = planets['jupiter'], planets['earth']
    astrometric = earth.at(t).observe(jupiter)
    ra, dec, distance = astrometric.radec()
    radeg = ra.to(u.deg)
    decdeg = dec.to(u.deg)
    skyfield_ephem = SkyCoord(radeg, decdeg)
    
    
    print()
    print('Horizons offset to jplephem: ', horizon_ephem.separation(jplephem))
    print()
    print('Horizons offset to "built in" ephemeris: ', horizon_ephem.separation(builtin_ephem))
    print()
    print('Horizons offset to Skyfield ephemeris: ', horizon_ephem.separation(skyfield_ephem))

    print()
    break









    



Time used for ephemeris:  2017-05-19T01:45:00

Horizons offset to jplephem:  0d00m13.0448s

Horizons offset to "built in" ephemeris:  0d00m15.7584s

Horizons offset to Skyfield ephemeris:  0d00m01.6335s

Conclusion: Use skyfield if you want to reproduce the JPL ephemerides

Use the jup310.bsp file for Jupiter. Need to confirm which of the avaiable .bsp files are approriate for inner solar system objects as well as the Sun/Moon



In [7]:

    
dt = 0.

for orbit in orbits:
    tstart = orbit[0]
    tend = orbit[1]
    print()
    on_time = (tend - tstart).total_seconds()
    
    point_time = tstart + 0.5*(tend - tstart)
    print('Time used for ephemeris: ', point_time.isoformat())
    
    astro_time = Time(point_time)
    solar_system_ephemeris.set('jpl')


    jupiter = get_body('Jupiter', astro_time)
    
    jplephem = SkyCoord(jupiter.ra.deg*u.deg, jupiter.dec.deg*u.deg)
    
    # Switch to the built in ephemris
    solar_system_ephemeris.set('builtin')
    jupiter = get_body('Jupiter', astro_time)
    
    builtin_ephem = SkyCoord(jupiter.ra.deg*u.deg, jupiter.dec.deg*u.deg)
    
    t = ts.from_astropy(astro_time)
    jupiter, earth = planets['jupiter'], planets['earth']
    astrometric = earth.at(t).observe(jupiter)
    ra, dec, distance = astrometric.radec()
    radeg = ra.to(u.deg)
    decdeg = dec.to(u.deg)
    skyfield_ephem = SkyCoord(radeg, decdeg)
    
    
    print()
    print('Skyfield offset to jplephem: ', skyfield_ephem.separation(jplephem))
    print()
    print('Skyfield offset to "built in" ephemeris: ', skyfield_ephem.separation(builtin_ephem))

    print()









    



Time used for ephemeris:  2017-05-19T01:45:00

Skyfield offset to jplephem:  0d00m14.6773s

Skyfield offset to "built in" ephemeris:  0d00m17.3375s


Time used for ephemeris:  2017-05-19T03:21:45

Skyfield offset to jplephem:  0d00m14.6587s

Skyfield offset to "built in" ephemeris:  0d00m17.315s


Time used for ephemeris:  2017-05-19T04:58:25

Skyfield offset to jplephem:  0d00m14.64s

Skyfield offset to "built in" ephemeris:  0d00m17.2925s


Time used for ephemeris:  2017-05-19T06:35:05

Skyfield offset to jplephem:  0d00m14.6213s

Skyfield offset to "built in" ephemeris:  0d00m17.2699s


Time used for ephemeris:  2017-05-19T08:11:50

Skyfield offset to jplephem:  0d00m14.6026s

Skyfield offset to "built in" ephemeris:  0d00m17.2474s


Time used for ephemeris:  2017-05-19T09:48:30

Skyfield offset to jplephem:  0d00m14.584s

Skyfield offset to "built in" ephemeris:  0d00m17.225s