[AstroPy] calculating Hour Angle to obtain Parallactic angle

Wed Sep 11 10:52:39 EDT 2013

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