title: consistency between attributes (ra,dec) and (hlon, hlat)
body: <p>I just start learning astronomy.</p>
  
  <p>I use pyephem to understand different coordinate systems such as geocentric equtorial,
  geocentric ecliptic, heliocentric ecliptic.</p>
  
  <p>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.</p>
  
  <p>My method(scripts attached below) is as follows:</p>
  
  <ol>
  <li>define function earthtilt to get earth axis tilt</li>
  <li>define functions radec2cart and cart2radec to perform coordinate transform
   between spherical coordinates and cartesian coordinates.</li>
  <li>define function R_x to get rotation matrix.</li>
  <li>get geocentric equtorial coordinates of sun/mars using (ra,dec) attributes.</li>
  <li>use R_x and earthtilt to rotate geocentric equtorial coordinates to geocentric ecliptic coordinates</li>
  <li>(heliocentric ecliptic cartesian coordinates of mars)=
   (geocentric ecliptic coordinates of mars) - (geocentric ecliptic coordinates of sun)</li>
  <li>change above heliocentric cartesian ecliptic coordinates of mars to 
   heliocentric spherical ecliptic coordinates of mars.</li>
  <li>compare above heliocentric spherical ecliptic coordinates of mars with 
   mars attributes hlon and hlat.</li>
  <li>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.</li>
  <li>Is my computation/concept wrong or ...?</li>
  </ol>
  
  <pre><code>
      ===================
      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]
      </code></pre>