In [8]:
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from astropy.io import fits
import astropy.table as table
import matplotlib as mpl
import matplotlib.image as mpimg
import SPaCT
import glob as glob
import pandas as pd
from scipy import ndimage
import ppxf_util as util
from scipy.signal import resample
import warnings
warnings.filterwarnings("ignore", message="using a non-integer number instead of an integer will result in an error in the future")
#plt.rc('text', usetex=True)
#plt.rc('font', family='serif')
#plt.rcParams['figure.figsize'] = 8, 6
SNR = 10.
max_nonzero = 5
In [4]:
#first, munge all the possible templates to get characteristics
template_files = glob.glob('miles_models/Mun*.fits')
ssp_rows = []
from astropy.table import Table
import astropy
for template in template_files:
template = template.rstrip('.fits').split('/')[1]
spectral_range = template[0]
IMF_type = template[1:3]
IMF_slope = float(template[3:7])
Z = SPaCT.plusminus(template[8])*float(template[9:13])
T = float(template[14:])
#print template + ':', spectral_range, IMF_type, IMF_slope, Z, T
ssp_i = [template, spectral_range, IMF_type, IMF_slope, Z, T]
ssp_rows.append(ssp_i)
ssps = Table(map(list, zip(*ssp_rows)), names = ['template name', 'spectral range', 'IMF type', 'IMF slope', 'Z', 't'])
In [5]:
#assign weights to templates
#first pick which templates will be nonzero
nonzero_templates = np.arange(len(template_files))
np.random.shuffle(nonzero_templates)
template_weights = np.empty(len(template_files))
for i in range(len(template_weights)):
if i in nonzero_templates[:np.random.randint(1, max_nonzero)]:
#print i, 'True'
template_weights[i] = np.random.random()
else:
template_weights[i] = 0.
template_weights = template_weights/template_weights.sum()
In [6]:
#now load in each template as a row in an array
all_templates = np.empty((len(template_files), fits.open(template_files[0])[0].header['NAXIS1']))
for i in range(len(template_files)):
all_templates[i] = fits.open(template_files[i])[0].data
clean_spectrum = (template_weights * all_templates.T).T.sum(axis = 0)
#now construct an original wavelength array, and redshift it randomly
template_header = fits.open(template_files[0])[0].header
l0 = np.arange(template_header['CRVAL1'], template_header['CRVAL1'] + template_header['NAXIS1'] * template_header['CDELT1'], template_header['CDELT1'])
z = np.random.uniform(0.001, 0.020)
#dz is the (random) redshift we'll be off by when running it through SPaCT
dz = np.sign(np.random.rand() - 1) * np.random.uniform(0.01, 0.05) * z
#finally, redshift the spectrum
l = l0 * (1 + z)
#downsample spectrum from MILES to SparsePak resolution
res_factor = 1.4 / template_header['CDELT1']
#cut out values outside SparsePak's resolution
l = l[(l >= 3907.) & (l <= 3907 + 1.4 * 1934)]
clean_spectrum = clean_spectrum[(l >= 3907.) & (l <= 3907 + 1.4 * 1934.)]
clean_spectrum = resample(clean_spectrum, 1934)
#need to actually construct the new l array, since resampling directly somehow shits on this
l = np.linspace(3907., 3907. + 1.4 * 1934, 1934)
lamRange1 = np.array([3907, 3907 + 1.4 * 1934.])
FWHM_gal = 4.877 # FWHM from SparsePak simulator
# Quadratic sigma difference in pixels SSP --> instrument FWHM
# The formula below is rigorously valid if the shapes of the
# instrumental spectral profiles are well approximated by Gaussians.
FWHM_dif = np.sqrt(4.877**2 - 1.5**2)
sigma = FWHM_dif/2.355/template_header['CDELT1'] # Sigma difference in pixels
clean_spectrum = ndimage.gaussian_filter1d(clean_spectrum, sigma)
#need to de-redshift spectrum in order to use pPXF
lamRange1 = lamRange1 / (1 + z)
FWHM_gal = FWHM_gal / (1 + z)
galaxy, logLam1, velscale = util.log_rebin(lamRange1, clean_spectrum)
#now add noise
galaxy /= np.median(galaxy)
noise = np.random.normal(loc = 0.0, scale = 1./SNR, size = len(galaxy))
galaxy += noise * galaxy
noise_samples = np.random.normal(loc = 0.0, scale = 1./SNR, size = len(galaxy) * 8).reshape(8, len(galaxy))
In [9]:
plt.rcParams['figure.figsize'] = 16, 12
ifu = np.vstack((galaxy, noise_samples))
pp = SPaCT.SP_pPXF(ifu, fiber = 0, l_summ = (3907., 1.4, 1934), z = z + dz, verbose = True, noise_plots = True, fit_plots = True, edgefibers = [1, 2, 3, 4, 5, 6, 7, 8], age_lim = 20.)
plt.plot(pp.weights)
plt.plot(template_weights)
Out[9]:
In [ ]:
In [ ]: