Skip to content
GitLab
Commit 7c2ce85a authored by Shuai Feng's avatar Shuai Feng
Browse files

fix PEP8 bugs

parent 1128112c
Pipeline #7131 failed with stage
in 0 seconds
Hide whitespace changes
Inline Side-by-side
csst_ifs_gehong/config.py
View file @ 7c2ce85a
......@@ -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
......
csst_ifs_gehong/cube3d.py
View file @ 7c2ce85a
......@@ -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)
csst_ifs_gehong/map2d.py
View file @ 7c2ce85a
......@@ -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
csst_ifs_gehong/spec1d.py
View file @ 7c2ce85a
......@@ -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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment