Often-times, the model spectrum is just not good enough to get an adequate fit using the normal Fitters.RVFitter. When that is the case, the automated flattening just doesn't work. Let's give this a shot without using models. The thing we want to minimize is the offset between the flux at one order and the flux in the adjacent order.
$$L = \sum_{i=1}^{N_{orders}-1} \left(\frac{O_i(\lambda)}{M(i, \lambda| \Theta)} - \frac{O_{i+1}(\lambda)}{M(i+1, \lambda | \theta)}\right)^2 $$Where $M(i, \lambda | \theta)$ is perhaps a 2d polynomial function of some order where $i$ is for the order number, $\lambda$ is the wavelength (or maybe pixel number?), and the $\theta$ are all the polynomial coefficients.
Try using Robust Linear Models from statsmodels. That might help it to ignore the big Balmer line. I will need to put the Chebyshev factors into a feature vector to use this...
Those might help a bit. How about a hybrid approach where I divide by a stellar model (or empirical spectrum of similar type as Adam suggests), and then do a full robust 2D fit.
I was able to sufficiently flatten a few of my spectra to get RV fits. I need to calibrate the RV though, so pretty much need to flatten every spectrum and fit the RV. This will take a while...
In [1]:
import SpecFlattener
import glob
import StarData
from astropy.io import fits
import SpectralTypeRelations
import logging
logger = logging.getLogger()
logger.setLevel(logging.INFO)
In [2]:
#hdf5_lib = '/media/ExtraSpace/Kurucz_FullGrid/CHIRON_grid_full.hdf5'
hdf5_lib = '/Volumes/DATADRIVE/Kurucz_Grid/HRS_grid_full.hdf5'
star_list = [f for f in glob.glob('../201*/H*-0.fits') if 'smoothed' not in f and
'uncorrected' not in f and 'flattened' not in f and 'B.fits' not in f and 'HRS' not in f]
print(len(star_list))
#star_list.index('../20131019/HIP_22913.fits')
for s in star_list:
print s
In [3]:
# Guess stellar properties
MS = SpectralTypeRelations.MainSequence()
def guess_teff_logg(fname):
header = fits.getheader(fname)
data = StarData.GetData(header['OBJECT'])
spt = data.spectype
teff = MS.Interpolate('Temperature', spt)
logg = 3.5 if 'I' in spt else 4.0
return teff, logg
In [4]:
teff, logg = guess_teff_logg(star_list[0])
print(teff, logg)
In [24]:
import HelperFunctions
orders = HelperFunctions.ReadExtensionFits(star_list[1])
import matplotlib.pyplot as plt
%matplotlib notebook
#nums = tuple(range(5, 16)) + tuple(range(18, 26))
nums = range(18, 35)
#for order in orders[4:25]:
good_orders = [o.copy() for i, o in enumerate(orders) if i in nums]
for order in good_orders:
plt.plot(order.x, order.y, 'k-', alpha=0.5)
In [25]:
reload(SpecFlattener)
print(len(nums))
output = SpecFlattener.flatten_spec(star_list[1], hdf5_lib, teff=teff, logg=logg, normalize_model=False,
first_order=4, last_order=15,
ordernums=nums)
final_orders, flattened, shifted_orders, mcf = output
In [27]:
%matplotlib notebook
for order in final_orders:
plt.plot(order.x, order.y, 'k-', alpha=0.5)
In [ ]: