<jupyter_start><jupyter_text>I study the apparent proper motion of radio sources to estimate the ICRF axis stability.<jupyter_code>import matplotlib.pyplot as plt
# from matplotlib.ticker import MultipleLocator
import numpy as np
import statsmodels.api as sm
from astropy.time import Time
from astropy.table import Table, join

from statsmodels.iolib.table import (SimpleTable, default_txt_fmt)
import seaborn as sns

from my_progs.stat_func.rms_calc import rms_calc
from my_progs.catalog.ga_func import glide_apex_calc
icrf3_def = Table.read("../data/icrf3sx-def-sou.txt", format="ascii")<jupyter_output><empty_output><jupyter_text>1 Distribution of the APM<jupyter_code>apm_tab = Table.read("../data/ts_nju_pm_fit_3sigma-10step.dat", format="ascii.csv")

# convert mas/yr into muas/yr
apm_tab["pmra"] = apm_tab["pmra"] * 1e3
apm_tab["pmra_err"] = apm_tab["pmra_err"] * 1e3
apm_tab["pmdec"] = apm_tab["pmdec"] * 1e3
apm_tab["pmdec_err"] = apm_tab["pmdec_err"] * 1e3<jupyter_output><empty_output><jupyter_text>Remove sources without apparant proper motion estimate.<jupyter_code>mask = apm_tab["num_cln"] >= 5
apm_tab = apm_tab[mask]<jupyter_output><empty_output><jupyter_text>Estimate the empirical cumulative probability function.<jupyter_code>ecdf_pmra = sm.distributions.ECDF(apm_tab["pmra"])
ecdf_pmdec = sm.distributions.ECDF(apm_tab["pmdec"])<jupyter_output><empty_output><jupyter_text>For the ICRF3 defining sources.<jupyter_code>apm_def = join(icrf3_def, apm_tab, keys="iers_name")
ecdf_pmra_def = sm.distributions.ECDF(apm_def["pmra"])
ecdf_pmdec_def = sm.distributions.ECDF(apm_def["pmdec"])<jupyter_output><empty_output><jupyter_text>Make the plot.<jupyter_code>bin_array = np.arange(-100, 100, 5)
fig, ax1 = plt.subplots()

color = "grey"
ax1.hist(apm_tab["pmra"],
         bins=bin_array,
         color=color,
         fill=False,
         label="All")
ax1.hist(apm_def["pmra"],
         bins=bin_array,
         color=color,
         label="Defining")
ax1.set_xlabel("$\\mu_{\\alpha*}$ (mas/yr)", fontsize=15)
ax1.set_ylabel("Nb Sources in  bins", color=color, fontsize=15)
ax1.tick_params(axis="y", labelcolor=color)
ax1.axis([-100, 100, 0, 400])
ax1.legend(loc="upper left")

ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis

color = "tab:red"
ax2.plot(ecdf_pmra.x,
         ecdf_pmra.y * len(apm_tab),
         color=color,
         lw=1,
         ls="--",
         label="All")
ax2.step(ecdf_pmra_def.x,
         ecdf_pmra_def.y * len(apm_def),
         color=color,
         lw=1,
         ls="-",
         label="Defining")
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_ylabel("Cumulative Nb sources", color=color,
               fontsize=15)  # we already handled the x-label with ax1
ax2.axis([-100, 100, 0, 3500])
ax2.legend(loc="center right")

plt.tight_layout()
plt.savefig("../plots/pmra-hist.eps")
per = [0, 5, 10, 25, 50, 75, 90, 95, 100] 
Q = np.percentile(apm_tab["pmra"], per)

line_fmt1 = "  {:8d}" * len(per)
line_fmt2 = "  {:8.2f}" * len(Q)

print("PM in R.A. (muas/yr)")
print("Percentile(%)", line_fmt1.format(*per))
print("             ", line_fmt2.format(*Q))

lim = 30
mask = (np.fabs(apm_tab["pmra"]) < lim)
print("[-{:2d}, {:2d}]: {:.2f}%".format(lim, lim, len(apm_tab[mask])/len(apm_tab)*100))
fig, ax1 = plt.subplots()

color = "grey"
ax1.hist(apm_tab["pmdec"],
         bins=bin_array,
         color=color,
         fill=False,
         label="All")
ax1.hist(apm_def["pmdec"],
         bins=bin_array,
         color=color,
         label="Defining")
ax1.set_xlabel("$\\mu_{\\delta}$ (mas/yr)", fontsize=15)
ax1.set_ylabel("Nb Sources in  bins", color=color, fontsize=15)
ax1.tick_params(axis="y", labelcolor=color)
ax1.axis([-100, 100, 0, 400])
ax1.legend(loc="upper left")

ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis

color = "tab:red"
ax2.step(ecdf_pmdec.x,
         ecdf_pmdec.y * len(apm_tab),
         color=color,
         lw=1,
         ls="--",
         label="All")
ax2.step(ecdf_pmdec_def.x,
         ecdf_pmdec_def.y * len(apm_def),
         color=color,
         lw=1,
         ls="-",
         label="Defining")
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_ylabel("Cumulative Nb sources", color=color,
               fontsize=15)  # we already handled the x-label with ax1
