[AstroPy] calculating Hour Angle to obtain Parallactic angle

Wed Sep 11 11:45:41 EDT 2013

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