[AstroPy] calculating Hour Angle to obtain Parallactic angle

Orfeo Colebatch orfeo.colebatch at utoronto.ca
Wed Sep 11 10:32:58 EDT 2013


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



More information about the AstroPy mailing list