Disclaimer: I updated my previous answer as I found a more accurate formula.
The rotation angle $\theta$ is the angle between the vernal point and the meridian at Greenwich. This also corresponds with the Sidereal Time at Greenwich, converted to radians.
![enter image description here](https://i.stack.imgur.com/xIhrH.png)
For a given UTC time and a given date, the corresponding Greenwich (Mean) Sidereal Time (in hours and decimal parts of an hour) can be calculated as
$$
\text{GMST} = 18.697\,374\,558 + 24.065\,709\,824\,419\,d,
$$
where $d$ is the interval, in UT1 days including any fraction of a day, since 2000 January 1, at 12h UT (source: US Naval Observatory). Of course, this
must be reduced to the range 0h to 24h. The rotation angle is then simply
$$
\theta = \frac{2\pi}{24}\text{GMST}.
$$
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
Best Answer
Some key reading if you want to understand this stuff is chapters two to five of the IERS Technical Note 36, the IERS Conventions (2010).
It's not just the J2000/FK5 frame (aka the EME2000 frame) that is associated with some epoch date. Every Earth-centered inertial frame has some epoch date. There are two fundamental reasons why this must be the case:
Note that the J2000/FK5 frame is now twice passé. The current best estimate of what constitutes an inertial frame is the International Celestial Reference Frame 2 (ICRF2). It's predecessor, the ICRF, represented a vast improvement over the J2000/FK5 frame. The ICRF2 is even better than the ICRF. The ICRF was supposed to be co-aligned with the J2000/FK5 frame on J2000.0 (12 noon Terrestrial Time on 2000 January 1). It turned out that this was not the case; there's a slight bias between the frames at the epoch. The J2000/FK5 frame also turns out to be rotating a tiny bit, about 3 milliarcseconds/year. Unless you are doing milliarcsecond astronomy, you can ignore that bias and rotation. For most applications, J2000/FK5=ICRF=ICRF2.
The first item, that the Earth's rotation axis is not constant, is important. The Earth's rotation axis precesses with a period of about 26,000 years. Accounting for the change in the precession between the epoch time and the time of interest (e.g., today) yields a transformation from the epoch frame (e.g. J2000) to the mean of date frame. In addition to this long-term precession, the Earth's axis also displays some shorter term variations in where it points. These short-term variations (from ~5.5 days to 18.6 years) are collectively called nutation. Accounting for the Earth's nutation on top of the precession yields the transformation to the true of date frame. Finally, the Earth rotates at about one revolution per sidereal day about this precessed and nutated axis. Applying this rotation of the true of date frame yields the Earth-centered, Earth-fixed frame. A somewhat widely-used name for the result of this process is the Earth RNP (rotation, nutation, and precession) matrix.
Well, almost. Precession and nutation are semi-analytical models, as is the concept of one revolution per sidereal day. There are some things those models just can't capture.
That one revolution per sidereal day is incorrect for two reasons. One is that the Earth's rotation rate is very gradually slowing down. Another is that when you look at the rotation rate very closely, the Earth sometimes rotates faster than nominal, other times slower. There are two key parameters that describe this, dUT1=UT1−UTC and ΔT=TT-UT1. If you care about this detail I suggest you use the latter as it's continuous. dUT1 has discontinuities at the leap seconds. This is a correction that you add when you compute the rotation part of the RNP matrix.
There are some things that the semi-analytical precession and nutation model just don't cover (yet). The Chandler wobble, for example. These are collectively called "polar motion", and can only be observed (and predicted to some extent). Polar motion needs to be applied after computing the RNP matrix. The full result is sometimes called the PRNP (polar motion, rotation, nutation, and precession) matrix. These fine scale variations in the Earth's orientation, along with dUT1 and ΔT, are called the "Earth Orientation Parameters". These are published on a regular basis as "IERS Bulletins A and B". I'll say more about this below.
You do not want to code this on your own. You can get code to do these calculations from a number of places. The best sites are:
International Astronomical Union (IAU) Standards of Fundamental Astronomy (SOFA)
The SOFA code is the "official" version of all of the concepts described above. You can get both Fortran and "Ctran" (Fortran converted to ugly C) versions of the SOFA code from the SOFA website. Also be sure to check out the cookbooks, particular the cookbook on "SOFA Tools for Earth Attitude".
Naval Observatory Vector Astrometry Software (NOVAS)
The US Naval Observatory is responsible for IERS Bulletins A and B. They have their own software, distinct from the SOFA code. The NOVAS software is available in Fortran, C, and Python. The difference in terms of results is negligible, in the microarcseconds. That's the kind of error one would expect from using double precision and performing the same computations a bit differently.
There are a number of others (e.g., JPL Spice, GSFC GMAT, Orekit), but I suggest you go to the source, and that would be either the IAU or the US Naval Observatory.
I mentioned IERS Bulletins A and B a couple of times above. The International Earth Orientation and Reference Systems Service (IERS) is the worldwide organization responsible for defining things such as the ICRF and for determining how the Earth is oriented. (Yes, the acronym doesn't match. It used to before they changed the name but not the acronym.) As far as those technical bulletins are concerned, they just contain numbers (and a tiny bit of text). These numbers are time-tabulated values for the Earth Orientation Parameters. These bulletins are updated monthly.
A couple of final points:
I put the link to IERS Technical Note 36 at the top of this answer. Read it.
Be very careful of time. There are a number of time scales involved in this modeling. A few of them that you will run into are:
TAI - International Atomic Time. Time according to an earthbound physicist who uses an atomic clock at sea level.
TT - Terrestrial Time. Time according to an earthbound astronomer. Physicists and astronomers disagree by 32.184 seconds.
UT1 - Universal Time. Conceptually, what a sundial says the time is, but smoothed to eliminate things such as the equation of time.
UTC - Coordinated Universal Time. That's what the clock on your computer shows if you are using network time protocol. Next time we have a leap second (probably the end of next year), you'll be able to see a minute with 61 seconds in it. UTC ticks at the same rate as TAI and TT but occasionally has leap seconds so as to stay within a second of UT1.
TCB - Barycentric Coordinate Time. A general relativistic timescale that on average ticks faster than clocks on the surface of the Earth.
TDB - Barycentric Dynamical Time. A general relativistic timescale that on average ticks at the same rate as clocks on the surface of the Earth.
GAST - Greenwich Apparent Sidereal Time. If you run into this time scale you are looking at an out-of-date concept for calculating Earth rotation. Use the newer Earth Rotation Angle concept, which relies on UT1.
Update
I didn't answer the title of the question,
This is the easy part. The only tricky aspect is that latitude is almost always geodetic latitude rather than geocentric latitude. I'll assume that latitude $\phi$, longitude $\lambda$, and altitude $h$ are in the WGS84 reference system. See [Department of Defense (2000), "World Geodetic System 1984: Its Definition and Relationships with Local Geodetic Systems" NIMA TR8350.2] (https://web.archive.org/web/20200101071826/https://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf)
Equations 4-14 and 4-15 in the above reference describe the transformation from latitude $\phi$, longitude $\lambda$, and altitude $h$ are in the WGS84 reference system to Cartesian Earth-centered, Earth-fixed (ECEF) coordinates. The equations below use two key parameters that describe the shape of the Earth (see tables 3.1 and 3.3 of the above reference):
$$ \begin{aligned} a &= 6378137\ \text{m} && \text{Earth equatorial radius} \\ e^2 &= 6.69437999014\times10^{-3} && \text{Square of Earth eccentricity} \end{aligned} $$
First you need to compute the "radius of curvature in the prime vertical" (equation 4-15 in the reference):
$$N = \frac a {\sqrt{1-e^2\sin^2\phi}}$$
Then simply compute the ECEF coordinates via equations 4-14: $$ \begin{aligned} x &= (N+h)\cos\phi \cos\lambda \\ y &= (N+h)\cos\phi \sin\lambda \\ z &= ((1-e^2)N+h) \sin\phi \end{aligned} $$