Commit 9daf7130 authored by Médéric Boquien's avatar Médéric Boquien

Initial implementation of the nebular lines.

parent fce2f0e4
......@@ -20,7 +20,7 @@ import numpy as np
from scipy import interpolate
import scipy.constants as cst
from pcigale.data import (Database, Filter, SspM2005, SspBC03, AgnFritz2006,
Dale2014, DL2007)
Dale2014, DL2007, Lines)
def read_bc03_ssp(filename):
......@@ -483,6 +483,39 @@ def build_fritz2006(base):
base.add_fritz2006_agn(agn)
def build_lines(base):
lines_dir = os.path.join(os.path.dirname(__file__), 'lines/')
for Z in ['0.0000001', '0.00001', '0.0004', '0.008', '0.02']:
filename = "{}lines_{}.dat".format(lines_dir, Z)
print("Importing {}...".format(filename))
wave, ratio1, ratio2, ratio3 = np.genfromtxt(filename, unpack=True,
usecols=(0, 3, 7, 11))
# Convert wavelength from Å to nm
wave *= 0.1
# Convert log(flux) into flux (arbitrary units)
ratio1 = 10**(ratio1-38.)
ratio2 = 10**(ratio2-38.)
ratio3 = 10**(ratio3-38.)
# Normalize all lines to Hβ
w = np.where(wave==486.1)
ratio1 = ratio1/ratio1[w]
ratio2 = ratio2/ratio2[w]
ratio3 = ratio3/ratio3[w]
lines = Lines(np.float(Z), -3., wave, ratio1)
base.add_lines(lines)
lines = Lines(np.float(Z), -2., wave, ratio2)
base.add_lines(lines)
lines = Lines(np.float(Z), -1., wave, ratio3)
base.add_lines(lines)
def build_base():
base = Database(writable=True)
base.upgrade_base()
......@@ -522,6 +555,11 @@ def build_base():
build_dale2014(base)
print("\nDONE\n")
print('#' * 78)
print("8- Build nebular lines\n")
build_lines(base)
print("\nDONE\n")
print('#' * 78)
base.session.close_all()
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -128,7 +128,7 @@ class CreationModule(object):
# We want to postfix the various keys of the SED with the same
# postfix as the module name, if any.
if '.' in name:
self.postfix = "." + name.split(".", 1)
self.postfix = "." + name.split(".", 1)[1]
else:
self.postfix = ""
......
# -*- coding: utf-8 -*-
# Copyright (C) 2014 University of Cambridge
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Médéric Boquien <mboquien@ast.cam.ac.uk>
from collections import OrderedDict
import numpy as np
from pcigale.data import Database
import scipy.constants as cst
from . import CreationModule
class Module(CreationModule):
"""
Module computing the nebular emission lines from the ultraviolet
to the near-infrared
Given the number of Lyman continuum photons, we compute the Hβ line
luminosity. We then compute the other lines using the
metallicity-dependent templates that provide the ratio between individual
lines and Hβ.
Information added to the SED: NAME_escape_fraction.
"""
parameter_list = OrderedDict([
('logU', (
'float',
"Radiation field intensity",
-2.
)),
('escape_fraction', (
'float',
"Fraction of Lyman continuum photons escaping the galaxy",
0.
)),
('f_dust', (
'float',
"Fraction of Lyman continuum photons absorbed by dust",
0.
)),
('lines_width', (
'float',
"Line width in km s¯¹",
300.
))
])
out_parameter_list = OrderedDict([
('logU', "Radiation field intensity"),
('escape_fraction', "Fraction of Lyman continuum photons escaping "
"the galaxy"),
('f_dust', "Fraction of Lyman continuum photons absorbed by dust"),
('lines_width', "Width of the lines")
])
def _init_code(self):
"""Get the nebular emission lines out of the database and resample
them to see the line profile.
"""
database = Database()
self.lines_template = {m:database.get_lines(m, self.parameters['logU'])
for m in database.get_lines_metallicities()}
database.session.close_all()
lines_width = self.parameters['lines_width'] * 1e3
for lines in self.lines_template.values():
new_wave = np.array([])
for line_wave in lines.wave:
width = line_wave * lines_width / cst.c
new_wave = np.concatenate((new_wave, np.linspace(line_wave -
3. * width,
line_wave
+ 3. *
width,
19)))
new_wave = np.sort(new_wave)
new_flux = np.zeros_like(new_wave)
for line_flux, line_wave in zip(lines.ratio, lines.wave):
width = line_wave * lines_width / cst.c
new_flux += (line_flux * np.exp(- 4. * np.log(2.) *
(new_wave - line_wave) ** 2. / (width * width)) /
(width * np.sqrt(np.pi / np.log(2.)) / 2.))
lines.wave = new_wave
lines.ratio = new_flux
# We compute the conversion coefficient to determine the fluxes using
# the formula of Krüger+95: h×ν(Hβ)*α(Hβ,10000K)/αB(10000K)
# Osterbrock 1989 gives:
# α(Hβ,10000K)=3.03×10¯¹⁴
# αB(10000K)2.59×10¯¹³
self.conv = (13.598 * cst.electron_volt * (1. / 4. - 1. / 16.) *
3.03e-14 / 2.59e-13 *
(1. - self.parameters['escape_fraction'] -
self.parameters['f_dust'])
)
print(self.conv)
def process(self, sed):
"""Add the nebular emission lines
Parameters
----------
sed : pcigale.sed.SED object
parameters : dictionary containing the parameters
"""
# Base name for adding information to the SED.
name = self.name or 'lines'
escape_fraction = self.parameters['escape_fraction']
f_dust = self.parameters['f_dust']
NLy_old = sed.info['ssp_n_ly_old']
NLy_young = sed.info['ssp_n_ly_young']
lines = self.lines_template[sed.info['ssp_metallicity']]
sed.add_module(name, self.parameters)
sed.add_info(name + '_escape_fraction', escape_fraction)
sed.add_info(name + '_f_dust', f_dust)
sed.add_contribution('lines_old', lines.wave, NLy_old * lines.ratio *
self.conv)
sed.add_contribution('lines_young', lines.wave, NLy_young * lines.ratio *
self.conv)
......@@ -27,6 +27,7 @@ from .ir_templates_dh2002 import IrTemplatesDH2002
from .ir_agn_templates_dale2014 import Dale2014
from .ir_models_dl2007 import DL2007
from .agn_fritz2006 import AgnFritz2006
from .lines import Lines
DATABASE_FILE = pkg_resources.resource_filename(__name__, 'data.db')
......@@ -200,6 +201,23 @@ class _Fritz2006AGN(BASE):
self.luminosity = agn.luminosity
class _Lines(BASE):
"""Storage for line templates
"""
__tablename__ = 'lines'
metallicity = Column(Float, primary_key=True)
logU = Column(Float, primary_key=True)
wave = Column(PickleType)
ratio = Column(PickleType)
def __init__(self, lines):
self.metallicity = lines.metallicity
self.logU = lines.logU
self.wave = lines.wave
self.ratio = lines.ratio
class Database(object):
"""Object giving access to pcigale database."""
......@@ -391,6 +409,20 @@ class Database(object):
else:
raise StandardError('The database is not writable.')
def add_lines(self, lines):
"""
Add ultraviolet and optical line templates to the database.
"""
if self.is_writable:
self.session.add(_Lines(lines))
try:
self.session.commit()
except exc.IntegrityError:
self.session.rollback()
raise StandardError('The line is already in the base')
else:
raise StandardError('The database is not writable')
def get_filter(self, name):
"""
Get a specific filter from the collection
......@@ -617,6 +649,39 @@ class Database(object):
"The DL2007 model for qpah <{0}>, umin <{1}>, and umax <{2}> "
"is not in the database.".format(qpah, umin, umax))
def get_lines(self, metallicity, logU):
"""
Get the line ratios corresponding to the given set of parameters.
Parameters
----------
metallicity: float
Gas phase metallicity
logU: float
Radiation field intensity
"""
result = (self.session.query(_Lines).
filter(_Lines.metallicity == metallicity).
filter(_Lines.logU == logU).
first())
if result:
return Lines(result.metallicity, result.logU, result.wave,
result.ratio)
else:
return None
def get_lines_metallicities(self):
"""Get the list of metallicities in the lines database
Returns
-------
metallicities: list
list of the line templates metallicities
"""
result = self.session.query(_Lines.metallicity).all()
return [_[0] for _ in result]
def get_filter_list(self):
"""Get the list of the filters in the database.
......
## -*- coding: utf-8 -*-
# Copyright (C) 2013 Centre de données Astrophysiques de Marseille
# Copyright (C) 2014 Institute of Astronomy, University of Cambridge
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Médéric Boquien
class Lines(object):
"""Nebular lines templates.
This class holds the data associated with the line templates
"""
def __init__(self, metallicity, logU, wave, ratio):
"""Create a new nebular lines template
Parameters
----------
metallicity: float
Gas phase metallicity
logU: float
Radiation field intensity
wave: array
Vector of the λ grid used in the templates [nm]
ratio: array
Line intensities relative to Hβ
"""
self.metallicity = metallicity
self.logU = logU
self.wave = wave
self.ratio = ratio
Markdown is supported
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