This question will likely get closed, but I think I am in a position to help you out, so I'll post some example code to get you going in the right direction.
Give the following a shot if you can install the SpiceyPy module (documentation here). The first function downloads all of the necessary CSPICE kernels from NASA's servers. If you end up using it, there is a special (under 'ea_latest') kernel for the latest high precision earth data, its updated somewhat regularly, so make sure to download that often (I set up an auto-download script to do this). Where ever you try to run this file from, create a sub-folder called "kernels", and these files will be downloaded automatically if you don't have them when you run the script. Thats what the first, long function does. The rest of the functions are pretty short, so you can see how easy it is to use this. I hope this helps.
import spiceypy.wrapper as spw
import numpy as np
import datetime
import os
import urllib
def get_ephem_kernels():
# The ephemeris data for the moon and the sun is now downloaded from JPL in the form of
# CSPICE kernels. We use CSPICE to load all of this data into our environment so that
# it can be called into our program whenever it is needed.
print 'Retrieving Ephemeris Data'
# This is the current site for the leapsecond data. This may change in the future as more leapseconds are added.
# check the path up to '.../lsk/' to find the most recent kernel
ls_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0011.tls'
ls_file = './kernels/ls.tls'
# This is the ephemeris data, using DE421. I use this because it goes continuously up to the year 2050.
# There are other versions that we could consider using at this point, but they will only provide very minor
# corrections to the sun and moon positions.
de_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/a_old_versions/de421.bsp'
de_file = './kernels/de.bsp'
# These three files introduce physical constants, and additional parameters related
# to the lunar reference frame. At the moment, they are not used, but may be needed in order
# to match the results of the VIIRS code.
pc_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/pck00010.tpc'
pc_file = './kernels/pck.tpc'
mn_pa_de_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/moon_pa_de421_1900-2050.bpc'
mn_pa_de_file = './kernels/mn_pa_de.bpc'
mn_pa_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/fk/satellites/moon_assoc_pa.tf'
mn_pa_file = './kernels/mn_pa.tf'
ea_pa_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/fk/planets/'
ea_pa_file = './kernels/ea_pa.tf'
ea_pck_in = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/earth_070425_370426_predict.bpc'
ea_pck_file = './kernels/ea_pck_predict.bpc'
ea_pck_hist = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/earth_720101_070426.bpc'
ea_pck_hist_file = './kernels/ea_pck_hist.bpc'
ea_latest = 'http://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/earth_latest_high_prec.bpc'
ea_latest_file = './kernels/ea_latest.bpc'
# The first thing that we figure out is whether or not the files exist on our machine.
# If they do, it will pass these statements. If not, it will download the necessary files.
if os.path.isfile(ls_file) == False:
urllib.urlretrieve(ls_in,ls_file)
if os.path.isfile(de_file) == False:
urllib.urlretrieve(de_in,de_file)
if os.path.isfile(pc_file) == False:
urllib.urlretrieve(pc_in,pc_file)
if os.path.isfile(mn_pa_de_file) == False:
urllib.urlretrieve(mn_pa_de_in,mn_pa_de_file)
if os.path.isfile(mn_pa_file) == False:
urllib.urlretrieve(mn_pa_in,mn_pa_file)
if os.path.isfile(ea_pa_file) == False:
urllib.urlretrieve(ea_pa_in,ea_pa_file)
if os.path.isfile(ea_pck_file) == False:
urllib.urlretrieve(ea_pck_in,ea_pck_file)
if os.path.isfile(ea_latest_file) == False:
urllib.urlretrieve(ea_latest,ea_latest_file)
if os.path.isfile(ea_pck_hist_file) == False:
urllib.urlretrieve(ea_pck_hist,ea_pck_hist_file)
# This uses CSPICE to load the data so that it can be called later. The data is loaded throughout
# the whole environment so it can be called within any function or class object.
spw.furnsh(ls_file)
spw.furnsh(de_file)
spw.furnsh(pc_file)
spw.furnsh(mn_pa_de_file)
spw.furnsh(mn_pa_file)
spw.furnsh(ea_pa_file)
spw.furnsh(ea_pck_file)
spw.furnsh(ea_latest_file)
spw.furnsh(ea_pck_hist_file)
def eq_to_cart(ra,dec,deg='yes'):
'''
This transforms equatorial coordinates to cartesian coordinates in whatever reference frame you are in.
This just gives a unit vector however, and any range will need to be supplied by you
'''
if deg == 'yes':
ra = np.deg2rad(ra)
dec = np.deg2rad(dec)
pv_out = spw.radrec(1.0,ra,dec)
return pv_out
def cart_to_eq(pv,deg='yes'):
'''
This transforms from cartesian to equatorial in your current reference frame. The range is in your
current units
'''
eq_data = spw.recrad(pv)
rang = eq_data[0]
ra = eq_data[1]
dec = eq_data[2]
if deg == 'yes':
ra = np.rad2deg(ra)
dec = np.rad2deg(dec)
return rang,ra,dec
def change_ref_frame(pv,time_ref,current='J2000',new='ITRF93'):
'''
This function will change the reference frame between two chosen frames using a rotation matrix
at a given time
'''
new_time = spw.str2et(str(time_ref))
pv = np.transpose(np.matrix(pv))
rot_mat = np.matrix(spw.pxform(current,new,new_time))
pv_rot = (rot_mat*pv).tolist()
for i in range(len(pv_rot)):
pv_rot[i] = pv_rot[i][0]
return pv_rot
if __name__ == '__main__':
get_ephem_kernels()
ra = 50.12
dec = -30.34
time_ref = datetime.datetime.strptime('2015-12-15 12:00:00','%Y-%m-%d %H:%M:%S')
pv = eq_to_cart(ra,dec)
rang,ra,dec = cart_to_eq(pv)
new_pv = change_ref_frame(pv,time_ref,current='J2000',new='ITRF93')
print new_pv