A suite of tests to explore interpolation amonst the stellar parameters


In [1]:
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline as IUS
from StellarSpectra.grid_tools import HDF5Interface
myInterface = HDF5Interface("../libraries/PHOENIX_submaster.hdf5")

In [2]:
def load_pixel(index, temp, logg, Z):
    '''
    Return the flux value at a specific pixel in a spectrum, defined by params
    '''
    params = {"temp":temp, "logg":logg, "Z":Z, "alpha":0.0}
    flux = myInterface.load_flux(params)
    return flux[index]

def get_wl(index):
    return myInterface.wl[index]

In [4]:
#Determine truncation indices for order of interest
wl = myInterface.wl
inds = np.argwhere((wl > 5120) & (wl < 5220))[np.array([0, -1])]
inds[0], inds[1]


Out[4]:
(array([345476]), array([355478]))

In [37]:
flux = myInterface.load_flux({"temp":6000, "logg":4.0, "Z":0.0, "alpha":0.0})
print(len(flux))
#plt.plot(flux[inds[0]:inds[1]])
#plt.show()


983561

In [6]:
temps = myInterface.points["temp"]

In [14]:
#Pixel indices for lines
cont = 5018
mg_bot = 5268
mg_side = 5242
fe_bot = 5053
fe_side = 5050

In [9]:
flux_cont = [load_pixel(cont, temp, 4.0, 0.0) for temp in temps]
flux_mg_bot = [load_pixel(mg_bot, temp, 4.0, 0.0) for temp in temps]
flux_mg_side = [load_pixel(mg_side, temp, 4.0, 0.0) for temp in temps]
flux_fe_bot = [load_pixel(fe_bot, temp, 4.0, 0.0) for temp in temps]
flux_fe_side = [load_pixel(fe_side, temp, 4.0, 0.0) for temp in temps]

In [11]:
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.plot(temps, flux_cont, label="continuum")
ax.plot(temps, flux_cont, "ko")
ax.plot(temps, flux_mg_bot, label="Mg bottom")
ax.plot(temps, flux_mg_bot, "ko")
ax.plot(temps, flux_mg_side, label="Mg side")
ax.plot(temps, flux_mg_side, "ko")
ax.plot(temps, flux_fe_bot, label="Fe bottom")
ax.plot(temps, flux_fe_bot, "ko")
ax.plot(temps, flux_fe_side, label="Fe side")
ax.plot(temps, flux_fe_side, "ko")

#try fitting a spline to fe_side
myspline = IUS(temps, flux_fe_side)
fine_temps = np.linspace(5000, 7000, num=300)
fine_fe_side = myspline(fine_temps)
ax.plot(fine_temps, fine_fe_side, "k", lw=0.5)


ax.set_xlabel(r"Temperature ($K$)")
ax.set_ylabel(r"$\propto f_\lambda$")
ax.legend(loc="upper left")
fig.savefig("../plots/interpolation_temp.png")
plt.show()

log g tests


In [62]:
loggs = np.arange(3.5, 5.6, 0.5)

In [63]:
flux_cont = [load_pixel(cont, 6000, logg, 0.0) for logg in loggs]
flux_mg_bot = [load_pixel(mg_bot, 6000, logg, 0.0) for logg in loggs]
flux_mg_side = [load_pixel(mg_side, 6000, logg, 0.0) for logg in loggs]
flux_fe_bot = [load_pixel(fe_bot, 6000, logg, 0.0) for logg in loggs]
flux_fe_side = [load_pixel(fe_side, 6000, logg, 0.0) for logg in loggs]

In [66]:
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.plot(loggs, flux_cont, label="continuum")
ax.plot(loggs, flux_cont, "ko")
ax.plot(loggs, flux_mg_bot, label="Mg bottom")
ax.plot(loggs, flux_mg_bot, "ko")
ax.plot(loggs, flux_mg_side, label="Mg side")
ax.plot(loggs, flux_mg_side, "ko")
ax.plot(loggs, flux_fe_bot, label="Fe bottom")
ax.plot(loggs, flux_fe_bot, "ko")
ax.plot(loggs, flux_fe_side, label="Fe side")
ax.plot(loggs, flux_fe_side, "ko")

ax.set_xlabel(r"Temperature ($K$)")
ax.set_ylabel(r"$\propto f_\lambda$")
ax.legend(loc="lower right")
plt.show()

Metallicity tests


In [16]:
Zs = myInterface.points["Z"]

In [21]:
flux_cont = [load_pixel(cont, 6000, 4.0, Z) for Z in Zs]
flux_mg_bot = [load_pixel(mg_bot, 6000, 4.0, Z) for Z in Zs]
flux_mg_side = [load_pixel(mg_side, 6000, 4.0, Z) for Z in Zs]
flux_fe_bot = [load_pixel(fe_bot, 6000, 4.0, Z) for Z in Zs]
flux_fe_side = [load_pixel(fe_side, 6000, 4.0, Z) for Z in Zs]

In [23]:
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.plot(Zs, flux_cont, label="continuum")
ax.plot(Zs, flux_cont, "ko")
ax.plot(Zs, flux_mg_bot, label="Mg bottom")
ax.plot(Zs, flux_mg_bot, "ko")
ax.plot(Zs, flux_mg_side, label="Mg side")
ax.plot(Zs, flux_mg_side, "ko")
ax.plot(Zs, flux_fe_bot, label="Fe bottom")
ax.plot(Zs, flux_fe_bot, "ko")
ax.plot(Zs, flux_fe_side, label="Fe side")
ax.plot(Zs, flux_fe_side, "ko")

ax.set_xlabel(r"Metallicity")
ax.set_ylabel(r"$\propto f_\lambda$")
ax.legend(loc="lower right")
plt.show()

In [22]:
flux_fe_bot


Out[22]:
[6378120.0, 5492418.5, 3951947.0]

Spline interpolation


In [13]:
from scipy.interpolate import RectBivariateSpline
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline as IUS

In [23]:
temps = myInterface.points["temp"]
print(temps)
loggs = myInterface.points["logg"]
print(loggs)


[5000 5100 5200 5300 5400 5500 5600 5700 5800 5900 6000 6100 6200 6300 6400
 6500 6600 6700 6800 6900 7000]
[ 3.5  4.   4.5  5.   5.5]

In [32]:
#This is the part that actually takes a while, loading all of the spectra. 
#If we could do it on a computer with a lot of memory, it might not
#be so bad
values = np.empty((temps.size, loggs.size))
for i,temp in enumerate(temps):
    for j, logg in enumerate(loggs):
        values[i,j] = load_pixel(0, temp, logg, Z=0.0)

In [25]:
def interpolate_pixel(index, mytemp, mylogg):
    
    #Create a 2D array of shape (temps.size, loggs.size) with the value of the flux at the pixel
    #values = np.empty((temps.size, loggs.size))
    
    #for i,temp in enumerate(temps):
    #    for j, logg in enumerate(loggs):
    #        values[i,j] = load_pixel(index, temp, logg, Z=0.0)
    
    #create a bivariate spline over this grid
    spl = RectBivariateSpline(temps, loggs, values)
    
    return spl(mytemp, mylogg)

In [36]:
timeit value = interpolate_pixel(0, 6010, 3.8)


1000 loops, best of 3: 787 µs per loop

In [35]:
print(value)


[[ 6129205.16129755]]

In [12]:
plt.imshow(values,origin="upper",interpolation="none")
plt.show()

In [ ]: