Change the way the Maraston models are built in order to capture the rapid...

Change the way the Maraston models are built in order to capture the rapid evolutionary phases, similarly to the changes done for the BC03 models.

database_builder/__init__.py

@@ -228,7 +228,8 @@ def build_m2005(base):
     m2005_dir = os.path.join(os.path.dirname(__file__), 'maraston2005/')
 
     # Age grid (1 Myr to 13.7 Gyr with 1 Myr step)
-    age_grid = np.arange(1, 13701)
+    time_grid = np.arange(1, 13701)
+    fine_time_grid = np.linspace(0.1, 13700, 137000)
 
     # Transpose the table to have access to each value vector on the first
     # axis
@@ -258,50 +259,57 @@ def build_m2005(base):
         # populations.
         mass_table = mass_table[1:7, mass_table[0] == metallicity]
 
-        # 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)
+        # Regrid the SSP data to the evenly spaced time grid. In doing so we
+        # assume 10 bursts every 0.1 Myr over a period of 1 Myr in order to
+        # capture short evolutionary phases.
+        # The time grid starts after 0.1 Myr, so we assume the value is the same
+        # as the first actual time step.
+        mass_table = interpolate.interp1d(mass_table[0] * 1e3, mass_table[1:],
+                                          assume_sorted=True)(fine_time_grid)
+        mass_table = np.mean(mass_table.reshape(5, -1, 10), axis=-1)
 
-        # Remove the age column from the mass table
-        mass_table = np.delete(mass_table, 0, 0)
+        # Extract the age and convert from Gyr to Myr
+        ssp_time = np.unique(spec_table[0]) * 1e3
+        spec_table = spec_table[1:]
 
         # Remove the metallicity column from the spec table
-        spec_table = np.delete(spec_table, 1, 0)
-
-        # Convert the wavelength from Å to nm
-        spec_table[1] = spec_table[1] * 0.1
+        spec_table = spec_table[1:]
 
-        # For all ages, the lambda grid is the same.
-        lambda_grid = np.unique(spec_table[1])
+        # Extract the wavelength and convert from Å to nm
+        ssp_wave = spec_table[0][:1221] * 0.1
+        spec_table = spec_table[1:]
 
-        # Creation of the age vs lambda flux table
-        tmp_list = []
-        for wavelength in lambda_grid:
-            [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 *= 1000  # Gyr to Myr
-            flux_regrid = interpolate.interp1d(age_grid_orig,
-                                               flux_orig)(age_grid)
+        # Extra the fluxes and convert from erg/s/Å to W/nm
+        ssp_lumin = spec_table[0].reshape(ssp_time.size, ssp_wave.size).T
+        ssp_lumin *= 10 * 1e-7
 
-            tmp_list.append(flux_regrid)
-        flux_age = np.array(tmp_list)
+        # We have to do the interpolation-averaging in several blocks as it is
+        # a bit RAM intensive
+        ssp_lumin_interp = np.empty((ssp_wave.size, time_grid.size))
+        for i in range(0, ssp_wave.size, 100):
+            fill_value = (ssp_lumin[i:i+100, 0], ssp_lumin[i:i+100, -1])
+            ssp_interp = interpolate.interp1d(ssp_time, ssp_lumin[i:i+100, :],
+                                              fill_value=fill_value,
+                                              bounds_error=False,
+                                              assume_sorted=True)(fine_time_grid)
+            ssp_interp = ssp_interp.reshape(ssp_interp.shape[0], -1, 10)
+            ssp_lumin_interp[i:i+100, :] = np.mean(ssp_interp, axis=-1)
 
         # To avoid the creation of waves when interpolating, we refine the grid
         # beyond 10 μm following a log scale in wavelength. The interpolation
         # is also done in log space as the spectrum is power-law-like
-        lambda_grid_resamp = np.around(np.logspace(np.log10(10000),
+        ssp_wave_resamp = np.around(np.logspace(np.log10(10000),
                                                    np.log10(160000), 50))
-        argmin = np.argmin(10000.-lambda_grid > 0)-1
-        flux_age_resamp = 10.**interpolate.interp1d(
-                                    np.log10(lambda_grid[argmin:]),
-                                    np.log10(flux_age[argmin:, :]),
+        argmin = np.argmin(10000.-ssp_wave > 0)-1
+        ssp_lumin_resamp = 10.**interpolate.interp1d(
+                                    np.log10(ssp_wave[argmin:]),
+                                    np.log10(ssp_lumin_interp[argmin:, :]),
                                     assume_sorted=True,
-                                    axis=0)(np.log10(lambda_grid_resamp))
+                                    axis=0)(np.log10(ssp_wave_resamp))
 
-        lambda_grid = np.hstack([lambda_grid[:argmin+1], lambda_grid_resamp])
-        flux_age = np.vstack([flux_age[:argmin+1, :], flux_age_resamp])
+        ssp_wave = np.hstack([ssp_wave[:argmin+1], ssp_wave_resamp])
+        ssp_lumin = np.vstack([ssp_lumin_interp[:argmin+1, :],
+                               ssp_lumin_resamp])
 
         # Use Z value for metallicity, not log([Z/H])
         metallicity = {-1.35: 0.001,
@@ -309,8 +317,8 @@ def build_m2005(base):
                        0.0: 0.02,
                        0.35: 0.04}[metallicity]
 
-        base.add_m2005(M2005(imf, metallicity, age_grid, lambda_grid,
-                             mass_table, flux_age))
+        base.add_m2005(M2005(imf, metallicity, time_grid, ssp_wave,
+                             mass_table, ssp_lumin))
 
 
 def build_bc2003(base, res):