Spectral Cube Modeling#
Use case: Extracting spatial-spectral features of interest from the cube of data and measuring their
attributes. Create line and PAH maps from the model parameters map.
Data: SL1 and SL2 spectral cubes generated using CUBISM, from Spitzer IRS spectral mapping observations of the center of nearby galaxy Messier 58.
Tools: astropy, specutils, dust-extinction, matplotlib.
Cross-intrument: NIRSpec, MIRI.
Documentation: This notebook is part of a STScI’s larger post-pipeline Data Analysis Tools Ecosystem.
Introduction#
Analysis of JWST spectral cubes requires extracting spatial-spectral features of interest and measuring their attributes. Here, we demonstrate 1-D spectral modeling applied to an individual spaxel, 900 spaxels in a cube, and spectra summed over spatial regions. Analysis of large JWST spectral data cubes can be computationally expensive, depending on number of spaxels modeled and number of free model parameters. This task can be made manageable by focusing on spatial sub-regions of interest or by minimizing model complexity. Sub-regions in the cube are selected by location and line ratios. The spaxels in these regions are then summed and modeled with a combination of continuum, emission line, polycyclic aromatic hydrocarbon (PAH), and silicate dust features. Best-fitting model spectra are fit using the lmfit package, utilizing a Levenberg-Marquardt least-squares method.
Imports#
time for timing
multiprocessing.Pool for parallel processing on cube spaxels
pickle for pickling and unpickling IO to multiprocessing pools
numpy for array processing and math
matplotlib.pyplot for plotting images and spectra
astropy.io for reading and writing FITS cubes and images
astropy.visualization to construct RGB images
lmfit Levenberg-Marquart fitting package for model-fitting
import time
import multiprocess as mp
from multiprocess import Pool
import pickle
import numpy as np
from numpy import exp, loadtxt, pi, sqrt, arange, interp
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from astropy.io import fits
from astropy.nddata import StdDevUncertainty
import astropy.units as u
from astropy.visualization import make_lupton_rgb
from astropy.visualization import quantity_support
from astropy.modeling import models
from specutils import Spectrum1D
import lmfit
from lmfit import Model, Minimizer, Parameters
from lmfit.model import save_modelresult
Definitions of functions#
Scrolling#
Limit scrolling in cell output.
%%javascript
IPython.OutputArea.auto_scroll_threshold=70;
Model Functions#
A variety of emission line, continuum and extinction models for spectral fitting.
def Modelcero(x):
"""Simple model of 0 to initialize some models"""
return 0*x
#powerlaw, zblackbody, gaussian and drude functions wrap existing astropy models.
#lmfit requires wavelength (x) to be passed as function parameter.
def powerlaw(x, c1, c2):
"""power law function (c1*nu^c2) normalized at 1 micron, with positive index convention"""
x0=1.0 #Normalize at 1 micron
c2=-c2
return models.PowerLaw1D.evaluate(x,c1,x0,c2)
def zblackbody(x, OnePlusZ, T, scale):
"""black body function"""
xs=x*u.um/OnePlusZ
bb=models.BlackBody(T*u.K, scale)
return bb(xs)
def gaussian(x, amplitude, xcen, std):
"""1-d gaussian"""
return models.Gaussian1D.evaluate(x,amplitude,xcen,std)
def drude(x, amplitude, peakx, frac_FWHM):
"""dust emission"""
FWHM=peakx*frac_FWHM
return models.Drude1D.evaluate(x,amplitude,peakx,FWHM)
#pahdust and sidust functions do not exist in astropy
def pahdust(x, OnePlusZ, amplitude_76, amplitude_113):
"""PAH dust emission"""
#Ionized PAH features
PAH_peakx= [ 5.27, 5.70, 6.22, 6.69, 7.42, 7.60, 7.85, 8.33, 8.61]
PAH_frac_FWHM= [0.034,0.035, 0.030,0.07, 0.126,0.044,0.053,0.050, 0.039]
PAH_rel_amplitude= [0.0, 0.0, 0.8, 0.0, 0.0, 1.0, 0.6, 0.0, 0.5]
PAH_amplitude=[]
for ampl in PAH_rel_amplitude: PAH_amplitude.append(amplitude_76*(ampl/PAH_rel_amplitude[5]))
#Neutral PAH features
PAH_peakx+= [10.68, 11.23,11.33]
PAH_peakx+= [11.99,12.62,12.69,13.48,14.04,14.19,15.90, 16.45,17.04,17.375,17.87,18.92]
PAH_frac_FWHM+=[0.020, 0.012,0.022]
PAH_frac_FWHM+=[0.045,0.042, 0.013,0.040,0.016,0.025,0.020,0.014,0.065, 0.012, 0.016,0.019]
PAH_rel_amplitude= [0.0, 0.6, 1.0]
PAH_rel_amplitude+=[0.2, 0.3, 0.1, 0.0, 0.1, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
for ampl in PAH_rel_amplitude: PAH_amplitude.append(amplitude_113*(ampl/PAH_rel_amplitude[2]))
pahflux=x-x
for peakx, frac_FWHM, ampl in zip(PAH_peakx,PAH_frac_FWHM,PAH_amplitude):
pahflux=pahflux+drude(x,ampl,peakx*OnePlusZ,frac_FWHM)
return pahflux #Amplitude unit
def sidust(x, T, amplitude):
"""Silicate dust emission"""
#Custom extinction curve built from weighted sum of three components:
#two Drude functions and an exponent 1.7 power-law.
d1=drude(x,0.80,10.0,0.25)
d2=drude(x,0.25,11.1,0.25)
#d3=drude(x,0.40,17.0,0.40) #Outside of wavelength range.
ext=d1+d2 #+d3
# Form linear combination of modified silicate and powerlaw.
beta=0.1
ext=(1. - beta)*ext + beta*(9.7/x)**1.7
si_em=amplitude*1.0E-6*3.97289E13/x**3/(exp(1.4387752E4/x/T)-1.)*ext
return si_em
# Enable tabulated extinction functions. Existing options in dust_extinction are insufficient
# at mid-IR wavelengths
def extinction_a(x, tau, a):
"""Extinction from table data"""
ext=np.interp(x, wave, a)
return exp(-tau*ext)
#Read extinction data from Chiar & Tielens (2005)
#Note, this is normalized to extinction in the K-band (1.004 at 2.14 um)
#agal is for Galactic Center, alocal is for local ISM extinction
Boxdata="https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/cube_fitting/"
wave, agal, alocal = np.loadtxt(Boxdata+'chiar+tielens_2005.dat', skiprows=14, usecols=(0, 1, 2), unpack=True)
#Compare to dust_extinction
#P92 is the only model that extends over the mid-IR range, but it is crude
#Normalized to extinction in the V-band
from dust_extinction.parameter_averages import F99
from dust_extinction.shapes import P92
a92 = P92()
#Extinction Comparison
tau=1.
xa=arange(5.,15.,0.01)
extin1=extinction_a(xa, tau, alocal)
extin2=extinction_a(xa, tau, agal)
lam=xa*u.um
scale=21.9 #Empirical scaling to adjust P2 curve to match Chiar & Tielens Local ISM curve
aa92=scale*a92(lam)
extin3 = exp(-tau*aa92)
with quantity_support():
plt.plot(lam,extin1,label='CT+05 Local ISM')
plt.plot(lam,extin2,label='CT+05 Gal. Cen.')
plt.plot(lam,extin3,label='P92 scaled')
plt.legend()
plt.show()
pl1=powerlaw(xa,1.0,-1.0)
Model Parameters#
Declare the parameters of model functions. Parameters are frozen by setting vary=’False’ in the lmfit setpar command. For now we are just fitting amplitudes of emission features. Fits with more variable parameters take longer.
def setpar(pars, name, value, vary, minus):
"""Set any parameter"""
pars[name].set(value,vary=vary,min=minus)
def gaussian_defpar(name, amplitude, mean, stddev, pars):
"""set the parameters for a gaussian"""
std=stddev
name=str(name)
xcen=mean*OneZ
setpar(pars,name+'_std',std,False,None)
#setpar(pars,name+'_std',std,True,0)
setpar(pars,name+'_xcen',xcen,False,None)
setpar(pars,name+'_amplitude',amplitude,True,0)
return #pars[name+'_std'].set(std,vary=False),pars[name+'_xcen'].set(xcen,vary=False),pars[name+'_amplitude'].set(amplitude)
def drude_defpar(name, amplitude, peakx, frac_FWHM, pars):
"""set the parameters for a drude"""
peakx=peakx*OneZ
setpar(pars,name+'_frac_FWHM',frac_FWHM,False,None)
#setpar(pars,name+'_frac_FWHM',frac_FWHM,True,0)
setpar(pars,name+'_peakx',peakx,False,None)
setpar(pars,name+'_amplitude',amplitude,True,0)
return #pars[name+'_frac_FWHM'].set(frac_FWHM,vary=False),pars[name+'_peakx'].set(peakx,vary=False),pars[name+'_amplitude'].set(amplitude)
def pahdust_defpar(name, OnePlusZ, amplitude_76, amplitude_113, pars):
"""set the parameters for pdr dust"""
setpar(pars,name+'_amplitude_76',amplitude_76,True,0)
setpar(pars,name+'_amplitude_113',amplitude_113,True,0)
setpar(pars,name+'_OnePlusZ',OnePlusZ,False,None)
return #pars[name+'_frac_FWHM'].set(frac_FWHM,vary=False),pars[name+'_peakx'].set(peakx,vary=False),pars[name+'_amplitude'].set(amplitude)
def sidust_defpar(name, T, amplitude, pars):
setpar(pars,name+'_amplitude',amplitude,True,0)
setpar(pars,name+'_T',T,False,0)
return
def gaussian_extractpars(prefix, result):
"""extract the parameters for a gaussian"""
amp = result.params[prefix+'amplitude'].value
cen = result.params[prefix+'xcen'].value
std = result.params[prefix+'std'].value
return amp, cen, std
def drude_extractpars(prefix, result):
"""extract the parameters for a drude"""
ampl = result.params[prefix+'amplitude'].value
peak = result.params[prefix+'peakx'].value
frac = result.params[prefix+'frac_FWHM'].value
return ampl, peak, frac
def pahdust_extractpars(prefix, result):
"""extract the parameters for pahdust model"""
ampl_76 = result.params[prefix+'amplitude_76'].value
ampl_113 = result.params[prefix+'amplitude_113'].value
return ampl_76, ampl_113
def sidust_extractpars(prefix, result):
"""extract the parameters for pahdust model"""
T = result.params[prefix+'T'].value
ampl = result.params[prefix+'amplitude'].value
return T,ampl
Line Flux from Model Parameters#
Helper functions to calculate the fluxes and widths of emission features from model parameters.
#NOTE that these are ANALYTIC estimates of line fluxes,
#NOT available in astropy modeling or specutils
def gaussianline_flux(amp,std,cen):
"""calculate the integrated flux for a gaussian line"""
#Units= amp:Jy cen:microns
c= 29979.2458*10**10 #c in microns/s
gaussianfactor=sqrt(2*pi)
return amp*gaussianfactor*c*std/(cen**2)*10**(-23) #erg s^{-1} cm^{-2}
def drudeline_flux(amp,frac,peak):
"""calculate the integrate flux for a drude line"""
#Units= amp:Jy peak:microns
c= 29979.2458*10**10 #c in microns/s
drudefactor=pi
return (amp*frac*drudefactor*c/((2*peak)))*10**(-23) #erg s^{-1} cm^{-2}
Extract and plot 1D spectrum and model components#
def extract_spec(a, b):
"""extract spec and errors from a,b coordinates"""
spec_pix=data_cube[:,a,b]
err_pix=error_cube[:,a,b]
return spec_pix, err_pix
def model_comps(model_result):
"""Evaluate model components"""
comps = model_result.eval_components()
plcomp=[]
h2comp=[]
pahcomp=[]
sidustcomp=[]
#print(comps.keys())
for key in comps.keys():
keyl=key.lower()
if keyl[0:3]=='pwl': plcomp.append(key)
if keyl[0]=='h': h2comp.append(key)
if keyl[0:3]=='pah': pahcomp.append(key)
if keyl[0:3]=='sid': sidustcomp.append(key)
plaw_model=x-x
h2_model=x-x
ion_model=x-x
pah_model=x-x
sidust_model=x-x
for comp in plcomp: plaw_model+=comps[comp]
for comp in h2comp: h2_model+=comps[comp]
for comp in pahcomp: pah_model+=comps[comp]
for comp in sidustcomp: sidust_model+=comps[comp]
atomiclines=['ArII','SIV','NeII']
for comp in atomiclines:
ion_model+=comps[comp+'_']
return([plaw_model,h2_model,ion_model,pah_model,sidust_model],['PL','H2','Ions','PAHs','Si Dust'])
def plot_fit(x, spec, specerr, model_result):
"""plot spectrum, model components, and residual"""
#Model Results
fit_model = model_result.best_fit
fit_residual = spec - fit_model
#Evaluate model components
mod_comps,mod_labels=model_comps(model_result)
#Spec-1D object
spec1 = Spectrum1D(spectral_axis=x* u.um, flux=spec*u.Jy, uncertainty=StdDevUncertainty(specerr))
with quantity_support():
f, ax = plt.subplots()
#ax.step, ax.plot, ax.scatter all work, but ax.scatter doesn't autscale correctly
ax.plot(spec1.spectral_axis, spec1.flux,label='Data', color='k')
#ax.set_xlabel(r"$Wavelength\ (\mu m)$")
ax.grid(linestyle=':')
#plt.errorbar only works on unitless data
#ax.errorbar(spec1.spectral_axis, spec1.flux, yerr=specerr,label='Data', color='k')
#ax.plot(spec1.spectral_axis, spec1.flux, yerr=specerr,label='Data', color='k')
plt.show()
#Plot results
fig1=plt.figure(figsize=(15, 6), dpi= 150)
frame1=fig1.add_axes((.1,.3,.8,.6))
frame1.set_xticklabels([])
plt.ylabel(r"$Flux\ (Jy)$")
plt.grid(linestyle=':')
plt.errorbar(x, spec, yerr=specerr,label='Data', color='k')
plt.plot(x, fit_model,label='Model',c='red')
mod_colors=['brown','g','b','magenta','orange']
for mcomp, mlabel, mcolr in zip(mod_comps,mod_labels,mod_colors):
plt.plot(x, mcomp, label=mlabel, c=mcolr)
plt.legend()
frame2=fig1.add_axes((.1,.1,.8,.2))
plt.plot(x, 0.*x, '-',c='r')
plt.xlabel(r"$Wavelength\ (\mu m)$")
plt.grid(linestyle=':')
plt.errorbar(x,fit_residual,yerr=specerr,c='k')
plt.show()
return
Input data#
With the present lack of JWST flight data, we utilize the SL1 and SL2 spectral cubes generated using CUBISM, from Spitzer IRS spectral mapping observations of the center of nearby galaxy Messier 58. This galaxy has a radio-loud Seyfert nucleus that is exciting strong H2 emission via shocks. There is also significant star-formation, revealed by PAH emission features.
# Spitzer IRS (CUBISM) cube loader does not exist in specutils
#Target
targname='M58'
#Redshift
z=0.005060
OneZ=1.+z
#Download and open the data cubes and their uncertainties
BoxPath="https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/cube_fitting/"
cubeSL1 = fits.open(BoxPath+targname+'_SL1_cube.fits')
errorsSL1 = fits.open(BoxPath+targname+'_SL1_cube_unc.fits')
cubeSL2 = fits.open(BoxPath+targname+'_SL2_cube.fits')
errorsSL2 = fits.open(BoxPath+targname+'_SL2_cube_unc.fits')
#Cube Info and Headers
cubeSL1.info()
cubeSL2.info()
hdr1 = cubeSL1[0].header
er_hdr1=errorsSL1[0].header
hdr2 = cubeSL2[0].header
er_hdr2=errorsSL2[0].header
#print(repr(hdr1))
#Flux Data
data_cube1 = cubeSL1[0].data
error_cube1 = errorsSL1[0].data
data_cube2 = cubeSL2[0].data
error_cube2 = errorsSL2[0].data
#Wavelength Data
xwave1 = cubeSL1[1].data
xwave1 = xwave1.field(0)[0]
xwave2 = cubeSL2[1].data
xwave2 = xwave2.field(0)[0]
x=[]
for line in xwave2:
x.append(float(line))
for line in xwave1:
x.append(float(line))
x=np.array(x)
#Change the units from MJy/sr to Jy/pix
correct_unit1=abs(hdr1['CDELT1'])*abs(hdr1['CDELT2'])*0.0003046118*(10**6)
data_cube1=data_cube1*correct_unit1
error_cube1=error_cube1*correct_unit1
correct_unit2=abs(hdr2['CDELT1'])*abs(hdr2['CDELT2'])*0.0003046118*(10**6)
data_cube2=data_cube2*correct_unit2
error_cube2=error_cube2*correct_unit2
#Concatenate SL1 and SL2
data_cube=np.concatenate((data_cube2,data_cube1),axis=0)
error_cube=np.concatenate((error_cube2,error_cube1),axis=0)
#Reorder x and data_cube
xcor = x.argsort()
data_cube = data_cube[xcor]
error_cube = error_cube[xcor]
x=np.sort(x)
#Cube dimensions and trimming
xsize, ysize, zsize = data_cube.shape
ytrim=0; ysize=ysize-ytrim
ztrim=3; zsize=zsize-ztrim
print('Trimmed cube dimensions:', xsize, ysize, zsize)
#Collapsed 2D image
cube_2dflux=np.sum(data_cube,axis=0)[0:ysize,0:zsize]
#Collapsed 1D spectrum
cube_1dflux=np.zeros(xsize,)
for a in arange(0,ysize):
for b in arange(0,zsize):
spec_pix,err_pix=extract_spec(a,b)
cube_1dflux=cube_1dflux+spec_pix
f,(ax1,ax2)=plt.subplots(1,2, figsize=(10,5))
ax1.set_title('log Flux (Sum over Wavelength)')
ax1.imshow(cube_2dflux, origin='lower', cmap='gray', norm=LogNorm())
ax2.set_title('Flux (Sum over Spaxels)')
ax2.plot(x,cube_1dflux)
ax2.set_xlabel(r"Wavelength $(\mu m)$")
ax2.set_ylabel('Flux (Jy)')
plt.show()
Model Parameter Starting Values#
Set starting parameters for spectral model components. Use a PAH dust model with only two free parameters (7.6 and 11.3 micron PAH flux). The rest of the PAH ratios are fixed by this ‘pahdust’ model, so that fitting the whole cube does not take so long. Later, we will fit summed region spectra with all PAH amplitudes free.
#Power Law
p_law = Model(powerlaw,prefix='pwl_')
generic_pars = p_law.make_params()
setpar(generic_pars,'pwl_c1',0.0002,True,None)
setpar(generic_pars,'pwl_c2',0.0015,True,None)
#PAH dust model with fixed PAH feature ratios
prefix='pahdust1_'
pahs = Model(pahdust,prefix=prefix)
generic_pars.update(pahs.make_params())
pahdust_defpar('pahdust1',1.00506,0.003,0.005,generic_pars)
#Silicate dust emission
prefix='sidust1_'
silicate= Model(sidust,prefix=prefix)
generic_pars.update(silicate.make_params())
sidust_defpar('sidust1', 190.,1.E-4,generic_pars)
#Gaussians with wavelengths from Smith et al 2007
#H2 pure-rotational lines
h2lines=['55','61','69','80','96','122']
h2wavel=[5.511,6.109,6.909,8.026,9.665,12.278]
amp_guess=0.05
std_guess=0.04
gaussH2=Model(Modelcero)
for line in h2lines:
prefix='h'+line+'_'
gaussH2+=Model(gaussian,prefix=prefix)
generic_pars.update(gaussH2.make_params())
for h2lin, h2w in zip(h2lines,h2wavel):
gaussian_defpar('h'+h2lin, amp_guess, h2w, std_guess,generic_pars)
#Atomic lines
atomiclines=['ArII','SIV','NeII']
atomicwavel=[6.985,10.511,12.813]
amp_guess=0.05
std_guess=0.05
gaussAt=Model(Modelcero)
for line in atomiclines:
prefix = line+'_'
gaussAt+=Model(gaussian,prefix=prefix)
generic_pars.update(gaussAt.make_params())
for atlin, atw in zip(atomiclines,atomicwavel):
gaussian_defpar(atlin,amp_guess, atw, std_guess,generic_pars)
Fit a Representative Spaxel#
Fit a single spaxel in the cube and use the fit parameters to adjust the generic model starting parameters.
#Select spaxel
spec=data_cube[:,8,16]
specerr=error_cube[:,8,16]
#Composite model (powerlaw + PAHs + H2 + Atomic lines, no extinction)
generic_mod = p_law + silicate + pahs + gaussH2 + gaussAt
print(generic_mod.param_names)
#Fit and plot
spax_result = generic_mod.fit(spec, generic_pars, weights=1./specerr, x=x)
plot_fit(x, spec, specerr, spax_result)
#Save results
save_modelresult(spax_result, 'OnePointResult.sav')
#Update generic model starting parameters
for par in spax_result.params:
value=spax_result.params[par].value
vary=spax_result.params[par].vary
minus=spax_result.params[par].min
if minus==float("-inf"):
minus=None
generic_pars[par].set(value,vary=vary,min=minus)
Plot of observed flux density (black), model fit (red), and residuals (bottom panel). This spaxel has strong H2 emission at 9.6 microns and 8 micron PAH and [Ne II] 12.8 micron emission from star-forming regions.
Fit the Entire Cube#
Identify NaNs, then create list of models (based on generic model), one entry per NaN-free spaxel. Next launch the multiprocessing Pool, limited to number of available cores minus one. Fitting the entire cube (900 spaxels) takes about 2 minutes to fit with 7 processes running in parallel.
cube_res = mp.Manager().dict()
#The fit method is redefined as a top-level function to make it pickle-able for multiprocessing.Pool
def modfit(spaxel, generic_pars, weights, x):
modf=generic_mod.fit(spaxel, generic_pars, weights=1./spaxerr, x=x)
return modf.dumps()
#Helper to either assign a NaN result or fit a single non-NaN spaxel
def fit_point(args):
a,b,spaxel,spaxerr,x,generic_pars,nancheck = args
# print(a, b, len(spaxel), mp.current_process().name)
if nancheck == False:
res_point=modfit(spaxel, generic_pars, weights=1./spaxerr, x=x)
cube_res[(a,b)] = res_point
elif nancheck == True:
cube_res[(a,b)] = float('nan')
return
#Extraction region dimensions (trimmed cube)
astart=0; aend=ysize
bstart=0; bend=zsize
#Check for NaNs, select data, and set model start parameters for each spaxel
pooldata=[]
for a in np.arange(astart,aend):
for b in np.arange(bstart,bend):
nancheck = False
for point in data_cube[:,a,b]:
if np.isnan(point) == True or point == 0:
nancheck = True
#Copy spaxel data vectors
spaxel = 1.0*np.array(data_cube[:,a,b])
spaxerr= 1.0*np.array(error_cube[:,a,b])
wavelen= 1.0*np.array(x)
pooldata.append([a,b,spaxel,spaxerr,wavelen,generic_pars,nancheck])
#Launch Multiprocessing Pool
start_time = time.time()
if __name__ == '__main__':
with Pool(mp.cpu_count() - 1) as pool:
pool.map(fit_point,pooldata)
print('Time count')
print("--- %s seconds ---" % (time.time() - start_time))
#Save model fit
with open('MPF_cube_result_model.sav', 'wb') as fp:
pickle.dump(list(cube_res.items()), fp)
#Load model fit to full cube
cube_res = mp.Manager().dict()
with open ('MPF_cube_result_model.sav', 'rb') as fp:
restore_cube = pickle.load(fp)
for line in restore_cube:
pos=line[0]
cube_res[pos] = line[1]
Total flux and Reduced Chi^2 Maps#
Display total flux image (observed flux cube collapsed along wavelength direction) and reduced Chi^2 image for model fit in square sub-region.
#Sum observed flux in each spaxel
cube_tflux=np.sum(data_cube,axis=0)[0:ysize,0:zsize]
#Sum model flux in each spaxel
tcube_chi=np.zeros((ysize, zsize))
tcube_modflux=np.zeros((ysize, zsize))
funcdefs = {}
modres=spax_result
for pos, dumpval in cube_res.items():
#for pos, dumpval in cube_res:
if str(cube_res[pos])=='nan':
tcube_modflux[pos]=float('nan')
tcube_chi[pos]=float('nan')
else:
result = modres.loads(dumpval, funcdefs=funcdefs)
tcube_modflux[pos]=np.sum(result.best_fit)
tcube_chi[pos]=result.redchi
f,(ax1,ax2,ax3)=plt.subplots(1,3, figsize=(10,15))
ax1.set_title('log Observed Flux')
ax1.imshow(cube_tflux, origin='lower', cmap='gray', norm=LogNorm())
ax2.set_title('log Model Flux')
ax2.imshow(tcube_modflux, origin='lower', cmap='gray', norm=LogNorm())
ax3.set_title('Reduced Chi Square')
ax3.imshow(cube_tflux-tcube_modflux, origin='lower', cmap='hot')
plt.show()
Collapsed cube maps, showing observed flux, model flux, and Chi^2/DF. White regions with NaNs in the input data cube also have NaNs in the output model cube along the left and top edges and in column 26. Residuals in the nucleus may either be the result of inadequate PSF sampling or mis-fit AGN model.
Line and PAH feature flux maps.#
Functions to extract line and PAH maps from model parameter cube#
def line_map(line, cube_res, modres):
"""Create total emission line flux map for Gaussian emission line model component"""
line_flux=np.zeros((35, 35))
for pos, dumpval in cube_res.items():
if str(cube_res[pos])=='nan':
line_flux[pos]=float('nan')
else:
result = modres.loads(dumpval, funcdefs=funcdefs)
prefix = line+'_'
amp, cen, std = gaussian_extractpars(prefix,result)
line_flux[pos]=gaussianline_flux(amp,std,cen)
return line_flux
def pahdust_maps(pahdust_comp, cube_res, modres):
"""Create total emission line flux map for Gaussian emission line model component"""
pah_ion_flux=np.zeros((35, 35))
pah_neutral_flux=np.zeros((35, 35))
for pos, dumpval in cube_res.items():
if str(cube_res[pos])=='nan':
pah_ion_flux[pos]=float('nan')
pah_neutral_flux[pos]=float('nan')
else:
result = modres.loads(dumpval, funcdefs=funcdefs)
prefix = pahdust_comp+'_'
ampl_76, ampl_113 = pahdust_extractpars(prefix,result)
pah_76_flux = drudeline_flux(ampl_76,0.044,7.60)
pah_113_flux = drudeline_flux(ampl_113,0.032,11.33)
pah_ion_flux[pos]=pah_76_flux
pah_neutral_flux[pos]=pah_113_flux
return pah_ion_flux, pah_neutral_flux
def sidust_map(sidust_comp, cube_res, modres):
"""Create silicate dust emission map"""
sidust_flux=np.zeros((35, 35))
for pos, dumpval in cube_res.items():
if str(cube_res[pos])=='nan':
sidust_flux[pos]=float('nan')
else:
result = modres.loads(dumpval, funcdefs=funcdefs)
prefix = sidust_comp+'_'
t,si_ampl = sidust_extractpars(prefix,result)
si_norm=1.0
sidust_flux[pos]=si_ampl*si_norm
return sidust_flux
H2 9.6 micron, [Ne II] 12.8 micron, PAH 7.6, PAH 11.3, and silicate dust emission maps from model#
h2_flux=line_map('h96', cube_res, modres)
ne_flux=line_map('NeII', cube_res, modres)
pah_ion_flux, pah_neutral_flux = pahdust_maps('pahdust1', cube_res, modres)
sidust_flux=sidust_map('sidust1', cube_res, modres)
f,(ax1,ax2,ax3,ax4,ax5)=plt.subplots(1,5, figsize=(10,20))
ax1.set_title('[Ne II] 12.8 um')
ax1.imshow(ne_flux, origin='lower')
ax2.set_title('H2 9.6 um')
ax2.imshow(h2_flux, origin='lower')
ax3.set_title('PAH 7.6 um')
ax3.imshow(pah_ion_flux, origin='lower')
ax4.set_title('PAH 11.3 um')
ax4.imshow(pah_neutral_flux, origin='lower')
ax5.set_title('Si dust em.')
ax5.imshow(sidust_flux, origin='lower')
plt.show()
#Make RGB images
stretch1=0.000000000000004
stretch2=0.000000000000005
m58_neii_pah0_h2 = make_lupton_rgb(2.5*ne_flux, h2_flux, pah_neutral_flux, minimum=0, Q=5, stretch=stretch1, filename="m58_NeII_H96_PAH113.png")
m58_pah_h2 = make_lupton_rgb(pah_neutral_flux, h2_flux, 1.2*pah_ion_flux, minimum=0, Q=5, stretch=stretch2, filename="m58_PAH113_H96_PAH76.png")
#Two-panel figure
f,(ax1,ax2)=plt.subplots(1,2, figsize=(10,15))
ax1.axis('off')
ax1.imshow(m58_neii_pah0_h2,origin='lower')
ax2.axis('off')
ax2.imshow(m58_pah_h2,origin='lower')
plt.show()
3-color RGB images. Top: individual feature maps. Bottom left: r = [Ne II] 12.8 micron, g = H2 S(3) 9.6 micron, b = PAH 11.3 micron. Ionized atomic gas in the active galactic nucleus shows up as red, shocked regions are green, and dust emission is blue. Bottom right: r = PAH 11.3 micron, g = H2 S(3) 9.6 micron, b = PAH 7.6 micron. Neutral PAH emission shows up as red, ionized PAH emission from PDRs in star-forming regions is blue, and shocked regions are green.
Extract and model cube sub-regions#
Make spectral region masks#
Define spectral extraction regions based on spatial location, line flux, and line or feature ratios.
#Nucleus extraction region, defined by 5x5 square
nuc_mask=np.zeros((35, 35))
for a in arange(-2,3):
for b in arange(-2,3):
nuc_mask[14-a,17-b]=int(1)
#Spiral arm extraction region
arm_mask=np.where(((pah_neutral_flux>=0.155e-14) &((pah_ion_flux/pah_neutral_flux)>0.8) & (nuc_mask==0)),1,0)
#Inner disk extraction region
pah_mask=np.where(((pah_neutral_flux>=0.155e-14) & (h2_flux<1e-15) & (nuc_mask==0)),1,0) - arm_mask
pah_mask[:,0:3]=int(0) #Trim edges of cube
pah_mask[28:,:]=int(0)
#Regions of strongest H2 emission
h2s_mask=np.where(((h2_flux>=1e-15)&((h2_flux/pah_neutral_flux)>0.7) & (nuc_mask==0)),1,0)
#H2 extraction region
h2_mask=np.where(((h2_flux>=1e-15) & (nuc_mask==0)),1,0)-h2s_mask
#M58 combined central regions mask (excluding spiral arm)
m58_mask=nuc_mask + h2_mask*2 + h2s_mask*3 + pah_mask*4 +arm_mask*5
f,((ax1,ax2,ax3,ax4,ax5,ax6))=plt.subplots(1,6,figsize=(15,10))
ax1.set_title('Nucleus')
ax1.imshow(nuc_mask,origin='lower')
ax2.set_title('Spiral Arm')
ax2.imshow(arm_mask,origin='lower')
ax3.set_title('Inner Disk')
ax3.imshow(pah_mask,origin='lower')
ax4.set_title('H2')
ax4.imshow(h2_mask,origin='lower')
ax5.set_title('H2-strongest')
ax5.imshow(h2s_mask,origin='lower')
ax6.set_title('All regions')
ax6.imshow(m58_mask,origin='lower')
plt.show()
Extraction region masks based on spatial and emission line or PAH feature strengths.
Extract spectra from mask regions#
Spectra summed over five different regions: Spiral Arm, Nucleus, PAH, H2, and Strongest H2
nuc_spec=np.zeros(xsize,)
nuc_err=np.zeros(xsize,)
pah_spec=np.zeros(xsize,)
pah_err=np.zeros(xsize,)
h2_spec=np.zeros(xsize,)
h2_err=np.zeros(xsize,)
arm_spec=np.zeros(xsize,)
arm_err=np.zeros(xsize,)
h2s_spec=np.zeros(xsize,)
h2s_err=np.zeros(xsize,)
for a in arange(0,ysize):
for b in arange(0,zsize):
if nuc_mask[a,b]==1:
spec_pix,err_pix=extract_spec(a,b)
nuc_spec=nuc_spec+spec_pix
nuc_err=nuc_err+err_pix**2
if pah_mask[a,b]==1:
spec_pix,err_pix=extract_spec(a,b)
pah_spec=pah_spec+spec_pix
pah_err=pah_err+err_pix**2
if h2_mask[a,b]==1:
spec_pix,err_pix=extract_spec(a,b)
h2_spec=h2_spec+spec_pix
h2_err=h2_err+err_pix**2
if arm_mask[a,b]==1:
spec_pix,err_pix=extract_spec(a,b)
arm_spec=arm_spec+spec_pix
arm_err=arm_err+err_pix**2
if h2s_mask[a,b]==1:
spec_pix,err_pix=extract_spec(a,b)
h2s_spec=h2s_spec+spec_pix
h2s_err=h2s_err+err_pix**2
nuc_err=sqrt(nuc_err)
pah_err=sqrt(pah_err)
h2_err=sqrt(h2_err)
arm_err=sqrt(arm_err)
h2s_err=sqrt(h2s_err)
f,((ax1,ax2,ax3),(ax4,ax5,ax6))=plt.subplots(2,3,figsize=(15,10))
ax1.set_title('Spiral Arm')
ax1.errorbar(x,arm_spec,yerr=arm_err, c='gold')
ax1.set_ylabel('Flux (Jy)')
ax2.set_title('Nucleus')
ax2.errorbar(x,nuc_spec,yerr=nuc_err, c='navy')
ax3.set_title('Inner Disk')
ax3.errorbar(x,pah_spec,yerr=pah_err, c='limegreen')
ax4.set_title('H2')
ax4.errorbar(x,h2_spec,yerr=h2_err, c='steelblue')
ax4.set_xlabel(r"Wavelength $(\mu m)$")
ax4.set_ylabel('Flux (Jy)')
ax5.set_title('All regions')
ax5.imshow(m58_mask,origin='lower')
ax5.set_xlabel('Pixel')
ax6.set_title('H2-Strongest')
ax6.errorbar(x,h2s_spec,yerr=h2s_err,c='teal')
ax6.set_xlabel(r"Wavelength $(\mu m)$")
plt.show()
Save spectra to files#
#np.savetxt(targname+'_spec_nuc.dat', np.column_stack((x, nuc_spec, nuc_err)),header='Wavelength[microns] Flux_Nucleus[Jy] Err_Flux_Nucleus')
#np.savetxt(targname+'_spec_h2.dat', np.column_stack((x, h2_spec, h2_err)), header='Wavelength[microns] Flux_H2[Jy] Err_H2')
#np.savetxt(targname+'_spec_h2s.dat', np.column_stack((x, h2_spec, h2_err)), header='Wavelength[microns] Flux_H2s[Jy] Err_H2s')
#np.savetxt(targname+'_spec_pah.dat', np.column_stack((x, pah_spec, pah_err)), header='Wavelength[microns] Flux_PAH[Jy] Err_PAH')
#np.savetxt(targname+'_spec_arm.dat', np.column_stack((x, arm_spec, arm_err)), header='Wavelength[microns] Flux_Arm[Jy] Err_Arm')
Fit spectra of regions#
Generic model components#
Model 14 individual PAH components with all amplitudes free. It would be computationally expensive to do this for the full cube, but only takes a couple of seconds for the summed region spectra.
#Power Law
p_law = Model(powerlaw,prefix='pwl_')
generic_pars = p_law.make_params() #Generate the parameters
setpar(generic_pars,'pwl_c1',0.005,True,None)
setpar(generic_pars,'pwl_c2',0.0015,True,None)
#Individual PAHs
drudes=Model(Modelcero)
allpahlines=['52','62','74','76','78','83','86','112','113','119','126','134','140','141']
for line in allpahlines:
prefix='PAH'+line+'_'
drudes+=Model(drude,prefix=prefix)
generic_pars.update(drudes.make_params())
drude_defpar('PAH52', 0.05, 5.27, 0.034,generic_pars)
drude_defpar('PAH62', 0.05, 6.22, 0.030,generic_pars)
drude_defpar('PAH74', 0.05, 7.42, 0.126,generic_pars)
drude_defpar('PAH76', 0.05, 7.60, 0.044,generic_pars)
drude_defpar('PAH78', 0.5, 7.85, 0.053,generic_pars)
drude_defpar('PAH83', 0.5, 8.33, 0.05,generic_pars)
drude_defpar('PAH86', 0.5, 8.61, 0.039,generic_pars)
drude_defpar('PAH112', 0.5, 11.23, 0.012,generic_pars)
drude_defpar('PAH113', 0.5, 11.33, 0.032,generic_pars)
drude_defpar('PAH119', 0.5, 11.99, 0.045,generic_pars)
drude_defpar('PAH126', 0.5, 12.62, 0.042,generic_pars)
drude_defpar('PAH134', 0.5, 13.48, 0.04,generic_pars)
drude_defpar('PAH140', 0.5, 14.04, 0.016,generic_pars)
drude_defpar('PAH141', 0.5, 14.19, 0.025,generic_pars)
#Silicate dust emission
prefix='sidust1_'
silicate= Model(sidust,prefix=prefix)
generic_pars.update(silicate.make_params())
sidust_defpar('sidust1', 190.,1.E-4,generic_pars)
#Gaussians with FWHM and X central from Smith et al 2007
#H2 molecular
h2lines=['55','61','69','80','96','122']
gaussH2=Model(Modelcero)
for line in h2lines:
prefix='h'+line+'_'
gaussH2+=Model(gaussian,prefix=prefix)
generic_pars.update(gaussH2.make_params())
gaussian_defpar('h55',0.05, 5.511, 0.02,generic_pars)
gaussian_defpar('h61',0.05, 6.109, 0.02,generic_pars)
gaussian_defpar('h69',0.05, 6.909, 0.02,generic_pars)
gaussian_defpar('h80',0.05, 8.026, 0.04,generic_pars)
gaussian_defpar('h96',0.05, 9.665, 0.04,generic_pars)
gaussian_defpar('h122',0.05, 12.278, 0.04,generic_pars)
#Atomic lines
atomiclines=['ArII','SIV','NeII']
gaussAt=Model(Modelcero)
for line in atomiclines:
prefix = line+'_'
gaussAt+=Model(gaussian,prefix=prefix)
generic_pars.update(gaussAt.make_params())
gaussian_defpar('ArII',0.05, 6.985, 0.02,generic_pars)
gaussian_defpar('SIV',0.05, 10.511, 0.04,generic_pars)
gaussian_defpar('NeII',0.05, 12.813, 0.04,generic_pars)
Spiral arm spectral fit#
#Spiral arm spectrum
spec=arm_spec
specerr=arm_err
#Model (powerlaw + individual PAHs + H2 and Atomic lines, no extinction)
model = p_law + drudes + gaussAt + gaussH2
#Adjust power law starting parameters
setpar(generic_pars,'pwl_c1',0.01,True,None)
setpar(generic_pars,'pwl_c2',0.001,True,None)
#Model Results
arm_result = model.fit(spec, generic_pars, weights=1./specerr, x=x)
arm_model = arm_result.best_fit
arm_residual = spec - arm_model
arm_modelcomps,arm_modlabels=model_comps(arm_result)
#Plot
#plot_fit(x, spec, specerr, arm_result)
Model fit to integrated spectrum of spiral arm region.
H2-strongest region spectral fit#
#H2 region spectrum
spec=h2s_spec
specerr=h2s_err
#Model (powerlaw + individual PAHs + H2 and Atomic lines, no extinction)
model = p_law + drudes + gaussAt + gaussH2
#Adjust power law starting parameters
setpar(generic_pars,'pwl_c1',0.002,True,None)
setpar(generic_pars,'pwl_c2',0.0015,True,None)
#Model Results
h2s_result = model.fit(spec, generic_pars, weights=1./specerr, x=x)
h2s_model = h2s_result.best_fit
h2s_residual = spec - h2s_model
h2s_modelcomps,h2s_modlabels=model_comps(h2s_result)
#Plot
#plot_fit(x, spec, specerr, h2s_result)
Model fit to integrated spectrum of H2-strongest region.
Nucleus spectral fit#
#Nucleus spectrum
spec=nuc_spec
specerr=nuc_err
#Model (powerlaw + AGN silicate emission + individual PAHs + H2 and Atomic lines, no extinction)
model = p_law + silicate + drudes + gaussH2 + gaussAt
#Adjust power law and silicate dust starting parameters
setpar(generic_pars,'pwl_c1',0.005,True,None)
setpar(generic_pars,'pwl_c2',0.01,True,None)
sidust_defpar('sidust1', 190.,2.E-4,generic_pars)
#Model Results
nuc_result = model.fit(spec, generic_pars, weights=1./specerr, x=x)
nuc_model = nuc_result.best_fit
nuc_residual = spec - nuc_model
nuc_modelcomps,nuc_modlabels=model_comps(nuc_result)
#Plot
plot_fit(x, spec, specerr, nuc_result)
Model fit to nucleus, including fixed-temperature AGN silicate emission.
nuc_result
#result.fit_report(show_correl=False)
Fit parameters for spectrum of nucleus
Summary Plot#
f,((ax1,ax2,ax3),(ax4,ax5,ax6))=plt.subplots(2,3,figsize=(15,10))
ax1.set_title('Spiral Arm')
ax1.set_ylabel('Flux (Jy)')
ax1.errorbar(x,arm_spec,yerr=arm_err, c='gold')
mod_colors=['brown','g','b','magenta','orange']
for mcomp,mlabl,mcolr in zip(arm_modelcomps,arm_modlabels,mod_colors):
ax1.plot(x,mcomp,label=mlabl,c=mcolr)
ax2.set_title('Nucleus')
ax2.errorbar(x,nuc_spec,yerr=nuc_err, c='navy')
#ax2.plot(x, fit_model,label='Model',c='red')
for mcomp,mlabl,mcolr in zip(nuc_modelcomps,nuc_modlabels,mod_colors):
ax2.plot(x,mcomp,label=mlabl,c=mcolr)
ax3.set_title('Inner Disk')
ax3.errorbar(x,pah_spec,yerr=pah_err, c='limegreen')
ax4.set_title('H2')
ax4.set_xlabel(r"Wavelength $(\mu m)$")
ax4.set_ylabel('Flux (Jy)')
ax4.errorbar(x,h2_spec,yerr=h2_err, c='steelblue')
ax5.set_title('All regions')
ax5.set_xlabel('Pixel')
ax5.imshow(m58_mask,origin='lower')
ax6.set_title('H2-Strongest')
ax6.set_xlabel(r"Wavelength $(\mu m)$")
ax6.errorbar(x,h2s_spec,yerr=h2s_err,c='teal')
for mcomp,mlabl,mcolr in zip(h2s_modelcomps,arm_modlabels,mod_colors):
ax6.plot(x,mcomp,label=mlabl,c=mcolr)
plt.show()
Notebook created by Patrick Ogle and Ivan Lopez.