testing IR luminosity based on SED models in .txt

DaisyInteLumTest.py

'''

Just integrating a model SED in observed frame from mbb & RP's output

Last Modified: 05 May 2016

'''

from astropy.cosmology import WMAP9
import numpy as np

RPsedOBS = np.loadtxt('forDaisy.txt')     # mJy
DaisyFIRobs = np.loadtxt('rangeFIR.txt')   # mJy

z = 0.6537
# Some constants
c = 299792458e6  # in microns
h = 6.6260693e-34  # J/s
k = 1.3806505e-23  # J/K
um_to_GHz = 299792458e-3
Jy_to_SI = 1e-26
Mpc_to_m = 3.086e+22

# micron, what we used to output the SED parameters
waverest_start = 1
waverest_end = 1000+1   # for df
# waverest_start = 42
# waverest_end = 123+1   # for df
wave_rest_um = range(waverest_start, waverest_end)
freq_rest_Hz = c / np.array(wave_rest_um)
df_Hz = abs(np.ediff1d(freq_rest_Hz))


def differential_lum(dSnu_SI):
    # differential

    DL = WMAP9.luminosity_distance(z).value*Mpc_to_m

    dLrest = 4*np.pi*DL**2/(1+z)*dSnu_SI
    return dLrest

lnu = differential_lum(RPsedOBS*1e-3*Jy_to_SI)
Lrest_SI = abs(np.trapz(lnu, df_Hz))
Lrest_SI = np.sum(lnu*df_Hz)
Lrest_Lsun = Lrest_SI / 3.826e26
print("L IR from {}-{} um: {} 1e12 Lsun").format(waverest_start, waverest_end-1, Lrest_Lsun/1e12)