Import subversion code

README.rst

+Python Code Investigating Galaxy Emission

database_builder/__init__.py

+# -*- coding: utf-8 -*-
+"""
+Copyright (C) 2012 Centre de données Astrophysiques de Marseille
+Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
+
+@author: Yannick Roehlly <yannick.roehlly@oamp.fr>
+
+This script is used the build pcigale internal database containing:
+- The various filter transmission tables;
+- The Maraston 2005 single stellar population (SSP) data;
+- The Dale and Helou 2002 infra-red templates.
+
+"""
+import sys
+import os
+sys.path.append(os.path.join(os.path.dirname(__file__), '../'))
+import glob
+import numpy as np
+from scipy import interpolate
+from pcigale.data import Database, Filter, SspM2005
+
+
+filters_dir = os.path.join(os.path.dirname(__file__), 'filters/')
+m2005_dir = os.path.join(os.path.dirname(__file__), 'maraston2005/')
+dh2002_dir = os.path.join(os.path.dirname(__file__), 'dh2002/')
+
+
+def build_base():
+    base = Database(writable=True)
+    base.upgrade_base()
+
+    print('#' * 78)
+    ########################################################################
+    # Filter transmission tables insertion                                 #
+    ########################################################################
+    print("1- Importing filters...\n")
+    for filter_file in glob.glob(filters_dir + '*.dat'):
+        with open(filter_file, 'r') as filter_file_read:
+            filter_name = filter_file_read.readline().strip('# \n\t')
+            filter_type = filter_file_read.readline().strip('# \n\t')
+            filter_description = filter_file_read.readline().strip('# \n\t')
+        filter_table = np.genfromtxt(filter_file)
+        # The table is transposed to have table[0] containing the wavelength
+        # and table[1] containing the transmission.
+        filter_table = filter_table.transpose()
+        # We convert the wavelength from Å to nm.
+        filter_table[0] *= 0.1
+
+        print("Importing %s... (%s points)" % (filter_name,
+                                               filter_table.shape[1]))
+
+        new_filter = Filter(filter_name, filter_description,
+                            filter_type, filter_table)
+
+        # We normalise the filter and compute the effective wavelength.
+        new_filter.normalise()
+
+        base.add_filter(new_filter)
+    print("\nDONE\n")
+    print('#' * 78)
+
+    ########################################################################
+    # Maraston 2005 SSP insertion                                          #
+    ########################################################################
+    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)
+
+    # Transpose the table to have access to each value vector on the first
+    # axis
+    kroupa_mass = np.genfromtxt(m2005_dir + 'stellarmass.kroupa').transpose()
+    salpeter_mass = \
+        np.genfromtxt(m2005_dir + '/stellarmass.salpeter').transpose()
+
+    for spec_file in glob.glob(m2005_dir + '*.rhb'):
+
+        print("Importing %s..." % spec_file)
+
+        spec_table = np.genfromtxt(spec_file).transpose()
+        metallicity = spec_table[1, 0]
+
+        if 'krz' in spec_file:
+            imf = 'kr'
+            mass_table = np.copy(kroupa_mass)
+        elif 'ssz' in spec_file:
+            imf = 'ss'
+            mass_table = np.copy(salpeter_mass)
+        else:
+            raise ValueError('Unknown IMF!!!')
+
+        # Keep only the actual metallicity values in the mass table
+        # 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)
+
+        # Remove the age column from the mass table
+        mass_table = np.delete(mass_table, 0, 0)
+
+        # 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
+
+        # For all ages, the lambda grid is the same.
+        lambda_grid = np.unique(spec_table[1])
+
+        # Creation of the age vs lambda flux table
+        tmpList = []
+        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
+            flux_regrid = interpolate.interp1d(age_grid_orig,
+                                               flux_orig)(age_grid)
+
+            tmpList.append(flux_regrid)
+        flux_age = np.array(tmpList)
+
+        base.add_ssp_m2005(SspM2005(imf, metallicity, age_grid,
+                                    lambda_grid, mass_table, flux_age))
+
+    print("\nDONE\n")
+    print('#' * 78)
+
+    ########################################################################
+    # Dale and Helou 2002 templates insertion                              #
+    ########################################################################
+    print("3- Importing Dale and Helou 2002 templates\n")
+
+    # Getting the alpha grid for the templates
+    dhcal = np.genfromtxt(dh2002_dir + 'dhcal.dat')
+    alpha_grid = dhcal[:, 1]
+
+    # Getting the lambda grid for the templates (we checked that all share the
+    # same grid).
+    first_template = np.genfromtxt(dh2002_dir + 'irdh01.spec', skip_header=1)
+    lambda_grid = first_template[:, 0] * 0.1  # Convert Å to nm
+
+    templates = []
+
+    for i in range(len(alpha_grid)):
+        filename = dh2002_dir + 'irdh' + ("%02d" % (i + 1)) + '.spec'
+        print("Importing %s..." % filename)
+        table = np.genfromtxt(filename, skip_header=1)[:, 1]  # Luminosity
+                                                              # column
+        # The table give the luminosity density in Lsun/Å normalised to 1 Lsun
+        # over the full spectrum. As we converted the wavelengths to nm, we
+        # must multiply the density per 10 to keep the normalisation.
+        table = table * 10
+        templates.append(table)
+
+    templates = np.array(templates)
+
+    data = (alpha_grid, lambda_grid, templates)
+
+    base.add_dh2002_infrared_templates(data)
+
+    print("\nDONE\n")
+    print('#' * 78)
+
+    base.session.close_all()
+
+if __name__ == '__main__':
+    build_base()

