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

"""
Explanation: Download the list of occultation periods from the MOC at Berkeley.
Note that the occultation periods typically only are stored at Berkeley for the future and not for the past. So this is only really useful for observation planning.
End of explanation
"""


tlefile = io.download_tle(outdir='../data')
print(tlefile)
times, line1, line2 = io.read_tle_file(tlefile)

"""
Explanation: Download the NuSTAR TLE archive.
This contains every two-line element (TLE) that we've received for the whole mission. We'll expand on how to use this later.
The times, line1, and line2 elements are now the TLE elements for each epoch.
End of explanation
"""


tstart = '2017-07-18T12:00:00'
tend = '2017-07-18T20:00:00'
orbits = planning.sunlight_periods(fname, tstart, tend)

"""
Explanation: Here is where we define the observing window that we want to use.
Note that tstart and tend must be in the future otherwise you won't find any occultation times and sunlight_periods will return an error.
End of explanation
"""


pa = planning.get_nustar_roll(tstart, 0)
print("NuSTAR Roll angle for Det0 in NE quadrant: {}".format(pa))

"""
Explanation: We want to know how to orient NuSTAR for the Sun.
We can more or less pick any angle that we want. But this angle has to be specified a little in advance so that the NuSTAR SOC can plan the "slew in" maneuvers. Below puts DET0 in the top left corner (north-east with respect to RA/Dec coordinates).
This is what you tell the SOC you want the "Sky PA angle" to be.
End of explanation
"""


offset = [-190., -47.]*u.arcsec

"""
Explanation: Set up the offset you want to use here:
The first element is the direction +WEST of the center of the Sun, the second is the offset +NORTH of the center of the Sun.
If you want multiple pointing locations you can either specify an array of offsets or do this "by hand" below.
End of explanation
"""


from astropy.coordinates import SkyCoord
for ind, orbit in enumerate(orbits):
    midTime = (0.5*(orbit[1] - orbit[0]) + orbit[0])
    sky_pos = planning.get_sky_position(midTime, offset)
    print("Orbit: {}".format(ind))
    print("Orbit start: {} Orbit end: {}".format(orbit[0].isoformat(), orbit[1].isoformat()))
    print('Aim time: {} RA (deg): {} Dec (deg): {}'.format(midTime.isoformat(), sky_pos[0], sky_pos[1]))
    skyfield_pos = planning.get_skyfield_position(midTime, offset, load_path='../data')
    print('SkyField Aim time: {} RA (deg): {} Dec (deg): {}'.format(midTime.isoformat(), skyfield_pos[0], skyfield_pos[1]))
    skyfield_ephem = SkyCoord(skyfield_pos[0], skyfield_pos[1])
    sunpy_ephem = SkyCoord(sky_pos[0], sky_pos[1])
    print("")
    print("Offset between SkyField and Astropy: {} arcsec".format(skyfield_ephem.separation(sunpy_ephem).arcsec))
    print("")
         

"""
Explanation: Loop over each orbit and see what the difference between the two methods is
Note that you may want to update the pointing for solar rotation. That's up to the user to decide and is not done here.
Looks like a fixed shift...probably some time-ephemeris issue.
End of explanation
"""


from astropy.coordinates import SkyCoord
from datetime import timedelta


for ind, orbit in enumerate(orbits):
    midTime = orbit[0]
    while(midTime < orbit[1]):
        
        sky_pos = planning.get_sky_position(midTime, offset)

        skyfield_pos = planning.get_skyfield_position(midTime, offset, load_path='../data', parallax_correction=True)
        skyfield_geo = planning.get_skyfield_position(midTime, offset, load_path='../data', parallax_correction=False)

        
        skyfield_ephem = SkyCoord(skyfield_pos[0], skyfield_pos[1])
        skyfield_geo_ephem = SkyCoord(skyfield_geo[0], skyfield_geo[1])
#        sunpy_ephem = SkyCoord(sky_pos[0], sky_pos[1])
        print('Offset between parallax-corrected positions and geoenctric is {} arcsec'.format(
            skyfield_geo_ephem.separation(skyfield_ephem).arcsec)
             )
        dra, ddec = skyfield_geo_ephem.spherical_offsets_to(skyfield_ephem)
        print('{0} delta-RA, {1} delta-Dec'.format(dra.to(u.arcsec), ddec.to(u.arcsec)))
        print('')
        midTime += timedelta(seconds=100)
        
    break

"""
Explanation: Okay, now check to see what the parallax does in each orbit.
Compare Astropy/Sunpy to what you get when you correct for the orbital parallax. Every step below is 100 seconds.
End of explanation
"""


