In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import astropy.io.fits as pf # Pyfits is depracated
import astropy.constants as ct
import ppxf.miles_util as lib # MILES library embedded in ppxf
import ppxf.ppxf_util as util
import ppxf as ppxf_module
import os
import scipy.interpolate as sci
from ppxf.ppxf import ppxf
from ppxf.ppxf_util import log_rebin
In [2]:
sed_example = '../../GAMAII_SPEC/WHAN_NA_UVUP/G09_Y1_CX1_374_modified.fit'
In [3]:
sed = pf.open(sed_example)
In [4]:
header = sed[0].header
In [5]:
header
Out[5]:
In [6]:
flux_all = pf.getdata(sed_example)[0, 3:-3]
flux_err_all = pf.getdata(sed_example)[1, 3:-3]
In [7]:
idx_clean = np.where(np.logical_not(np.isnan(flux_all)))
In [8]:
flux = flux_all[idx_clean]
flux_err = flux_err_all[idx_clean]
In [9]:
wl_i = header['WMIN']
wl_step = header['CD1_1']
wl = np.arange(flux.size)*wl_step + wl_i
In [10]:
redshift = header['Z']
In [11]:
# %matplotlib notebook
In [12]:
sns.set_style("white")
plt.rcParams["axes.edgecolor"] = "0.15"
plt.rcParams["axes.linewidth"] = 1.
plt.subplots(1,1, figsize=(7,4))
plt.plot(wl,flux,'-')
plt.xlabel("Wavelength ($\mathrm{\AA}$)", fontsize=12)
plt.ylabel("Flux (10$^{-17}$ erg s$^{-1}$ cm$^2$ $\mathrm{\AA^{-1}}$)", fontsize=12)
plt.tick_params('both', labelsize=12)
plt.show()
In [13]:
wl_range_index = (wl>3540) & (wl<7409) # Limiting the range of wl given the library -- in this case MILES
In [14]:
ppxf_dir = os.path.dirname(os.path.realpath(ppxf_module.__file__))
In [15]:
ppxf_dir
Out[15]:
In [16]:
flux_pp = flux[wl_range_index]
wl_pp = wl[wl_range_index]
In [17]:
print (flux_pp)
In [18]:
print(ct.c)
In [19]:
print(ct.c.to('km/s'))
In [20]:
light_kms = float(str(ct.c.to('km/s')).split(' ')[0]) # selecting the numerical values only
In [21]:
vel_scale = light_kms*np.log(wl_pp[1]/wl_pp[0]) # velocity scale - Eq.8 Capellari 2017
In [22]:
fwhm_gama = 2. # MUST BE CHECKED!!!!! This is just a temporary value!!!
In [23]:
# TODO
# check the following information for GAMA:
# 1) Are the spectra corrected by the air? -- Do I need to correct them before runing ppxf?
# 2) What is the FWHM (espectral resolution) of GAMA spectra? - IDEM
# 3) Do we need gas templates?
In [24]:
miles_templates = ppxf_dir+'/miles_models/Mun1.30*.fits'
In [35]:
miles = lib.miles(miles_templates, vel_scale, fwhm_gama)
In [26]:
regularization_dimensions = miles.templates.shape[1:]
stellar_templates = miles.templates.reshape(miles.templates.shape[0], -1)
print (stellar_templates.shape)
In [27]:
regularization_err = 0.013 # Desired regularization error
In [28]:
wl_range = np.array([np.min(wl_pp), np.max(wl_pp)])/(1 + redshift)
In [29]:
dv = light_kms*(miles.log_lam_temp[0] - np.log(wl_pp[0])) # eq.(8) of Cappellari (2017)
vel = light_kms*np.log(1 + redshift) # eq.(8) of Cappellari (2017)
start = [vel, 180.] # (km/s), starting guess for [V, sigma]
In [30]:
# TODO: REBIN ? See error below.
In [31]:
new_shape = stellar_templates[:,0].shape
rebinning_function = sci.interp1d(wl_pp, flux_pp)
wl_pp_rebin = np.linspace(start=wl_pp.min(), stop=wl_pp.max(), num=int(new_shape[0]))
flux_pp_rebin = rebinning_function(wl_pp_rebin)
In [32]:
galaxy = flux_pp_rebin/np.median(flux_pp_rebin) # Normalize spectrum to avoid numerical issues
In [33]:
noise = np.full_like(galaxy, 0.01635) # Assume constant noise per pixel here
# MUST BE CHECKED!!!!! This is just a temporary value!!!
In [37]:
%matplotlib notebook
In [85]:
pp = ppxf(templates=stellar_templates, galaxy=galaxy, noise=noise, velscale=vel_scale, start=start, plot=True,
lam=wl_pp_rebin, moments=4, degree=250)
In [86]:
pp.bestfit
Out[86]:
In [87]:
pp.chi2
Out[87]:
In [88]:
pp.bias
Out[88]:
In [89]:
pp.clean
Out[89]:
In [90]:
pp.apoly
Out[90]:
In [91]:
pp.component
Out[91]:
In [92]:
pp.sol
Out[92]:
In [93]:
pp.ftol
Out[93]:
In [94]:
pp??
In [103]:
plt.subplot(2,1,1)
plt.plot(wl_pp_rebin, pp.galaxy, '-', linewidth=1., label="Spectrum")
plt.plot(wl_pp_rebin, pp.bestfit, '-', linewidth=1., label="Best Fit")
plt.legend(loc='best')
plt.xlabel("Wavelength ($\mathrm{\AA}$)")
plt.ylabel("Norm. flux")
plt.subplot(2,1,2)
plt.plot(wl_pp_rebin, (galaxy-pp.bestfit), '.', linewidth=1., label="Residues", alpha=0.8)
plt.axhline(y=0, c='black')
plt.legend(loc='best')
plt.xlabel("Wavelength ($\mathrm{\AA}$)")
plt.ylabel("Residues")
plt.savefig('./../Figs/ppxf_teste.png')
plt.show()
In [ ]: