I just start learning astronomy.
I use pyephem to understand different coordinate systems such as geocentric equtorial, geocentric ecliptic, heliocentric ecliptic.
In phephem, planet body such as mars has attributes (ra,dec) in geocentric equtorial, (hlon, hlat) in heliocentric ecliptic. I try to check the consistency between them.
My method(scripts attached below) is as follows:
- define function earthtilt to get earth axis tilt
- define functions radec2cart and cart2radec to perform coordinate transform between spherical coordinates and cartesian coordinates.
- define function R_x to get rotation matrix.
- get geocentric equtorial coordinates of sun/mars using (ra,dec) attributes.
- use R_x and earthtilt to rotate geocentric equtorial coordinates to geocentric ecliptic coordinates
- (heliocentric ecliptic cartesian coordinates of mars)= (geocentric ecliptic coordinates of mars) - (geocentric ecliptic coordinates of sun)
- change above heliocentric cartesian ecliptic coordinates of mars to heliocentric spherical ecliptic coordinates of mars.
- compare above heliocentric spherical ecliptic coordinates of mars with mars attributes hlon and hlat.
- I got : difference between ecliptic longitude is about 16.9". difference between ecliptic latitude is about 9.2". Such a defference is a little too large.
- Is my computation/concept wrong or ...?
===================
import ephem
import numpy as np
import math
def R_x(theta):
return np.array([[1,0,0],[0,math.cos(theta),math.sin(theta)],[0,-math.sin(theta),math.cos(theta)]])
def radec2cart(rho,ra,dec):
return np.array([rho*math.cos(dec)*math.cos(ra),rho*math.cos(dec)*math.sin(ra),rho*math.sin(dec)])
def cart2radec(X):
r=math.sqrt(X[0]*X[0]+X[1]*X[1])
ra=ephem.hours(math.atan2(X[1],X[0])).norm
dec=ephem.degrees(math.atan2(X[2],r))
rho=math.sqrt(X[0]*X[0]+X[1]*X[1]+X[2]*X[2])
return (rho,ra,dec)
def earthtilt(date):
pos=ephem.Ecliptic('90','0',epoch=date)
return pos.to_radec()[1]
def verifyhlonhlat(star,date):
star.compute(date)
sun=ephem.Sun(date)
starx=radec2cart(star.earth_distance,star.g_ra,star.g_dec)
eps=earthtilt(date)
starx=R_x(eps).dot(starx)
sunx=radec2cart(sun.earth_distance,sun.g_ra,sun.g_dec)
sunx=R_x(eps).dot(sunx)
starx=starx-sunx
distance,hlon,hlat=cart2radec(starx)
hlon=ephem.degrees(hlon)
print([distance-star.sun_distance,ephem.degrees(hlon-star.hlon),ephem.degrees(hlat-star.hlat)])
verifyhlonhlat(ephem.Mars(),ephem.Date('2016/6/7'))
===========
>>> verifyhlonhlat(ephem.Mars(),ephem.Date('2016/6/7'))
[1.3654961393161358e-05, -0:00:16.9, 0:00:09.2]
