In [1]:

    

#Import all of the relevant python modules
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter as FSF
from matplotlib.ticker import MaxNLocator
from matplotlib.ticker import MultipleLocator


    

In [9]:

    

from StellarSpectra.model import Model
from StellarSpectra.spectrum import DataSpectrum
from StellarSpectra.grid_tools import SPEX, HDF5Interface
from StellarSpectra.utils import saveall
from StellarSpectra import utils
import scipy.sparse as sp
import numpy as np

myDataSpectrum = DataSpectrum.open("../../data/Gl51/Gl51RA.hdf5", orders=np.array([0]))
myInstrument = SPEX()
myHDF5Interface = HDF5Interface("../../libraries/PHOENIX_SPEX_M.hdf5")
myErrorHDF5Interface = HDF5Interface("../../libraries/PHOENIX_SPEX_M.hdf5")


    

In [4]:

    

#Load a model using the JSON file
myModel = Model.from_json("Gl51_model0_final.json", myDataSpectrum, myInstrument, myHDF5Interface, myErrorHDF5Interface)
myOrderModel = myModel.OrderModels[0]
spec = myModel.get_data()

wl = spec.wls[0]
fl = spec.fls[0]

model_fl = myOrderModel.get_spectrum()

residuals = fl - model_fl

cheb = myOrderModel.get_Cheb()


    

    




Determine Chunk Log: Wl is 2048
Determine Chunk Log: Wl is 2048
Creating OrderModel 0
creating region 





    




INFO:RegionCovarianceMatrix 22056.92494864523:Created region
INFO:CovarianceMatrix:Added region 0 to self.RegionList
INFO:RegionCovarianceMatrix 22652.978100190612:Created region
INFO:CovarianceMatrix:Added region 1 to self.RegionList
INFO:RegionCovarianceMatrix 22084.78296470615:Created region
INFO:CovarianceMatrix:Added region 2 to self.RegionList
INFO:RegionCovarianceMatrix 22625.371116290782:Created region
INFO:CovarianceMatrix:Added region 3 to self.RegionList
INFO:RegionCovarianceMatrix 22609.26887539992:Created region
INFO:CovarianceMatrix:Added region 4 to self.RegionList






    




 0 0 {'sigma': 122.01839684048089, 'loga': -14.7541719333674, 'mu': 22056.92494864523}
creating region  1 1 {'sigma': 99.9072991296247, 'loga': -14.73332879603619, 'mu': 22652.978100190612}
creating region  2 2 {'sigma': 120.99653620176579, 'loga': -14.83385277442167, 'mu': 22084.78296470615}
creating region  3 3 {'sigma': 94.83258590913002, 'loga': -14.744408866651673, 'mu': 22625.371116290782}
creating region  4 4 {'sigma': 85.14706896796783, 'loga': -14.81901727721987, 'mu': 22609.26887539992}
creating region 





    




INFO:RegionCovarianceMatrix 23380.627852924943:Created region
INFO:CovarianceMatrix:Added region 5 to self.RegionList






    




 5 5 {'sigma': 73.25433304613068, 'loga': -14.81368835095075, 'mu': 23380.627852924943}

In [5]:

    

draws = np.load("Gl51_draws.npy").T


    

In [6]:

    

#Colorbrewer bands
s3 = '#fee6ce'
s2 = '#fdae6b'
s1 = '#e6550d'


    

In [14]:

    

fig, ax = plt.subplots(nrows=2, figsize=(7,3), sharex=True)

wlm = wl * 1e-4 #convert to microns

l0, = ax[0].plot(wlm, fl*1e14, "b")
l1, = ax[0].plot(wlm, model_fl*1e14, "r")

ax[0].annotate("data", (0.86, 0.75), xycoords="axes fraction", ha="right", color="b", size=9)
ax[0].annotate("model", (0.97, 0.75), xycoords="axes fraction", ha="right", color="r", size=9)

#ax[0].legend([l0, l1], ["data", "model"], loc="upper right", ncol=2, prop={'size':10})
ax[0].yaxis.set_major_locator(MaxNLocator(nbins=7, prune='lower'))

#Jump again and now plot the envelope of the draws
min_spec, max_spec = utils.std_envelope(draws)
ax[1].fill_between(wlm, 3e14*min_spec, 3e14*max_spec, zorder=0, color=s3)
ax[1].fill_between(wlm, 2e14*min_spec, 2e14*max_spec, zorder=0, color=s2)
ax[1].fill_between(wlm, 1e14*min_spec, 1e14*max_spec, zorder=0, color=s1)

ax[1].plot(wl, residuals*1e14, "w", lw=1.8, zorder=1, ls="steps-mid")
l2, = ax[1].plot(wlm, residuals*1e14, "k", zorder=1)

ax[1].xaxis.set_major_formatter(FSF("%.1f"))
ax[1].set_xlabel(r"$\lambda$ [$\mu$m]")

scale = 1.5 * np.max(np.abs(residuals)) * 1e14

ax[0].yaxis.set_major_locator(MaxNLocator(nbins=6, prune='lower'))   
ax[1].yaxis.set_major_locator(MaxNLocator(nbins=6, prune='upper'))  

ax[1].set_ylim(-scale, scale)
ax[1].set_xlim(np.min(wlm), np.max(wlm))
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", prop={'size':10})

fig.subplots_adjust(hspace=0, left=0.1, right=0.9, bottom=0.17, top=0.93)
fig.text(0.01, 0.83, r"$f_\lambda\, \times 10^{-14}\, \quad [\textrm{erg}\;\textrm{cm}^{-2}\;\textrm{s}^{-1}\;\textrm{\AA}^{-1}]$", rotation="vertical")

saveall(fig, "../../plots/residuals_Gl51_logg")
plt.show()


    

In [ ]: