Add the possibility to have an exponential burst with the sfhdelayed module in...

Add the possibility to have an exponential burst with the sfhdelayed module in order to model the latest episode of star formation.

CHANGELOG.md

@@ -6,6 +6,8 @@
 - It is now possible to fit any physical property indicated by the code (e.g. equivalent width, dust luminosity, etc.). For this the physical property needs to be given in the input file and the properties to be fitted must be given in the properties filed in pcigale.ini. (Héctor Salas & Médéric Boquien)
 
 ### Changed
+- The `sfhdelayed` module has been extended to optionally include an exponential burst to model the latest episode of star formation. (Médéric Boquien & Barbara Lo Faro)
+
 ### Fixed
 - The histogram bin width was not computed optimally when some models were invalid. (David Corre & Médéric Boquien)
 - Missing import in the `m2005` module. (Médéric Boquien, reported by Dominika Wylezalek)

pcigale/sed_modules/sfhdelayed.py

+# -*- coding: utf-8 -*-
 # Copyright (C) 2013 Centre de données Astrophysiques de Marseille
-# Copyright (C) 2014 Laboratoire d'Astrophysique de Marseille
+# Copyright (C) 2014, 2016 Laboratoire d'Astrophysique de Marseille
 # Copyright (C) 2014 University of Cambridge
+# Copyright (C) 2018 Universidad de Antofagasta
 # Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
-# Author: Yannick Roehlly, Véronique Buat & Médéric Boquien
 
 """
-Delayed tau model for star formation history
-============================================
+Delayed tau model for star formation history with an optional exponential burst
+===============================================================================
 
 This module implements a star formation history (SFH) described as a delayed
-rise of the SFR up to a maximum, followed by an exponential decrease.
+rise of the SFR up to a maximum, followed by an exponential decrease. Optionally
+a decreasing exponential burst can be added to model a recent episode of star
+formation.
 
 """
 
@@ -21,10 +24,13 @@ from . import SedModule
 
 
 class SFHDelayed(SedModule):
-    """Delayed tau model for Star Formation History
+    """Delayed tau model for Star Formation History with an optionally
+    exponential burst.
 
     This module sets the SED star formation history (SFH) proportional to time,
-    with a declining exponential parametrised with a time-scale (tau_main).
+    with a declining exponential parametrised with a time-scale τ. Optionally
+    an exp(-t_/τ_burst) component can be added to model the latest episode of
+    star formation.
 
     """
 
@@ -34,37 +40,66 @@ class SFHDelayed(SedModule):
             "e-folding time of the main stellar population model in Myr.",
             2000.
         )),
-        ("age", (
+        ("age_main", (
             "cigale_list(dtype=int, minvalue=0.)",
-            "Age of the oldest stars in the galaxy in Myr. The precision "
-            "is 1 Myr.",
+            "Age of the main stellar population in the galaxy in Myr. The "
+            "precision is 1 Myr.",
             5000
         )),
+        ("tau_burst", (
+            "cigale_list()",
+            "e-folding time of the late starburst population model in Myr.",
+            50.
+        )),
+        ("age_burst", (
+            "cigale_list(dtype=int, minvalue=1.)",
+            "Age of the late burst in Myr. The precision is 1 Myr.",
+            20
+        )),
+        ("f_burst", (
+            "cigale_list(minvalue=0., maxvalue=0.9999)",
+            "Mass fraction of the late burst population.",
+            0.
+        )),
         ("sfr_A", (
             "cigale_list(minvalue=0.)",
-            "Multiplicative factor controlling the amplitude of SFR.",
+            "Value of SFR at t = 0 in M_sun/yr.",
             1.
         )),
         ("normalise", (
             "boolean()",
             "Normalise the SFH to produce one solar mass.",
             True
-        ))
+        )),
     ])
 
     def _init_code(self):
         self.tau_main = float(self.parameters["tau_main"])
-        age = int(self.parameters["age"])
+        age_main = int(self.parameters["age_main"])
+        self.tau_burst = float(self.parameters["tau_burst"])
+        self.age_burst = int(self.parameters["age_burst"])
+        self.f_burst = float(self.parameters["f_burst"])
         sfr_A = float(self.parameters["sfr_A"])
         normalise = bool(self.parameters["normalise"])
 
-        # Time grid and SFR
-        time_grid = np.arange(age)
-        self.sfr = (time_grid * np.exp(-time_grid / self.tau_main) /
-                    self.tau_main**2)
+        # Time grid for each component
+        t = np.arange(age_main)
+        t_burst = np.arange(self.age_burst)
 
-        # Compute the galaxy mass and normalise the SFH to 1 solar mass
-        # produced if asked to.
+        # SFR for each component
+        self.sfr = t * np.exp(-t / self.tau_main) / self.tau_main**2
+        sfr_burst = np.exp(-t_burst / self.tau_burst)
+
+        # Height of the late burst to have the desired produced mass fraction
+        sfr_burst *= (self.f_burst / (1.-self.f_burst) * np.sum(self.sfr) /
+                      np.sum(sfr_burst))
+
+        # We add the age burst exponential for ages superior to age_main -
+        # age_burst
+        self.sfr[-(t_burst[-1]+1):] += sfr_burst
+
+        # Compute the integral of the SFH and normalise it to 1 solar mass
+        # if asked to.
         self.sfr_integrated = np.sum(self.sfr) * 1e6
         if normalise:
             self.sfr /= self.sfr_integrated
@@ -80,13 +115,16 @@ class SFHDelayed(SedModule):
         sed : pcigale.sed.SED object
 
         """
+
         sed.add_module(self.name, self.parameters)
 
         # Add the sfh and the output parameters to the SED.
         sed.sfh = self.sfr
         sed.add_info("sfh.integrated", self.sfr_integrated, True)
         sed.add_info("sfh.tau_main", self.tau_main)
-
+        sed.add_info("sfh.age_burst", self.age_burst)
+        sed.add_info("sfh.tau_burst", self.tau_burst)
+        sed.add_info("sfh.f_burst", self.f_burst)
 
 # SedModule to be returned by get_module
 Module = SFHDelayed