database_builder/dh2002/Readme.rst

+ The files in this directory are the Dale and Helou (2002) templates as used
+ by the Fortran Cigale code. They are different from the pristine Dale and
+ Helou templates as the stellar contribution has been removed (see Noll et
+ al., 2009).

database_builder/dh2002/dhcal.dat

+# Ndhmod
+# 64
+# Ndh alpha lgfnu60_fnu100 lgLIR_M[Lsol] lgLIR_C[Lsol]
+01 0.0625  0.2135 14.350 14.168                                                 
+02 0.1250  0.2104 14.326 14.148                                                 
+03 0.1875  0.2069 14.300 14.125                                                 
+04 0.2500  0.2032 14.272 14.102                                                 
+05 0.3125  0.1991 14.241 14.077                                                 
+06 0.3750  0.1947 14.208 14.048                                                 
+07 0.4375  0.1898 14.171 14.016                                                 
+08 0.5000  0.1846 14.132 13.984                                                 
+09 0.5625  0.1787 14.089 13.946                                                 
+10 0.6250  0.1723 14.042 13.905                                                 
+11 0.6875  0.1653 13.990 13.860                                                 
+12 0.7500  0.1576 13.934 13.810                                                 
+13 0.8125  0.1490 13.872 13.754                                                 
+14 0.8750  0.1395 13.803 13.693                                                 
+15 0.9375  0.1291 13.729 13.626                                                 
+16 1.0000  0.1175 13.647 13.550                                                 
+17 1.0625  0.1046 13.556 13.466                                                 
+18 1.1250  0.0903 13.456 13.372                                                 
+19 1.1875  0.0744 13.346 13.268                                                 
+20 1.2500  0.0567 13.224 13.151                                                 
+21 1.3125  0.0371 13.091 13.021                                                 
+22 1.3750  0.0154 12.944 12.877                                                 
+23 1.4375 -0.0084 12.784 12.717                                                 
+24 1.5000 -0.0345 12.611 12.542                                                 
+25 1.5625 -0.0629 12.423 12.350                                                 
+26 1.6250 -0.0933 12.222 12.143                                                 
+27 1.6875 -0.1257 12.008 11.921                                                 
+28 1.7500 -0.1596 11.785 11.685                                                 
+29 1.8125 -0.1944 11.554 11.434                                                 
+30 1.8750 -0.2298 11.318 11.167                                                 
+31 1.9375 -0.2649 11.081 10.879                                                 
+32 2.0000 -0.2990 10.849 10.551                                                 
+33 2.0625 -0.3316 10.625 10.130                                                 
+34 2.1250 -0.3621 10.412  9.452                                                 
+35 2.1875 -0.3900 10.214  8.321                                                 
+36 2.2500 -0.4151 10.035  7.106                                                 
+37 2.3125 -0.4372  9.876  6.024                                                 
+38 2.3750 -0.4565  9.735  5.079                                                 
+39 2.4375 -0.4731  9.614  4.268                                                 
+40 2.5000 -0.4871  9.511  3.584                                                 
+41 2.5625 -0.4989  9.424  3.009                                                 
+42 2.6250 -0.5086  9.353  2.536                                                 
+43 2.6875 -0.5166  9.294  2.146                                                 
+44 2.7500 -0.5232  9.246  1.825                                                 
+45 2.8125 -0.5284  9.209  1.572                                                 
+46 2.8750 -0.5327  9.178  1.363                                                 
+47 2.9375 -0.5361  9.154  1.197                                                 
+48 3.0000 -0.5388  9.136  1.066                                                 
+49 3.0625 -0.5408  9.123  0.969                                                 
+50 3.1250 -0.5424  9.113  0.891                                                 
+51 3.1875 -0.5435  9.107  0.838                                                 
+52 3.2500 -0.5444  9.103  0.794                                                 
+53 3.3125 -0.5450  9.100  0.765                                                 
+54 3.3750 -0.5453  9.100  0.750                                                 
+55 3.4375 -0.5454  9.102  0.745                                                 
+56 3.5000 -0.5455  9.103  0.740                                                 
+57 3.5625 -0.5454  9.106  0.745                                                 
+58 3.6250 -0.5451  9.110  0.760                                                 
+59 3.6875 -0.5448  9.115  0.774                                                 
+60 3.7500 -0.5445  9.119  0.789                                                 
+61 3.8125 -0.5440  9.125  0.813                                                 
+62 3.8750 -0.5436  9.130  0.832                                                 
+63 3.9375 -0.5430  9.136  0.861                                                 
+64 4.0000 -0.5425  9.142  0.885