from astropy.time import Time
from astropy.coordinates import solar_system_ephemeris, EarthLocation, get_body, get_body_barycentric, AltAz, SkyCoord, TETE
from astropy.coordinates.solar_system import _apparent_position_in_true_coordinates
from astropy.coordinates.builtin_frames.utils import get_polar_motion
import astropy.units as u
import numpy as np
from datetime import datetime, timedelta
import cProfile

APR_6_2017 = 1491436800.0
RAD2ASEC = (180.0 * 3600.0 / np.pi)

# Site Location
class_lat = -22.959748*u.deg
class_lon = -67.787260*u.deg
class_elev= 5186*u.m
class_loc =  EarthLocation.from_geodetic(class_lon,class_lat,class_elev)

# For comparison with JPL Horizons as this is what they use
solar_system_ephemeris.set('de430')

class body_path(object):
    
    def __init__(self, body_name, ctimes, loc=class_loc):
        """
        This is a class to compute the locations of a celestial body
        as a function of [ctime]. Specify [body_name], such as "sun", 
        "moon", "jupiter", etc. You can specify an EarthLocation (astropy)
        or leave it at it's default of the CLASS location. 
        
        """
        self.body_name=body_name
        self.ctimes = ctimes
        self.loc = loc
        self.ast_times = Time(ctimes, format='unix', location=self.loc)
        self.lst = self.ast_times.sidereal_time('apparent')
        self.dut1 = self.ast_times.delta_ut1_utc
        self.xp, self.yp = get_polar_motion(self.ast_times)
        self.radec = None
        self.altaz = None
        self.apparent = None
        self.astrometric = None
    
    def compute_radec(self):
        if self.radec is not None:
            return
        self.radec = get_body(self.body_name, self.ast_times)
        self.apparent = self.radec.transform_to(TETE(obstime=self.ast_times,location=self.loc))
        earth_loc = get_body_barycentric('earth', self.ast_times)
        obsgeoloc, obsgeovel = self.loc.get_gcrs_posvel(self.ast_times)
        earth_loc += obsgeoloc
        self.astrometric = SkyCoord(self.radec.transform_to('icrs').cartesian - earth_loc)
    
    def compute_altaz(self):
        if self.altaz is not None:
            return
        if self.radec is None:
            self.compute_radec()
        self.altaz = self.radec.transform_to(AltAz(obstime=self.ast_times,location=self.loc))
    
start = APR_6_2017
nsamples = 1441
delta = 60.0
times = []
ctimes = []
for i in range(nsamples):
    ctime = start + i * delta
    ctimes.append(ctime)
    times.append(datetime.utcfromtimestamp(ctime))

obj = body_path('moon', ctimes)
obj.compute_radec()
obj.compute_altaz()


print("UTC,LST,ICRS_RA,ICRS_DEC,RA,DEC,AZ,EL,DIST,DUT1,XP,YP")
for i in range(nsamples):
    times[i] = times[i].strftime("%Y-%m-%d %H:%M:%S.%f")
    print ('%s,%.12lf,%.12lf,%.12lf,%.12lf,%.12lf,%.12lf,%.12lf,%.6lf,%.7lf,%.7lf,%.7lf' % (times[i], obj.lst[i].value, obj.astrometric.ra[i].value , obj.astrometric.dec[i].value, 
                                                      obj.apparent.ra[i].value, obj.apparent.dec[i].value, 
                                                      obj.altaz.az[i].value, obj.altaz.alt[i].value, obj.astrometric.distance[i].to(u.meter).value,
                                                      obj.dut1[i], obj.xp[i] * RAD2ASEC, obj.yp[i] * RAD2ASEC ))