EDIT: Just found out that I need to convert latitude, longitude and elevation of a location on earth to J2000 coordinates and nothing to do with ra/dec or the moon. Sorry for this. Your answers did give me a lot of insights. Please see the edited question below.
Question: how do i convert latitude, longitude and elevation to J2000 coordinates (XYZ). Is there a conversion present in ephem? I checked the docs but I couldnt find something I need (or mightve overlooked something due to my lack of knowledge in this field). Thanks
***************** OLD (Disregard) ******************
I have the moon position in Right Ascension (RA) and Declination (Dec) and I want to convert them into X Y Z coordinates. Is there a built-in PyEphem function for this? Also, what is the math behind it? Thanks.
EDIT: I am using the J2000 coordinate system (which is equatorial i think, this is my first time working with astronomy). I have the distance to moon available. The ra/dec values are already in J2000 (equatorial) coordinates.
X points North
Y points West
Z points towards the sky
Best answer:
It has just come to my attention that, in June 2011, the Naval Observatory released a Python interface to the powerful NOVAS reference software with which the highest-precision astronomical computations are performed:
http://aa.usno.navy.mil/software/novas/novas_py/novaspy_intro.php
With this library you can get the answer you are seeking, at far higher precision than PyEphem has ever offered:
from novas import compat as novas
jd_tt = novas.julian_date(2012, 9, 8, 12.00)
delta_t = 66.603 # from http://maia.usno.navy.mil/ser7/deltat.preds
lat = 42.3583 # positive is north
lon = -71.0603 # negative is west
observer = novas.make_observer_on_surface(lat, lon, 0, 0, 0)
print novas.geo_posvel(jd_tt, delta_t, observer)
On my machine this gives the answer:
((-3.5081406460494928e-06, 3.135277713623258e-05, 2.858239912567112e-05), (-0.00019753847060164693, -2.2329942271278055e-05, 2.4885824275734915e-07))
Try it yourself and see if this gives you the kind of results that you need!
Newer answer:
It appears that the answer is “no” — PyEphem, to my surprise, gives no easy way to get the answer to the question "where, in x, y, z coordinates, is (say) Boston at time t ?”
This is a surprise because “libastro”, the library behind PyEphem, of course has to compute this internally in order to figure out where other objects are relative to an observer. It seems to do so in two places. In parallax.c it defines ta_par() which talks only about angles on the outside, but on its inside you can see that it temporarily computes the x, y, z of the observer. You can even see the important constant 298.257 hidden inside there, which measures how flat the earth is, since it is not a perfect sphere.
The other place is in earthsat.c which looks like a completely different code base from the rest of “libastro”, and so it duplicates some of the logic. Its EarthFlat constant of 298.25 is a bit less precise, but is doing the same job. And its function, GetSitPosition(), actually exposes x,y,z coordinates instead of keeping them hidden. But it is declared static so there is no way to call in to this useful function from outside!
So for the moment, PyEphem gives you no way to compute your x,y,z directly. But it does provide one important piece of information: the current sidereal time, which you will (I think) be able to use to figure how far around the earth Boston (or wherever) has traveled by time t, which will be important in figuring out your coordinates.
I will see if I can work up a quick solution in Python that combines the hour angle from PyEphem with some explicit trigonometry to get you an answer. But, for the moment, no: PyEphem does not expose this information directly, sadly enough; I will put it on the list of things for a future version!
Older answer, from when the question was about the x,y,z position of the Moon:
PyEphem does not, alas, have built-in functions for converting from the polar coordinates used in amateur astronomy to the x/y/z coordinates which will let you map out how objects are distributed in space around the Earth. But the conversion is easy to do yourself:
import ephem
import math
m = ephem.Mars('2012/8/1')
print m.ra, m.dec
x = math.cos(m.dec) * math.cos(m.ra)
y = math.cos(m.dec) * math.sin(m.ra)
z = math.sin(m.dec)
print x, y, z
print 'sanity check: vector length =', math.sqrt(x*x + y*y + z*z)
The output of this script is:
12:58:51.20 -6:24:05.6
-0.961178016954 -0.252399543786 -0.111495695074
sanity check: vector length = 1.0
The position of Mars for the random date that I used here are quite reasonable values: an RA that is almost one hour more than halfway around the great circle (since 12h would be exactly halfway), and a declination that pushes the position a bit south. Thus the x, y, and z that we get out: the z is a slightly negative number since -6° is indeed south of the equator, and x and y are both negative since going 13h around a 24h circle puts you down in the negative/negative quadrant of a normal unit circle.
Note that although J2000 has a north and south — so that we can truthfully say that the slightly negative z is a southward direction — it does not have an east and west, since the earth turning below it is constantly swinging east and west in all directions. Instead, RA measures from “the first point of Ares” which is the direction in which the sun lies during the spring Equinox. So x and y are not east or west; they are coordinates pointing out into the solar system on a fixed axis defined by the direction that the Earth sits in every Spring.
This x y z vector I have created is a “unit vector” — a little vector that has the magnitude 1.0, as I verified in the script to make sure I had the formulae correct. If you were computing x y and z coordinates for objects whose distance from the earth you knew, then you could get a real vector — whose magnitude were distances, instead of fractions of 1 — by multiplying each of the three x y and z by the distance to the object.
Does that help you out? From your description — and your question about east and west — I could not tell if you wanted RA and dec turned into x y z or whether you are actually wanting the azimuth and altitude converted (but the math is the same either way). That would look something like:
x = math.cos(m.alt) * math.cos(m.az)
y = math.cos(m.alt) * math.sin(m.az)
z = math.sin(m.alt)
What are you trying to accomplish with these coordinates? That could help us make sure that we are giving them to you in a useful format.