Commit 52af6e0f authored by Yannick Roehlly's avatar Yannick Roehlly

Use Myr for times

Using Gyr leads to some errors when comparing times (because they are no
really equals when they shoud). So we switch to the use of Myr as it is
the precision for our times.

Maybe later, we'll go further and use integers for these values.
parent 82092706
......@@ -46,7 +46,7 @@ def read_bc03_ssp(filename):
Returns
-------
time_grid: numpy 1D array of floats
Vector of the time grid of the SSP in Gyr.
Vector of the time grid of the SSP in Myr.
wavelength: numpy 1D array of floats
Vector of the wavelength grid of the SSP in nm.
spectra: numpy 2D array of floats
......@@ -119,9 +119,9 @@ def read_bc03_ssp(filename):
if counter == 0:
what_line = file_structure.next()
# The time grid is in year, we want Gyr.
# The time grid is in year, we want Myr.
time_grid = np.array(time_grid, dtype=float)
time_grid = time_grid * 1.e-9
time_grid = time_grid * 1.e-6
# The first "long" vector encountered is the wavelength grid. The value
# are in Ångström, we convert it to nano-meter.
......@@ -179,8 +179,8 @@ def build_base():
########################################################################
print("2- Importing Maraston 2005 SSP\n")
# Age grid (1My to 13.7Gy with 1My step)
age_grid = np.arange(1e-3, 13.701, 1e-3)
# Age grid (1 Myr to 13.7 Gyr with 1 Myr step)
age_grid = np.arange(1, 13701)
# Transpose the table to have access to each value vector on the first
# axis
......@@ -208,8 +208,10 @@ def build_base():
# we don't take the first column which contains metallicity
mass_table = mass_table[1:, mass_table[0] == metallicity]
# Interpolate the mass table over the new age grid
mass_table = interpolate.interp1d(mass_table[0], mass_table)(age_grid)
# Interpolate the mass table over the new age grid. We multiply per
# 1000 because the time in Maraston files is given in Gyr.
mass_table = interpolate.interp1d(mass_table[0] * 1000,
mass_table)(age_grid)
# Remove the age column from the mass table
mass_table = np.delete(mass_table, 0, 0)
......@@ -229,6 +231,7 @@ def build_base():
[age_grid_orig, lambda_grid_orig, flux_orig] = \
spec_table[:, spec_table[1, :] == wavelength]
flux_orig = flux_orig * 10 * 1.e-7 # From erg/s^-1/Å to W/nm
age_grid_orig = age_grid_orig * 1000 # Gyr to Myr
flux_regrid = interpolate.interp1d(age_grid_orig,
flux_orig)(age_grid)
......@@ -246,8 +249,8 @@ def build_base():
########################################################################
print("3- Importing Bruzual and Charlot 2003 SSP\n")
# Time grid (1My to 20Gy with 1My step)
time_grid = np.arange(1e-3, 20., 1e-3)
# Time grid (1 Myr to 20 Gyr with 1 Myr step)
time_grid = np.arange(1, 20000)
# Metallicities associated to each key
metallicity = {
......
......@@ -33,7 +33,7 @@ class SspBC03(object):
The metallicity. Possible values are 0.0001, 0.0004, 0.004, 0.008,
0.02 (solar metallicity) and 0.05.
time_grid : array of floats
The time [Gyr] grid used in the colors_table and the lumin_table.
The time [Myr] grid used in the colors_table and the lumin_table.
wavelength_grid : array of floats
The wavelength [nm] grid used in the lumin_table.
color_table: 2 axis array of floats
......@@ -78,9 +78,9 @@ class SspBC03(object):
Parameters
----------
sfh_time : array of floats
Time grid of the star formation history. It must be increasing and
not run beyond 20 Gyr. The SFH will be regrided to the SSP time
grid.
Time grid [Myr] of the star formation history. It must be
increasing and not run beyond 20 Gyr. The SFH will be regrided to
the SSP time.
sfh_sfr: array of floats
Star Formation Rates in Msun/yr at each time of the SFH time grid.
norm: boolean
......@@ -102,6 +102,7 @@ class SspBC03(object):
- "b4_vn": Amplitude of 4000 Å narrow break (Balogh et al. 1999)
- "b4_sdss" : Amplitude of 4000 Å break (Stoughton et al. 2002)
- "b_912": Amplitude of Lyman discontinuity
"""
# We work on a copy of SFH (as we change it)
sfh_time, sfh_sfr = np.copy((sfh_time, sfh_sfr))
......@@ -117,12 +118,12 @@ class SspBC03(object):
sfh_time, sfh_sfr,
left=0., right=0.)
# Step between two item in the time grid in Gyr
# Step between two item in the time grid in Myr
step = self.time_grid[1] - self.time_grid[0]
# If needed, we normalise the SFH to 1 solar mass produced.
if norm:
sfh_sfr = sfh_sfr / np.trapz(sfh_sfr * 1e9,
sfh_sfr = sfh_sfr / np.trapz(sfh_sfr * 1.e6,
self.time_grid[:idx + 1])
# As both the SFH and the SSP (limited to the age of the SFH) data now
......@@ -132,9 +133,9 @@ class SspBC03(object):
color_table = self.color_table[:, :idx + 1]
lumin_table = self.lumin_table[:, :idx + 1]
# The 1.e9 * step is because the SFH is in solar mass per year.
color_info = 1.e9 * step * np.dot(color_table, sfh_sfr[::-1])
luminosity = 1.e9 * step * np.dot(lumin_table, sfh_sfr[::-1])
# The 1.e6 * step is because the SFH is in solar mass per year.
color_info = 1.e6 * step * np.dot(color_table, sfh_sfr[::-1])
luminosity = 1.e6 * step * np.dot(lumin_table, sfh_sfr[::-1])
bc03_info = dict(zip(
["m_star", "m_gas", "n_ly", "b_4000", "b4_vn", "b4_sdss", "b_912"],
......
......@@ -48,7 +48,7 @@ class SspM2005(object):
* -0.33 (corresponding to 0.5 Zsun)
* -1.35 (corresponding to 1/50 Zsun)
time_grid : array of floats
The time [Gyr] grid used in the mass_table and the spec_table.
The time [Myr] grid used in the mass_table and the spec_table.
wavelength_grid : array of floats
The wavelength [nm] grid used in the spec_table.
mass_table : (6, n) array of floats
......@@ -92,10 +92,10 @@ class SspM2005(object):
Parameters
----------
sfh_time : array of floats
Time grid of the star formation history. It must be increasing and
not run beyond 13.7 Gyr. As the SFH will be regrided to the SSP
time grid, it is better to have a SFH age grid compatible, i.e.
with a precision limited to 1e-3 Gyr.
Time grid [Myr)] of the star formation history. It must be
increasing and not run beyond 13.7 Gyr. As the SFH will be
regrided to the SSP time grid, it is better to have a SFH age grid
compatible, i.e. with a precision limited to 1 Myr.
sfh_sfr: array of floats
Star Formation Rates in Msun/yr at each time of the SFH time grid.
norm: boolean
......@@ -120,34 +120,32 @@ class SspM2005(object):
# We work on a copy of SFH (as we change it)
sfh_time, sfh_sfr = np.copy((sfh_time, sfh_sfr))
# Index, in the SSP time grid, of the time nearest to the age of
# the SFH.
idx = np.abs(self.time_grid - np.max(sfh_time)).argmin()
# Step between two item in the time grid in Myr
step = self.time_grid[1] - self.time_grid[0]
# Number of step to go to the age of the SFH on the SSP age grid.
nb_steps = 1 + np.round((np.max(sfh_time) - self.time_grid[0]) / step)
# We regrid the SFH to the time grid of the SSP using a linear
# interpolation. If the SFH does no start at 0, the first SFR values
# will be set to 0.
sfh_sfr = np.interp(self.time_grid[:idx + 1],
sfh_sfr = np.interp(self.time_grid[:nb_steps],
sfh_time, sfh_sfr,
left=0., right=0.)
# Step between two item in the time grid in Gyr
step = self.time_grid[1] - self.time_grid[0]
# If needed, we normalise the SFH to 1 solar mass produced.
if norm:
sfh_sfr = sfh_sfr / np.trapz(sfh_sfr * 1e9,
self.time_grid[:idx + 1])
sfh_sfr = sfh_sfr / np.trapz(sfh_sfr * 1.e6,
self.time_grid[:nb_steps])
# As both the SFH and the SSP (limited to the age of the SFH) data now
# share the same time grid, the convolution is just a matter of
# reverting one and computing the sum of the one to one product; this
# is done using the dot product.
mass_table = self.mass_table[:, :idx + 1]
spec_table = self.spec_table[:, :idx + 1]
mass_table = self.mass_table[:, :nb_steps]
spec_table = self.spec_table[:, :nb_steps]
# The 1.e9 * step is because the SFH is in solar mass per year.
masses = 1.e9 * step * np.dot(mass_table, sfh_sfr[::-1])
spectra = 1.e9 * step * np.dot(spec_table, sfh_sfr[::-1])
# The 1.e6 * step is because the SFH is in solar mass per year.
masses = 1.e6 * step * np.dot(mass_table, sfh_sfr[::-1])
spectra = 1.e6 * step * np.dot(spec_table, sfh_sfr[::-1])
return masses, spectra
......@@ -14,8 +14,8 @@ from pcigale.data import Database
# Time lapse used to compute the average star formation rate. We use a
# constant to keep it easily changeable for advanced user while limiting the
# number of parameters. The value is in Gyr.
AV_LAPSE = 0.1
# number of parameters. The value is in Myr.
AV_LAPSE = 100
class Module(common.SEDCreationModule):
......@@ -43,17 +43,17 @@ class Module(common.SEDCreationModule):
),
"separation_age": (
"float",
"Age [Gyr] of the separation between the young and the old star "
"populations. The default value in 10^7 years (0.01 Gyr). Set "
"Age [Myr] of the separation between the young and the old star "
"populations. The default value in 10^7 years (10 Myr). Set "
"to 0 not to differentiate ages (only an old population).",
0.01
10
)
}
out_parameter_list = {
"sfr": "Instantaneous Star Formation Rate in solar mass per year, "
"at the age of the galaxy.",
"average_sfr": "Average SFR in the last 0.1 Gyr (default) of the "
"average_sfr": "Average SFR in the last 100 Myr (default) of the "
"galaxy history.",
"m_star": "Total mass in stars in Solar mass.",
"m_gas": "Mass returned to the ISM by evolved stars in Solar mass.",
......@@ -105,7 +105,7 @@ class Module(common.SEDCreationModule):
# SFR of the galaxy
sfr = sfh_sfr[len(sfh_sfr) - 1]
# Average SFR on the last AV_LAPSE Gyr of its history
# Average SFR on the last AV_LAPSE Myr of its history
average_sfr = np.mean(sfh_sfr[sfh_age <= AV_LAPSE])
# Base name for adding information to the SED.
......
......@@ -14,8 +14,8 @@ from pcigale.data import Database
# Time lapse used to compute the average star formation rate. We use a
# constant to keep it easily changeable for advanced user while limiting the
# number of parameters. The value is in Gyr.
AV_LAPSE = 0.1
# number of parameters. The value is in Myr.
AV_LAPSE = 100
class Module(common.SEDCreationModule):
......@@ -29,11 +29,11 @@ class Module(common.SEDCreationModule):
- imf, metallicity, galaxy_age
- sfr: star formation rate normalised to 1 solar mass formed at the
age of the galaxy.
- average_sfr: SFR averaged on the last 0.1 Gyr of the galaxy history
- average_sfr: SFR averaged on the last 100 Myr of the galaxy history
- mass_total, mass_alive, mass_white_dwarf,mass_neutrino,
mass_black_hole, mass_turn_off : stellar masses in solar mass.
- age: age of the oldest stars in the galaxy.
- old_young_separation_age: age (in Gyr) separating the young and the
- old_young_separation_age: age (in Myr) separating the young and the
old star populations (if 0, there is only one population)
- mass_total_old, mass_alive_old, mass_white_dwarf_old,
mass_neutrino_old, mass_black_hole_old, mass_turn_off_old: old
......@@ -65,17 +65,17 @@ class Module(common.SEDCreationModule):
),
'separation_age': (
'float',
"Age [Gyr] of the separation between the young and the old star "
"populations. The default value in 10^7 years (0.01 Gyr). Set "
"to 0 not to differentiate ages (only an old population).",
0.01
"Age [Myr] of the separation between the young and the old star "
"populations. The default value in 10^7 years (10 Myr). Set to "
"0 not to differentiate ages (only an old population).",
10
)
}
out_parameter_list = {
'sfr': 'Instantaneous Star Formation Rate in solar mass per year, '
'at the age of the galaxy.',
'average_sfr': 'Average SFR in the last 0.1 Gyr (default) of the '
'average_sfr': 'Average SFR in the last 100 Myr (default) of the '
'galaxy history.',
'mass_total': 'Total stellar mass of the galaxy in solar mass.',
'mass_alive': 'Mass of alive stars in solar mass.',
......@@ -83,7 +83,7 @@ class Module(common.SEDCreationModule):
'mass_neutrino': 'Mass of neutrino stars in solar mass.',
'mass_black_hole': 'Mass of black holes in solar mass.',
'mass_turn_off': 'Mass in the turn-off in solar mass.',
'old_young_separation_age': 'Age (in Gyr) separating the old and '
'old_young_separation_age': 'Age (in Myr) separating the old and '
'the young star populations (0 if there '
'is only one population).',
'mass_total_old': 'Total stellar mass of the old population in solar '
......@@ -152,7 +152,7 @@ class Module(common.SEDCreationModule):
# SFR of the galaxy
sfr = sfh_sfr[len(sfh_sfr) - 1]
# Average SFR on the last AV_LAPSE Gyr of its history
# Average SFR on the last AV_LAPSE Myr of its history
average_sfr = np.mean(sfh_sfr[sfh_age <= AV_LAPSE])
# Base name for adding information to the SED.
......
......@@ -10,9 +10,9 @@ Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
import numpy as np
from . import common
# Time lapse used in the age grid in Gyr. If should be consistent with the
# Time lapse used in the age grid in Myr. If should be consistent with the
# time lapse in the SSP modules.
AGE_LAPSE = 1e-3
AGE_LAPSE = 1
class Module(common.SEDCreationModule):
......@@ -27,12 +27,12 @@ class Module(common.SEDCreationModule):
parameter_list = {
"tau_main": (
"float",
"e-folding time of the main stellar population model in Gyr.",
"e-folding time of the main stellar population model in Myr.",
None
),
"tau_burst": (
"float",
"e-folding time of the late starburst population model in Gyr.",
"e-folding time of the late starburst population model in Myr.",
None
),
"f_burst": (
......@@ -42,25 +42,25 @@ class Module(common.SEDCreationModule):
),
"age": (
"float",
"Age of the oldest stars in the galaxy in Gyr. The precision "
"Age of the oldest stars in the galaxy in Myr. The precision "
"is 1 Myr.",
None
),
"burst_age": (
"float",
"Age of the late burst in Gyr. Precision is 1 Myr.",
"Age of the late burst in Myr. Precision is 1 Myr.",
None
)
}
out_parameter_list = {
"tau_main": "e-folding time of the main stellar population model "
"in Gyr.",
"in Myr.",
"tau_burst": "e-folding time of the late starburst population model "
"in Gyr.",
"in Myr.",
"f_burst": "Produced mass fraction of the late burst population.",
"age": "Age of the oldest stars in the galaxy in Gyr.",
"burst_age": "Age of the late burst in Gyr."
"age": "Age of the oldest stars in the galaxy in Myr.",
"burst_age": "Age of the late burst in Myr."
}
def _process(self, sed, parameters):
......@@ -101,7 +101,7 @@ class Module(common.SEDCreationModule):
(-time_grid[mask] + age - burst_age) / tau_burst)
# We normalise the SFH to have one solar mass produced.
sfr = sfr / np.trapz(sfr * 1e9, time_grid)
sfr = sfr / np.trapz(sfr * 1.e6, time_grid)
# Base name for adding information to the SED.
name = self.name or "sfh2exp"
......
......@@ -16,7 +16,7 @@ class Module(common.SEDCreationModule):
"""Module reading the SFH from a file
This module is used to read the Star Formation Histories from a FITS or
VO-Table file. The first column must contain the time values (in Gyr) and
VO-Table file. The first column must contain the time values (in Myr) and
each other column may contain the Star Formation Rates (in solar mass per
year) corresponding. Each SFR may be cut and normalised to 1 solar mass
produced at the desired age.
......@@ -27,7 +27,7 @@ class Module(common.SEDCreationModule):
"filename": (
"str",
"Name of the file containing the SFH. The first column must be "
"the time [Gyr] and the other column must contain the SFR "
"the time [Myr] and the other column must contain the SFR "
"[Msun/yr].",
None
),
......@@ -39,7 +39,7 @@ class Module(common.SEDCreationModule):
),
"age": (
"float",
"List of ages where each SFH will be looked at.",
"List of ages [Myr] where each SFH will be looked at.",
None
)
}
......@@ -71,7 +71,7 @@ class Module(common.SEDCreationModule):
time_grid = time_grid[time_grid <= age]
# The we normalise it to 1 solar mass produced.
sfr = sfr / np.trapz(sfr * 1e9, time_grid)
sfr = sfr / np.trapz(sfr * 1.e6, time_grid)
# Base name for adding information to the SED.
name = self.name or 'loadfile'
......
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