ax2.axis([-100, 100, 0, 3500])
ax2.legend(loc="center right")

plt.tight_layout()
plt.savefig("../plots/pmdec-hist.eps")
per = [0, 5, 10, 25, 50, 75, 90, 95, 100] 
Q = np.percentile(apm_tab["pmdec"], per)

line_fmt1 = "  {:8d}" * len(per)
line_fmt2 = "  {:8.2f}" * len(Q)

print("PM in decl. (muas/yr)")
print("Percentile(%)", line_fmt1.format(*per))
print("             ", line_fmt2.format(*Q))

lim = 30
mask = (np.fabs(apm_tab["pmdec"]) < lim)
print("[-{:2d}, {:2d}]: {:.2f}%".format(lim, lim, len(apm_tab[mask])/len(apm_tab)*100))
mean_pmra = np.mean(apm_tab["pmra"])
median_pmra = np.median(apm_tab["pmra"])
std_pmra = np.std(apm_tab["pmra"])

mean_pmdec = np.mean(apm_tab["pmdec"])
median_pmdec = np.median(apm_tab["pmdec"])
std_pmdec = np.std(apm_tab["pmdec"])
mean_pmra_def = np.mean(apm_def["pmra"])
median_pmra_def = np.median(apm_def["pmra"])
std_pmra_def = np.std(apm_def["pmra"])

mean_pmdec_def = np.mean(apm_def["pmdec"])
median_pmdec_def = np.median(apm_def["pmdec"])
std_pmdec_def = np.std(apm_def["pmdec"])
se = np.vstack([[len(apm_tab), mean_pmra, median_pmra, std_pmra],
                [len(apm_tab), mean_pmdec, median_pmdec, std_pmdec],
                [len(apm_def), mean_pmra_def, median_pmra_def, std_pmra_def],
                [len(apm_def), mean_pmdec_def, median_pmdec_def, std_pmdec_def]])

se = np.round(se, 2)

colnames = ["#sou", "Mean", "Median", "Std"]
rownames = ["All-RA", "All-Dec", "Def-RA", "Def-Dec"]

tabl = SimpleTable(se, colnames, rownames, txt_fmt=default_txt_fmt)
print(tabl)
print("Unit: muas/yr")<jupyter_output>==================================
         #sou   Mean Median  Std  
----------------------------------
All-RA  2831.0 -5.86  -0.48 256.23
All-Dec 2831.0 -6.47   0.15 467.04
Def-RA   290.0  0.81  -0.52   19.1
Def-Dec  290.0 -3.55   -0.4  31.94
----------------------------------
Unit: muas/yr<jupyter_text>2. Global spin from the ICRF3-defining source subset<jupyter_code>from tool_func import vsh_fit_for_pm, vsh_fit_for_pm2
pmt, sig, output = vsh_fit_for_pm(apm_def)
pmt1, sig1, output1 = vsh_fit_for_pm2(apm_def)
se = np.vstack((pmt[:4], sig[:4], pmt1, sig1)).T
se = np.round(se, 2)

colnames = ["Full", "+/-", "Rotation", "+/-"]
rownames = ["Rx", "Ry", "Rz", "R"]

tabl = SimpleTable(se, colnames, rownames, txt_fmt=default_txt_fmt)
print(tabl)

se = np.vstack((pmt[4:], sig[4:])).T
se = np.round(se, 2)

colnames = ["Glide", "+/-"]
rownames = ["Gx", "Gy", "Gz", "G"]

tabl = SimpleTable(se, colnames, rownames, txt_fmt=default_txt_fmt)
print(tabl)

print("Unit: muas/yr")
print("Number of sources used in the fitting is", len(apm_def))
print("Apex for rotation (deg)")
print("      Full      Rotation")
print("RA : {:3.0f} +/- {:3.0f}  {:3.0f} +/- {:3.0f}".format(
    output["R_ra"], output["R_ra_err"], output1["R_ra"], output1["R_ra_err"]))
print("Dec: {:3.0f} +/- {:3.0f}  {:3.0f} +/- {:3.0f}".format(
    output["R_dec"], output["R_dec_err"], output1["R_dec"],
    output1["R_dec_err"]))

print("Apex for glide (deg)")
print("      Glide")
print("RA : {:3.0f} +/- {:3.0f}  ".format(
    output["G_ra"], output["G_ra_err"]))
print("Dec: {:3.0f} +/- {:3.0f}".format(
    output["G_dec"], output["G_dec_err"]))<jupyter_output>Apex for rotation (deg)
      Full      Rotation
RA : 183 +/-  28  107 +/-  27
Dec: -22 +/-  18  -32 +/-  25
Apex for glide (deg)
      Glide
RA : 307 +/-  12  
Dec: -35 +/-  11
{'fraction_non_alphanumeric': 0.19406095337327428, 'fraction_numerical': 0.04024485543110185, 'mean_word_length': 3.9670116429495472, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '310622', 'license': 'MIT', 'n_tokens_mistral': 3424, 'n_tokens_neox': 3175, 'n_words': 576, 'path': 'progs/.ipynb_checkpoints/apparent-proper-motion-analysis-checkpoint.ipynb', 'repo_name': 'Niu-Liu/icrf-axis-stability', 'size': '7678'}