In [1]:
%matplotlib
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter as FSF
from matplotlib.ticker import MaxNLocator
from matplotlib.ticker import MultipleLocator
import numpy as np
from IPython.html import widgets # Widget definitions
from IPython.display import display, clear_output, HTML # Used to display widgets in the notebook
from IPython.html.widgets import interact, interactive, fixed
[ -13.37354867, 5188.75488375, 5.79597804], [ -13.58059031, 5151.82322042, 6.15930413], [ -13.38718559, 5142.85381703, 7.06802524], [ -13.76186598, 5202.25227836, 7.45325824], [ -13.64347323, 5150.74527037, 6.67246923], [ -13.70379713, 5169.07086766, 6.47508788], [ -13.7420562, 5198.61903481, 5.48643032], [ -13.7557276, 5170.68109, 5.07069537], [ -13.4095891, 5235.30766, 4.82242126], [ -13.70503714, 5187.83335919, 5.30911125], [ -13.74135283, 5234.57274185, 5.49249104]
In [2]:
from StellarSpectra.model import Model
from StellarSpectra.spectrum import DataSpectrum
from StellarSpectra.grid_tools import TRES, SPEX, HDF5Interface
from StellarSpectra.utils import saveall
from StellarSpectra import utils
import scipy.sparse as sp
myDataSpectrum = DataSpectrum.open("../../data/WASP14/WASP14-2009-06-14.hdf5", orders=np.array([21, 22, 23]))
myInstrument = TRES()
# myInstrument = SPEX()
myHDF5Interface = HDF5Interface("../../libraries/PHOENIX_TRES_F.hdf5")
myErrorHDF5Interface = HDF5Interface("../../libraries/PHOENIX_TRES_F.hdf5")
#Load a model using the JSON file
#Update this:
myModel = Model.from_json("WASP14_PHOENIX_model0_final.json", myDataSpectrum, myInstrument, myHDF5Interface, myErrorHDF5Interface)
#old
#myModel = Model.from_json("../WASP14_22_model_final_region.json", myDataSpectrum, myInstrument, myHDF5Interface)
In [4]:
myOrderModel = myModel.OrderModels[1]
model_flux = myOrderModel.get_spectrum()
In [5]:
spec = myModel.get_data()
wl = spec.wls[1]
fl = spec.fls[1]
model_fl = myOrderModel.get_spectrum()
residuals = fl - model_fl
mask = spec.masks[1]
cov = myModel.OrderModels[1].get_Cov().todense()
In [6]:
draws = np.load("PHOENIX_draws.npy").T
In [9]:
utils.visualize_draws(draws, num=50)
In [7]:
#Colorbrewer bands
s3 = '#fee6ce'
s2 = '#fdae6b'
s1 = '#e6550d'
In [9]:
fig, ax = plt.subplots(nrows=2, figsize=(7,3), sharex=True)
l0, = ax[0].plot(wl, fl*1e13, "b")
l1, = ax[0].plot(wl, model_flux*1e13, "r")
ax[0].annotate("data", (0.86, 0.1), xycoords="axes fraction", ha="right", color="b", size=9)
ax[0].annotate("model", (0.97, 0.1), xycoords="axes fraction", ha="right", color="r", size=9)
#ax[0].legend([l0, l1], ["data", "model"], loc="lower right", ncol=2, prop={'size':10})
#Jump again and now plot the envelope of the draws
min_spec, max_spec = utils.std_envelope(draws)
ax[1].fill_between(wl, 3e13*min_spec, 3e13*max_spec, zorder=0, color=s3)
ax[1].fill_between(wl, 2e13*min_spec, 2e13*max_spec, zorder=0, color=s2)
ax[1].fill_between(wl, 1e13*min_spec, 1e13*max_spec, zorder=0, color=s1)
ax[1].plot(wl, residuals*1e13, "w", lw=1.8, zorder=1, ls="steps-mid")
l2, = ax[1].plot(wl, residuals*1e13, "k", zorder=1, ls="steps-mid")
ax[1].set_xlabel(r"$\lambda$ [\AA]")
ax[1].xaxis.set_major_formatter(FSF("%.0f"))
ax[1].set_xlim(np.min(wl), np.max(wl))
scale = 1.2 * np.max(np.abs(residuals)) * 1e13
ax[1].set_ylim(-scale, scale)
ax[1].annotate("residuals", (0.97, 0.8), xycoords="axes fraction", ha="right", color="k", size=9)
#ax[1].legend([l2], ["residuals"], loc="upper right", ncol=2, prop={'size':10})
ax[0].yaxis.set_major_locator(MaxNLocator(nbins=6, prune='lower'))
ax[1].yaxis.set_major_locator(MaxNLocator(nbins=6, prune='upper'))
fig.text(0.01, 0.83, r"$f_\lambda\, \times 10^{-13}\, \quad [\textrm{erg}\;\textrm{cm}^{-2}\;\textrm{s}^{-1}\;\textrm{\AA}^{-1}]$", rotation="vertical")
fig.subplots_adjust(hspace=0, left=0.1, right=0.9, bottom=0.17, top=0.93)
saveall(fig, "../../plots/residuals_PHOENIX_logg")
plt.show()
In [ ]: