Libraries


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

FIT or FITS file's


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]:
SIMPLE  =                    T / file does conform to FITS standard             
BITPIX  =                  -32 / number of bits per data pixel                  
NAXIS   =                    2 / number of data axes                            
NAXIS1  =                 4944 / length of data axis 1                          
NAXIS2  =                    5 / length of data axis 2                          
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
BSCALE  =                   1. / True_value = BSCALE * FITS_value + BZERO       
BZERO   =                   0. / True_value = BSCALE * FITS_value + BZERO       
DATE    = '2015-03-24T14:06:21' / file creation date (YYYY-MM-DDThh:mm:ss UT)   
ORIGIN  = 'GAMA    '           / Source of the data                             
OBSERVAT= 'SSO     '           / Name of observatory                            
LATITUDE=            -31.27704 / Telescope latitude in deg                      
LONGITUD=             149.0661 / Telescope longitude in deg                     
ALTITUDE=                1164. / Telescope height above sea level in m          
TELESCOP= 'AAT     '           / 3.9m Anglo-Australian Telescope                
INSTRUME= 'AAOMEGA-2dF'        / AAOmega/2dF multi-fibre spectrograph           
SPECTID = 'RD      '           / Spectrograph ID                                
GRATID  = '385R 580V'          / Disperser ID                                   
DICHROIC= 'X5700   '           / Dichroic name                                  
PLATE   =                    0 / 2dF plate number                               
RDNOISE =                 3.49 / Readout noise in e-                            
GAIN    =                1.799 / Readout gain in e-/ADU                         
FIBRE   =                  374 / Fibre number                                   
PIVOT   =                  374 / Pivot number                                   
X       =                79130 / Fibre button x position on plate               
Y       =               185274 / Fibre button y position on plate               
THETA   =               334.53 / Fibre angle in deg                             
N_EXP   =                    2 / Number of exposures                            
T_EXP   =                3000. / Total exposure time in s                       
N_NIGHTS=                    1 / Number of nights over which exposures were spre
DATE-OBS= '2008-03-31'         / UT date of observation                         
UTMJD   =     54556.3916837732 / Mean MJD of start of exposures                 
UTSTART = '09:10:29'           / UT of start of first exposure                  
UTEND   = '10:04:40'           / UT of end of last exposure                     
ZDSTART =                37.59 / Zenith distance at start of first exposure in d
ZDEND   =                34.15 / Zenith distance at end of last exposure in deg 
HASTART =               -17.49 / Hour angle at start of first exposure in deg   
HAEND   =                -5.49 / Hour angle at end of last exposure in deg      
HEL_VC  =               -23.65 / Heliocentric velocity correction in km/s       
CTYPE1  = 'Wavelength'         / Type of coordinate on axis 1                   
CUNIT1  = 'Angstrom'           / Units for axis 1                               
CRPIX1  =                2472. / WCS reference point                            
CDELT1  =       1.035900920889 / Wavelength at the ref pixel 2476 for STECKMAP  
CRVAL1  =              3727.79 / WCS reference value                            
CD1_1   =       1.035900920889 / WCS transformation matrix element              
CD1_2   =                   0. / WCS transformation matrix element              
CD2_1   =                   0. / WCS transformation matrix element              
CD2_2   =                   1. / WCS transformation matrix element              
ROW1    = 'Spectrum'           / Flux-calibrated spectrum in 10^-17 erg/s/cm^2/A
ROW2    = 'Error   '           / 1 sigma error spectrum                         
ROW3    = 'Spectrum_nocalib'   / Spectrum without flux calibration              
ROW4    = 'Error_nocalib'      / 1 sigma error spectrum (no flux calibration)   
ROW5    = 'Sky     '           / Sky spectrum (no flux calibration)             
FIELDID = 'G09_Y1_CX1'         / ID of GAMA AAT field                           
SPECID  = 'G09_Y1_CX1_374'     / ID of GAMA AAT spectrum                        
RA      =            133.96354 / RA of spectrum in deg (J2000)                  
DEC     =              2.79225 / DEC of spectrum in deg (J2000)                 
WMIN    =                3736. / Minimum wavelength in A                        
WMAX    =              8856.46 / Maximum wavelength in A                        
Z       =              0.21647 / Heliocentric redshift                          
NQ      =                    4 / Normalised redshift quality (use nQ>2 for scien
PROB    =                   1. / Probability that redshift is correct           
SN      =                12.81 / S/N as measured by runz                        
SN_ALT  =                 11.3 / Median S/N per pixel in 4500-6500 A band       
RFLUXFAC=           7.0767E-19 / Ratio of r-band petrosian to fibre flux        
CATAID  =               423115 / ID of matched GAMA object                      
OBJECT  =               423115 / ID of matched GAMA object                      
GAMANAME= 'GAMAJ085551.24+024732.1' / Name of matched GAMA object               
IC_FLAG =                 4104 / Catalogue source(s) of matched GAMA object     
DIST    =                   0. / Distance between spec and object in arcsec     
IS_SBEST=                    T / GAMA AAT spectrum with best redshift for this o

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()


pPXF config

Getting the Libraries' path -- those that already came with ppxf


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]:
'/home/mlldantas/anaconda3/envs/py3k6/lib/python3.6/site-packages/ppxf'

In [16]:
flux_pp = flux[wl_range_index]
wl_pp   = wl[wl_range_index]

In [17]:
print (flux_pp)


[-15.374364   -4.797919   -4.4835305 ...  20.24269    20.094873
  19.964918 ]

In [18]:
print(ct.c)


  Name   = Speed of light in vacuum
  Value  = 299792458.0
  Uncertainty  = 0.0
  Unit  = m / s
  Reference = CODATA 2014

In [19]:
print(ct.c.to('km/s'))


299792.458 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)


(2664, 150)

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)


Best Fit:       Vel     sigma        h3        h4
 comp. 0:     58976       583    -0.080    -0.013
chi2/DOF: 29.05
method = capfit ; Jac calls: 7 ; Func calls: 40 ; Status: 2
Nonzero Templates:  1  /  150

In [86]:
pp.bestfit


Out[86]:
array([-1.22329505, -0.37344859, -0.5042784 , ...,  1.59928294,
        1.60279836,  1.58855085])

In [87]:
pp.chi2


Out[87]:
29.05012921059821

In [88]:
pp.bias


Out[88]:
0

In [89]:
pp.clean


Out[89]:
False

In [90]:
pp.apoly


Out[90]:
array([-1.57303319, -0.7229255 , -0.8534935 , ...,  1.23876014,
        1.24252058,  1.22852125])

In [91]:
pp.component


Out[91]:
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

In [92]:
pp.sol


Out[92]:
array([ 5.89764502e+04,  5.83136952e+02, -8.03152732e-02, -1.33812228e-02])

In [93]:
pp.ftol


Out[93]:
0.0001

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 [ ]: