Observer.target_hour_angle() and parallactic_angle() do not precess coordinates

[bjweiner]

I am doing something that needs fairly accurate parallactic angles, and discovered that the astroplan calculation of parallactic angle does not agree with the skycalc C program, for an ICRS SkyCoord and obstime in 2025. It does agree with skycalc in year 2000. The attached script sets up a calculation for a few times from 2000 to 2099 (same date and time at 25 year intervals) and prints the hour angle and parallactic angle, and it shows only minuscule change over 99 years. The hour angle should change by almost 5 minutes of time from 2000-04-02 to 2099-04-02, with a corresponding change in parallactic angle.

I tried using the plot_airmass and plot_parallactic functions to make plots for the same object in year 2000 and year 3000, and found that the plot_airmass correctly shows transit changing by about 50 minutes from 2000 to 3000, but plot_parallactic shows no change.

I believe this is because the Observer.target_hour_angle() and .parallactic_angle() methods are using the target.ra and .dec in the native frame (like ICRS) without precessing them to the time of observation.

astroplan_parang_test.py

[bjweiner]

It occurred to me that the parallactic angle should have an extra step. The methods should transform the target RA and Dec from the target native frame to the current RA and Dec (Tim Pickering noted to me that it should be apparent RA and Dec, at the observer location) to calculate the hour angle and parallactic angle at the current time. For hour angle, the hour angle at current time can be returned, but the parallactic angle should be transformed back to the frame that the user input. In other words, if the user calls this method on a target with equinox 2000 (or ICRS) coordinates, for an obstime in 2026, we need to precess to 2026 to calculate HA and parallactic angle. Now the parallactic angle is the sky PA in equinox 2026 that points to the zenith. So the PA should be transformed back to equinox 2000 and returned, because the user is probably expecting a PA in the 2000 system (the system of their input coordinates).

[tepickering]

to be clear, the definition of parallactic_angle at

astroplan/astroplan/observer.py

Line 610 in a4d89a0

def parallactic_angle(self, time, target, grid_times_targets=False,

is using the provided target.ra and target.dec as-is without checking what frame those coordinates are defined in (e.g. ICRS or FK5). the correct way to calculate parallactic_angle would be to transform the target coordinates to an AltAz frame and then calculate it using azimuth, elevation, and latitude. this will do all the proper accounting for precession and earth position/time. https://astronomy.stackexchange.com/questions/34086/how-to-calculate-parallactic-angle-from-a-fixed-alt-az-position provides a nice summary of how to do this using equations from jean meeus's book.

[bjweiner]

To give an example, here is a standalone function I wrote to calculate the hour angle and parallactic angle for a target given an AltAz frame. For an Observer method, one would write it a little differently calling observer.altaz() to create that frame. The one thing this routine does not do, is transform the parallactic angle back from FK5(equinox=obstime) to the frame of the input coordinates if wanted (SkyCoord provides easy transformation of coordinates between frames, but transforming position angles is not built in as far as I know).

# arguments are a SkyCoord and AltAz
def parangle(field, altaz_frame):
    field_altaz = field.transform_to(altaz_frame)
    lon = field_altaz.location.lon
    lat = field_altaz.location.lat
    field_now = field_altaz.transform_to(FK5(equinox=field_altaz.obstime))
    ranow = field_now.ra
    decnow = field_now.dec
    lmst = field_altaz.obstime.sidereal_time('mean', longitude=lon)
    hourang = lmst - ranow
    parang = Angle(np.arctan2(np.sin(hourang), (np.tan(lat) * np.cos(decnow) - np.sin(decnow) * np.cos(hourang))) )
    return hourang, parang