diff --git a/csst_ifs_gehong/config.py b/csst_ifs_gehong/config.py
index d5c06857a35b17dc15b7d923373161660d482862..9f8b85382a526c70636857e62b3c8d2894b2161c 100644
--- a/csst_ifs_gehong/config.py
+++ b/csst_ifs_gehong/config.py
@@ -22,13 +22,12 @@ class config():
dpix : float, optional
Pixel size in the spatial direction, by default 0.2arcsec
"""
- def __init__(self, wave_min=3500.0, wave_max=10000.0,
+ def __init__(self, wave_min=3500.0, wave_max=10000.0,
dlam=2.0, inst_fwhm=0.1,
nx=30, ny=30, dpix=0.2):
self.dlam = dlam
self.wave = np.arange(wave_min, wave_max, dlam)
self.wave_min = wave_min
-
self.inst_fwhm = inst_fwhm
self.nx = nx
self.ny = ny
diff --git a/csst_ifs_gehong/cube3d.py b/csst_ifs_gehong/cube3d.py
index 0bbf67b1959af6ab342559a893776024eb12b891..ae7811d291463725525c4aa18163c06fc5613c0a 100644
--- a/csst_ifs_gehong/cube3d.py
+++ b/csst_ifs_gehong/cube3d.py
@@ -9,86 +9,78 @@ class Cube3D():
Class of 3-dimentional spectral cube
"""
def __init__(self, config, stellar_map=None, gas_map=None):
-
self.config = config
-
self.nx = config.nx
self.ny = config.ny
self.dpix = config.dpix
self.fov_x = config.fov_x
self.fov_y = config.fov_y
-
self.wave = config.wave
self.nz = len(self.wave)
self.wave0 = np.min(self.wave)
self.inst_fwhm = config.inst_fwhm
-
self.flux = np.zeros((self.nx, self.ny, self.nz))
-
self.stellar_map = stellar_map
self.gas_map = gas_map
-
+
def make_cube(self, stellar_tem=None, hii_tem=None):
-
for i in range(self.nx):
for j in range(self.ny):
if self.stellar_map is not None:
- ss = StellarContinuum(self.config, stellar_tem, mag=self.stellar_map.mag[i, j],
- age=self.stellar_map.age[i, j], feh=self.stellar_map.feh[i, j],
- vel=self.stellar_map.vel[i, j], vdisp=self.stellar_map.vdisp[i, j],
+ ss = StellarContinuum(self.config, stellar_tem, mag=self.stellar_map.mag[i, j],
+ age=self.stellar_map.age[i, j], feh=self.stellar_map.feh[i, j],
+ vel=self.stellar_map.vel[i, j], vdisp=self.stellar_map.vdisp[i, j],
ebv=self.stellar_map.ebv[i, j])
if self.gas_map is not None:
- gg = HII_Region(self.config, hii_tem, halpha=self.gas_map.halpha[i, j],
- logz=self.gas_map.zh[i, j], vel=self.gas_map.vel[i, j],
+ gg = HII_Region(self.config, hii_tem, halpha=self.gas_map.halpha[i, j],
+ logz=self.gas_map.zh[i, j], vel=self.gas_map.vel[i, j],
vdisp=self.gas_map.vdisp[i, j], ebv=self.gas_map.ebv[i, j])
self.flux[i, j, :] = ss.flux + gg.flux
else:
self.flux[i, j, :] = ss.flux
-
+
def wcs_info(self):
# wcs = fits.Header()
- wcs_dict = {'CTYPE1': 'WAVE ',
- 'CUNIT1': 'Angstrom',
- 'CDELT1': self.config.dlam,
- 'CRPIX1': 1,
- 'CRVAL1': np.min(self.wave),
- 'CTYPE2': 'RA---TAN',
- 'CUNIT2': 'deg',
- 'CDELT2': self.dpix / 3600.,
- 'CRPIX2': np.round(self.ny / 2.),
- 'CRVAL2': 0.5,
- 'CTYPE3': 'DEC--TAN',
- 'CUNIT3': 'deg',
- 'CDELT3': self.dpix / 3600.,
- 'CRPIX3': np.round(self.nx / 2.),
- 'CRVAL3': 1,
+ wcs_dict = {'CTYPE1': 'WAVE ',
+ 'CUNIT1': 'Angstrom',
+ 'CDELT1': self.config.dlam,
+ 'CRPIX1': 1,
+ 'CRVAL1': np.min(self.wave),
+ 'CTYPE2': 'RA---TAN',
+ 'CUNIT2': 'deg',
+ 'CDELT2': self.dpix / 3600.,
+ 'CRPIX2': np.round(self.ny / 2.),
+ 'CRVAL2': 0.5,
+ 'CTYPE3': 'DEC--TAN',
+ 'CUNIT3': 'deg',
+ 'CDELT3': self.dpix / 3600.,
+ 'CRPIX3': np.round(self.nx / 2.),
+ 'CRVAL3': 1,
'BUNIT': '10**(-17)*erg/s/cm**2/Angstrom'}
input_wcs = astropy.wcs.WCS(wcs_dict)
self.wcs_header = input_wcs.to_header()
-
+
def insert_spec(self, spec, dx=0, dy=0):
x0 = np.int(np.round(self.config.nx / 2.))
y0 = np.int(np.round(self.config.ny / 2.))
self.flux[x0 + dx, y0 + dy, :] = self.flux[x0 + dx, y0 + dy, :] + spec.flux
-
- def savefits(self, filename, path='./'):
-
+
+ def savefits(self, filename, path='./'):
hdr = fits.Header()
-
hdr['FILETYPE'] = 'SCICUBE'
hdr['CODE'] = 'CSST-IFS-GEHONG'
hdr['VERSION'] = '0.0.1'
-
+
hdr['OBJECT'] = 'NGC1234'
hdr['RA'] = 0.0
hdr['DEC'] = 0.0
-
+
hdu0 = fits.PrimaryHDU(header=hdr)
-
+
self.wcs_info()
hdr = self.wcs_header
hdu1 = fits.ImageHDU(self.flux, header=hdr)
-
+
# Output
hdulist = fits.HDUList([hdu0, hdu1])
- hdulist.writeto(path + filename, overwrite=True)
\ No newline at end of file
+ hdulist.writeto(path + filename, overwrite=True)
diff --git a/csst_ifs_gehong/map2d.py b/csst_ifs_gehong/map2d.py
index 5ffd512a9207b50653be39876a0841c69d3ba7cd..db5e1a2d49b6c169177d1ff3f01ef63800d87eed 100644
--- a/csst_ifs_gehong/map2d.py
+++ b/csst_ifs_gehong/map2d.py
@@ -49,23 +49,23 @@ def Sersic2D(x, y, mag=12.0, r_eff=1.0, n=2.0, ellip=0.5,
x_min = -(x - x_0) * sin_theta + (y - y_0) * cos_theta
z = (abs(x_maj / a) ** 2 + abs(x_min / b) ** 2) ** (1 / 2)
profile = np.exp(-bn * (z ** (1 / n) - 1))
-
+
# Normalization
integral = a * b * 2 * np.pi * n * np.exp(bn) / (bn ** (2 * n)) * sp.gamma(2 * n)
prof_norm = profile / integral * pixelscale
-
+
# Calibration
total_flux = 10. ** ((22.5 - mag) * 0.4)
sb_mag = 22.5 - 2.5 * np.log10(prof_norm * total_flux / pixelscale)
-
+
return sb_mag
-def VelMap2D(x, y, vmax=200.0, rt=1.0, ellip=0.5,
+def VelMap2D(x, y, vmax=200.0, rt=1.0, ellip=0.5,
theta=0.0, x_0=0.0, y_0=0.0):
"""
- VelMap2D: Caculate the velocity at given spatial position base on the rotating disk model.
- The rotation curve is adpot as the tanh model.
+ VelMap2D: Caculate the velocity at given spatial position base on the rotating disk model.
+ The rotation curve is adpot as the tanh model.
Parameters
----------
@@ -102,15 +102,15 @@ def VelMap2D(x, y, vmax=200.0, rt=1.0, ellip=0.5,
x_min = -(x - x_0) * sin_theta + (y - y_0) * cos_theta
z = (abs(x_maj / a) ** 2 + abs(x_min / b) ** 2) ** (1 / 2)
profile = vmax * np.tanh(z) * ((x_maj / a) / z)
-
+
return profile
-def GradMap2D(x, y, a0=10.0, r_eff=1.0, gred=-1.0, ellip=0.5,
+def GradMap2D(x, y, a0=10.0, r_eff=1.0, gred=-1.0, ellip=0.5,
theta=0.0, x_0=0.0, y_0=0.0):
"""
- GradMap2D: Caculate the intensity at given spatial position base on the disk model.
- The radial profile is adpot as the gradient model.
+ GradMap2D: Caculate the intensity at given spatial position base on the disk model.
+ The radial profile is adpot as the gradient model.
Parameters
----------
@@ -149,7 +149,7 @@ def GradMap2D(x, y, a0=10.0, r_eff=1.0, gred=-1.0, ellip=0.5,
x_min = -(x - x_0) * sin_theta + (y - y_0) * cos_theta
z = (abs(x_maj / a) ** 2 + abs(x_min / b) ** 2) ** (1 / 2)
profile = a0 + z * gred
-
+
return profile
@@ -207,7 +207,7 @@ class Map2d(object):
xsh_rot = xsh * np.cos(pa_radians) + ysh * np.sin(pa_radians)
ysh_rot = -xsh * np.sin(pa_radians) + ysh * np.cos(pa_radians)
return ysh_rot, xsh_rot
-
+
def sersic_map(self, mag=12.0, r_eff=2.0, n=2.5, ellip=0.5, theta=-50.0):
"""
Generate 2D map of Sersic model
@@ -246,11 +246,11 @@ class Map2d(object):
self.n = n
self.ellip = ellip
self.theta = theta
- self.map = Sersic2D(self.x, self.y, mag=self.mag,
- r_eff=self.reff, n=self.n,
- ellip=self.ellip, theta=self.theta,
+ self.map = Sersic2D(self.x, self.y, mag=self.mag,
+ r_eff=self.reff, n=self.n,
+ ellip=self.ellip, theta=self.theta,
pixelscale=self.xsamp * self.ysamp)
-
+
def tanh_map(self, vmax=200.0, rt=2.0, ellip=0.5, theta=-50.0):
"""
Generate 2D velocity map of rotating disk according to tanh rotation curve
@@ -284,9 +284,9 @@ class Map2d(object):
self.rt = rt / self.xsamp
self.ellip = ellip
self.theta = theta
- self.map = VelMap2D(self.x, self.y, vmax=self.vmax, rt=self.rt,
+ self.map = VelMap2D(self.x, self.y, vmax=self.vmax, rt=self.rt,
ellip=self.ellip, theta=self.theta)
-
+
def gred_map(self, a0=10, r_eff=1, gred=-1, ellip=0.5, theta=0):
"""
Generate 2D maps according to the radial gradient form
@@ -319,9 +319,9 @@ class Map2d(object):
self.gred = gred
self.ellip = ellip
self.theta = theta
- self.map = GradMap2D(self.x, self.y, a0=self.a0, r_eff=self.reff,
+ self.map = GradMap2D(self.x, self.y, a0=self.a0, r_eff=self.reff,
gred=self.gred, ellip=self.ellip, theta=self.theta)
-
+
def load_map(self, image):
"""
Generate 2D map according to input image
@@ -339,10 +339,10 @@ class Map2d(object):
class StellarPopulationMap():
"""
- Class of 2D maps for the parameters of stellar population, such as
- surface brightness, median age and metallicity of stellar population,
+ Class of 2D maps for the parameters of stellar population, such as
+ surface brightness, median age and metallicity of stellar population,
velocity and velocity dispersion maps, and dust extinction.
-
+
Parameters
----------
config : class
@@ -360,7 +360,7 @@ class StellarPopulationMap():
ebv : class, optional
Class of the map of dust extinction, by default None
"""
- def __init__(self, config, sbright=None, logage=None,
+ def __init__(self, config, sbright=None, logage=None,
feh=None, vel=None, vdisp=None, ebv=None):
self.nx = config.nx
@@ -368,29 +368,29 @@ class StellarPopulationMap():
self.dpix = config.dpix
self.fov_x = config.fov_x
self.fov_y = config.fov_y
-
+
if (sbright is None):
print('Input SurfaceBrightness Map is empty!')
else:
self.sbright = sbright.map
self.mag = self.sbright - 2.5 * np.log10(self.dpix * self.dpix)
-
+
if (logage is None):
print('Input Age Map is empty!')
else:
self.logage = logage.map
self.age = 10 ** self.logage / 1e9
-
+
if (feh is None):
print('Input Metallicity Map is empty!')
else:
self.feh = feh.map
-
+
if (vel is None):
print('Input Velocity Map is empty!')
else:
self.vel = vel.map
-
+
if (vdisp is None):
print('Input VelocityDispersion Map is empty!')
else:
@@ -399,7 +399,7 @@ class StellarPopulationMap():
if len(self.vdisp[ind_overrange]) > 0:
print("Notice: Spaxel with <10km/s in the input vdisp map will be automatically adjusted to 10km/s.")
self.vdisp[ind_overrange] = 10
-
+
if (ebv is None):
print('Input EBV Map is empty!')
else:
@@ -408,12 +408,12 @@ class StellarPopulationMap():
if len(self.ebv[ind_overrange]) > 0:
print("Notice: Spaxel with < 0 mag in the input ebv map will be automatically adjusted to 0 mag.")
self.ebv[ind_overrange] = 0
-
+
class IonizedGasMap():
"""
- Class of 2D maps for the parameters of ionized gas, such as
- Halpha flux map, gas-phase metallicity map,
+ Class of 2D maps for the parameters of ionized gas, such as
+ Halpha flux map, gas-phase metallicity map,
velocity and velocity dispersion maps, and dust extinction.
Parameters
@@ -432,13 +432,13 @@ class IonizedGasMap():
Class of the map of dust extinction, by default None
"""
def __init__(self, config, halpha=None, zh=None, vel=None, vdisp=None, ebv=None):
-
+
self.nx = config.nx
self.ny = config.ny
self.dpix = config.dpix
self.fov_x = config.fov_x
self.fov_y = config.fov_y
-
+
if (halpha is None):
print('Input Halpha Map is empty!')
else:
@@ -470,4 +470,4 @@ class IonizedGasMap():
ind_overrange = (self.ebv < 0)
if len(self.ebv[ind_overrange]) > 0:
print("Notice: Spaxel with < 0 mag in the input ebv map will be automatically adjusted to 0 mag.")
- self.ebv[ind_overrange] = 0
\ No newline at end of file
+ self.ebv[ind_overrange] = 0
diff --git a/csst_ifs_gehong/spec1d.py b/csst_ifs_gehong/spec1d.py
index c15cd429447715cf6c5463ef08dcb6a979e3c794..9b226c02caa4f9e3abb674563153884b23829416 100644
--- a/csst_ifs_gehong/spec1d.py
+++ b/csst_ifs_gehong/spec1d.py
@@ -47,8 +47,8 @@ def readcol(filename, **kwargs):
def log_rebin(lamRange, spec, oversample=False, velscale=None, flux=False):
"""
- Logarithmically rebin a spectrum, while rigorously conserving the flux.
- This function is taken from ppxf.
+ Logarithmically rebin a spectrum, while rigorously conserving the flux.
+ This function is taken from ppxf.
Parameters
----------
@@ -175,7 +175,7 @@ def calibrate(wave, flux, mag, filtername='SLOAN_SDSS.r'):
else:
filter_file = data_path + '/data/filter/' + filtername + '.filter'
wave0, response0 = readcol(filter_file)
-
+
# Setting the response
func = interp1d(wave0, response0)
response = np.copy(wave)
@@ -183,17 +183,17 @@ def calibrate(wave, flux, mag, filtername='SLOAN_SDSS.r'):
response[ind_extra] = 0
ind_inside = (wave < max(wave0)) & (wave > min(wave0))
response[ind_inside] = func(wave[ind_inside])
-
+
# Flux map of datacube for given filter band
preflux = np.sum(flux * response * np.mean(np.diff(wave))) / np.sum(response * np.mean(np.diff(wave)))
-
+
# Real flux from magnitude for given filter
realflux = (mag * u.STmag).to(u.erg / u.s / u.cm**2 / u.AA).value
# Normalization
flux_ratio = realflux / preflux
flux_calibrate = flux * flux_ratio * 1e17 # Units: 10^-17 erg/s/A/cm^2
-
+
return flux_calibrate
# ----------------
@@ -218,11 +218,10 @@ def Calzetti_Law(wave, Rv=4.05):
"""
wave_number = 1. / (wave * 1e-4)
reddening_curve = np.zeros(len(wave))
-
+
idx = (wave >= 1200) & (wave < 6300)
- reddening_curve[idx] = 2.659 * (-2.156 + 1.509 * wave_number[idx] - 0.198 *
- (wave_number[idx] ** 2)) + 0.011 * (wave_number[idx] ** 3) + Rv
-
+ reddening_curve[idx] = 2.659 * (-2.156 + 1.509 * wave_number[idx] - 0.198 * (wave_number[idx] ** 2)) + 0.011 * (wave_number[idx] ** 3) + Rv
+
idx = (wave >= 6300) & (wave <= 22000)
reddening_curve[idx] = 2.659 * (-1.857 + 1.040 * wave_number[idx]) + Rv
return reddening_curve
@@ -277,7 +276,7 @@ def SingleEmissinoLine(wave, line_wave, FWHM_inst):
float array
Spectra of single emission line
"""
- sigma = FWHM_inst / 2.355
+ sigma = FWHM_inst / 2.355
flux = norm.pdf(wave, line_wave, sigma)
return flux
@@ -295,20 +294,20 @@ class EmissionLineTemplate():
dlam : float, optional
Wavelength width per pixel, by default 0.1A
model : str, optional
- Emission line model, including 'hii' for HII region and 'nlr' for narrow line region of AGN,
+ Emission line model, including 'hii' for HII region and 'nlr' for narrow line region of AGN,
by default 'hii'
"""
-
+
def __init__(self, config, lam_range=[500, 15000], dlam=0.1, model='hii'):
-
+
self.lam_range = lam_range
self.wave = np.arange(lam_range[0], lam_range[1], 0.1)
self.FWHM_inst = config.inst_fwhm
self.model = model
-
+
# HII region model of fsps-cloudy
if model == 'hii':
-
+
# Loading emission line flux table
flux_table_file = data_path + '/data/fsps.nebular.fits'
line_table = fits.open(flux_table_file)
@@ -317,18 +316,18 @@ class EmissionLineTemplate():
line_list = line_table[1].data
line_wave = line_list['Wave']
line_names = line_list['Name']
-
+
w = (line_wave > lam_range[0]) & (line_wave < lam_range[1])
self.line_names = line_names[w]
self.line_wave = line_wave[w]
-
+
# Make parameter grid
grid = line_table[2].data
self.logz_grid = grid['logZ']
-
- # Narrow line region model of
+
+ # Narrow line region model of
if model == 'nlr':
-
+
# Loading emission line flux table
flux_table_file = data_path + '/data/AGN.NLR.fits'
line_table = fits.open(flux_table_file)
@@ -337,19 +336,19 @@ class EmissionLineTemplate():
line_list = line_table[1].data
line_wave = line_list['Wave']
line_names = line_list['Name']
-
+
w = (line_wave > lam_range[0]) & (line_wave < lam_range[1])
self.line_names = line_names[w]
self.line_wave = line_wave[w]
-
+
# Make parameter grid
grid = line_table[2].data
self.logz_grid = grid['logZ']
-
+
# Flux ratio
flux_ratio = line_table[3].data
self.flux_ratio = flux_ratio
-
+
# Make emission line
nline = len(line_wave)
for i in range(nline):
@@ -389,9 +388,9 @@ class HII_Region():
ValueError
The value of logZ should be between -2 and 0.5.
"""
-
+
def __init__(self, config, temp, halpha=100.0, logz=0.0, vel=100.0, vdisp=120.0, ebv=0.1):
-
+
# Check Input Parameters
# if (halpha < 0):
# print("Notice: Your input Halpha flux (halpha) is < 0 erg/s/A/cm^2. ")
@@ -403,25 +402,25 @@ class HII_Region():
# if (ebv < 0):
# ebv = 0
# print('Notice: Your input dust extinction (ebv) is < 0 mag, which will be automatically adjusted to 0 mag. ')
-
- # Make emission line spectra through adding emission lines
+
+ # Make emission line spectra through adding emission lines
emlines = temp.emission_lines * flux_ratio
flux_combine = np.sum(emlines, axis=1)
flux_calibrate = flux_combine * halpha # Units: erg/s/A/cm^2
-
+
# Dust attenuation
if np.isscalar(ebv):
flux_dust = reddening(temp.wave, flux_calibrate, ebv=ebv)
-
+
# Broadening caused by Velocity Dispersion
velscale = 10
lam_range = [np.min(temp.wave), np.max(temp.wave)]
flux_logwave, logLam = log_rebin(lam_range, flux_dust, velscale=velscale)[:2]
-
+
sigma_gas = vdisp / velscale # in pixel
sigma_LSF = temp.FWHM_inst / (np.exp(logLam)) * 3e5 / velscale # in pixel
-
- if sigma_gas > 0:
+
+ if sigma_gas > 0:
sigma_dif = np.zeros(len(flux_logwave))
idx = (sigma_gas > sigma_LSF)
sigma_dif[idx] = np.sqrt(sigma_gas ** 2. - sigma_LSF[idx] ** 2.)
@@ -431,12 +430,12 @@ class HII_Region():
else:
print('Notice: Your input velocity dispersion (vdisp) is <= 0 km/s. The nagtive value will be automatically adjusted to 0 km/s. ')
flux_broad = flux_logwave
-
+
# Redshift
redshift = vel / 3e5
wave_r = np.exp(logLam) * (1 + redshift)
flux_red = np.interp(config.wave, wave_r, flux_broad)
-
+
self.wave = config.wave
self.flux = flux_red
@@ -472,45 +471,45 @@ class AGN_NLR():
ValueError
The value of logZ should be between -2 and 0.5.
"""
-
+
def __init__(self, config, temp, halpha=100.0, logz=0.0, vel=100.0, vdisp=120.0, ebv=0.1):
-
+
if (logz >= -2.3) & (logz <= 0.54):
indz = np.argmin(np.abs(logz - temp.logz_grid))
flux_ratio = temp.flux_ratio[indz, :]
else:
raise ValueError('The value of logZ is not in the range of [-2.3, 0.54]')
-
- # Make emission line spectra through adding emission lines
+
+ # Make emission line spectra through adding emission lines
emlines = temp.emission_lines * (flux_ratio / flux_ratio[6])
flux_combine = np.sum(emlines, axis=1)
flux_calibrate = flux_combine * halpha # Units: 1e-17 erg/s/A/cm^2
-
+
# Dust attenuation
if np.isscalar(ebv):
flux_dust = reddening(temp.wave, flux_calibrate, ebv=ebv)
-
+
# Broadening caused by Velocity Dispersion
velscale = 10
lam_range = [np.min(temp.wave), np.max(temp.wave)]
flux_logwave, logLam = log_rebin(lam_range, flux_dust, velscale=velscale)[:2]
-
+
sigma_gas = vdisp / velscale # in pixel
sigma_LSF = temp.FWHM_inst / (np.exp(logLam)) * 3e5 / velscale # in pixel
-
- if sigma_gas > 0:
+
+ if sigma_gas > 0:
sigma_dif = np.zeros(len(flux_logwave))
idx = (sigma_gas > sigma_LSF)
sigma_dif[idx] = np.sqrt(sigma_gas ** 2. - sigma_LSF[idx] ** 2.)
idx = (sigma_gas <= sigma_LSF)
sigma_dif[idx] = 0.1
flux_broad = gaussian_filter1d(flux_logwave, sigma_dif)
-
+
# Redshift
redshift = vel / 3e5
wave_r = np.exp(logLam) * (1 + redshift)
flux_red = np.interp(config.wave, wave_r, flux_broad)
-
+
self.wave = config.wave
self.flux = flux_red
@@ -536,42 +535,42 @@ class AGN_BLR():
Wavelength range, by default [500, 15000]
"""
- def __init__(self, config, hbeta_flux=100.0, hbeta_fwhm=2000.0, ebv=0.1,
+ def __init__(self, config, hbeta_flux=100.0, hbeta_fwhm=2000.0, ebv=0.1,
vel=0., lam_range=[500, 15000]):
wave_rest = np.arange(lam_range[0], lam_range[1], 0.1)
-
+
line_names = ['Hepsilon', 'Hdelta', 'Hgamma', 'Hbeta', 'Halpha']
line_waves = [3970.079, 4101.742, 4340.471, 4861.333, 6562.819]
line_ratio = [0.101, 0.208, 0.405, 1.000, 2.579] # From Ilic et al. (2006)
-
+
# Make emission lines
for i in range(len(line_names)):
if i == 0:
- emission_line = SingleEmissinoLine(wave_rest, line_waves[i],
+ emission_line = SingleEmissinoLine(wave_rest, line_waves[i],
hbeta_fwhm / 3e5 * line_waves[i])
emission_lines = emission_line
else:
- emission_line = SingleEmissinoLine(wave_rest, line_waves[i],
+ emission_line = SingleEmissinoLine(wave_rest, line_waves[i],
hbeta_fwhm / 3e5 * line_waves[i])
emission_lines = np.vstack((emission_lines, emission_line))
emlines = emission_lines.T * line_ratio
flux_combine = np.sum(emlines, axis=1)
-
+
# Flux callibration
flux_calibrate = flux_combine * hbeta_flux # Units: 1e-17 erg/s/A/cm^2
-
+
# Dust attenuation
if np.isscalar(ebv):
flux_dust = reddening(wave_rest, flux_calibrate, ebv=ebv)
else:
flux_dust = flux_calibrate
-
+
# Redshift
redshift = vel / 3e5
wave_r = wave_rest * (1 + redshift)
flux_red = np.interp(config.wave, wave_r, flux_dust)
-
+
self.wave = config.wave
self.flux = flux_red
@@ -593,39 +592,39 @@ class AGN_FeII():
ebv : float, optional
Dust extinction, by default 0.1
"""
-
+
def __init__(self, config, hbeta_broad=100.0, r4570=0.4, ebv=0.1, vel=100.0):
filename = data_path + '/data/FeII.AGN.fits'
-
+
# Loading FeII template
hdulist = fits.open(filename)
data = hdulist[1].data
wave_rest = data['WAVE']
flux_model = data['FLUX']
-
+
# Determine the flux of FeII
Fe4570_temp = 100
Fe4570_model = hbeta_broad * r4570
Ratio_Fe4570 = Fe4570_model / Fe4570_temp
-
+
# Flux calibration
flux_calibrate = flux_model * Ratio_Fe4570
-
+
# Dust attenuation
if np.isscalar(ebv):
flux_dust = reddening(wave_rest, flux_calibrate, ebv=ebv)
else:
flux_dust = flux_calibrate
-
+
# Redshift
redshift = vel / 3e5
wave_r = wave_rest * (1 + redshift)
flux_red = np.interp(config.wave, wave_r, flux_dust)
-
+
self.wave = config.wave
self.flux = flux_red
-
+
class AGN_Powerlaw():
@@ -645,29 +644,29 @@ class AGN_Powerlaw():
Ebv : float, optional
Dust extinction, by default 0.1
"""
-
+
def __init__(self, config, m5100=1000.0, alpha=-1.5, vel=100.0, ebv=0.1):
wave_rest = np.linspace(1000, 20000, 10000)
flux = wave_rest ** alpha
-
+
# Flux calibration
flux_calibrate = calibrate(wave_rest, flux, m5100, filtername='5100')
-
+
# Dust attenuation
if np.isscalar(ebv):
flux_dust = reddening(wave_rest, flux_calibrate, ebv=ebv)
else:
flux_dust = flux_calibrate
-
+
# Redshift
redshift = vel / 3e5
wave_r = wave_rest * (1 + redshift)
flux_red = np.interp(config.wave, wave_r, flux_dust)
-
+
self.wave = config.wave
self.flux = flux_red
-
+
class AGN():
@@ -681,7 +680,7 @@ class AGN():
nlr_template : class
Class of emission line template
bhmass : float, optional
- Black hole mass used for calculating the luminosity of power law spectrum at 5100A,
+ Black hole mass used for calculating the luminosity of power law spectrum at 5100A,
by default 1e6 solar mass
edd_ratio : float, optional
Eddinton ratio used for calculating the luminosity of power law spectrum at 5100A, by default 0.05
@@ -702,10 +701,10 @@ class AGN():
dist : float, optional
Luminosity distance of AGN, by default 20.0Mpc
"""
- def __init__(self, config, nlr_template, bhmass=1e6, edd_ratio=0.05,
- halpha_broad=100.0, halpha_narrow=100.0, vdisp_broad=5000.0, vdisp_narrow=500.0,
+ def __init__(self, config, nlr_template, bhmass=1e6, edd_ratio=0.05,
+ halpha_broad=100.0, halpha_narrow=100.0, vdisp_broad=5000.0, vdisp_narrow=500.0,
vel=1000.0, logz=0.0, ebv=0.1, dist=20.0):
-
+
# Check Input Parameters
if (bhmass < 0):
raise Exception("Input black hole mass (bhmass) should be >= 0 solarmass.")
@@ -724,17 +723,17 @@ class AGN():
raise Exception('Input metallicity of narrow lines (logz) should be < -2.3')
if (dist <= 0):
raise Exception('Input luminosity distance (dist) should be > 0 Mpc')
-
+
self.wave = config.wave
self.flux = self.wave * 0
-
- if (vdisp_narrow > 0) & (halpha_narrow > 0):
+
+ if (vdisp_narrow > 0) & (halpha_narrow > 0):
NLR = AGN_NLR(config, nlr_template, halpha=halpha_narrow, logz=logz,
vel=vel, vdisp=vdisp_narrow, ebv=ebv)
self.flux = self.flux + NLR.flux
if (halpha_broad > 0) & (vdisp_broad > 0):
- BLR = AGN_BLR(config, hbeta_flux=halpha_broad / 2.579,
+ BLR = AGN_BLR(config, hbeta_flux=halpha_broad / 2.579,
hbeta_fwhm=vdisp_broad / 2.355, ebv=ebv, vel=vel)
self.flux = self.flux + BLR.flux
@@ -743,7 +742,7 @@ class AGN():
PL = AGN_Powerlaw(config, m5100=m5100, ebv=ebv, vel=vel)
Fe = AGN_FeII(config, hbeta_broad=halpha_broad / 2.579, ebv=ebv, vel=vel)
self.flux = self.flux + PL.flux + Fe.flux
-
+
def BHmass_to_M5100(bhmass, edd_ratio=0.05, dist=21.0):
"""
@@ -763,7 +762,7 @@ def BHmass_to_M5100(bhmass, edd_ratio=0.05, dist=21.0):
float
Magnitude at 5100A
"""
-
+
# Calculate bolometric luminosity
Ledd = 3e4 * bhmass
Lbol = Ledd * edd_ratio
@@ -817,7 +816,7 @@ def age_metal(filename):
class StellarContinuumTemplate(object):
"""
- Class of single stellar population template.
+ Class of single stellar population template.
Parameters
----------
@@ -826,17 +825,17 @@ class StellarContinuumTemplate(object):
velscale : array
velocity scale in km/s per pixels, by default 50.0km/s
pathname : string, optional
- path with wildcards returning the list files to use,
+ path with wildcards returning the list files to use,
by default data_path+'/data/EMILES/Ech*_baseFe.fits'
normalize : bool, optional
Set to True to normalize each template to mean=1, by default False
"""
def __init__(self, config, velscale=50,
- pathname=data_path + '/data/EMILES/Ech*_baseFe.fits',
+ pathname=data_path + '/data/EMILES/Ech*_baseFe.fits',
normalize=False):
FWHM_inst = config.inst_fwhm
-
+
files = glob.glob(pathname)
assert len(files) > 0, "Files not found %s" % pathname
@@ -857,7 +856,7 @@ class StellarContinuumTemplate(object):
lam_range_temp = h2['CRVAL1'] + np.array([0, h2['CDELT1'] * (h2['NAXIS1'] - 1)])
sspNew, log_lam_temp = log_rebin(lam_range_temp, ssp, velscale=velscale)[:2]
# wave=((np.arange(hdr['NAXIS1'])+1.0)-hdr['CRPIX1'])*hdr['CDELT1']+hdr['CRVAL1']
-
+
templates = np.empty((sspNew.size, n_ages, n_metal))
age_grid = np.empty((n_ages, n_metal))
metal_grid = np.empty((n_ages, n_metal))
@@ -870,7 +869,7 @@ class StellarContinuumTemplate(object):
# Quadratic sigma difference in pixels Vazdekis --> galaxy
# The formula below is rigorously valid if the shapes of the
# instrumental spectral profiles are well approximated by Gaussians.
-
+
# FWHM of Emiles templates
Emile_wave = np.exp(log_lam_temp)
Emile_FWHM = np.zeros(h2['NAXIS1'])
@@ -879,9 +878,9 @@ class StellarContinuumTemplate(object):
Emile_FWHM[np.where((Emile_wave >= 3540) & (Emile_wave < 8950))] = 2.5
Lwave = Emile_wave[np.where(Emile_wave >= 8950)]
Emile_FWHM[np.where(Emile_wave >= 8950)] = 60 * 2.35 / 3.e5 * Lwave # sigma=60km/s at lambda > 8950
-
+
LSF = Emile_FWHM
-
+
FWHM_eff = Emile_FWHM.copy() # combined FWHM from stellar library and instrument(input)
if np.isscalar(FWHM_inst):
FWHM_eff[Emile_FWHM < FWHM_inst] = FWHM_inst
@@ -918,8 +917,8 @@ class StellarContinuumTemplate(object):
self.n_metal = n_metal
self.LSF = log_rebin(lam_range_temp, LSF, velscale=velscale)[0]
self.velscale = velscale
-
- def fmass_ssp(self):
+
+ def fmass_ssp(self):
isedpath = data_path + '/data/EMILES/model/'
massfile = isedpath + 'out_mass_CH_PADOVA00'
@@ -928,7 +927,7 @@ class StellarContinuumTemplate(object):
# fage = self.age_grid[:,0]
fMs = np.zeros((n_ages, n_metal))
-
+
Metal, Age, Ms = readcol(massfile, usecols=(2, 3, 6))
for i in range(n_metal):
for j in range(self.n_ages):
@@ -936,12 +935,12 @@ class StellarContinuumTemplate(object):
fMs[j, i] = Ms[locmin]
return fMs
-
+
class StellarContinuum():
"""
The class of stellar continuum
-
+
Parameters
----------
config : class
@@ -961,9 +960,9 @@ class StellarContinuum():
ebv : float, optional
Dust extinction, by default 0.1
"""
- def __init__(self, config, template, mag=15.0, age=1.0, feh=0.0,
+ def __init__(self, config, template, mag=15.0, age=1.0, feh=0.0,
vel=100.0, vdisp=100.0, ebv=0.1):
-
+
# Check Input Parameters
if (mag > 26) or (mag < 8):
print("Notice: Your input magnitude (mag) is > 26 mag or < 8 mag.")
@@ -986,19 +985,19 @@ class StellarContinuum():
minloc = np.argmin(abs(feh - metals))
tpls = SSP_temp[:, :, minloc]
# fmass = template.fmass_ssp()[:, minloc]
-
+
# Select age bins
Ages = template.age_grid[:, 0]
minloc = np.argmin(abs(age - Ages))
Stellar = tpls[:, minloc]
-
+
wave = np.exp(template.log_lam_temp)
-
+
# Broadening caused by Velocity Dispersion
sigma_gal = vdisp / template.velscale # in pixel
sigma_LSF = template.LSF / template.velscale # in pixel
-
- if sigma_gal > 0:
+
+ if sigma_gal > 0:
sigma_dif = np.zeros(len(Stellar))
idx = (sigma_gal > sigma_LSF)
sigma_dif[idx] = np.sqrt(sigma_gal ** 2. - sigma_LSF[idx] ** 2.)
@@ -1013,21 +1012,21 @@ class StellarContinuum():
# Dust Reddening
if np.isscalar(ebv):
flux0 = reddening(wave, flux0, ebv=ebv)
-
+
# Redshift
redshift = vel / 3e5
wave_r = wave * (1 + redshift)
-
+
flux = np.interp(config.wave, wave_r, flux0)
-
+
# Calibration
if np.isscalar(mag):
flux = calibrate(config.wave, flux, mag, filtername='SLOAN_SDSS.r')
-
+
# Convert to input wavelength
self.wave = config.wave
self.flux = flux
-
+
#####################
# Single Star Spectra
#####################
@@ -1048,15 +1047,15 @@ class SingleStarTemplate():
FWHM_inst = config.inst_fwhm
filename = data_path + '/data/Starlib.XSL.fits'
-
+
hdulist = fits.open(filename)
lam = hdulist[1].data['Wave']
flux = hdulist[2].data
par = hdulist[3].data
-
+
lam_range_temp = np.array([3500, 12000])
TemNew, log_lam_temp = log_rebin(lam_range_temp, flux[1, :], velscale=velscale)[:2]
-
+
# FWHM of XLS templates
Temp_wave = np.exp(log_lam_temp)
Temp_FWHM = np.zeros(len(log_lam_temp))
@@ -1064,9 +1063,9 @@ class SingleStarTemplate():
Temp_FWHM[(Temp_wave >= 5330) & (Temp_wave < 9440)] = 11 * 2.35 / 3e5 * Temp_wave[(Temp_wave >= 5330) & (Temp_wave < 9440)]
# sigma = 13km/s at 5330 < lambda < 9440
Temp_FWHM[(Temp_wave >= 9440)] = 16 * 2.35 / 3e5 * Temp_wave[(Temp_wave >= 9440)] # sigma=16km/s at lambda > 9440
-
+
LSF = Temp_FWHM
-
+
FWHM_eff = Temp_FWHM.copy() # combined FWHM from stellar library and instrument(input)
if np.isscalar(FWHM_inst):
FWHM_eff[Temp_FWHM < FWHM_inst] = FWHM_inst
@@ -1094,7 +1093,7 @@ class SingleStarTemplate():
self.logg_grid = par['logg']
self.LSF = Temp_FWHM
self.velscale = velscale
-
+
class SingleStar():
@@ -1135,34 +1134,34 @@ class SingleStar():
feh = 0.81
if (feh < -2.69):
raise Exception('Your input metallicity (feh) is beyond the range of stellar template [-2.69, 0.81]')
-
+
StarTemp = template.templates
-
+
# Select metal bins
idx_FeH = (np.abs(template.feh_grid - feh) < 0.5)
tpls = StarTemp[:, idx_FeH]
-
+
# Select Teff bins
Teff_FeH = template.teff_grid[idx_FeH]
minloc = np.argmin(abs(teff - Teff_FeH))
starspec = tpls[:, minloc]
-
+
wave = np.exp(template.log_lam_temp)
-
+
# Dust Reddening
if np.isscalar(ebv):
starspec = reddening(wave, starspec, ebv=ebv)
-
+
# Redshift
redshift = vel / 3e5
wave_r = wave * (1 + redshift)
-
+
flux = np.interp(config.wave, wave_r, starspec)
-
+
# Calibration
if np.isscalar(mag):
flux = calibrate(config.wave, flux, mag, filtername='SLOAN_SDSS.r')
-
+
# Convert to input wavelength
self.wave = config.wave
self.flux = flux