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
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/Kurucz_TRES.hdf5")
myErrorHDF5Interface = HDF5Interface("../../libraries/Kurucz_TRES.hdf5")
#Load a model using the JSON file
#Taken from:
#/n/home07/iczekala/StellarSpectra/output/WASP14/Kurucz/21_22_23/logg/cov/2014-08-06/run18
myModel = Model.from_json("WASP14_Kurucz_logg_model_final.json", myDataSpectrum, myInstrument, myHDF5Interface, myErrorHDF5Interface)
In [3]:
myOrderModel = myModel.OrderModels[1]
model_flux = myOrderModel.get_spectrum()
In [4]:
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 [5]:
#draws = utils.random_draws(cov, num=4)
draws = np.load("krucuz_residual_draws.npy")
In [6]:
utils.visualize_draws(draws, num=200)
In [6]:
#Colorbrewer bands
s3 = '#fee6ce'
s2 = '#fdae6b'
s1 = '#e6550d'
In [25]:
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})
#plot 2 draws, then
sep = np.std(draws[0])
#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.8 * 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_Kurucz_logg")
#ax[1].set_ylim(-9e-14, 9e-14)
#saveall(fig, "../../plots/residuals_Kurucz_logg_scale")
plt.show()
In [ ]: