Computing Parallactic Angle

[JSKenyon]

Hello! Apologies if this is not the correct place for posting questions. I am looking to compute parallactic angles but I have managed to get myself completely muddled. Basically, given the location of a telescope on the Earth's surface and the ra and dec of a source/field, what is the easiest (and most accurate) way to compute the parallactic angle?

[JSKenyon]

Just an update on what I have managed so far. There are still many gaps in my understanding so I would really appreciate some assistance. The following almost works:

from datetime import datetime, timezone
import numpy as np
import skyfield.api as sky
import skyfield.toposlib as skytop
import skyfield.trigonometry as skytrig

datetimes = np.array(
    [
       datetime(2013, 8, 25, 10, 21, 39, 499798, tzinfo=timezone.utc),
       datetime(2013, 8, 25, 10, 21, 44, 500294, tzinfo=timezone.utc),
       datetime(2013, 8, 25, 10, 21, 49, 499502, tzinfo=timezone.utc),
       datetime(2013, 8, 25, 10, 21, 54, 499998, tzinfo=timezone.utc),
       datetime(2013, 8, 25, 10, 21, 59, 500494, tzinfo=timezone.utc)
    ]
)

timescale = sky.load.timescale()
times = timescale.from_datetimes(datetimes)

field_centre = sky.Star(
    ra=sky.Angle(degrees=85.65057465185124),
    dec=sky.Angle(degrees=49.85200932178193)
)

earth = sky.load('de421.bsp')['earth']

# Input antenna positions in an ECEF frame (ITRF).
ant_positions_ecef = np.array(
    [
        [-1.60108719e+06, -5.04133984e+06,  3.55581586e+06],
        [-1.60116259e+06, -5.04182900e+06,  3.55509589e+06],
        [-1.60062293e+06, -5.04233178e+06,  3.55462705e+06],
        [-1.59992611e+06, -5.04277298e+06,  3.55431979e+06]
    ]
)

# Specify input antenna positions as ITRS positions. Is there an easy way to
# go from these to geodetic values (lat, lon, dist)?
ant_positions_itrf = [
    skytop.ITRSPosition(sky.Distance(m=pos)) for pos in ant_positions_ecef
]

apparent_positions = [
    (earth + pos).at(times).observe(field_centre).apparent()
    for pos in ant_positions_itrf
]

# Get the zenith for each antenna at each time. This is the part that
# doesn't currently work and which I may be misunderstanding.
zeniths = [
    (earth + pos).at(times).from_altaz(alt_degrees=90, az_degrees=0)
    for pos in ant_positions_itrf
]

# The parallactic angle can then be computed as follows (ordering of aguments
# may be incorrect).
parallactic_angles = [
    skytrig.position_angle_of(a.radec(), z.radec())
    for a, z in zip(apparent_positions, zeniths)
]

print(parallactic_angles)

Edit: Meant to give a little context. I am looking at computing parallactic angles for a radio interferometer, hence multiple antenna (telescope) locations.

[Bernmeister]

Although this is NOT python, the algorithms/calculations are still valid (I believe). Look for the function parallacticAngle and that'll show you how to calculate.

EDIT:

If you download the source for my indicator, find astrobase.py which has an implementation you might be able to use as a basis called getZenithAngleOfBrightLimb(). Note that I mistakenly called the parameters bodyLat and bodyLong whereas they should be called observerLat and observerLon, respectively.

[JSKenyon]

@Bernmeister thanks for that - I will definitely take a look! I was sort of hoping to get the above code snippet working as I know that I can get skyfield to give decent results if I feed it geodetic coordinates computed by another package. I wanted to see if it was possible to do the entire calculation using skyfield as it will allow me to avoid a dependency on another package which has thread safety issues.

[Bernmeister]

@JSKenyon I implemented my python version without actually understanding the underlying detail (and so it is above by pay grade to sensibly comment on your code). Perhaps also look at the references in my code as that could help you determine if you are calculating the constituent parts correctly (hour angle or sidereal time). I made sure I got the same answer as the frink example along each step and that helped me to get to the final answer.

[brandon-rhodes]

Good morning, @JSKenyon!

With the approach of the holidays, I'm now finding time to dig into Skyfield issues like this one; thanks for your patience, and I hope my answer is in time to at least be of some use to you.

You're right that there doesn't seem to be a simple way to convert a plain ITRF $x,y,z$ position into a geographic coordinate with latitude, longitude, and elevation. I should probably add one, as the current way of doing so is a little inefficient: because Skyfield is built around the fixed ICRS reference frame, an ITRF position first has to be rotated into the ICRS using a specific time t, then the geographic coordinate has to be generated by rotating back the other way. Before we close this issue, I'll see if I can stare at Skyfield's code and see where we might fit in a shortcut to move directly from $x,y,z$ to geography.

But in the meantime, here's a working version of your script—try it out and let me know if the resulting angles look reasonable!

from datetime import datetime, timezone
import numpy as np
import skyfield.api as sky
import skyfield.toposlib as skytop
import skyfield.trigonometry as skytrig

datetimes = np.array(
    [
       datetime(2013, 8, 25, 10, 21, 39, 499798, tzinfo=timezone.utc),
       datetime(2013, 8, 25, 10, 21, 44, 500294, tzinfo=timezone.utc),
       datetime(2013, 8, 25, 10, 21, 49, 499502, tzinfo=timezone.utc),
       datetime(2013, 8, 25, 10, 21, 54, 499998, tzinfo=timezone.utc),
       datetime(2013, 8, 25, 10, 21, 59, 500494, tzinfo=timezone.utc)
    ]
)

timescale = sky.load.timescale()
times = timescale.from_datetimes(datetimes)

field_centre = sky.Star(
    ra=sky.Angle(degrees=85.65057465185124),
    dec=sky.Angle(degrees=49.85200932178193)
)

earth = sky.load('de421.bsp')['earth']

# Input antenna positions in an ECEF frame (ITRF).
ant_positions_ecef = np.array(
    [
        [-1.60108719e+06, -5.04133984e+06,  3.55581586e+06],
        [-1.60116259e+06, -5.04182900e+06,  3.55509589e+06],
        [-1.60062293e+06, -5.04233178e+06,  3.55462705e+06],
        [-1.59992611e+06, -5.04277298e+06,  3.55431979e+06]
    ]
)

# Specify input antenna positions as ITRS positions.
ant_positions_itrf = [
    skytop.ITRSPosition(sky.Distance(m=pos)) for pos in ant_positions_ecef
]

# Convert ITRS positions into geographic positions.
t = times[0]  # doesn't matter what time
ant_positions_geo = [
    sky.wgs84.geographic_position_of(pos.at(t)) for pos in ant_positions_itrf
]

apparent_positions = [
    (earth + pos).at(times).observe(field_centre).apparent()
    for pos in ant_positions_geo
]

# Get the zenith for each antenna at each time. This is the part that
# doesn't currently work and which I may be misunderstanding.
zeniths = [
    (earth + pos).at(times).from_altaz(alt_degrees=90, az_degrees=0)
    for pos in ant_positions_geo
]

# The parallactic angle can then be computed as follows (ordering of aguments
# may be incorrect).
parallactic_angles = [
    skytrig.position_angle_of(a.radec(), z.radec())
    for a, z in zip(apparent_positions, zeniths)
]

print(parallactic_angles)

[JSKenyon]

Thank you for the response @brandon-rhodes. I have been on leave so I have not yet had a chance to test it out. I will do so during the course of the day and get back to you!

[JSKenyon]

This does work (although I am still struggling to get adequate agreement between skyfield and python-casacore, for any radio astronomer who stumbles over this). Thanks for the assist @brandon-rhodes!

[JSKenyon]

Hi @brandon-rhodes - this is still something I am interested in. The legacy package we currently rely on is a difficult dependency. Did you end up adding this functionality? No worries if not - just want to make sure that I am not longing for something that already exists.

[brandon-rhodes]

Oh, perhaps I misunderstood — I thought you were going to be able to switch from your old dependency to the calculation in the script that I offered above. Did it wind up giving results different from the ones you need?

No, I have not had time yet to sit down and see if Skyfield can be tweaked to make the computation any more convenient. But in any case the calculation would be the same as in the script I gave above, just packaged a little differently.

[JSKenyon]

Oh, perhaps I misunderstood — I thought you were going to be able to switch from your old dependency to the calculation in the script that I offered above. Did it wind up giving results different from the ones you need?

I think that I could never quite get them to agree but it has also been a long time since I did those experiments. I should probably revisit them.

No, I have not had time yet to sit down and see if Skyfield can be tweaked to make the computation any more convenient. But in any case the calculation would be the same as in the script I gave above, just packaged a little differently.

All good! Thanks for the response. I think the onus is on me to revisit my experiments and determine if the results are similar enough - convincing astronomers to move away from older software is always tough in my experience!

[brandon-rhodes]

Thanks for the update. I'll be interested to learn whether the script now returns accurate results, since I wouldn't want to turn it into a real routine called "parallactic_angle" or something if it really isn't accurately computing that quantity.

If you do find disagreement, and can't yourself see how to tweak the Python code, I'd need to know what exact angles the other tool printed for the same circumstances, so that I would have a goal to work towards. Thanks!

[Bernmeister]

@brandon-rhodes I've taken a look at my code and have noticed a difference in the naming of skyfield.trigonometry.position_angle_of( anglepair1, anglepair2 ) and what is calculated for return, versus what I have done in the past.

Note: this is not a difference in values returned (a calculation error) but I suspect what is calculated in skyfield is not actually the position angle.

Assuming I've not misunderstood the reading of the references and then my implementation, here is my take on computing the zenith angle of the bright limb (how much is the moon tilted). At the end of this function, to compute the zenith angle I need the parallactic angle. The function below returns the same value as does skyfield.trigonometry.position_angle_of( anglepair1, anglepair2 ) and so I believe there is a naming issue here (not the calculation/value itself).

    @staticmethod
    def getZenithAngleOfBrightLimb( utcNow, sunRA, sunDec, bodyRA, bodyDec, observerLat, observerLon ):
        # Astronomical Algorithms by Jean Meeus, Second Edition, Equation 48.5
        y = math.cos( sunDec ) * math.sin( sunRA - bodyRA )
        x = math.sin( sunDec ) * math.cos( bodyDec ) - math.cos( sunDec ) * math.sin( bodyDec ) * math.cos( sunRA - bodyRA )
        positionAngleOfBrightLimb = math.atan2( y, x )

        # Multiply by 15 to convert from decimal time to decimal degrees; section 22 of Practical Astronomy with Your Calculator.
        localSiderealTime = math.radians( getSiderealTime( utcNow, observerLon ) * 15 )

        # Astronomical Algorithms by Jean Meeus, Second Edition, page 92.
        # https://tycho.usno.navy.mil/sidereal.html
        # http://www.wwu.edu/skywise/skymobile/skywatch.html
        hourAngle = localSiderealTime - bodyRA

        # Astronomical Algorithms by Jean Meeus, Second Edition, Equation 14.1
        y = math.sin( hourAngle )
        x = math.tan( observerLat ) * math.cos( bodyDec ) - math.sin( bodyDec ) * math.cos( hourAngle )
        parallacticAngle = math.atan2( y, x )

        return ( positionAngleOfBrightLimb - parallacticAngle ) % ( 2.0 * math.pi )

I can include the function getSiderealTime( utcNow, longitude ) if need be.

Bottom line: I suspect skyfield.trigonometry.position_angle_of( anglepair1, anglepair2 ) is computing the parallactic angle rather than the position angle. The usual caveat of please check my working applies here more than ever!

[brandon-rhodes]

I can include the function getSiderealTime( utcNow, longitude ) if need be.

We might need that at some point — but could you instead try sharing one particular moment in time, and what the result of your function is for, say, the Moon and Sun? That way I could plot their positions and make physical sense of the angle your routine is returning. Then, I could compare that with Skyfield to address your concern. Thanks!

[Bernmeister]

Here is the sidereal time function:

    # Compute the sidereal decimal time for the given longitude (in floating point radians).
    #
    # Reference:
    #  'Practical Astronomy with Your Calculator' by Peter Duffett-Smith.
    @staticmethod
    def getSiderealTime( utcNow, longitude ):
        # Find the Julian date.  Section 4 of the reference.
        # Assume the date is always later than 15th October, 1582.
        y = utcNow.year
        m = utcNow.month
        d = \
            utcNow.day + \
            ( utcNow.hour / 24 ) + \
            ( utcNow.minute / ( 60 * 24 ) ) + \
            ( utcNow.second / ( 60 * 60 * 24 ) ) + \
            ( utcNow.microsecond / ( 60 * 60 * 24 * 1000 ) )

        if m == 1 or m == 2:
            yPrime = y - 1
            mPrime = m + 12

        else:
            yPrime = y
            mPrime = m

        A = int( yPrime / 100 )
        B = 2 - A + int( A / 4 )
        C = int( 365.25 * yPrime )
        D = int( 30.6001 * ( mPrime + 1 ) )
        julianDate = B + C + D + d + 1720994.5

        # Find universal time.  Section 12 of the reference.
        S = julianDate - 2451545.0
        T = S / 36525.0
        T0 = ( 6.697374558 + ( 2400.051336 * T ) + ( 0.000025862 * T * T ) ) % 24
        universalTimeDecimal = ( ( ( utcNow.second / 60 ) + utcNow.minute ) / 60 ) + utcNow.hour
        A = universalTimeDecimal * 1.002737909
        greenwichSiderealTimeDecimal = ( A + T0 ) % 24

        # Find local sidereal time.  Section 14 of the reference.
        longitudeInHours = math.degrees( longitude ) / 15

        return ( greenwichSiderealTimeDecimal + longitudeInHours ) % 24 # Local sidereal time as a decimal time.

[Bernmeister]

We might need that at some point

See above.

— but could you instead try sharing one particular moment in time, and what the result of your function is for, say, the Moon and Sun? That way I could plot their positions and make physical sense of the angle your routine is returning. Then, I could compare that with Skyfield to address your concern. Thanks!

See below...

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# Test computing parallactic angle / postiion angle.

import datetime
import math

import ephem

from skyfield.api import load, wgs84
from skyfield.trigonometry import position_angle_of

latitude = -33.0
longitude = 151.0
elevation = 10.0

utcNow = datetime.datetime.now( datetime.timezone.utc )

print( f"Current date/time is { utcNow } at latitude/longitude { latitude } / { longitude }" )

# Functions

def getSiderealTime( utcNow, longitude ):
    y = utcNow.year
    m = utcNow.month
    d = \
        utcNow.day + \
        ( utcNow.hour / 24 ) + \
        ( utcNow.minute / ( 60 * 24 ) ) + \
        ( utcNow.second / ( 60 * 60 * 24 ) ) + \
        ( utcNow.microsecond / ( 60 * 60 * 24 * 1000 ) )

    if m == 1 or m == 2:
        yPrime = y - 1
        mPrime = m + 12

    else:
        yPrime = y
        mPrime = m

    A = int( yPrime / 100 )
    B = 2 - A + int( A / 4 )
    C = int( 365.25 * yPrime )
    D = int( 30.6001 * ( mPrime + 1 ) )
    julianDate = B + C + D + d + 1720994.5

    # Find universal time.  Section 12 of the reference.
    S = julianDate - 2451545.0
    T = S / 36525.0
    T0 = ( 6.697374558 + ( 2400.051336 * T ) + ( 0.000025862 * T * T ) ) % 24
    universalTimeDecimal = ( ( ( utcNow.second / 60 ) + utcNow.minute ) / 60 ) + utcNow.hour
    A = universalTimeDecimal * 1.002737909
    greenwichSiderealTimeDecimal = ( A + T0 ) % 24

    # Find local sidereal time.  Section 14 of the reference.
    longitudeInHours = math.degrees( longitude ) / 15

    return ( greenwichSiderealTimeDecimal + longitudeInHours ) % 24 # Local sidereal time as a decimal time.

def getZenithAngleOfBrightLimb( utcNow, sunRA, sunDec, bodyRA, bodyDec, observerLat, observerLon ):
    # Astronomical Algorithms by Jean Meeus, Second Edition, Equation 48.5
    y = math.cos( sunDec ) * math.sin( sunRA - bodyRA )
    x = math.sin( sunDec ) * math.cos( bodyDec ) - math.cos( sunDec ) * math.sin( bodyDec ) * math.cos( sunRA - bodyRA )
    positionAngleOfBrightLimb = math.atan2( y, x )

    # Multiply by 15 to convert from decimal time to decimal degrees; section 22 of Practical Astronomy with Your Calculator.
    localSiderealTime = math.radians( getSiderealTime( utcNow, observerLon ) * 15 )
    print( f"localSiderealTime = { localSiderealTime }" )

    hourAngle = localSiderealTime - bodyRA

    y = math.sin( hourAngle )
    x = math.tan( observerLat ) * math.cos( bodyDec ) - math.sin( bodyDec ) * math.cos( hourAngle )
    parallacticAngle = math.atan2( y, x )
    print( f"parallacticAngle = { parallacticAngle }" )

    return ( positionAngleOfBrightLimb - parallacticAngle ) % ( 2.0 * math.pi )

# PyEphem
ephemNow = ephem.Date( utcNow )

observer = ephem.city( "Sydney" )
observer.lat = str( latitude )
observer.lon = str( longitude )
observer.elev = elevation
observer.date = ephemNow

moon = ephem.Moon( observer )
sun = ephem.Sun( observer )

brightLimb = getZenithAngleOfBrightLimb(
    ephemNow.datetime().replace( tzinfo = datetime.timezone.utc ),
    sun.ra, sun.dec,
    moon.ra, moon.dec,
    float( observer.lat ), float( observer.lon ) )

print( f"PyEphem bright limb angle of moon { brightLimb }" )

# Skyfield

EPHEMERIS_PLANETS = load( "de421.bsp" )

timeScale = load.timescale( builtin = True )
now = timeScale.utc( utcNow.year, utcNow.month, utcNow.day, utcNow.hour, utcNow.minute, utcNow.second )

location = wgs84.latlon( latitude, longitude, elevation )
locationAtNow = ( EPHEMERIS_PLANETS[ "EARTH" ] + location ).at( now )

moonAtNow = locationAtNow.observe( EPHEMERIS_PLANETS[ "MOON" ] )
moonAtNowApparent = moonAtNow.apparent()

sunAltAz = locationAtNow.observe( EPHEMERIS_PLANETS[ "SUN" ] ).apparent().altaz()
positionAngle = position_angle_of( moonAtNowApparent.altaz(), sunAltAz ).radians

print( f"Skyfield parallactic angle of moon { brightLimb }" )

Results:

Current date/time is 2024-03-10 10:42:30.278240+00:00 at latitude/longitude -33.0 / 151.0
localSiderealTime = 2.1060892046556803
parallacticAngle = 2.2995619490748647
PyEphem bright limb angle of moon 3.1457734275624
Skyfield parallactic angle of moon 3.1457734275624

[brandon-rhodes]

The function below returns the same value as does skyfield.trigonometry.position_angle_of( anglepair1, anglepair2 ) and so I believe there is a naming issue here (not the calculation/value itself).

While their return values are quite close, they are not actually the same — you accidentally printed the same value twice; both of your print statements say { brightLimb }. 🙂

Your script didn't print the same value when I ran it as it did for you, but that was easy to fix; the only problem was that you left the date as now rather than fixing the date to a specific value. I fixed it like this (which was easy thanks to your very detailed printout, that included the date and time, thank you!):

#utcNow = datetime.datetime.now( datetime.timezone.utc )
utcNow = datetime.datetime(2024, 3, 10, 10, 42, 30, 278240)

I was then able to run the script and get the same output as you. To make things a bit clearer for myself, I put these lines at the end:

print('Moon:')
print(moonAtNowApparent.altaz())
print('Sun:')
print(sunAltAz)
print(f"PyEphem bright limb angle of moon {brightLimb / math.tau * 360}°")
print(f"Skyfield parallactic angle of moon {positionAngle / math.tau * 360}°")

And got:

Moon:
(<Angle -27deg 56' 32.8">, <Angle 242deg 45' 05.6">, <Distance 0.00240557 au>)
Sun:
(<Angle -29deg 25' 53.6">, <Angle 242deg 45' 54.0">, <Distance 0.993294 au>)
PyEphem bright limb angle of moon 180.23954070372852°
Skyfield parallactic angle of moon 180.45088519461905°

Remarks:

• You chose a very interesting moment to run the script, a moment at which the Sun and Moon had almost exactly the same azimuth! They differ by less than a minute of arc. Was that deliberate, or were you just very lucky? 🙂
• Both are below the horizon, with the Sun being something like 1.5° farther below the horizon along the same line of azimuth.
• We would therefore expect the lower edge of the Moon to be the most illuminated, since the Sun is directly below it.
• And both routines succeed in detecting this! They both give very close to 180°, meaning that the Moon is illuminated from directly below. So all the math seems to be working great.

We can turn, then, to the question of terminology. Strictly speaking, a position angle is not any old ‘relative angle between one point on a sphere and another', but is specifically ‘the angle between two astronomical positions, measured relative to the North Celestial Pole.’ So Skyfield's position_angle_of() only returns a real position angle if it's supplied coordinates in RA/dec.

If instead position_angle_of() is given alt/az, then — I'm not sure there's a uniform name for a relative angle in the sky, measured against the zenith? It's not ‘parallactic angle’, I don't think, because that would require the other point to be the zenith, and position_angle_of() can take any point as the second argument, not only the zenith? In this paper:

https://www.seas.upenn.edu/%7Eamyers/MoonPaperOnline.pdf

—on page 4, it defines an angle α that it calls the ‘moon pointer’, that's definitely oriented with respect to azimuth and altitude, but it requires the two objects to be the Moon and Sun. If you pass two other objects, you get some other ‘pointer’ angle measured with respect to up / left / down / right. (I note that, for their investigation, they define 0° as right, rather than up.) And I don't know an official name for that concept.

Are you thinking that position_angle_of() should be something like relative_angle_in_sky() or something that would sound more agnostic about which coordinate system was in use?

[Bernmeister]

While their return values are quite close, they are not actually the same — you accidentally printed the same value twice; both of your print statements say { brightLimb }. 🙂

Mea culpa! Put it down to a (not quite) late night scramble to cobble together a test script and over-zealous cut 'n' paste.

• You chose a very interesting moment to run the script, a moment at which the Sun and Moon had almost exactly the same azimuth! They differ by less than a minute of arc. Was that deliberate, or were you just very lucky? 🙂

Fluke!

So all the math seems to be working great.

I'll take the credit for printing out the current date/time; all else is on the shoulders of giants!

We can turn, then, to the question of terminology. ... Are you thinking that position_angle_of() should be something like relative_angle_in_sky() or something that would sound more agnostic about which coordinate system was in use?

This is where I humbly say I just implemented what others created. I'll read the reference you found and see if it makes sense to me (comparing what I implemented from Meeus) and how the terminology aligns.

Originally I wanted to create an icon that represented the moon as it appeared in the sky in real time (using PyEphem). So I stumbled on the terminology of bright limb and managed to implement enough code to give me an angle to create a dynamic SVG icon! When I started migrating over to Skyfield and noticed your position_angle_of() function, I compared the two results and was glad to see we are in agreement. I did notice the terminology difference but felt it was above my pay grade to say something...

[JSKenyon]

I have managed to dig out my old test and have tried to simplify it to demonstrate the problem. I have included the code for calling the reference implementation (python-casacore, see function casa_parangles) but it is not called (spares you installing dependencies). Instead, I simply tabulated its output and compare the skyfield output to those results. Note that there is a commented block in the skyfield code which can be enabled to produce results with smaller errors. Please forgive the length of this test - tricky to make it any smaller while remaining representative.

import numpy as np

import skyfield.api as sky
import skyfield.toposlib as skytop
import skyfield.trigonometry as skytrig

from datetime import datetime, timezone

DATETIMES = np.array(
    [
       datetime(2013, 8, 25, 10, 21, 39, 499798, tzinfo=timezone.utc),
       datetime(2013, 8, 25, 10, 21, 44, 500294, tzinfo=timezone.utc),
       datetime(2013, 8, 25, 10, 21, 49, 499502, tzinfo=timezone.utc),
       datetime(2013, 8, 25, 10, 21, 54, 499998, tzinfo=timezone.utc),
       datetime(2013, 8, 25, 10, 21, 59, 500494, tzinfo=timezone.utc)
    ]
)

ANTPOS_ECEF = np.array(
    [
        [-1.60108719e+06, -5.04133984e+06,  3.55581586e+06],
        [-1.60116259e+06, -5.04182900e+06,  3.55509589e+06],
        [-1.60062293e+06, -5.04233178e+06,  3.55462705e+06],
        [-1.59992611e+06, -5.04277298e+06,  3.55431979e+06]
    ]
)

FIELD_CENTRE = [85.65057465185124, 49.85200932178193]  # RA, DEC in deg.

REFERENCE_PA = np.array(
    [
        [-1.45528121, -1.45543223, -1.45561187, -1.45577755],
        [-1.45553615, -1.4556872 , -1.45586687, -1.45603258],
        [-1.45579106, -1.45594214, -1.45612185, -1.45628759],
        [-1.45604607, -1.45619719, -1.45637694, -1.45654271],
        [-1.45630113, -1.45645229, -1.45663207, -1.45679787]
    ]
)

def parallactic_angle(ha, dec, lat):
    numer = np.sin(ha)
    denom = np.cos(dec) * np.tan(lat) - np.sin(dec) * np.cos(ha)

    return np.arctan2(numer, denom)

def compare_parangles(parangs0, parangs1):

    parangle_diff = np.abs(parangs0 - parangs1)

    max_diff = sky.Angle(radians=parangle_diff.max()).arcseconds()
    min_diff = sky.Angle(radians=parangle_diff.min()).arcseconds()
    mean_diff = sky.Angle(radians=parangle_diff.mean()).arcseconds()

    print(f"Max abs difference in arcsec: {max_diff}")
    print(f"Min abs difference in arcsec: {min_diff}")
    print(f"Mean abs difference in arcsec: {mean_diff}")

def casa_parangles(
    times,
    ant_positions_ecef,
    field_centre
):

    epoch = "J2000"
    n_ant = ant_positions_ecef.shape[0]
    n_time = np.unique(times).size

    times = np.array([t.isoformat() for t in times])

    cms = measures()

    # Assume all antenna are pointed in the same direction.
    field_centre_meas = \
        cms.direction(epoch, *(pq.quantity(fi, 'deg') for fi in field_centre))

    # Zenith is above the observer. AZELGEO is probably the correct choice,
    # as it will use geodetic coordinates i.e. based on a referece ellipsoid.
    # AZEL will assume spherical geometery which is incorrect (skyfield and
    # casacore are loosely in agreement for AZELGEO).

    zenith_azel = cms.direction(
        "AZELGEO", *(pq.quantity(fi, 'deg') for fi in (0, 90))
    )

    # Express antenna positions in ITRF (International Terrestrial Reference
    # Frame).

    ant_positions_itrf = [
        cms.position(
            'itrf', *(pq.quantity(p, 'm') for p in pos)
        ) for pos in ant_positions_ecef
    ]

    cm_angles = np.zeros((n_time, n_ant), dtype=np.float64)

    for ti, t in enumerate(np.unique(times)):
        cms.do_frame(cms.epoch("UTC", pq.quantity(t)))
        for rpi, rp in enumerate(ant_positions_itrf):
            cms.do_frame(rp)
            app_field_centre = cms.measure(field_centre_meas, "APP")
            cm_angles[ti, rpi] = \
                cms.posangle(app_field_centre, zenith_azel).get_value("rad")

    return cm_angles

def skyfield_parangles(
    times,
    ant_positions_ecef,
    field_centre
):

    n_ant = ant_positions_ecef.shape[0]
    n_time = np.unique(times).size

    # Create a "Star" object to be the field centre.
    field_centre = sky.Star(
        ra=sky.Angle(degrees=field_centre[0]),
        dec=sky.Angle(degrees=field_centre[1])
    )

    # Our time values are stored as MJD in seconds (which is weird).
    # This example avoids the problem by working with datetimes.
    ts = sky.load.timescale()
    times = ts.from_datetimes(times)

    planets = sky.load('de421.bsp')
    earth = planets['earth']

    # Specify input antenna positions as ITRS positions.
    ant_positions_itrf_sf = [
        skytop.ITRSPosition(sky.Distance(m=pos)) for pos in ant_positions_ecef
    ]

    # Convert ITRS positions into geographic positions.
    t = times[0]  # Time is irrelevant here, but is a required input.
    ant_positions_geo = [
        sky.wgs84.geographic_position_of(pos.at(t))
        for pos in ant_positions_itrf_sf
    ]

    # This is a hands-on approach that seems to work.
    sf_angles = np.zeros((n_time, n_ant), dtype=np.float64)

    # Apparent ra and dec of source relative to earth at each time.
    app_ra, app_dec, _ = \
        earth.at(times).observe(field_centre).apparent().radec(epoch='date')

    # Local apparent sidereal time, per antenna, per time.
    last = [sky.Angle(hours=pos.lst_hours_at(times)) for pos in ant_positions_geo]

    for ai in range(n_ant):
        app_ha = sky.Angle(radians=(last[ai].radians - app_ra.radians))

        sf_angles[:, ai] = parallactic_angle(
            app_ha.radians,
            app_dec.radians,
            ant_positions_geo[ai].latitude.radians
        )

    return sf_angles

if __name__ == "__main__":

    parangles_sf = skyfield_parangles(
        DATETIMES,
        ANTPOS_ECEF,
        FIELD_CENTRE,
    )

    parangles_sf = np.arctan2(np.sin(parangles_sf), np.cos(parangles_sf))

    print("CASACORE (REFERENCE) VS SKYFIELD")
    compare_parangles(REFERENCE_PA, parangles_sf)

Edit: After staring at the above some more, I saw that there was a simpler way of accomplishing my objective. Does the above seem correct? The residual errors are now at the level of 0.12 arcseconds and I don't know if skyfield or the reference is more correct.

[brandon-rhodes]

Does the above seem correct? The residual errors are now at the level of 0.12 arcseconds and I don't know if skyfield or the reference is more correct.

Your approach does look sensible! And the very close agreement with the other library is encouraging — and perhaps close enough for any practical application? I have never tried to apply a correction for atmospheric dispersion (is that what the angles are used for here?), so I don't know what error bars are permissible on the parallactic angle.

This evening I have had some time to sit back down with this problem, and try to figure out how to solve it in Skyfield itself. From which there is good and bad news!

• Good: I was able to perform the calculation two different ways, without adding any new features to Skyfield, and getting exactly the same answer each time.
• Bad: it disagrees with your answer by 1.3 arcseconds, making the agreement worse!

One key was re-reading the code and realizing that Star is entirely the wrong way to create an ra/dec position for this kind of problem (which is what I had been trying previously), because a Skyfield Star object must be viewed from a barycenter-referenced position in the Solar System, because if it has any proper motion then Skyfield needs to know whether the observer is closer to the star than the barycenter — which will very slightly advance the star's apparent position along its direction of motion — or is farther from the star, which introduces a light-time delay that puts the star slightly backwards in the sky along the course of its proper motion.

So here's what I came up with, and added down at the bottom of your script:

    from skyfield import api as sky
    from skyfield.api import load, tau
    from skyfield.positionlib import position_of_radec
    from skyfield.toposlib import ITRSPosition

    dt = datetime(2013, 8, 25, 10, 21, 39, 499798, tzinfo=timezone.utc)
    ts = load.timescale()
    t = ts.from_datetime(dt)

    antenna_icrf = [-1.60108719e+06, -5.04133984e+06,  3.55581586e+06]
    antenna_pos = ITRSPosition(sky.Distance(m=antenna_icrf))
    antenna_geo = sky.wgs84.geographic_position_of(antenna_pos.at(t))

    # Approach 1: figure out the altitude and azimuth of the field
    # centre and the Celestial Pole, and compute the bearing from the
    # centre to the pole.

    field_centre = position_of_radec(
        ra_hours=FIELD_CENTRE[0] / 15.0, dec_degrees=FIELD_CENTRE[1],
        center=antenna_geo, t=t,
    )
    north_pole = position_of_radec(
        ra_hours=0.0, dec_degrees=90.0, epoch=t,
        center=antenna_geo, t=t,
    )
    pa = skytrig.position_angle_of(field_centre.altaz(), north_pole.altaz())

    print(pa.radians)
    print('Difference (arcsec):', abs(pa.radians + REFERENCE_PA[0][0])
          / tau * 360 * 60 * 60)

    # Approach 2: figure out the equinox-of-date right ascension and
    # declination of the field centre and zenith, and compute the
    # bearing from the centre to the zenith.

    zenith = antenna_geo.at(t).from_altaz(alt_degrees=90.0, az_degrees=0.0)
    pa = skytrig.position_angle_of(
        field_centre.radec(epoch=t),
        zenith.radec(epoch=t),
    )
    pa_radians = pa.radians - tau  # since REFERENCE_PA is a negative angle
    print(pa_radians)
    print('difference (arcsec):', abs(pa_radians - REFERENCE_PA[0][0])
          / tau * 360 * 60 * 60)

I still don't love the awkward ITRSPosition-then-geographic_position_of maneuver, but I'll have to tackle that some other evening.

The output:

1.4552746776107097
Difference (arcsec): 1.3474020113110952
-1.4552746776107135
difference (arcsec): 1.3474020105324953

Maybe in a few days I'll sit down with this again, and try to puzzle out why this approach doesn't agree as closely as yours.

Oh, and, terminology: I'm thinking that maybe what I called position_angle_of() is really simply the universal concept of the "bearing" from one position to another. That's what you'd call it on a Mercator map. So maybe I should rename that routine bearing_to or something, and then if you do it relative to the celestial pole it becomes the "position angle", or relative to the zenith it becomes "parallactic angle".

Anyway, goodnight!

[brandon-rhodes]

Waking this morning and looking back over the code, I realize why Star is needed for the field center after all: because its position is shifted by the aberration of light, so we need to know the speed of the observation platform through space (relative to the Solar System Barycenter), and so we need to bring earth as a planet in anyway — it's not an extra complexity that can be ignored, as I had thought.

from skyfield import api as sky
from skyfield.api import load, tau
from skyfield.positionlib import position_of_radec
from skyfield.toposlib import ITRSPosition

dt = datetime(2013, 8, 25, 10, 21, 39, 499798, tzinfo=timezone.utc)
ts = load.timescale()

# from skyfield.data import iers
# url = load.build_url('finals2000A.all')
# with load.open(url) as f:
#     finals_data = iers.parse_x_y_dut1_from_finals_all(f)
# iers.install_polar_motion_table(ts, finals_data)

t = ts.from_datetime(dt)

antenna_icrf = [-1.60108719e+06, -5.04133984e+06,  3.55581586e+06]
antenna_pos = ITRSPosition(sky.Distance(m=antenna_icrf))
antenna_geo = sky.wgs84.geographic_position_of(antenna_pos.at(t))

planets = sky.load('de421.bsp')
earth = planets['earth']

star = sky.Star(
    ra=sky.Angle(degrees=FIELD_CENTRE[0]),
    dec=sky.Angle(degrees=FIELD_CENTRE[1])
)
field_centre = (earth + antenna_geo).at(t).observe(star).apparent()

# Approach 1: figure out the altitude and azimuth of the field
# centre and the Celestial Pole, and compute the bearing from the
# centre to the pole.

north_pole = position_of_radec(
    ra_hours=0.0, dec_degrees=90.0, epoch=t,
    center=antenna_geo, t=t,
)
pa = skytrig.position_angle_of(field_centre.altaz(), north_pole.altaz())

print(pa.radians)
print('Difference (arcsec):', abs(pa.radians + REFERENCE_PA[0][0])
      / tau * 360 * 60 * 60)

# Approach 2: figure out the equinox-of-date right ascension and
# declination of the field centre and zenith, and compute the
# bearing from the centre to the zenith.

zenith = antenna_geo.at(t).from_altaz(alt_degrees=90.0, az_degrees=0.0)
pa = skytrig.position_angle_of(
    field_centre.radec(epoch=t),
    zenith.radec(epoch=t),
)
pa_radians = pa.radians - tau  # since REFERENCE_PA is a negative angle
print(pa_radians)
print('difference (arcsec):', abs(pa_radians - REFERENCE_PA[0][0])
      / tau * 360 * 60 * 60)

The result agrees to around 0.4 arcseconds with your other tool, which is a happier result than last night's.

1.4552792249858868
Difference (arcsec): 0.4094385514783524
-1.4552792249858904
difference (arcsec): 0.4094385507455525

[JSKenyon]

Hi @brandon-rhodes and apologies for the long silence - a combination of vacation and other urgent work. Thank you for all your help. I will try out the various approaches.

Your approach does look sensible! And the very close agreement with the other library is encouraging — and perhaps close enough for any practical application? I have never tried to apply a correction for atmospheric dispersion (is that what the angles are used for here?), so I don't know what error bars are permissible on the parallactic angle.

The particular use-case is correcting for the apparent rotation of the sky relative to an interferometer (which manifests as errors in polarisation). The effect of ~0.12 arcsecond differences may well be below the noise. I think that I am happy to move over to skyfield, I was just hoping that we could pinpoint the source of the discrepancy so that I could document it e.g. skyfield includes the effect of X, unlike package Y.