This computes ASPQ1 for a subsample of the 1p7 and 1p8 catalogs, using the offset tables in 1p7 and 1p8. These computed values are compared to the ASPQ1 values in 1p7.

If the algorithm were the same as the algorithm used in previous versions, then the computed ASPQ1 values using the 1p7 catalog and the 1p7 offsets would match the ASPQ1 values in 1p7. They do not match, which means the algorithm is not exactly the same. The correlation is decent, but not exact.

Using the 1p8 catalog to compute ASPQ1, but still using the 1p7 offsets, just introduces some spread.

Using the 1p7 catalog to compute ASPQ1, while using the 1p8 offsets, gives a very different result, since the shape of the offset distribution as a function of $(\delta r, \delta mag)$ differs from 1p7 to 1p8.

On average, the 1p8 offsets are smaller than the 1p7 ones

import sys
sys.path.insert(0, '..')

import os
import tables
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

from pathlib import Path
from scipy.stats import pearsonr
from astropy.table import Table, hstack
from astropy import units as u
from astropy.coordinates import SkyCoord

from agasc_gaia import agasc_update

# set the sky coordinates
def set_coord(agasc):
    a_ra = np.array(agasc['RA'])
    a_dec = np.array(agasc['DEC'])
    a_pm = (agasc['PM_DEC'] != -9999) & (agasc['PM_RA'] != -9999)
    a_ra[a_pm] += agasc['PM_RA'][a_pm] * (2023 - agasc['EPOCH'][a_pm]) / 1000 / 3600 / np.cos(np.deg2rad(agasc['DEC'][a_pm]))
    a_dec[a_pm] += agasc['PM_DEC'][a_pm] * (2023 - agasc['EPOCH'][a_pm]) / 1000 / 3600
    agasc['coord'] = SkyCoord(ra=a_ra, dec=a_dec, unit="deg")

# this function will be used to plot offsets

def plot_offsets(
    spoiler,
    levels=None,
    title=None,
):
    if title is None:
        title = "Centroid offset (arcsec) with spoiler star"
    if levels is None:
        levels = np.array([0.0005, 0.005, 0.05, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 10.0, 12.5, 15.0])
    # plt.figure(figsize=(6, 6))
    contour = plt.contour(
        spoiler.d_mag_bins,
        spoiler.separation_bins,
        spoiler.offset_lookup_data,
        levels=levels,
        # norm="log",
    )
    plt.clabel(contour, inline=1, fontsize=10)
    plt.xlabel("Spoiler $\Delta$Mag")
    plt.ylabel("Spoiler radial offset (arcsec)")
    plt.title(title)

spoiler_1p8 = agasc_update.OffsetLookup('data/offset_lookup_1p8rc10.h5')
spoiler_1p7 = agasc_update.OffsetLookup('data/offset_lookup_1p7.fits')

SKA = Path(os.environ['SKA'])
with tables.open_file('data/agasc1p8rc10.h5') as h5:
    agasc_1p8 = Table(h5.root.data[:])

with tables.open_file(SKA/ 'data' / 'agasc' / 'agasc1p7.h5') as h5:
    agasc_1p7 = Table(h5.root.data[:])

set_coord(agasc_1p8)
set_coord(agasc_1p7)

# remove extra entry in 1p7
if np.count_nonzero(agasc_1p7['AGASC_ID'] == 154534513) > 1:
    print("Warning: duplicate AGASC_ID=154534513")
    agasc_1p7 = agasc_1p7[~((agasc_1p7['AGASC_ID'] == 154534513) & (agasc_1p7['MAG_CATID'] == 1))]
    assert np.count_nonzero(agasc_1p7['AGASC_ID'] == 154534513) == 1

Warning: duplicate AGASC_ID=154534513

assert len(np.unique(agasc_1p7["AGASC_ID"])) == len(agasc_1p7)
assert len(np.unique(agasc_1p8["AGASC_ID"])) == len(agasc_1p8)

agasc_1p7.sort('AGASC_ID')
agasc_1p8.sort('AGASC_ID')

assert np.all(agasc_1p7["AGASC_ID"] == agasc_1p8["AGASC_ID"])

idx = np.arange(len(agasc_1p7))[(np.abs(agasc_1p7["RA"]) < 20) & (np.abs(agasc_1p7["DEC"]) < 20)]
agasc_1p7 = agasc_1p7[idx]
agasc_1p8 = agasc_1p8[idx]

# this is a copy of the function from agasc_update.py, which does not take the spoiler as an argument
def get_aspq_1(agasc, spoiler):
    assert "coord" in agasc.colnames, "agasc must have 'coord' column"
    idxc, idxcatalog, d2d, d3d = agasc["coord"].search_around_sky(
        agasc["coord"], 80 * u.arcsec
    )
    pairs = hstack([agasc[idxc], agasc[idxcatalog]])
    pairs["idx1"] = idxc
    pairs["idx2"] = idxcatalog
    pairs["d2d"] = d2d.to(u.arcsec)
    pairs["d3d"] = d3d
    pairs["d_mag"] = pairs["MAG_ACA_2"] - pairs["MAG_ACA_1"]

    pairs = pairs[
        (pairs["AGASC_ID_1"] != pairs["AGASC_ID_2"]) & (np.abs(pairs["d_mag"]) < 10)
    ]
    # ASPQ1 is in units of 50 milli-arcsec
    pairs["ASPQ1_3"] = spoiler(pairs["d_mag"], pairs["d2d"]) / 0.05

    pairs = pairs.group_by("idx1")
    agasc_2 = pairs[["AGASC_ID_1"]][pairs.groups.indices[:-1]].copy()
    agasc_2["ASPQ1_3"] = pairs["ASPQ1_3"].groups.aggregate(np.max)
    agasc_2.rename_column("AGASC_ID_1", "AGASC_ID")
    agasc_2.add_index("AGASC_ID")

    aspq1 = np.zeros(len(agasc))
    sel = np.in1d(agasc["AGASC_ID"], agasc_2["AGASC_ID"])
    aspq1[sel] = agasc_2.loc[agasc["AGASC_ID"][sel]]["ASPQ1_3"]
    return aspq1  # , agasc_2

agasc_1p7["aspq1_1p7"] = get_aspq_1(agasc_1p7, spoiler_1p7)
agasc_1p7["aspq1_1p8"] = get_aspq_1(agasc_1p7, spoiler_1p8)
agasc_1p8["aspq1_1p7"] = get_aspq_1(agasc_1p8, spoiler_1p7)
agasc_1p8["aspq1_1p8"] = get_aspq_1(agasc_1p8, spoiler_1p8)

Plots

aspq1_bins = np.linspace(0, 400, 50)

plt.plot(
    agasc_1p7["ASPQ1"],
    agasc_1p7["aspq1_1p7"],
    ".",
    alpha=0.3
)
plt.xlabel("ASPQ1 1p7")
plt.ylabel("Calculated ASPQ1 (1p7/1p7)")
plt.title("Calculated ASPQ1 with 1p7 catalog and 1p7 offsets")

Text(0.5, 1.0, 'Calculated ASPQ1 with 1p7 catalog and 1p7 offsets')

plt.plot(
    agasc_1p7["ASPQ1"],
    agasc_1p8["aspq1_1p7"],
    ".",
    alpha=0.3
)
plt.xlabel("ASPQ1 1p7")
plt.ylabel("Calculated ASPQ1 (1p8/1p7)")
plt.title("Calculated ASPQ1 with 1p8 catalog and 1p7 offsets")

Text(0.5, 1.0, 'Calculated ASPQ1 with 1p8 catalog and 1p7 offsets')

plt.plot(
    agasc_1p7["aspq1_1p7"],
    agasc_1p7["aspq1_1p8"],
    "."
)
plt.xlabel("ASPQ1 1p7")
plt.ylabel("Calculated ASPQ1 (1p7/1p8)")
plt.title("Calculated ASPQ1 with 1p7 catalog and 1p8 offsets")

Text(0.5, 1.0, 'Calculated ASPQ1 with 1p7 catalog and 1p8 offsets')

sel = (agasc_1p7["ASPQ1"] > 0) & (agasc_1p7["aspq1_1p7"] > 0)
vals = np.histogram2d(
    agasc_1p7["ASPQ1"][sel],
    agasc_1p7["aspq1_1p7"][sel],
    bins=[aspq1_bins, aspq1_bins]
)
vals, bins1, bins2 = vals
log_vals = np.log(vals+1)
x, y = np.meshgrid(bins1, bins2)

correlation = pearsonr(agasc_1p7["ASPQ1"][sel], agasc_1p7["aspq1_1p7"][sel])
correlation

PearsonRResult(statistic=0.9852674001118998, pvalue=0.0)

norm = np.where(np.ones_like(vals) > np.sum(vals, axis=0, keepdims=True), np.ones_like(vals), np.sum(vals, axis=0, keepdims=True))
plt.pcolor(x, y, vals.T / norm)

<matplotlib.collections.PolyQuadMesh at 0x137978c90>

plt.pcolor(x, y, log_vals.T)

<matplotlib.collections.PolyQuadMesh at 0x138bac2d0>

aspq1_bins = np.linspace(0, 400, 10)
agasc_1p7['spq1_bin'] = np.digitize(agasc_1p7['ASPQ1'], aspq1_bins)
agasc_1p8['spq1_bin'] = np.digitize(agasc_1p7['ASPQ1'], aspq1_bins)

for bin in range(len(aspq1_bins)):
    sel = (agasc_1p7['spq1_bin'] == bin)
    if np.count_nonzero(sel) == 0:
        continue
    m_spq1 = np.mean(agasc_1p7['ASPQ1'][sel])
    plt.hist(
        agasc_1p7['aspq1_1p7'][sel] - agasc_1p7['ASPQ1'][sel],
        bins=np.linspace(-100, 100, 100),
        histtype='step',
        density=True,
        label=f'{m_spq1:.1f}'
    )
# plt.yscale('log')
plt.legend()
plt.ylim(ymax=0.1)

(0.0, 0.1)

Plotting the Offset Tables

On average, the 1p8 offsets are smaller than the 1p7 ones, but this depends very much on where in the ($\delta r$, $\delta mag$) plane one stands

fig, axes = plt.subplots(1, 2, figsize=(12, 6), sharey=True, sharex=True)
plt.sca(axes[0])
plot_offsets(spoiler_1p7, title="1p7")
plt.sca(axes[1])
plot_offsets(spoiler_1p8, title="1p8")

plt.suptitle("Centroid offset (arcsec) with spoiler star")
plt.xlim((-6, 6))
plt.ylim((0, 50))

(0.0, 50.0)

Note the the 1p8 offsets are defined for negative angular separation, which is not physical and will never happen. This should have no effect, but can be fixed.

fig, axes = plt.subplots(1, 2, figsize=(12, 6), sharey=True, sharex=True)
plt.sca(axes[0])
plot_offsets(spoiler_1p7, title="1p7")
plt.sca(axes[1])
plot_offsets(spoiler_1p8, title="1p8")

plt.suptitle("Centroid offset (arcsec) with spoiler star")

Text(0.5, 0.98, 'Centroid offset (arcsec) with spoiler star')