[AstroPy] calculating Hour Angle to obtain Parallactic angle

Wed Sep 11 12:19:18 EDT 2013

Hi Orfeo,

I wrote a Python module a few years back that you might find useful. It computes the airmass for a target throughout a give night (along with the HA, Parallactic angle, …) for a give site/observatory.

http://www.ucolick.org/~magee/observer/

It's dependent on the PyEphem (http://rhodesmill.org/pyephem/) module

Cheers,
Dan

On Sep 11, 2013, at 8:45 AM, Orfeo Colebatch <orfeo.colebatch at utoronto.ca> wrote:

> Thanks Erik, 
> 
> You're spot on, thanks for the tips on converting the radians of ephem to the coord.Angle of astropy. The hour angle being returned is now matching the website. I just need to put it in decimal hours by dividing by 15 and i should be set to put it into Ian's Astro-Python Code ( http://www.mpia-hd.mpg.de/homes/ianc/python/_modules/spec.html#parangle ) in order to calculate the hour angle. 
> 
> 
> Regrads
> orfeo
> 
> ps here's my code
> 
> 
> import ephem
> 
> from astropy import coordinates as coord
> 
> from astropy.time import Time
> 
> from astropy import units as u
> 
> from datetime import datetime as dt
> 
> def HA():
> 
>    obs = ephem.Observer()
> 
>    time = dt.strptime('2013/10/21' +' ' + '13:00:00', '%Y/%m/%d %H:%M:%S')
> 
>    obs.date = time
> 
>    obs.lon =ephem.degrees('-79.4')
> 
>    obs.lat = ephem.degrees('43.66')
> 
>    obs.elevation = 174.0
> 
>    sun = ephem.Sun()
>    sun.compute(obs)
> 
>    ra = sun.g_ra
> 
>    c= coord.ICRSCoordinates(ra=ra, dec = sun.dec, unit=(u.radian, u.radian))
> 
>    t = coord.Angle(obs.sidereal_time(),u.radian)
> 
>    t.lat = 43.66
> 
>    t.lon =-79.4
> 
>    ha = coord.angles.RA.hour_angle(c.ra,t)
> 
>    return ha
> 
> if __name__ == '__main__':
>    print ' Calculated hour angle is ', HA()
> 
> 
> 
> 
> 
> ________________________________________
> From: etollerud at gmail.com [etollerud at gmail.com] on behalf of Erik Tollerud [erik.tollerud at gmail.com]
> Sent: 11 September 2013 10:52
> To: Orfeo Colebatch
> Cc: astropy at scipy.org
> Subject: Re: [AstroPy] calculating Hour Angle to obtain Parallactic angle
> 
> Ah, that's your problem.  As I think you've now surmised, the time
> that goes into `hour_angle` is the Local Sidereal Time, *not* the
> local time.  Given that you're already using pyephem, the easiest
> thing is to just use the Observer.sidereal_time method.  (E.g.,
> obs.sidereal_time() ) in your example.
> 
> One other thing you should know, though: pyephem uses radians for
> essentially everything. So in you're example, when you do e.g.,
> ``obs.lon =-79.4``, pyephem is treating that is -79.4 radians, not
> degrees.  I think you can use strings somehow to tell it to use
> degrees, but you'll have to look more closely at the documentation, as
> I don't remember exactly how that works.  Regardless, though,
> `obs.sidereal_time` will probably give you radians out, so you'll need
> to do something like ``t = coord.Angle(obs.sidereal_time(),
> u.radian)`` to make sure the units work out.
> 
> On Wed, Sep 11, 2013 at 10:32 AM, Orfeo Colebatch
> <orfeo.colebatch at utoronto.ca> wrote:
>> Thanks Erik and Tim,
>> 
>> 
>> I forgot to include the t object definition. I used the astropy.time module to set
>> 
>> t = Time('2013-06-21 13:00:00', format = 'iso', scale = 'utc')
>> 
>> I wasn't sure if i needed to set this time (t) to Toronto Canada time (UTC-4 currently) but i thought as long as i used the actual UTC time which corresponds to 09:00:00 Toronto time that would be ok. I corrected my code below (also changed the ephem time to 13:00:00). Still no change in the Hour Angle.
>> 
>> 
>> I'm looking into sidereal time now as that looks like it's probably the cause, my vague understanding of LST is that i need to correct my Local time (or UTC) by the longitude i sit at (Toronto -79.4 W).
>> 
>> Hopefully i can figure this out but if you have any further pointers i'm all ears.
>> 
>> Regards
>> orfeo
>> 
>> 
>> 
>> 
>> 
>> 
>> import ephem
>> 
>> from astropy import coordinates as coord
>> 
>> from astropy.time import Time
>> 
>> from astropy import units as u
>> 
>> from datetime import datetime as dt
>> 
>> def HA():
>>    obs = ephem.Observer()
>> 
>>    time = dt.strptime('2013/6/21' +' ' + '13:00:00', '%Y/%m/%d %H:%M:%S')
>> 
>>    obs.date = time
>> 
>>    obs.lon =-79.4
>> 
>>    lat = 43.66
>> 
>>    obs.elevation = 174.0
>> 
>>    sun = ephem.Sun()
>> 
>>    sun.compute(obs)
>> 
>>    ra = sun.g_ra
>> 
>>    c= coord.ICRSCoordinates(ra=ra, dec = sun.dec, unit=(u.radian, u.radian))
>> 
>>    t = Time('2013-06-21 13:00:00', format = 'iso', scale = 'utc')
>> 
>>    t.lat = 43.66
>> 
>>    t.lon =-79.4
>> 
>>    ha = coord.angles.RA.hour_angle(c.ra,t)
>> 
>>    return ha
>> 
>> 
>> if __name__== '__main__':
>>    print 'Calculating Hour Angle as ', HA()
>> 
>> 
>> --------------------------------------------
>> 
>> In [1]: import ha_astropy
>> 
>> In [2]: ha_astropy.HA()
>> Out[2]: <Angle 277.78141 deg>
>> 
>> 
>> 
>> 
>> 
>> 
>> 
>> 
>> 
>> ________________________________________
>> From: etollerud at gmail.com [etollerud at gmail.com] on behalf of Erik Tollerud [erik.tollerud at gmail.com]
>> Sent: 10 September 2013 21:00
>> To: Orfeo Colebatch; astropy at scipy.org
>> Subject: Re: [AstroPy] calculating Hour Angle to obtain Parallactic angle
>> 
>> 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
> 
> 
> 
> --
> Erik
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy
> 
> !DSPAM:424,5230903518667127246057!
> 