Spectral Cube Modeling

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)
../../_images/4d53835931bf49df5ba8577dc57c4a59f0ab7d934be56439bdeb6f2eebb879dc.png

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

Filename: /home/runner/.astropy/cache/download/url/06853bf85036af1a6d55c1fb3a08260f/contents
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     286   (35, 30, 117)   float32   
  1  WCS-TAB       1 BinTableHDU     13   1R x 1C   [117E]   
Filename: /home/runner/.astropy/cache/download/url/bec801dcb9ee22f478e4a56b91820bd4/contents
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     286   (35, 30, 77)   float32   
  1  WCS-TAB       1 BinTableHDU     13   1R x 1C   [77E]   
Trimmed cube dimensions: 194 30 32

/tmp/ipykernel_2096/4105009505.py:39: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
  x.append(float(line))
/tmp/ipykernel_2096/4105009505.py:41: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
  x.append(float(line))
../../_images/42129ab91ef9d911cef7eb6ddaae58bda6fe5e274caeddd4bb28494897f50d20.png

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)
    

['pwl_c1', 'pwl_c2', 'sidust1_T', 'sidust1_amplitude', 'pahdust1_OnePlusZ', 'pahdust1_amplitude_76', 'pahdust1_amplitude_113', 'h55_amplitude', 'h55_xcen', 'h55_std', 'h61_amplitude', 'h61_xcen', 'h61_std', 'h69_amplitude', 'h69_xcen', 'h69_std', 'h80_amplitude', 'h80_xcen', 'h80_std', 'h96_amplitude', 'h96_xcen', 'h96_std', 'h122_amplitude', 'h122_xcen', 'h122_std', 'ArII_amplitude', 'ArII_xcen', 'ArII_std', 'SIV_amplitude', 'SIV_xcen', 'SIV_std', 'NeII_amplitude', 'NeII_xcen', 'NeII_std']
../../_images/87899aef14759b9110c2ed14f39aeadbf4320746729a4252f9b9451dc94fad1c.png ../../_images/c27e8dd7c39ba7a0ae525737ea8e6386e6638b75e4d13118f2919e8faa8721ba.png

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]

Time count
--- 819.6926140785217 seconds ---

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()
../../_images/7934afe29a21dba43f394c5cdd7727bd11eeb468d41a223feb7c81fc12e1e19b.png

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()
../../_images/49bc12c0c40767c884bb46f4a711ddf25b4d8d6baf075359c2ca1c53dab7f739.png ../../_images/cfefdbe42fda5149ed2a83483f8c43d8f7538fed95a0ddf801eff70d4fe18704.png

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

/tmp/ipykernel_2096/486002270.py:8: RuntimeWarning: invalid value encountered in divide
  arm_mask=np.where(((pah_neutral_flux>=0.155e-14) &((pah_ion_flux/pah_neutral_flux)>0.8) & (nuc_mask==0)),1,0)
/tmp/ipykernel_2096/486002270.py:16: RuntimeWarning: invalid value encountered in divide
  h2s_mask=np.where(((h2_flux>=1e-15)&((h2_flux/pah_neutral_flux)>0.7) & (nuc_mask==0)),1,0)
../../_images/fd62a55585f957c0203eeb499ce31270b9eb27b58f5de42cbb4aa4d3d0da3756.png

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()
../../_images/3e04721f7648da9b7616001a8b5f11edaddab3a712f26e551ef7e5b332dbc217.png

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)
../../_images/8baac928440f6c7c2b82107e009ca41a1a6726311d05df9f4632b728ecc0eae7.png ../../_images/d3d3d39b52da89329e32f7beec30ff78f265468f1ecdd85479b1f3a20cd9ee81.png

Model fit to nucleus, including fixed-temperature AGN silicate emission.

nuc_result
#result.fit_report(show_correl=False)

Fit Result

Model: ((((Model(powerlaw, prefix='pwl_') + Model(sidust, prefix='sidust1_')) + ((((((((((((((Model(Modelcero) + Model(drude, prefix='PAH52_')) + Model(drude, prefix='PAH62_')) + Model(drude, prefix='PAH74_')) + Model(drude, prefix='PAH76_')) + Model(drude, prefix='PAH78_')) + Model(drude, prefix='PAH83_')) + Model(drude, prefix='PAH86_')) + Model(drude, prefix='PAH112_')) + Model(drude, prefix='PAH113_')) + Model(drude, prefix='PAH119_')) + Model(drude, prefix='PAH126_')) + Model(drude, prefix='PAH134_')) + Model(drude, prefix='PAH140_')) + Model(drude, prefix='PAH141_'))) + ((((((Model(Modelcero) + Model(gaussian, prefix='h55_')) + Model(gaussian, prefix='h61_')) + Model(gaussian, prefix='h69_')) + Model(gaussian, prefix='h80_')) + Model(gaussian, prefix='h96_')) + Model(gaussian, prefix='h122_'))) + (((Model(Modelcero) + Model(gaussian, prefix='ArII_')) + Model(gaussian, prefix='SIV_')) + Model(gaussian, prefix='NeII_')))

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()
../../_images/a20ee6e1e5063435258a2cd8ccc00a5ab836fc3b37f503c4581708c4bdefcb5a.png
Space Telescope Logo

Notebook created by Patrick Ogle and Ivan Lopez.

Contents