In [1]:
import matplotlib.pyplot as plt
import numpy as np
from StellarSpectra.grid_tools import HDF5Interface
import StellarSpectra.constants as C
from StellarSpectra.utils import saveall
In [2]:
from matplotlib.ticker import FormatStrFormatter as FSF
from matplotlib.ticker import MaxNLocator
from matplotlib.ticker import MultipleLocator
In [3]:
#interface = HDF5Interface("../libraries/PHOENIX_F.hdf5")
interface = HDF5Interface("../../libraries/PHOENIX_TRES_F.hdf5")
In [4]:
wl_raw = interface.wl
#ind = (wl_raw > 5162.01) & (wl_raw < 5162.54) #for hires
#ind = (wl_raw > 5161.9) & (wl_raw < 5162.69) #for lowres
ind = (wl_raw > 5160.9) & (wl_raw < 5163.3) #for lowres, slightly wider wl range
#ind = (wl_raw > 3000) & (wl_raw < 9000) #for statistics
#ind = (wl_raw > 5150) & (wl_raw < 5200) #for lowres, all
wl_raw = wl_raw[ind]
Fake interpolated spectrum will be T=6000, L=4.0, Z=-0.5
In [5]:
flt60g40m05 = interface.load_flux({"temp":6000,"logg":4.0, "Z":-0.5, "alpha":0.0})[ind]
flt59g40m05 = interface.load_flux({"temp":5900,"logg":4.0, "Z":-0.5, "alpha":0.0})[ind]
flt61g40m05 = interface.load_flux({"temp":6100,"logg":4.0, "Z":-0.5, "alpha":0.0})[ind]
flt60g45m05 = interface.load_flux({"temp":6000,"logg":4.5, "Z":-0.5, "alpha":0.0})[ind]
flt60g35m05 = interface.load_flux({"temp":6000,"logg":3.5, "Z":-0.5, "alpha":0.0})[ind]
flt60g40m00 = interface.load_flux({"temp":6000,"logg":4.0, "Z":0.0, "alpha":0.0})[ind]
flt60g40m10 = interface.load_flux({"temp":6000,"logg":4.0, "Z":-1.0, "alpha":0.0})[ind]
In [6]:
spec_temp = (flt61g40m05 + flt59g40m05)/2.
spec_logg = (flt60g35m05 + flt60g45m05)/2.
spec_Z = (flt60g40m10 + flt60g40m00)/2.
In [7]:
def get_error(spec):
return (spec - flt60g40m05)/flt60g40m05
In [8]:
err_temp = get_error(spec_temp)
err_logg = get_error(spec_logg)
err_Z = get_error(spec_Z)
In [9]:
change_temp = (flt61g40m05 - flt59g40m05)/flt60g40m05
change_logg = (flt60g45m05 - flt60g35m05)/flt60g40m05
change_Z = (flt60g40m00 - flt60g40m10)/flt60g40m05
In [35]:
specs = [spec_temp, spec_logg, spec_Z]
errs = [err_temp, err_logg, err_Z]
change_spec = [change_temp, change_logg, change_Z]
labels = [r"$T_{\rm eff}$", r"$\log g$", r"[${\rm Fe}/{\rm H}$]"]
colors = ["b", "g", "r"]
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(6.5, 4.0), sharex=True)
aa = ax[0] #row 0, all cols
bb = ax[1] #row 1, all cols
for a in aa:
a.plot(wl_raw, flt60g40m05*1e-7, "k", label="True")
for a,spec,color,label in zip(aa, specs, colors, labels):
a.plot(wl_raw, spec*1e-7, color=color, label=label)
for a,label,color in zip(aa, labels, colors):
a.annotate(label, (0.1, 0.1), xycoords="axes fraction", ha="left", color=color, size=9)
aa[0].annotate("True", (0.1, 0.2), xycoords="axes fraction", ha="left", color="k", size=9)
#ax[0].legend(loc="lower left", prop={'size':8})
for b,change,color in zip(bb, change_spec, colors):
env = 100*np.abs(change)
b.fill_between(wl_raw, -env, env, color=color, alpha=0.2)
b.axhline(0, color="k", ls=":")
for b,err,color in zip(bb, errs, colors):
b.plot(wl_raw, 100*err, color=color)
scale = 200*np.max(np.abs(err_Z))
ax[0,0].set_xlim(wl_raw[0], wl_raw[-1])
ax[0,0].set_ylabel(r"$f_\lambda\, \times 10^{7}\, [\textrm{erg}\;\textrm{cm}^{-2}\;\textrm{s}^{-1}\;\textrm{\AA}^{-1}]$")
ax[1,0].set_ylabel(r"$e_\lambda / f_\lambda \quad$" + "[\%]")
for a in aa:
a.set_ylim(0.4, 1.0)
for a in aa[1:]:
a.yaxis.set_ticklabels([])
for b in bb:
b.set_ylim(-5.5, 5.5)
b.xaxis.set_major_formatter(FSF("%.0f"))
b.xaxis.set_major_locator(MultipleLocator(1.))
for b in bb[1:]:
b.yaxis.set_ticklabels([])
fig.text(0.48, 0.02, r"$\lambda$ [\AA]")
fig.subplots_adjust(left=0.1, right=0.9, hspace=0.1, wspace=0.1, top=0.95, bottom=0.12)
saveall(fig, "../../plots/interp_errors")
plt.show()
In contrast, what is the actual change in 100K, 0.5 dex surface gravity or metallicity?
In [33]:
for change in change_spec:
plt.plot(wl_raw, 100*change)
plt.show()
In [5]:
#errspec = interface.get_error_spectrum({"temp":6000,"logg":4.5, "Z":-0.5, "alpha":0.0}, axis="temp")[ind]
errspec1 = interface.get_error_spectrum({"temp":5900,"logg":3.5, "Z":-1.0, "alpha":0.0}, axis="temp")[ind]
errspec2 = interface.get_error_spectrum({"temp":5900,"logg":3.5, "Z":0.0, "alpha":0.0}, axis="temp")[ind]
errspec3 = interface.get_error_spectrum({"temp":5900,"logg":4.5, "Z":-1.0, "alpha":0.0}, axis="temp")[ind]
errspec4 = interface.get_error_spectrum({"temp":5900,"logg":4.5, "Z":0.0, "alpha":0.0}, axis="temp")[ind]
errspec5 = interface.get_error_spectrum({"temp":6100,"logg":3.5, "Z":-1.0, "alpha":0.0}, axis="temp")[ind]
errspec6 = interface.get_error_spectrum({"temp":6100,"logg":3.5, "Z":0.0, "alpha":0.0}, axis="temp")[ind]
errspec7 = interface.get_error_spectrum({"temp":6100,"logg":4.5, "Z":-1.0, "alpha":0.0}, axis="temp")[ind]
errspec8 = interface.get_error_spectrum({"temp":6100,"logg":4.5, "Z":0.0, "alpha":0.0}, axis="temp")[ind]
In [6]:
#estimate the covariance between these
errspec = np.array([errspec1, errspec2, errspec3, errspec4, errspec5, errspec6, errspec7, errspec8])
#cov_mat = np.cov(errspec,rowvar=0)
cov_mat = np.dot(errspec.T, errspec)/8.
plt.imshow(cov_mat, origin="upper", interpolation="none")
plt.show()
In [10]:
def force_to_zero(mat, dist=10):
'''
If a row and column are greater than 10 apart, set to 0.
'''
nmat = mat.copy()
nrows, ncols = nmat.shape
for i in range(nrows):
for j in range(ncols):
if np.abs(i - j) > dist:
nmat[i,j] = 0.0
return nmat
In [35]:
mat2 = force_to_zero(cov_mat, dist=50)
In [36]:
plt.imshow(mat2, origin="upper", interpolation="none")
plt.show()
In [20]:
draws = [np.random.multivariate_normal(np.zeros_like(wl_raw), cov_mat) for i in range(25)]
In [37]:
draws2 = [np.random.multivariate_normal(np.zeros_like(wl_raw), mat2) for i in range(2)]
In [38]:
for draw in draws2:
plt.plot(wl_raw, draw, lw=0.3)
for spec in errspec:
plt.plot(wl_raw, spec)
plt.show()
In [8]:
#Over metallicity
flt61g45m05_int = (flt61g45m10 + flt61g45p00)/2
flt61g35m05_int = (flt61g35m10 + flt61g35p00)/2
flt59g45m05_int = (flt59g45m10 + flt59g45p00)/2
flt59g35m05_int = (flt59g35m10 + flt59g35p00)/2
#Plot all of these to see what they look like
errspecm0 = (flt61g45m05 - flt61g45m05_int)/flt61g45m05
errspecm1 = (flt61g35m05 - flt61g35m05_int)/flt61g35m05
errspecm2 = (flt59g45m05 - flt59g45m05_int)/flt59g45m05
errspecm3 = (flt59g35m05 - flt59g35m05_int)/flt59g35m05
errspecm = (errspecm0 + errspecm1 + errspecm2 + errspecm3)/4
In [9]:
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)
errspecs = [errspecm0, errspecm1, errspecm2, errspecm3, errspecm]
labels = ["0", "1", "2", "3", "avg"]
for errspec, label in zip(errspecs, labels):
ax.plot(wl_raw, errspec*100, label=label)
ax.legend()
plt.show()
In [89]:
fig, ax = plt.subplots(nrows=2, figsize=(8,8), sharex=True)
colors = ["b", "c", "r", "k"]
#labels = ["0.0", "-0.5", "-1.0", "-0.5 int"]
labels = ["-1.0", "-0.5", "0.0", "0.0 int"]
fls = [flt60g45m10, flt60g45m05, flt60g45p00, flt60g45m05_int]
for fl, color, label in zip(fls, colors, labels):
ax[0].plot(wl_raw, fl/fl[0], "o", color=color, label=label)
for wl in wl_raw:
ax[0].axvline(wl, color="0.2", lw=0.1)
ax[0].legend(loc="lower right")
ax[1].axhline(0.0, color="0.2")
ax[1].plot(wl_raw, errspec*100, "k-")
ax[1].plot(wl_raw, errspec*100, "ko")
ax[-1].set_xlim(wl_raw[0], wl_raw[-1])
ax[1].set_ylabel("Error [\%]")
ax[1].xaxis.set_major_formatter(FSF("%.0f"))
ax[-1].set_xlabel(r"$\lambda$ [\AA]")
fig.savefig("../plots/interp_tests/metallicity.png")
plt.show()
In [3]:
from StellarSpectra.model import Model
from StellarSpectra.spectrum import DataSpectrum
from StellarSpectra.grid_tools import TRES, HDF5Interface
import scipy.sparse as sp
import matplotlib
myDataSpectrum = DataSpectrum.open("../../data/WASP14/WASP14-2009-06-14.hdf5", orders=np.array([22]))
myInstrument = TRES()
myHDF5Interface = HDF5Interface("../../libraries/PHOENIX_TRES_F.hdf5")
myErrorHDF5Interface = HDF5Interface("../../libraries/PHOENIX_TRES_F.hdf5")
In [4]:
#Load a model using the JSON file
myModel = Model.from_json("../WASP14_22_model_final_region.json", myDataSpectrum, myInstrument,
myHDF5Interface, myErrorHDF5Interface)
myOrderModel = myModel.OrderModels[0]
spec = myModel.get_data()
wl = spec.wls[0]
fl = spec.fls[0]
sigma = spec.sigmas[0]
model_fl = myOrderModel.get_spectrum()
residuals = fl - model_fl
cheb = myOrderModel.get_Cheb()
In [5]:
params = {'vz': -4.8456827988388413, 'logOmega': -19.658942686778765, 'vsini': 5.0328046854058668,
'temp': 6220., 'logg': 4.2, 'Z': -0.4}
myModel.update_Model(params)
In [6]:
cov = np.asarray(myModel.OrderModels[0].get_Cov_interp().todense())
In [7]:
ind2 = (wl > 5160.9) & (wl < 5163.3)
wlt = wl[ind2]
In [8]:
start = np.argwhere(ind2 == True)[0][0]
end = np.argwhere(ind2 == True)[-1][0]
In [9]:
covt = cov[start:end, start:end]
In [35]:
#norm = matplotlib.colors.Normalize(vmin=-31, vmax=-27.5, clip=True)
#norm = matplotlib.colors.Normalize(vmin=-6.9981851210520075e-32, vmax=1.9066636993100171e-30, clip=True)
scale = np.max(np.abs(covt))
norm = matplotlib.colors.SymLogNorm(7e-32, vmin=-scale, vmax=scale)
#vmin=-31, vmax=-27.5
In [57]:
def plot_matrix(wl, cov, title=False, axl=False):
'''
Construct the plot!
'''
fig = plt.figure(figsize=(2.95,2.35))
ax = fig.add_subplot(111)
left, right = wl[0], wl[-1]
#cmap=plt.cm.Blues,
img = ax.imshow(cov, interpolation="none", origin="upper", extent=[left, right, right, left],
norm=norm, cmap=plt.cm.RdBu)
ax.xaxis.set_major_formatter(FSF("%.0f"))
ax.xaxis.set_major_locator(MultipleLocator(1.))
ax.yaxis.set_major_formatter(FSF("%.0f"))
ax.yaxis.set_major_locator(MultipleLocator(1.))
if axl:
ax.set_xlabel(r"$\lambda$ [\AA]")
ax.set_ylabel(r"$\lambda$ [\AA]")
else:
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
if title:
cax = fig.add_axes([0.82, 0.22, 0.03, 0.65])
cb = fig.colorbar(img, cax=cax)
cb.set_ticks(matplotlib.ticker.FixedLocator([-1e-30, -1e-31, 0, 1e-31, 1e-30]))
#fig.text(0.25, 0.91,
# r"$\mathsf{C}_{ij}^\textrm{int}\,[\textrm{erg}^2\;\textrm{cm}^{-4}\;\textrm{s}^{-2}\;\textrm{\AA}^{-2}]$")
#img.set_clim(-6.9981851210520075e-32, 1.9066636993100171e-30)
fig.subplots_adjust(left=0.2, bottom=0.12, top=0.95, right=0.79)
return fig,ax
In [58]:
fig, ax = plot_matrix(wlt, covt, title=True, axl=True)
ax.annotate(r"$\mathcal{K}^{\rm int}$", (5163.08, 5161.3), ha="right", size=9)
saveall(fig, "../../plots/interp_matrix")
plt.show()
In [ ]: