pa_header.py

from astropy.io import fits
from astropy.time import Time
from astropy.coordinates import SkyCoord, EarthLocation
import astropy.units as u


subaru_loc = EarthLocation(lat=19.825504 * u.deg, lon=-155.4760187 * u.deg)

def get_pa_from_header(filename):
    with fits.open(filename) as hdul:
        ## Step 1: get exposure time
        date = hdul[0].header["DATE"]
        time = Time(date, format="isot", scale="ut1", location=subaru_loc)

        ## Step 2: Parse coordinate
        # (does not do any proper motion correction!)
        coord = SkyCoord(
            ra=hdul[1].header["RA"],
            dec=hdul[1].header["DEC"],
            unit=(u.hourangle, u.deg),
            frame="ICRS",
            equinox="J2000",
            obstime=time,
            location=subaru_loc
        )

        ## Step 3: Get local hour angle
        ha = time.sidereal_time("apparent").hourangle - coord.ra.hourangle

        ## Step 4: Calculate PA from HA/DEC
        _ha = ha * np.pi / 12  # hour angle to radian
        sin_ha, cos_ha = np.sin(_ha), np.cos(_ha)
        sin_dec, cos_dec = np.sin(coord.dec.rad), np.cos(coord.dec.rad)
        pa = np.arctan2(sin_ha, np.tan(subaru_loc.lat.rad) * cos_dec - sin_dec * cos_ha)
        return np.rad2deg(pa)