[AstroPy] calculating Hour Angle to obtain Parallactic angle
Erik Tollerud
erik.tollerud at gmail.com
Tue Sep 10 21:00:42 EDT 2013
This example seems to be incomplete, because I can't see where `t` is
defined. That's important, because `hour_angle` doesn't compute
sidereal time based on location, it just accepts and LST and gives you
out an hour angle. Are you perhaps missing a line somewhere above?
On Tue, Sep 10, 2013 at 2:14 PM, Orfeo Colebatch
<orfeo.colebatch at utoronto.ca> wrote:
> 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
>
>
>
>
>
>
>
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy
>
--
Erik
More information about the AstroPy
mailing list