In [6]:
#Odyssey
from StellarSpectra import grid_tools
from StellarSpectra import constants as C
import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline
from StellarSpectra.grid_tools import TRES
raw_library_path = "/home/ian/Grad/Research/Disks/StellarSpectra/libraries/raw/PHOENIX/"
mygrid = grid_tools.PHOENIXGridInterface(base=raw_library_path)
In [2]:
spec = mygrid.load_file({"temp":5000, "logg":3.5, "Z":0.0,"alpha":0.0})
wl_dict_t = spec.calculate_log_lam_grid(min=True)
wl_t = wl_dict_t["wl"]
wl_dict_f = spec.calculate_log_lam_grid(min=False)
wl_f = wl_dict_f["wl"]
In [5]:
def calculate_min_v(wl_dict):
CDELT1 = wl_dict["CDELT1"]
min_v = C.c_kms * (10**CDELT1 - 1)
return min_v
In [7]:
print(calculate_min_v(wl_dict_t), "km/s")
print(calculate_min_v(wl_dict_f), "km/s")
What is the smallest velocity spacing of the spectra? Presumably this would be set by the macroturbulence (is it measured as FWHM? or as sigma). Look this up in the paper.
Let's compare the different kinds of interpolation to make sure that we are properly capturing all of the "information" contained in the spectrum. Try spline interpolation on three different grids to a single common grid, and then use np.allclose to determine if we lost anything.
log_lam_grid(min=True)log_lam_grid(min=False)log_lam_grid(min=True) grid to log_lam_grid(min=False) gridlog_lam_grid(min=False) grid to log_lam_grid(min=True) gridCompare all using np.allclose to see if any information has been lost.
In [60]:
#Create a finely spaced grid, say 0.1kms
wl_dict_fine = grid_tools.create_log_lam_grid(wl_start=3000, wl_end=13000, min_vc=0.06/C.c_kms)
wl_fine = wl_dict_fine['wl']
print(calculate_min_v(wl_dict_fine), "km/s") #note, this grid will be a power of 2, so min_vc might be what you specify
How can we be sure that using a 0.2 km/s grid and k=5 spline that we are not loosing any information? Basically, what if we interpolate to this grid, then use it to interpolate back to the original, linearly spaced grid? If this test passes, then we know that we are doing OK. However, how does using a larger wavelength grid like this slow down the FFT?
In [59]:
interp = InterpolatedUnivariateSpline(spec.wl, spec.fl, k=5) #This higher power is needed to go from the raw grid to evenly spaced
fl_fine = interp(wl_fine)
#now try interpolating *back* to the original grid
interp_back = InterpolatedUnivariateSpline(wl_fine, fl_fine, k=5)
fl_back = interp_back(spec.wl)
np.allclose(spec.fl, fl_back)
Out[59]:
In [58]:
interp = InterpolatedUnivariateSpline(spec.wl, spec.fl, k=5) #This higher power is needed to go from the raw grid to evenly spaced
fl_fine = interp(wl_fine)
fl_t = interp(wl_t)
interp_fine = InterpolatedUnivariateSpline(wl_fine,fl_fine, k=3)
fl_fine_t = interp_fine(wl_t)
np.allclose(fl_t, fl_fine_t)
Out[58]:
The closest we can get away with is min_v = 0.5 km/s and spline k=4 that still satisfies np.allclose = True. So, because the grid actually contains information in the irregularly spaced pixels, we need to go to a very high oversampling (X4) to preserve this.
A grid that is spaced at 0.1 km/s will not work.
In [46]:
timeit fl_fine_t = interp_fine(wl_t)
In [27]:
fig, ax = plt.subplots(nrows=2, sharex=True)
ax[0].plot(wl_t, fl_t)
ax[0].plot(wl_t, fl_fine_t)
ax[1].plot(wl_t, fl_t - fl_fine_t)
plt.show()
This is only bad at the 1e-6 level.
In [29]:
interp = InterpolatedUnivariateSpline(spec.wl, spec.fl, k=5)
fl_t = interp(wl_t)
fl_f = interp(wl_f)
fl_t_f = InterpolatedUnivariateSpline(wl_t,fl_t, k=5)(wl_f)
np.allclose(fl_f, fl_t_f)
Out[29]:
In [15]:
import matplotlib.pyplot as plt
In [19]:
fig, ax = plt.subplots(nrows=2, sharex=True)
ax[0].plot(wl_f, fl_f)
ax[0].plot(wl_f, fl_t_f)
ax[1].plot(wl_f, fl_f - fl_t_f)
plt.show()
The result is that some informaiton around spectra lines is lost at the 1e-3 level. I think this isn't the best...
By eye, it looks like the lines that have a semi-circular bottom are the worst performing (at the edges), while the Gaussian lines don't do so badly.
What are some other tests that we could try? What about interpolating to a fine (0.1 km/s) grid, and THEN interpolating to a minimum spaced grid (0.2 km/s), and then comparing this to interpolating directly from the natural grid to 0.2 km/s?
In [16]:
wl_dict = grid_tools.create_log_lam_grid(wl_start=3000, wl_end=13000, min_vc=0.06/C.c_kms)
wl = wl_dict.pop("wl")
wl_dict.update({"type":"test spectrum"})
spec = grid_tools.LogLambdaSpectrum(wl, np.ones_like(wl), metadata=wl_dict)
instrument = TRES()
In [17]:
timeit -n 1 -r 1 spec.instrument_and_stellar_convolve(instrument, 10.)
In [ ]: