In this notebook we assess the change in handling the weighted mask for Condition #2 of Figueria et al 2016.
For that paper the spectrum is broken into several small spectra on regions masked telluric lines (>2%).
The updated version just applies a boolen mask after pixel weights are calculated.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from eniric.atmosphere import Atmosphere
from eniric.legacy import mask_clumping, RVprec_calc_masked
from scripts.phoenix_precision import convolve_and_resample
from eniric.snr_normalization import snr_constant_band
from eniric.precision import pixel_weights, rv_precision
from eniric.utilities import band_limits, load_aces_spectrum, wav_selector
In [2]:
wav_, flux_ = load_aces_spectrum([3900, 4.5, 0, 0])
# Small section in K bands to experiment with
wav, flux = wav_selector(wav_, flux_, 2.1025, 2.1046)
In [3]:
# Telluric mask
atm_ = Atmosphere.from_band("K", bary=True)
atm = atm_.at(wav)
mask = atm.mask
In [4]:
# Clumping method
wclump, fclump = mask_clumping(wav, flux, mask)
print("# Number of clumps = ", len(wclump))
# print(wclump, fclump)
print(len(wclump))
In [5]:
wis_0 = pixel_weights(wav, flux, grad=False)
wis_1 = pixel_weights(wav, flux, grad=True)
wis_0 *= mask[:-1]
wis_1 *= mask
wis_0[wis_0 == 0] = np.nan
wis_1[wis_1 == 0] = np.nan
In [6]:
plt_setting = {"figsize": (15, 6)}
plt.figure(**plt_setting)
plt.plot(wav, flux / np.max(flux), label="Star")
plt.plot(atm.wl, atm.transmission, label="Telluric")
plt.plot(atm.wl, atm.mask, "--", label="Mask")
plt.axhline(0.98)
plt.legend()
Out[6]:
In [7]:
plt.figure(**plt_setting)
plt.plot(wav[:-1], wis_0, "bs-", label="Mask Grad False")
plt.plot(wav, wis_1, "ko--", label="Mask Grad true")
w, f = (wclump[0], fclump[0])
wis1 = pixel_weights(w, f, grad=True)
wis0 = pixel_weights(w, f, grad=False)
plt.plot(w[:-1], wis0 * 1.05, "g+:", label="Clump, Grad False")
plt.plot(w, wis1 * 1.05, "r.-.", label="Clump, Grad True")
plt.legend()
plt.xlim(wclump[0][0] * 0.99999, wclump[0][-1] * 1.00001)
plt.show()
In [8]:
plt.figure(**plt_setting)
plt.plot(wav[:-1], wis_0, "bs-", label="grad False")
plt.plot(wav, wis_1, "ko--", label="grad true")
w, f = (wclump[1], fclump[1])
wis1 = pixel_weights(w, f, grad=True)
wis0 = pixel_weights(w, f, grad=False)
plt.plot(w[:-1], wis0 * 1.05, "g+:", label="Clump grad False")
plt.plot(w, wis1 * 1.05, "r.-.", label="Clump grad True")
plt.legend()
plt.xlim(wclump[-1][0] * 0.999999, wclump[-1][-1] * 1.00001)
plt.show()
Ffrom these two examples the calculations with the same gradients produces the same pixel weights The clumped version produces less weight though.
The masked version produces slightly different value for the last pixel due to how it is calculated. With the new graident all pixels are kept except in the clumped version the last pixel s the end and not in the middle so is not calculated with central difference but finite backward difference.
In [9]:
# Old and new indicate the split method.
print("Old with gradient {:0.06f}".format(RVprec_calc_masked(wav, flux, atm.mask, grad=True)))
print("New with gradient {:0.06f}".format(rv_precision(wav, flux, atm.mask, grad=True)))
print("Old without finite diff{:0.06f}".format(RVprec_calc_masked(wav, flux, atm.mask, grad=False)))
print("New with finite diff{:0.06f}".format(rv_precision(wav, flux, atm.mask, grad=False)))
Differences between versions with same gradient is at the 4th sf. These are not the correct scale, this will be addressed in the next section.
Assessing the changes to actual Figueira et al. values by not splitting on telluric lines first then applying mask to the weights, or slitting then calculating weights.
Convolving the spectra to vsini=1 an R=100k of a M0 spectra in the Z,Y,J,H,K bands to be consistent with paper. Using Old gradient and new gradient methods to access difference. The old gradient drops the last pixel, which drops many pixels when spectra is split between telluric lines.
In [10]:
# Explore relative difference of different bands
wav_, flux_ = load_aces_spectrum([3900, 4.5, 0, 0])
wav, flux = wav_selector(wav_, flux_, 0.7, 2.5)
table = []
table.append("Band, Cond#1, Split, Masked, ratio, Cond#1, Split, Masked, ratio")
table.append("Grad, False , True ")
# Get J band SNR normalization value
wav_j, flux_j = convolve_and_resample(
wav, flux, vsini=1, R=100000, band="J", sampling=3
)
snr_norm = snr_constant_band(wav_j, flux_j, snr=100, band="J")
for band in ["Z", "Y", "J", "H", "K"]:
atm = Atmosphere.from_band(band, bary=True)
w, f = convolve_and_resample(wav, flux, vsini=1, R=100000, band=band, sampling=3)
f /= snr_norm
atm = atm.at(w)
a = RVprec_calc_masked(w, f, atm.mask, grad=True)
b = RVprec_calc_masked(w, f, atm.mask, grad=False)
c = rv_precision(w, f, atm.mask, grad=True)
d = rv_precision(w, f, atm.mask, grad=False)
e = rv_precision(w, f, grad=True)
f = rv_precision(w, f, grad=False)
false_ratio = (d - b) / b
true_ratio = (c - a) / a
table.append(
"{0:5}, {1:4.02f}, {2:6.02f}, {3:6.02f}, {4:5.04f}, {5:6.02f}, {6:6.02f}, {7:6.02f}, {8:5.04f}".format(
band,
f.value,
b.value,
d.value,
false_ratio,
e.value,
a.value,
c.value,
true_ratio,
)
)
In [11]:
for line in table:
print(line)
This table shows that the methods produce very similar results. With new gradient it is 0.01 % difference.
This confirms that the new weighted mask is suitable.