[AstroPy] calculating Hour Angle to obtain Parallactic angle

[Orfeo Colebatch]

Hi All,

We have a alt/az sun tracker whose coding requires the calculation of the parallactic angle in order to do some coordinate transformations between the alt/az mount and a camera monitoring the image it produces.

I have found Ian's Astro Python code (http://www.mpia-hd.mpg.de/homes/ianc/python/_modules/spec.html#parangle) which works when i use it conjunction with the coordinates output by this site (http://www.jgiesen.de/astro/suncalc/) when the Hour Angle is converted to decimal hours.

I have tried to calculate the same hour angle using coordinates.angles.RA.hour_angle (among other codes) bu have not been able to reproduce that Hour Angle (HA) from the above site.

E.g for 21st June 2013, 13:00 UTC (9:00 GMT-4), 43.66 deg latitude, -79.4 W longitude I get a HA of -64.73 degreees  (or +295 degrees)


By doing the following in python with ephem and astropy i get the following


In [1]: import ephem

In [2]: from astropy import coordinates as coord

In [3]: from astropy.time import Time

In [4]: from astropy import units as u

In [5]: obs = ephem.Observer()

In [6]: from datetime import datetime as dt

In [7]: time = dt.strptime('2013/6/21' +' ' + '09:00:00', '%Y/%m/%d %H:%M:%S')

In [8]: obs.date = time

In [9]: obs.lon =-79.4

In [10]: obs.lat = 43.66

In [11]: obs.elevation = 174.0

In [12]: sun = ephem.Sun()

In [13]: sun.compute(obs)

In [14]: ra = sun.g_ra

In [15]: ra
Out[15]: 1.573768747713852

In [16]: coord.ICRSCoordinates(ra=ra, dec = sun.dec, unit=(u.radian, u.radian))
Out[16]: <ICRSCoordinates RA=90.17031 deg, Dec=23.43624 deg>

In [17]: c= coord.ICRSCoordinates(ra=ra, dec = sun.dec, unit=(u.radian, u.radian
))

In [18]: c
Out[18]: <ICRSCoordinates RA=90.17031 deg, Dec=23.43624 deg>

In [19]: c.ra
Out[19]: <RA 90.17031 deg>

In [20]: sun.ra
Out[20]: 1.573725506944724

In [21: c.ra * ephem.pi/180
Out[21]: <Angle 1.57377 deg>

In [22]: c.ra
Out[22]: <RA 90.17031 deg>

In [23]: coord.angles.RA.hour_angle(c.ra,t)
Out[23]: <Angle 277.95469 deg>

In [24]: t.lat = 43.66

In [25]: t.lon =-79.4

In [26]: coord.angles.RA.hour_angle(c.ra,t)
Out[26]: <Angle 277.95469 deg>



As you can see the angle is off by 20 degrees. I know this must be a unit problem, but i'm not sure what as of yet. If anyone has any idea's i'd be interested. Perhaps i should be calculating RA inside astropy? and not switching between ephem and astropy?


I do have access to this IDL code (http://idlastro.gsfc.nasa.gov/) and the code for the above sun web page( http://www.esrl.noaa.gov/gmd/grad/solcalc/calcdetails.html) but i would like to avoid coding that all up in python.

If anyone has any advise i'm all ears.

kind regards


Orfeo






-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/astropy/attachments/20130910/f8754f16/attachment.html>