Parallactic angle calculation and plot in Julia

parallactic.jl

using AstroLib, AstroAngles, Dates

"""
Given a source position in RA and DEC, a time in Julian days, and an observatory,
return the parallactic angle of the target.

You can also pass optional keyword arguments that are forwarded to `eq2hor` to
adjust assumptions around refraction correction, temperature, nutation, etc.
"""
function postime2pa(ra_deg, dec_deg, time_jd, observatory; kwargs...)
    obs = AstroLib.observatories["keck"]
    az_east_from_north, ha_deg = eq2hor(ra_deg, dec_deg, time_jd, observatory; kwargs...);
    # pa = -atand(
    #     -sind(ha_deg),
    #     cosd(dec_deg)*tand(obs.latitude) - sind(dec_deg)*cosd(ha_deg)
    # )
    pa = asind(
        sind(az_east_from_north) * cosd(obs.latitude) / cosd(dec_deg)
    )
    return pa
end


# Calculate the PA at one moment
pa = postime2pa(
    hms"03 42 31.80",
    dms"12 16 22.5",
    jdcnv(now(Dates.UTC)),
    "keck"
)

# Calculate the PA over time and plot with Makie

using Plots

ts = range(
    start=DateTime(2023,2,8,18,0,0,0),
    step=Minute(10),
    stop=DateTime(2023,2,9,7,0,0,0),
)

pas = postime2pa.(
    hms"03 42 31.80",
    dms"12 16 22.5",
    jdcnv.(ts),
    "keck"
)


plot(ts, pas, xlabel="time", ylabel="PA - deg")