Made by Yoonsoo P. Bach
The radiance factor is derivable from the shape model (i.e., the absolute value of the projected cross-sectional area as a function of time). The original code was written in IDL
using POLYSHADE
, developed by one of our co-authors (Sunho Jin). Here, I will just utilize the projected area produced by that code, and show you how we found the 1-sigma models and draw the RADF figure.
We did a cross-check of the ranges of projected area values fron this IDL
code with ParaView
values (see 04_cross_section
notebook).
%load_ext version_information
import time
now = time.strftime("%Y-%m-%d %H:%M:%S (%Z = GMT%z)")
print(f"This notebook was generated at {now} ")
vv = %version_information scipy, numpy, matplotlib, pandas, numba, version_information
for i, pkg in enumerate(vv.packages):
print(f"{i} {pkg[0]:10s} {pkg[1]:s}")
This notebook was generated at 2019-04-08 11:39:27 (KST = GMT+0900) 0 Python 3.6.8 64bit [GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)] 1 IPython 6.5.0 2 OS Darwin 18.2.0 x86_64 i386 64bit 3 scipy 1.2.1 4 numpy 1.16.2 5 matplotlib 3.0.3 6 pandas 0.24.2 7 numba 0.42.0 8 version_information 1.0.3
As described in the appendix of the original publication, we multiplied a linear factor to the total projected area, and matched it with the observed magnitude (converted magnitude to a linear scale, since the error in magnitude is not Gaussian and thus the chi-square minimization process should not be used).
First, import and definitions:
from pathlib import Path
import numpy as np
import pandas as pd
from scipy.interpolate import UnivariateSpline as usp
from scipy.optimize import curve_fit
from matplotlib import pyplot as plt
from matplotlib import rcParams
from astropy.time import Time
from astropy import units as u
plt.style.use('default')
rcParams.update({'font.size':13})
V_SUN = -26.76
D2R = np.pi / 180
DATAPATH = Path("data")
def arr(x): return np.array(x)
def linscale(x, a):
return a * x
def iauHG(alpha, Gpar, degree=True):
""" IAU H-G system. See ``iauHa2H0``.
"""
if degree:
alpha = np.deg2rad(alpha)
def exp(x):
return np.exp(x)
tanp2 = np.tan(alpha / 2) # tangent per 2 value
term1 = Gpar * exp(-1.87 * tanp2**1.22)
term2 = (1 - Gpar) * exp(-3.33 * tanp2**0.63)
return term1 + term2
def iauHa2H0(Halpha, alpha, Gpar, degree=True):
''' Converts H(alpha) = V(1, 1, alpha) to H(0) using IAU HG model.
Parameters
----------
Halpha: float
The reduced magnitude in mag unit.
alpha: float
The phase angle.
Gpar: float
The slope parameter. Usually 0.1 to 0.3.
degree: bool
Whether ``alpha`` is in degrees.
'''
H = Halpha + 2.5 * np.log10(iauHG(alpha=alpha, Gpar=Gpar, degree=degree))
return H
def Ha2RADF(Halpha, area):
''' Converts H(alpha) = V(1, 1, alpha) to I/F (radiacne factor RADF)
Parameters
Halpha: float
The reduced magnitude in mag unit.
area: float or Quantity
The cross sectional area at the observed epoch in square meter or given
as a `Quantity`.
'''
if isinstance(area, u.Quantity):
area = area.to(u.m^2)
m_c = 5 * np.log10(1.4960e11)
term1 = Halpha - V_SUN - m_c - 2.5 * np.log10(np.pi / area)
RADF = 10**(-term1 / 2.5)
return RADF
def powf(x, a, b):
return a * 10**(-b * x)
def get_minmax(a, axis):
maximum = np.nanmax(arr(a), axis=axis)
minimum = np.nanmin(arr(a), axis=axis)
return minimum, maximum
Then run the code to save a 4-column csv file, containing the two offsets, the linear factor, and the reduced chi-square statistic:
photall = pd.read_csv(DATAPATH / "phot_targ_4latex_new.csv")
photall = photall.sort_values(by="datetime_jd")
area_new = pd.read_csv(DATAPATH / "toutatis_new.csv")
dg_area = area_new.groupby(["off1", "off2"])
Gpars = [0.0, 0.1, 0.2, 0.5]
starts = {"all":0, "no1":1}
for usedata, i_init in starts.items():
phot = photall.copy()[i_init:]
jd_obs = phot["datetime_jd"]
dof = len(jd_obs) - 1
dm = phot["merr_red"]
for Gpar in Gpars:
H0 = arr(iauHa2H0(phot["m_red"], phot["alpha"], Gpar=Gpar, degree=True))
I0 = arr(10**(- H0 / 2.5)) # Convert to linear-scale
dI0 = arr(I0 * dm * np.log(10) / 2.5)
chisqdf_path = DATAPATH / f"Chisq_Gpar_{1000*Gpar:04.0f}_{usedata}.csv"
if not chisqdf_path.exists():
chisqsel = dict(off1=[], off2=[], factor=[], chisqred=[])
for name, group in dg_area:
off1, off2 = name
area_orig = group["area"]
ff = usp(group["jd"], area_orig)
area = ff(jd_obs)
popt, pcov = curve_fit(linscale, xdata=area, ydata=I0,
sigma=dI0, absolute_sigma=True)
chisqred = np.sum( ((I0 - linscale(area, *popt)) / dI0)**2 ) / dof
chisqsel["off1"].append(off1)
chisqsel["off2"].append(off2)
chisqsel["factor"].append(*popt)
chisqsel["chisqred"].append(chisqred)
chisqdf = pd.DataFrame.from_dict(chisqsel)
chisqdf = chisqdf.sort_values("chisqred")
chisqdf.to_csv(chisqdf_path, index=False)
Then the following code will generate the 1-sigma models as well as the photometry results obtained in data_reduction
:
# Gpars = [0.0, 0.10, 0.2, 0.5]
# sigmas = [1, 2]
Gpars = [0.10]
sigmas = [1]
for usedata, i_init in starts.items():
phot = photall.copy()[i_init:]
Ha = phot["m_red"]
jd_obs = phot["datetime_jd"]
dm = phot["merr_red"]
Ha = phot["m_red"]
for Gpar in Gpars:
for sigma in sigmas:
area_new = pd.read_csv(DATAPATH / "toutatis_new.csv")
dg_area = area_new.groupby(["off1", "off2"])
t_ref = Time('2018-04-01', format='isot').jd - 1
dof = len(phot) - 1
t_val = arr(jd_obs - t_ref)
H0 = arr(iauHa2H0(phot["m_red"], phot["alpha"], Gpar=Gpar, degree=True))
I0 = arr(10**(-H0/2.5)) # Convert to linear-scale
dI0 = arr(I0 * dm * np.log(10) / 2.5)
chisqdf_path = DATAPATH / f"Chisq_Gpar_{1000*Gpar:04.0f}_{usedata:s}.csv"
chisqdf = pd.read_csv(chisqdf_path)
# 1-sigma models
minchisq = chisqdf["chisqred"].min()
models_1sigma = chisqdf[chisqdf["chisqred"] < minchisq + sigma**2]
fig, axs = plt.subplots(2, 1, figsize=(6, 8))
for i, row in models_1sigma.iterrows():
areas = dg_area.get_group((row["off1"], row["off2"]))
modelmag = -2.5 * np.log10(arr(row["factor"] * areas["area"]))
xx = areas["jd"] - t_ref
xx2 = np.linspace(xx.min(), xx.max(), 100)
ff = usp(xx, modelmag, s=0)
if row["chisqred"] == minchisq:
axs[0].plot(xx2, ff(xx2), 'r-', zorder=999, label="Best-fit model")
axs[1].errorbar(jd_obs - t_ref, H0 - ff(jd_obs - t_ref), yerr=dm,
marker='x', ls='', elinewidth=0.5, capsize=3, color='r')
else:
axs[0].plot(xx2, ff(xx2), 'k:', lw=1, alpha=0.1)
# Just for legend
axs[0].plot([None], [None], 'k:', lw=1, alpha=1,
label=r"Models within {:d}-$\sigma$ CI".format(sigma))
# Obseved data (reduced mag)
axs[0].plot(t_val, Ha, 'b+', label=r"$H_\mathrm{V}(\alpha)$")
# Calculated H-mag
axs[0].errorbar(t_val, H0, dm, marker='x', ls='', elinewidth=0.5,
capsize=3, zorder=999, color='r',
label=r"$H_\mathrm{V}(0)$")
# Put the texts indicating phase angles
for i, a in enumerate(phot["alpha"]):
if i == 0:
axs[0].text(t_val[i]-0.1, 14.65, #H0[i] - dm[i] - 0.3,
r'$\alpha = {:.2f}^{{\circ}}$'.format(a), rotation=90)
else:
axs[0].text(t_val[i]-0.1, 14.65, #H0[i] - dm[i] - 0.3,
r'${:.2f}^{{\circ}}$'.format(a), rotation=90)
axs[0].set_ylim(16.5, 14.5)
axs[1].set_ylim(-0.5, 0.5)
axs[1].axhline(0, color='r', ls='-')
axs[0].set_title(f"G = {Gpar}")
axs[0].legend(loc=8)
axs[0].set_ylabel("[mag]")
axs[1].set_ylabel("difference [mag]")
for ax in [axs[0], axs[1]]:
ax.set_xlabel("Date of April, 2018")
plt.tight_layout()
# plt.savefig(f"fig_Halpha_Gpar_{1000*Gpar:04.0f}_{sigma:d}_sigma_{usedata:s}.pdf")
And the following code will generate the RADF plot as in the original publication.
numba
, see below.First, define numba
functions:
import numba as nb
@nb.njit(parallel=True, fastmath=True)
def brute_force_powf(x, y, yerr, arr_a, arr_pow, arr_chisq):
for i in nb.prange(arr_a.shape[0]):
for j in nb.prange(arr_pow.shape[0]):
a = arr_a[i]
b = arr_pow[j]
arr_chisq[i, j] = np.sum(( (y - a * 10**(-b * x)) / yerr )**2)
@nb.njit(parallel=True)
def analysis_powf(x, y, yerr, arr_amp, arr_pow, arr_chisqred, arr_vals,
arr_slopes, arr_sOEs, amax_sOE, amin_sOE):
da_sOE = amax_sOE - amin_sOE
dof = y.shape[0] - 2
xx = np.arange(0, 5.1, 0.1)
for i in nb.prange(arr_vals.shape[0]):
for j in nb.prange(arr_vals.shape[1]):
a = arr_amp[i]
b = arr_pow[j]
chisqred = np.sum(( (y - a * 10**(-b * x)) / yerr )**2) / dof
for k in nb.prange(arr_vals.shape[2]):
arr_vals[i, j, k] = a * 10**(-b * xx[k])
arr_chisqred[i, j] = chisqred
# Save slope and sOE
slope = 10**(-b * (0.3 - 5.0))
sOE = (a / da_sOE * (10**(-b * amin_sOE) - 10**(-b * amax_sOE)))
arr_slopes[i, j] = slope
arr_sOEs[i, j] = sOE
photall = pd.read_csv(DATAPATH / "phot_targ_4latex_new.csv")
photall = photall.sort_values(by="datetime_jd")
Gpars = [0.1]
starts = {"all":0, "no1":1}
sigmas = [1]
amax_sOE = 2.5
amin_sOE = 0.7
da_sOE = amax_sOE - amin_sOE
aa = np.arange(0, 5.1, 0.1)
title_str = (r"Power-law ($ I/F = p_\mathrm{{V}} \cdot 10^{{-b \alpha}}$) fit"
+ "\n"
+ r"$p_\mathrm{{V}} = {:.3f}_{{-{:.3f}}}^{{+{:.3f}}}$ ; "
+ r"$ f_1 $ "
+ r"$ = {:.3f}_{{-{:.3f}}}^{{+{:.3f}}}$ ; "
+ r"$ f_2 $ "
+ r"$ = {:.3f}_{{-{:.3f}}}^{{+{:.3f}}}$")
for usedata, i_init in starts.items():
phot = photall.copy()[i_init:]
alpha = arr(phot["alpha"])
Ha = arr(phot["m_red"])
jd_obs = arr(phot["datetime_jd"])
dm = arr(phot["merr_red"])
Ha = arr(phot["m_red"])
dof = len(alpha) - 2
for Gpar in Gpars:
# The best-fit shape model
area_new = pd.read_csv(DATAPATH / "toutatis_new.csv")
dg_area = area_new.groupby(["off1", "off2"])
chisqdf = pd.read_csv(DATAPATH / f"Chisq_Gpar_{Gpar:04.0f}_{usedata:s}.csv")
bestshapemodel = dg_area.get_group((chisqdf["off1"][0], chisqdf["off2"][1]))
ff_area = usp(bestshapemodel["jd"], bestshapemodel["area"])
bestarea = ff_area(jd_obs)
RADF = Ha2RADF(Ha, bestarea)
RADF_err = RADF * np.log(10) / 2.5 * dm
popt, pcov = curve_fit(powf, alpha, RADF, sigma=RADF_err, absolute_sigma=True)
bestpV = popt[0]
bestslope = powf(0.3, *popt) / powf(5, *popt)
bestA5 = powf(5, *popt)
bestsOE = (powf(amin_sOE, *popt) - powf(amax_sOE, *popt)) / da_sOE
dpV, dpow = np.diag(np.sqrt(pcov))
chisqredmin = np.sum(((RADF - powf(alpha, *popt)) / RADF_err)**2) / dof
pVs = np.arange(bestpV - 5 * dpV, bestpV + 5 * dpV, 0.0001)
pows = np.arange(popt[1] - 5 * dpow, popt[1] + 5 * dpow, 0.0001)
chisqred = np.empty(shape=(len(pVs), len(pows)))
vals_ksig = np.empty(shape=(len(pVs), len(pows), len(aa))) # 3-d arr
slopes = np.empty_like(chisqred)
sOEs = np.empty_like(chisqred)
analysis_powf(x=alpha, y=RADF, yerr=RADF_err, arr_amp=pVs, arr_pow=pows,
arr_chisqred=chisqred, arr_vals=vals_ksig, arr_slopes=slopes,
arr_sOEs=sOEs, amax_sOE=amax_sOE, amin_sOE=amin_sOE)
for ksig in sigmas:
mask1sig = chisqred < (chisqredmin + ksig**2)
# Find the min/max of RADF at each alpha for 1-sigma models
RADF_ksig_min = []
RADF_ksig_max = []
for i in range(len(aa)):
RADF_ksig_max.append(np.nanmax(vals_ksig[:,:,i][mask1sig]))
RADF_ksig_min.append(np.nanmin(vals_ksig[:,:,i][mask1sig]))
#mask1sig = np.tile(chisqred < (chisqredmin + 1), len(aa))
#mask1sig = np.reshape(mask1sig, (len(pVs), len(pows), len(aa)))
#vals_ksig[~mask1sig] = np.nan
#RADF_min = np.nanmin(vals_ksig, axis=(0,1))
#RADF_max = np.nanmax(vals_ksig, axis=(0,1))
sl_ksig_min, sl_ksig_max = get_minmax(slopes[mask1sig], axis=None)
OE_ksig_min, OE_ksig_max = get_minmax(sOEs[mask1sig], axis=None)
dpV_neg = bestpV - RADF_ksig_min[0]
dpV_pos = RADF_ksig_max[0] - bestpV
dA5_neg = bestA5 - RADF_ksig_min[-1]
dA5_pos = RADF_ksig_max[-1] - bestA5
dsl_neg = bestslope - sl_ksig_min
dsl_pos = sl_ksig_max - bestslope
dOE_neg = bestsOE - OE_ksig_min
dOE_pos = OE_ksig_max - bestsOE
fig, axs = plt.subplots(2, 1, figsize=(7, 9))
axf, axc = axs[0], axs[1] # axis for "f"it and "c"hi-square
axf.set_title(title_str.format(bestpV, dpV_neg, dpV_pos,
bestslope, dsl_neg, dsl_pos,
bestsOE, dOE_neg, dOE_pos))
axf.errorbar(alpha, RADF, RADF_err, marker='x', ls='', elinewidth=0.5,
capsize=3, label="$ I/F $")
axf.fill_between(aa, RADF_ksig_min, RADF_ksig_max, color='k', alpha=0.2,
label=r"{:d}-$\sigma$ CI of model".format(ksig))
axf.plot(aa, powf(aa, *popt), color='r', label="Best-fit model")
axf.legend()
axf.errorbar(5, bestA5, marker='o', color='k', capsize=3, elinewidth=1,
yerr=[[dA5_neg], [dA5_pos]])
axf.text(4.4, bestA5 - 0.01, f"{bestA5:.3f}")
axf.text(4.4, bestA5 + dA5_pos - 0.01, f"{RADF_ksig_max[-1]:.3f}")
axf.text(4.4, bestA5 - dA5_neg, f"{RADF_ksig_min[-1]:.3f}")
axf.set_xlim(0, 5.05)
axf.set_ylim(0.06, 0.22)
axf.set_xlabel(r"Phase angle, $\alpha$ $[{^{\circ}}]$")
axf.set_ylabel(r"Radiance factor, $I/F$")
axf.grid(ls=':', which='both')
# xy order is differnt from numpy axis order!
axc.set_title(r"$\chi^2_\mathrm{red} $ map of the fitting (color in log scale)")
axc.set_xlabel("Power $b$")
axc.set_ylabel("Albedo $p_\mathrm{{V}}$")
c = axc.contourf(pows, pVs, np.log10(chisqred), levels=30)
cs = axc.contour(pows, pVs, chisqred, levels=[chisqredmin+1], colors='r')
c2 = axc.contour(pows, pVs, chisqred, levels=[5, 10, 20, 50, 100], colors='k')
axc.axhline(RADF_ksig_min[0], color='k', ls=':')
axc.axhline(RADF_ksig_max[0], color='k', ls=':')
#axc.axhline(bestpV, color='r', ls=':')
axc.clabel(cs, fmt=r"1-$\sigma$", inline=1, fontsize=15, colors='r')
axc.clabel(c2, fmt="$%.0f $", inline=1, fontsize=12, colors='k')
axc.errorbar(pows[len(pows)//2], pVs[len(pVs)//2],
yerr=[[dpV_neg], [dpV_pos]],
marker='x', ls='', capsize=0, elinewidth=2, color='r', ms=8,
label=r"$ \chi^2_\mathrm{{red, min}} = {:.3f} $".format(chisqredmin))
axc.plot(np.nan, np.nan, 'r-',
label=r"1-$\sigma$ contour ($\chi^2_\mathrm{{red}} = {:.3f}$)".format(chisqredmin+1))
axc.legend(loc=4, fontsize=12)
#plt.colorbar(c, label=r"$\log_{10} (\chi^2_\mathrm{red}) $")
plt.tight_layout()
plt.savefig(Path("figs") / f"fig_RADF_new_{usedata:s}_{ksig:d}_sigma.pdf")
A better approach might be using MC technique, but I believed the primitive brute-force grid search will suffice our goal accuracy for this work.
The following code will generate similar result without numba, but much slower. Due to the speed issue, I increased the grid size by factor of 5 and reduced the searching space from 5-sigma to 3-sigma.
def powf(x, a, b):
return a * 10**(-b * x)
def get_minmax(a, axis):
maximum = np.max(arr(a), axis=axis)
minimum = np.min(arr(a), axis=axis)
return minimum, maximum
amax_sOE = 2.5
amin_sOE = 0.7
da_sOE = amax_sOE - amin_sOE
title_str = (r"Power-law ($\mathrm{{RADF}} = p_\mathrm{{V}} \cdot 10^{{-b \alpha}}$) fit"
+ "\n"
+ r"$p_\mathrm{{V}} = {:.3f}_{{-{:.3f}}}^{{+{:.3f}}}$ ; "
+ r"$ f_1 $ "
+ r"$ = {:.3f}_{{-{:.3f}}}^{{+{:.3f}}}$ ; "
+ r"$ f_2 $ "
+ r"$ = {:.3f}_{{-{:.3f}}}^{{+{:.3f}}}$")
# The best-fit shape model
area_new = pd.read_csv("toutatis_new.csv")
dg_area = area_new.groupby(["off1", "off2"])
chisqdf = pd.read_csv(Path(f"Chisq_Gpar_0100.csv"))
bestshapemodel = dg_area.get_group((chisqdf["off1"][0], chisqdf["off2"][1]))
ff_area = usp(bestshapemodel["jd"], bestshapemodel["area"])
ksig = 1
phot = pd.read_csv("phot_targ_4latex_new.csv")
phot = phot.sort_values(by="datetime_jd")
alpha = phot["alpha"]
Ha = phot["m_red"]
jd_obs = phot["datetime_jd"]
dm = phot["merr_red"]
Ha = phot["m_red"]
dof = len(alpha) - 2
bestarea = ff_area(jd_obs)
RADF = Ha2RADF(Ha, bestarea)
RADF_err = RADF * np.log(10) / 2.5 * dm
popt, pcov = curve_fit(powf, alpha, RADF, sigma=RADF_err, absolute_sigma=True)
bestpV = popt[0]
bestslope = powf(0.3, *popt) / powf(5, *popt)
bestA5 = powf(5, *popt)
bestsOE = (powf(amin_sOE, *popt) - powf(amax_sOE, *popt)) / da_sOE
dpV, dpow = np.diag(np.sqrt(pcov))
aa = np.arange(0, 5.1, 0.1)
chisqredmin = np.sum(((RADF - powf(alpha, *popt)) / RADF_err)**2) / dof
pVs = np.arange(bestpV - 3 * dpV, bestpV + 3 * dpV, 0.0005)
pows = np.arange(popt[1] - 3 * dpow, popt[1] + 3 * dpow, 0.0005)
vals_ksig = []
slopes = []
sOEs = []
for p1 in pVs:
for p2 in pows:
chisqred = np.sum(((RADF - powf(alpha, p1, p2)) / RADF_err)**2) / dof
if chisqred < (chisqredmin + ksig**2):
slope = powf(0.3, p1, p2) / powf(5.0, p1, p2)
sOE = (powf(amin_sOE, p1, p2) - powf(amax_sOE, p1, p2)) / da_sOE
vals_ksig.append(powf(aa, p1, p2))
slopes.append(slope)
sOEs.append(sOE)
RADF_ksig_min, RADF_ksig_max = get_minmax(vals_ksig, axis=0)
sl_ksig_min, sl_ksig_max = get_minmax(slopes, axis=None)
OE_ksig_min, OE_ksig_max = get_minmax(sOEs, axis=None)
dpV_neg = bestpV - RADF_ksig_min[0]
dpV_pos = RADF_ksig_max[0] - bestpV
dA5_neg = bestA5 - RADF_ksig_min[-1]
dA5_pos = RADF_ksig_max[-1] - bestA5
dsl_neg = bestslope - sl_ksig_min
dsl_pos = sl_ksig_max - bestslope
dOE_neg = bestsOE - OE_ksig_min
dOE_pos = OE_ksig_max - bestsOE
fig, ax = plt.subplots(figsize=(6, 4.5))
ax.set_title(title_str.format(bestpV, dpV_neg, dpV_pos,
bestslope, dsl_neg, dsl_pos,
bestsOE, dOE_neg, dOE_pos))
ax.errorbar(alpha, RADF, RADF_err, marker='x', ls='', elinewidth=0.5,
capsize=3, label="RADF")
ax.fill_between(aa, RADF_ksig_min, RADF_ksig_max, color='k', alpha=0.2,
label=r"{:d}-$\sigma$ CI of model".format(ksig))
ax.plot(aa, powf(aa, *popt), color='r', label="Best-fit model")
ax.legend()
ax.errorbar(5, bestA5, marker='o', color='k', capsize=3, elinewidth=1,
yerr=[[dA5_neg], [dA5_pos]])
ax.text(4.4, bestA5 - 0.01, f"{bestA5:.3f}")
ax.text(4.4, bestA5 + dA5_pos - 0.01, f"{RADF_ksig_max[-1]:.3f}")
ax.text(4.4, bestA5 - dA5_neg, f"{RADF_ksig_min[-1]:.3f}")
ax.set_xlim(0, 5.05)
ax.set_ylim(0.06, 0.22)
ax.set_xlabel(r"Phase angle, $\alpha$ $[{^{\circ}}]$")
ax.set_ylabel(r"Radiance factor, $I/F$")
ax.grid(ls=':', which='both')
plt.tight_layout()
plt.show()