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 [2]:

    

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

myDataSpectrum = DataSpectrum.open("../../data/WASP14/WASP14-2009-06-14.hdf5", orders=np.array([22]))
myInstrument = TRES()
myHDF5Interface = HDF5Interface("../../libraries/PHOENIX_TRES_F.hdf5")


    

In [3]:

    

#Load a model using the JSON file
myModel = Model.from_json("../WASP14/model_final.json", myDataSpectrum, myInstrument, myHDF5Interface)
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 8192
Creating OrderModel 0

In [22]:

    

fig, ax = plt.subplots(nrows=2, figsize=(3.5,2.7), sharex=True)

ax[0].plot(wl, fl*1e13, "w", lw=1.8)
ax[0].plot(wl, model_fl*1e13, "w", lw=1.8)

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

#ax[0].legend([l0, l1], ["data", "model"], loc="lower right", ncol=2, prop={'size':10})
ax[0].yaxis.set_major_locator(MaxNLocator(nbins=5, prune='lower'))
ax[0].fill_betweenx(np.array([0.1, 2.7]), 5144, 5150, color="0.5", alpha=0.5)
ax[0].set_ylim(0.1, 2.7)

ax[1].fill_betweenx(np.array([-0.9, 0.9]), 5144, 5150, color="0.5", alpha=0.5)
ax[1].plot(wl, residuals*1e13, "w", lw=1.8)
l2, = ax[1].plot(wl, residuals*1e13, "k")
ax[1].xaxis.set_major_formatter(FSF("%.0f"))
ax[1].set_xlabel(r"$\lambda$ [\AA]")
ax[1].yaxis.set_major_locator(MaxNLocator(nbins=5))
ax[1].set_xlim(5135, 5200)
ax[1].set_ylim(-0.9, 0.9)
#legend1 = ax[1].legend([l0, l1], ["data", "model", ], loc="upper right", prop={'size':10}, ncol=2)
#ax[1].legend([l2], ["residuals"], loc="lower center", prop={'size':10})
#ax[1].add_artist(legend1)

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

ax[1].annotate("residuals", (0.97, 0.82), xycoords="axes fraction", ha="right", color="k", size=9)

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

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


    

In [21]:

    

fig, ax = plt.subplots(nrows=2, sharex=True, figsize=(3.5, 2.7))

#Use the model and residuals from the previos plot.
ind = (wl > 5144) & (wl < 5150)

l0, = ax[0].plot(wl[ind], fl[ind]*1e13, "b", label="data", drawstyle="steps-mid")
l1, = ax[0].plot(wl[ind], model_fl[ind]*1e13, "r", label="model", drawstyle="steps-mid")


l2, = ax[1].plot(wl[ind], residuals[ind]*1e13, "k", label="residuals", drawstyle="steps-mid")

for a in ax:
    a.yaxis.set_major_locator(MaxNLocator(nbins=6, prune='upper'))

#ax[0].legend([l0, l1], ["data", "model"], loc="lower left", prop={'size':10})
#ax[1].legend([l2], ["residuals"], loc="upper left", prop={'size':10})
ax[1].set_ylim(-0.2, 0.2)

ax[1].set_xlabel(r"$\lambda$ [\AA]")
ax[1].xaxis.set_major_formatter(FSF("%.0f"))
fig.text(0.01, 0.8, r"$f_\lambda\, \times 10^{-13}\, [\textrm{erg}\;\textrm{cm}^{-2}\;\textrm{s}^{-1}\;\textrm{\AA}^{-1}]$", rotation="vertical")
fig.subplots_adjust(bottom=0.16, top=0.93, hspace=0, left=0.17, right=0.83)
saveall(fig, "../../plots/class0_residuals")
plt.show()


    

In [23]:

    

def estimated_autocorrelation(x):
    '''From http://stackoverflow.com/questions/14297012/estimate-autocorrelation-using-python'''
    n = len(x)
    variance = np.var(x)
    x = x - np.mean(x)
    r = np.correlate(x, x, mode = 'full')[-n:] #takes the positive half of the correlation
    #assert N.allclose(r, N.array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)]))
    result = r/(variance*(np.arange(n, 0, -1))) #divides by the number pixels used in the estimation
    return result


    

In [24]:

    

#Now estimate the autocorrelation of the residuals 
fig = plt.figure(figsize=(3.5,1.7))
ax = fig.add_subplot(111)
ax.axhline(0, ls=":", color="k")
ax.plot(estimated_autocorrelation(residuals), color="0.3", drawstyle="steps-mid")
ax.set_xlabel("pixel offset")
ax.set_ylabel("autocorrelation")
ax.set_xlim(-0.3,35)
ax.set_ylim(-0.2,1.1)

fig.subplots_adjust(bottom=0.22, top=0.97, left=0.17, right=0.83)
saveall(fig, "../../plots/class0_autocorrelation")
plt.show()


    

In [ ]: