Update the nebular models from Akio Inoue. Now log U is sampled in steps of...

Update the nebular models from Akio Inoue. Now log U is sampled in steps of 0.1 dex down to -4. Also deviations from the 10000K case B assumption are taken into account.

CHANGELOG.md

@@ -3,12 +3,15 @@
 ## Unreleased
 ### Added
 - The stellar mass-weighted age is now provided. This is a much more usual measure of the age than the age of the oldest star. This is accessible with the `stellar.age_m_star` keyword in the `bc03` module with with the `stellar.age_mass` keyword in the `m2005` module. (Médéric Boquien)
+- The nebular models have been expanded from log U=-3 to log U=-4. (Médéric Boquien & Akio Inoue)
+- The nebular models are now sampled in steps of 0.1 dex in log U rather than 1.0 dex steps. (Médéric Boquien & Akio Inoue)
 
 ### Changed
 - We do not output the break strength from the `bc03` module anymore as these were not computed properly. (Médéric Boquien)
 
 ### Fixed
 - When the pcigale.ini file was missing, pcigale would crash and display a fairly cryptic backtrace. Now it explicitly states that the file could not be found. (Médéric Boquien)
+- The nebular emission now takes into account deviations from the 10000K case B assumption. In practice this yields fluxes about 10% fainter. (Médéric Boquien & Akio Inoue)
 
 ### Optimised
 - By default the MKL library created many threads for each for the parallel processes. Not only was this not necessary as a high-level parallelisation already exists, but it generated a strong oversubscription on the CPU and on the RAM. The slowdown was over a factor of ~2 in some cases. Now we mandate KML to use only 1 thread fo each process. (Médéric Boquien & Yannick Roehlly)

database_builder/__init__.py

@@ -612,58 +612,55 @@ def build_fritz2006(base):
 def build_nebular(base):
     models_lines = []
     models_cont = []
-    lines_dir = os.path.join(os.path.dirname(__file__), 'nebular/')
 
-    # Number of Lyman continuum photon to normalize the nebular continuum
-    # templates
-    nlyc_continuum = {'0.0001': 2.68786E+53, '0.0004': 2.00964E+53,
-                      '0.004': 1.79593E+53, '0.008': 1.58843E+53,
-                      '0.02': 1.24713E+53, '0.05': 8.46718E+52}
+    nebular_dir = os.path.join(os.path.dirname(__file__), 'nebular/')
+    print("Importing {}...".format(nebular_dir + 'lines.dat'))
+    lines = np.genfromtxt(nebular_dir + 'lines.dat')
 
-    for Z in ['0.0001', '0.0004', '0.004', '0.008', '0.02', '0.05']:
-        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))
+    wave_lines = np.genfromtxt(nebular_dir + 'line_wavelengths.dat')
+    print("Importing {}...".format(nebular_dir + 'continuum.dat'))
+    cont = np.genfromtxt(nebular_dir + 'continuum.dat')
 
-        # Convert wavelength from Å to nm
-        wave *= 0.1
+    # Convert wavelength from Å to nm
+    wave_lines *= 0.1
+    wave_cont = cont[:1600, 0] * 0.1
 
-        # Convert log(flux) into flux (arbitrary units)
-        ratio1 = 10**(ratio1-38.)
-        ratio2 = 10**(ratio2-38.)
-        ratio3 = 10**(ratio3-38.)
+    # Get the list of metallicities
+    metallicities = np.unique(lines[:, 1])
 
-        # Normalize all lines to Hβ
-        w = np.where(wave == 486.1)
-        ratio1 = ratio1/ratio1[w]
-        ratio2 = ratio2/ratio2[w]
-        ratio3 = ratio3/ratio3[w]
+    # Keep only the fluxes
+    lines = lines[:, 2:]
+    cont = cont[:, 1:]
 
-        models_lines.append(NebularLines(np.float(Z), -3., wave, ratio1))
-        models_lines.append(NebularLines(np.float(Z), -2., wave, ratio2))
-        models_lines.append(NebularLines(np.float(Z), -1., wave, ratio3))
+    # We select only models with ne=100. Other values could be included later
+    lines = lines[:, 1::3]
+    cont = cont[:, 1::3]
 
-        filename = "{}continuum_{}.dat".format(lines_dir, Z)
-        print("Importing {}...".format(filename))
-        wave, cont1, cont2, cont3 = np.genfromtxt(filename, unpack=True,
-                                                  usecols=(0, 3, 7, 11))
+    # Convert lines to W and to a linear scale
+    lines = 10**(lines-7)
 
-        # Convert wavelength from Å to nm
-        wave *= 0.1
+    # Convert continuum to W/nm
+    cont *= np.tile(1e-7 * cst.c * 1e9 / wave_cont**2,
+                    metallicities.size)[:, np.newaxis]
 
-        # Normalize flux from erg s¯¹ Hz¯¹ (Msun/yr)¯¹ to W nm¯¹ photon¯¹ s¯¹
-        conv = 1e-7 * cst.c * 1e9 / (wave * wave) / nlyc_continuum[Z]
-        cont1 *= conv
-        cont2 *= conv
-        cont3 *= conv
+    # Import lines
+    for idx, metallicity in enumerate(metallicities):
+        spectra = lines[idx::6, :]
+        for logU, spectrum in zip(np.around(np.arange(-4., -.9, .1), 1),
+                                  spectra.T):
+            models_lines.append(NebularLines(metallicity, logU, wave_lines,
+                                             spectrum))
 
-        models_cont.append(NebularContinuum(np.float(Z), -3., wave, cont1))
-        models_cont.append(NebularContinuum(np.float(Z), -2., wave, cont2))
-        models_cont.append(NebularContinuum(np.float(Z), -1., wave, cont3))
+    # Import continuum
+    for idx, metallicity in enumerate(metallicities):
+        spectra = cont[1600 * idx: 1600 * (idx+1), :]
+        for logU, spectrum in zip(np.around(np.arange(-4., -.9, .1), 1),
+                                  spectra.T):
+            models_cont.append(NebularContinuum(metallicity, logU, wave_cont,
+                                                spectrum))
 
-    base.add_nebular_continuum(models_cont)
     base.add_nebular_lines(models_lines)
+    base.add_nebular_continuum(models_cont)
 
 
 def build_schreiber2016(base):

database_builder/nebular/line_wavelengths.dat

+303.8   # He 2 303.8A
+584.3   # He 1 584.3A
+625.6   # He 1 625.6A
+972.5   # H  1 972.5A	Ly gamma (1-4)
+977.0   # C  3 977.0A
+1026    # H  1  1026A	Ly beta (1-3)
+1216    # H  1  1216A	Ly alpha (1-2)
+1335    # C  2  1335A	C II
+1397    # TOTL  1397A	Si IV 1394, 1397, 1398
+1486    # TOTL  1486A	N IV]
+1549    # TOTL  1549A	C IV 1548, 1551
+1640    # He 2  1640A
+1665    # TOTL  1665A	O III] 1661, 1666
+1750    # TOTL  1750A	N III]
+1888    # TOTL  1888A	Si III] 1883, 1892
+1909    # TOTL  1909A	C III]
+2326    # TOTL  2326A	C II]
+2400    # Fe 2  2400A
+2471    # O II  2471A
+2798    # TOTL  2798A	Mg II 2796, 2803
+2829    # He 1  2829A
+2836    # Fe 4  2836A
+2945    # He 1  2945A
+3096    # Fe 4  3096A
+3188    # He 1  3188A
+3203    # He 2  3203A
+3671    # Ca B  3671A	H 24 (Case B)
+3674    # Ca B  3674A	H 23 (Case B)
+3676    # Ca B  3676A	H 22 (Case B)
+3679    # Ca B  3679A	H 21 (Case B)
+3683    # Ca B  3683A	H 20 (Case B)
+3687    # Ca B  3687A	H 19 (Case B)
+3692    # Ca B  3692A	H 18 (Case B)
+3697    # Ca B  3697A	H 17 (Case B)
+3704    # Ca B  3704A	H 16 (Case B)
+3712    # Ca B  3712A	H 15 (Case B)
+3722    # Ca B  3722A	H 14 (Case B)
+3727    # TOTL  3727A	[OII] 3726, 3729
+3734    # Ca B  3734A	H 13 (Case B)
+3750    # Ca B  3750A	H 12 (Case B)
+3771    # Ca B  3771A	H 11 (Case B)
+3798    # H  1  3798A	H 10
+3820    # He 1  3820A
+3835    # H  1  3835A	H 9
+3869    # Ne 3  3869A
+3889    # He 1  3889A
+3889    # H  1  3889A	H 8
+3965    # He 1  3965A
+3968    # Ne 3  3968A
+3970    # H  1  3970A	H epsilon (2-7)
+4026    # He 1  4026A
+4070    # S II  4070A
+4074    # S  2  4074A
+4102    # H  1  4102A	H delta (2-6)
+4300    # Fe 2  4300A
+4340    # H  1  4340A	H gamma (2-5)
+4363    # TOTL  4363A
+4471    # He 1  4471A
+4659    # Fe 3  4659A
+4686    # He 2  4686A
+4861    # H  1  4861A	H beta (2-4)
+4922    # He 1  4922A
+4959    # O  3  4959A
+5007    # O  3  5007A
+5016    # He 1  5016A
+5876    # He 1  5876A
+6300    # O  1  6300A
+6312    # S  3  6312A
+6548    # N  2  6548A
+6563    # H  1  6563A	H alpha (2-3)
+6584    # N  2  6584A
+6678    # He 1  6678A
+6716    # S II  6716A
+6720    # S  2  6720A
+6731    # S II  6731A
+7065    # He 1  7065A
+7135    # Ar 3  7135A
+7325    # TOTL  7325A	[OII] 7323, 7332
+7751    # Ar 3  7751A
+8334    # Ca B  8334A	Pa 24 (Case B)
+8346    # Ca B  8346A	Pa 23 (Case B)
+8359    # Ca B  8359A	Pa 22 (Case B)
+8374    # Ca B  8374A	Pa 21 (Case B)
+8392    # Ca B  8392A	Pa 20 (Case B)
+8413    # Ca B  8413A	Pa 19 (Case B)
+8438    # Ca B  8438A	Pa 18 (Case B)
+8467    # Ca B  8467A	Pa 17 (Case B)
+8502    # Ca B  8502A	Pa 16 (Case B)
+8545    # Ca B  8545A	Pa 15 (Case B)
+8598    # Ca B  8598A	Pa 14 (Case B)
+8665    # Ca B  8665A	Pa 13 (Case B)
+8750    # Ca B  8750A	Pa 12
+8863    # Ca B  8863A	Pa 11
+9015    # H  1  9015A	Pa 10
+9069    # S  3  9069A
+9229    # H  1  9229A	Pa 9
+9532    # S  3  9532A
+9546    # H  1  9546A	Pa epsilon (3-8)
+10050   # H  1 1.005m	Pa delta (3-7)
+10830   # TOTL 1.083m	He I
+10940   # H  1 1.094m	Pa gamma (3-6)
+12820   # H  1 1.282m	Pa beta (3-5)
+18750   # H  1 1.875m	Pa alpha (3-4)
+21660   # H  1 2.166m	Br gamma (4-7)
+26250   # H  1 2.625m	Br beta (4-6)
+40510   # H  1 4.051m	Br alpha (4-5)
+69800   # Ar 2 6.980m
+74580   # H  1 7.458m	Pf alpha (5-6)
+90000   # Ar 3 9.000m
+105100  # S  4 10.51m
+128100  # Ne 2 12.81m
+155500  # Ne 3 15.55m
+186700  # S  3 18.67m
+334700  # S  3 33.47m
+348100  # Si 2 34.81m
+360100  # Ne 3 36.01m
+518000  # O  3 51.80m
+572100  # N  3 57.21m
+631700  # O  1 63.17m
+883300  # O  3 88.33m
+1217000 # N  2 121.7m
+1455000 # O  1 145.5m
+1576000 # C  2 157.6m
+2054000 # N  2 205.4m

pcigale/sed_modules/nebular.py

@@ -29,7 +29,10 @@ class NebularEmission(SedModule):
 
     parameter_list = OrderedDict([
         ('logU', (
-            'cigale_list(options=-3. & -2. & -1.)',
+            'cigale_list(options=-4.0 & -3.9 & -3.8 & -3.7 & -3.6 & -3.5 & '
+            '-3.4 & -3.3 & -3.2 & -3.1 & -3.0 & -2.9 & -2.8 & -2.7 & -2.6 & '
+            '-2.5 & -2.4 & -2.3 & -2.2 & -2.1 & -2.0 & -1.9 & -1.8 & -1.7 & '
+            '-1.6 & -1.5 & -1.4 & -1.3 & -1.2 & -1.1 & -1.0)',
             "Ionisation parameter",
             -2.
         )),
@@ -104,20 +107,15 @@ class NebularEmission(SedModule):
                 lines.wave = new_wave
                 lines.ratio = new_flux
 
-            # We compute the conversion coefficient to determine the fluxes using
-            # the formula of Inoue 2011: LHβ=Q(H)*γHβ(10000K)/αβ(10000K)
-            gamma_Hbeta = 1.23e-38  # Inoue 2011, W m³
-            alpha_B = 2.58e-19  # Ferland 1980, m³ s¯¹
-
             # To take into acount the escape fraction and the fraction of Lyman
             # continuum photons absorbed by dust we correct by a factor
             # k=(1-fesc-fdust)/(1+(α1/αβ)*(fesc+fdust))
+            alpha_B = 2.58e-19  # Ferland 1980, m³ s¯¹
             alpha_1 = 1.54e-19  # αA-αB, Ferland 1980, m³ s¯¹
             k = (1. - self.fesc - self.fdust) / (1. + alpha_1 / alpha_B * (
                 self.fesc + self.fdust))
 
-            self.conv_line = gamma_Hbeta / alpha_B * k
-            self.conv_cont = k
+            self.corr = k
         self.idx_Ly_break = None
         self.absorbed_old = None
         self.absorbed_young = None
@@ -164,14 +162,14 @@ class NebularEmission(SedModule):
             sed.add_info('nebular.logU', self.logU)
 
             sed.add_contribution('nebular.lines_old', lines.wave,
-                                 lines.ratio * NLy_old * self.conv_line)
+                                 lines.ratio * NLy_old * self.corr)
             sed.add_contribution('nebular.lines_young', lines.wave,
-                                 lines.ratio * NLy_young * self.conv_line)
+                                 lines.ratio * NLy_young * self.corr)
 
             sed.add_contribution('nebular.continuum_old', cont.wave,
-                                 cont.lumin * NLy_old * self.conv_cont)
+                                 cont.lumin * NLy_old * self.corr)
             sed.add_contribution('nebular.continuum_young', cont.wave,
-                                 cont.lumin * NLy_young * self.conv_cont)
+                                 cont.lumin * NLy_young * self.corr)
 
 # SedModule to be returned by get_module
 Module = NebularEmission