Make model computations much faster. On a test run we spent 20% of the time...

Make model computations much faster. On a test run we spent 20% of the time just computing the luminosity distance. Now the information is just provided by the redshifting module. If the information is not present we just assume a distance of 10 pc.

pcigale/creation_modules/redshifting.py

@@ -21,6 +21,7 @@ is changed, this module may need to be adapted.
 from collections import OrderedDict
 from ..creation_modules import CreationModule
 from pcigale.sed.cosmology import cosmology
+from scipy.constants import parsec
 
 
 class Redshifting(CreationModule):
@@ -44,6 +45,12 @@ class Redshifting(CreationModule):
         """
         self.redshift = float(self.parameters["redshift"])
         self.universe_age = cosmology.age(self.redshift).value * 1000.
+        if self.redshift == 0.:
+            self.luminosity_distance = 10. * parsec
+        else:
+            self.luminosity_distance = (
+                cosmology.luminosity_distance(self.redshift).value * 1e6 *
+                parsec)
 
     def process(self, sed):
         """Redshift the SED
@@ -71,6 +78,7 @@ class Redshifting(CreationModule):
         sed.luminosities /= 1. + self.redshift
 
         sed.add_info("redshift", self.redshift)
+        sed.add_info("universe.luminosity_distance", self.luminosity_distance)
         sed.add_info("universe.age", self.universe_age)
         sed.add_module(self.name, self.parameters)
 

pcigale/sed/__init__.py

@@ -40,7 +40,7 @@ import numpy as np
 from collections import OrderedDict
 from . import utils
 from .io.vo import save_sed_to_vo
-from scipy.constants import c
+from scipy.constants import c, parsec
 from scipy.interpolate import interp1d
 
 # Time lapse used to compute the average star formation rate. We use a
@@ -142,7 +142,8 @@ class SED(object):
         """
 
         # Fλ flux density in W/m²/nm
-        f_lambda = utils.luminosity_to_flux(self.luminosity, self.info['redshift'])
+        f_lambda = utils.luminosity_to_flux(self.luminosity,
+                                            self.info['universe.luminosity_distance'])
 
         # Fν flux density in mJy
         f_nu = utils.lambda_flambda_to_fnu(self.wavelength_grid, f_lambda)
@@ -382,17 +383,17 @@ class SED(object):
             transmission_r = np.interp(wavelength_r, transmission[0],
                                        transmission[1])
 
-            if 'redshift' in self.info.keys():
-                redshift = self.info['redshift']
+            if 'universe.luminosity_distance' in self.info.keys():
+                dist = self.info['universe.luminosity_distance']
             else:
-                redshift = 0.
+                dist = 10. * parsec
 
             # TODO: Can we avoid to normalise as the filter transmission is
             # already normalised?
             f_lambda = utils.luminosity_to_flux(
                 (np.trapz(transmission_r * l_lambda_r, wavelength_r) /
                  np.trapz(transmission_r, wavelength_r)),
-                redshift
+                dist
             )
 
             # Add the Fλ fluxes from the spectral lines.

pcigale/sed/utils.py

@@ -4,8 +4,7 @@
 # Authors: Yannick Roehlly, Médéric Boquien
 
 import numpy as np
-from scipy.constants import c, pi, parsec
-from .cosmology import cosmology
+from scipy.constants import c, pi
 
 
 def lambda_to_nu(wavelength):
@@ -68,7 +67,7 @@ def best_grid(wavelengths1, wavelengths2):
     return new_grid
 
 
-def luminosity_to_flux(luminosity, redshift=0):
+def luminosity_to_flux(luminosity, dist):
     """
     Convert a luminosity (or luminosity density) to a flux (or flux density).
 
@@ -76,11 +75,10 @@ def luminosity_to_flux(luminosity, redshift=0):
 
     Parameters
     ----------
-    luminosity : float or array of floats
+    luminosity: float or array of floats
         Luminosity (typically in W) or luminosity density (W/nm or W/Hz).
-    redshift :
-        Redshift. If redshift is 0 (the default) the flux at a luminosity
-        distance of 10 pc is returned.
+    dist: float
+        Luminosity distance of the object in metres
 
     Returns
     -------
@@ -88,20 +86,10 @@ def luminosity_to_flux(luminosity, redshift=0):
         The flux (typically in W/m²) of flux density (W/m²/nm or W/m²/Hz).
 
     """
-    if redshift == 0:
-        dist = 10 * parsec
-    else:
-        dist = cosmology.luminosity_distance(redshift) * 1.e6 * parsec
-
-    flux = luminosity / (4 * pi * np.square(dist))
 
-    # astropy 0.3 cosmology functions return quantities
-    try:
-        result = flux.value
-    except AttributeError:
-        result = flux
+    flux = luminosity / (4. * pi * dist * dist)
 
-    return result
+    return flux
 
 
 def lambda_flambda_to_fnu(wavelength, flambda):