__init__.py 41.5 KB
Newer Older
Yannick Roehlly's avatar
Yannick Roehlly committed
1
# -*- coding: utf-8 -*-
2
3
# Copyright (C) 2012, 2013 Centre de données Astrophysiques de Marseille
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
Yannick Roehlly's avatar
Yannick Roehlly committed
4
# Authors: Yannick Roehlly, Médéric Boquien, Laure Ciesla
Yannick Roehlly's avatar
Yannick Roehlly committed
5

6
"""
Yannick Roehlly's avatar
Yannick Roehlly committed
7
8
9
10
11
12
13
14
15
16
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
17
import io
18
import itertools
Yannick Roehlly's avatar
Yannick Roehlly committed
19
20
import numpy as np
from scipy import interpolate
21
import scipy.constants as cst
22
from astropy.table import Table
23
from pcigale.data import (Database, Filter, M2005, BC03, Fritz2006,
24
                          Dale2014, DL2007, DL2014, NebularLines,
25
                          NebularContinuum, SKIRTOR2016, Schreiber2016, THEMIS)
Yannick Roehlly's avatar
Yannick Roehlly committed
26
27


28
29
30
31
32
def read_bc03_ssp(filename):
    """Read a Bruzual and Charlot 2003 ASCII SSP file

    The ASCII SSP files of Bruzual and Charlot 2003 have se special structure.
    A vector is stored with the number of values followed by the values
Yannick Roehlly's avatar
Yannick Roehlly committed
33
    separated by a space (or a carriage return). There are the time vector, 5
34
35
36
37
38
39
40
41
42
43
44
45
    (for Chabrier IMF) or 6 lines (for Salpeter IMF) that we don't care of,
    then the wavelength vector, then the luminosity vectors, each followed by
    a 52 value table, then a bunch of other table of information that are also
    in the *colors files.

    Parameters
    ----------
    filename : string

    Returns
    -------
    time_grid: numpy 1D array of floats
Yannick Roehlly's avatar
Yannick Roehlly committed
46
              Vector of the time grid of the SSP in Myr.
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
    wavelength: numpy 1D array of floats
                Vector of the wavelength grid of the SSP in nm.
    spectra: numpy 2D array of floats
             Array containing the SSP spectra, first axis is the wavelength,
             second one is the time.

    """

    def file_structure_generator():
        """Generator used to identify table lines in the SSP file

        In the SSP file, the vectors are store one next to the other, but
        there are 5 informational lines after the time vector. We use this
        generator to the if we are on lines to read or not.
        """
        if "chab" in filename:
            bad_line_number = 5
        else:
            bad_line_number = 6
        yield("data")
        for i in range(bad_line_number):
            yield("bad")
        while True:
            yield("data")

    file_structure = file_structure_generator()
    # Are we in a data line or a bad one.
74
    what_line = next(file_structure)
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
    # Variable conting, in reverse order, the number of value still to
    # read for the read vector.
    counter = 0

    time_grid = []
    full_table = []
    tmp_table = []

    with open(filename) as file_:
        # We read the file line by line.
        for line in file_:
            if what_line == "data":
                # If we are in a "data" line, we analyse each number.
                for item in line.split():
                    if counter == 0:
                        # If counter is 0, then we are not reading a vector
                        # and the first number is the length of the next
                        # vector.
                        counter = int(item)
                    else:
                        # If counter > 0, we are currently reading a vector.
                        tmp_table.append(float(item))
                        counter -= 1
                        if counter == 0:
                            # We reached the end of the vector. If we have not
                            # yet store the time grid (the first table) we are
                            # currently reading it.
                            if time_grid == []:
                                time_grid = tmp_table[:]
                            # Else, we store the vector in the full table,
                            # only if its length is superior to 250 to get rid
                            # of the 52 item unknown vector and the 221 (time
                            # grid length) item vectors at the end of the
                            # file.
                            elif len(tmp_table) > 250:
                                full_table.append(tmp_table[:])

                            tmp_table = []

            # If at the end of a line, we have finished reading a vector, it's
            # time to change to the next structure context.
            if counter == 0:
117
                what_line = next(file_structure)
118

Yannick Roehlly's avatar
Yannick Roehlly committed
119
    # The time grid is in year, we want Myr.
120
    time_grid = np.array(time_grid, dtype=float)
121
    time_grid *= 1.e-6
122
123
124
125

    # The first "long" vector encountered is the wavelength grid. The value
    # are in Ångström, we convert it to nano-meter.
    wavelength = np.array(full_table.pop(0), dtype=float)
126
    wavelength *= 0.1
127
128
129
130

    # The luminosities are in Solar luminosity (3.826.10^33 ergs.s-1) per
    # Ångström, we convert it to W/nm.
    luminosity = np.array(full_table, dtype=float)
131
    luminosity *= 3.826e27
132
133
134
135
136
137
138
139
    # Transposition to have the time in the second axis.
    luminosity = luminosity.transpose()

    # In the SSP, the time grid begins at 0, but not in the *colors file, so
    # we remove t=0 from the SSP.
    return time_grid[1:], wavelength, luminosity[:, 1:]


140
def build_filters(base):
141
    filters = []
142
    filters_dir = os.path.join(os.path.dirname(__file__), 'filters/')
Yannick Roehlly's avatar
Yannick Roehlly committed
143
144
145
146
147
148
149
150
151
    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()
152

Yannick Roehlly's avatar
Yannick Roehlly committed
153
154
155
        # We convert the wavelength from Å to nm.
        filter_table[0] *= 0.1

156
157
158
159
160
161
162
        # We convert to energy if needed
        if filter_type == 'photon':
            filter_table[1] *= filter_table[0]
        elif filter_type != 'energy':
            raise ValueError("Filter transmission type can only be "
                             "'energy' or 'photon'.")

Yannick Roehlly's avatar
Yannick Roehlly committed
163
164
165
        print("Importing %s... (%s points)" % (filter_name,
                                               filter_table.shape[1]))

166
        new_filter = Filter(filter_name, filter_description, filter_table)
Yannick Roehlly's avatar
Yannick Roehlly committed
167

168
169
170
        # We normalise the filter and compute the pivot wavelength. If the
        # filter is a pseudo-filter used to compute line fluxes, it should not
        # be normalised.
171
172
        if not (filter_name.startswith('PSEUDO') or
                filter_name.startswith('linefilter')):
173
174
            new_filter.normalise()
        else:
175
            new_filter.pivot_wavelength = np.mean(
176
177
                filter_table[0][filter_table[1] > 0]
            )
178
        filters.append(new_filter)
Yannick Roehlly's avatar
Yannick Roehlly committed
179

180
    base.add_filters(filters)
Yannick Roehlly's avatar
Yannick Roehlly committed
181

Médéric Boquien's avatar
Médéric Boquien committed
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
def build_filters_gazpar(base):
    filters = []
    filters_dir = os.path.join(os.path.dirname(__file__), 'filters_gazpar/')
    for filter_file in glob.glob(filters_dir + '**/*.pb', recursive=True):
        with open(filter_file, 'r') as filter_file_read:
            _ = filter_file_read.readline() # We use the filename for the name
            filter_type = filter_file_read.readline().strip('# \n\t')
            _ = filter_file_read.readline() # We do not yet use the calib type
            filter_desc = filter_file_read.readline().strip('# \n\t')

        filter_name = filter_file.replace(filters_dir, '')[:-3]
        filter_name = filter_name.replace('/', '.')

        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()
199

Médéric Boquien's avatar
Médéric Boquien committed
200
201
202
        # We convert the wavelength from Å to nm.
        filter_table[0] *= 0.1

203
204
205
206
207
208
209
        # We convert to energy if needed
        if filter_type == 'photon':
            filter_table[1] *= filter_table[0]
        elif filter_type != 'energy':
            raise ValueError("Filter transmission type can only be "
                             "'energy' or 'photon'.")

Médéric Boquien's avatar
Médéric Boquien committed
210
211
212
        print("Importing %s... (%s points)" % (filter_name,
                                               filter_table.shape[1]))

213
        new_filter = Filter(filter_name, filter_desc, filter_table)
Médéric Boquien's avatar
Médéric Boquien committed
214

215
216
217
        # We normalise the filter and compute the pivot wavelength. If the
        # filter is a pseudo-filter used to compute line fluxes, it should not
        # be normalised.
Médéric Boquien's avatar
Médéric Boquien committed
218
219
220
        if not filter_name.startswith('PSEUDO'):
            new_filter.normalise()
        else:
221
            new_filter.pivot_wavelength = np.mean(
Médéric Boquien's avatar
Médéric Boquien committed
222
223
224
225
226
                filter_table[0][filter_table[1] > 0]
            )
        filters.append(new_filter)

    base.add_filters(filters)
227
228
229

def build_m2005(base):
    m2005_dir = os.path.join(os.path.dirname(__file__), 'maraston2005/')
Yannick Roehlly's avatar
Yannick Roehlly committed
230

Yannick Roehlly's avatar
Yannick Roehlly committed
231
    # Age grid (1 Myr to 13.7 Gyr with 1 Myr step)
232
233
    time_grid = np.arange(1, 13701)
    fine_time_grid = np.linspace(0.1, 13700, 137000)
Yannick Roehlly's avatar
Yannick Roehlly committed
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248

    # 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:
249
            imf = 'krou'
Yannick Roehlly's avatar
Yannick Roehlly committed
250
251
            mass_table = np.copy(kroupa_mass)
        elif 'ssz' in spec_file:
252
            imf = 'salp'
Yannick Roehlly's avatar
Yannick Roehlly committed
253
254
255
256
257
            mass_table = np.copy(salpeter_mass)
        else:
            raise ValueError('Unknown IMF!!!')

        # Keep only the actual metallicity values in the mass table
258
259
260
261
        # we don't take the first column which contains metallicity.
        # We also eliminate the turn-off mas which makes no send for composite
        # populations.
        mass_table = mass_table[1:7, mass_table[0] == metallicity]
Yannick Roehlly's avatar
Yannick Roehlly committed
262

263
264
265
266
267
268
269
270
        # 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)
Yannick Roehlly's avatar
Yannick Roehlly committed
271

272
273
274
        # Extract the age and convert from Gyr to Myr
        ssp_time = np.unique(spec_table[0]) * 1e3
        spec_table = spec_table[1:]
Yannick Roehlly's avatar
Yannick Roehlly committed
275
276

        # Remove the metallicity column from the spec table
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
        spec_table = spec_table[1:]

        # Extract the wavelength and convert from Å to nm
        ssp_wave = spec_table[0][:1221] * 0.1
        spec_table = spec_table[1:]

        # 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

        # 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)
298

299
300
301
        # 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
302
        ssp_wave_resamp = np.around(np.logspace(np.log10(10000),
303
                                                   np.log10(160000), 50))
304
305
306
307
        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:, :]),
308
                                    assume_sorted=True,
309
                                    axis=0)(np.log10(ssp_wave_resamp))
310

311
312
313
        ssp_wave = np.hstack([ssp_wave[:argmin+1], ssp_wave_resamp])
        ssp_lumin = np.vstack([ssp_lumin_interp[:argmin+1, :],
                               ssp_lumin_resamp])
314

315
316
317
318
319
        # Use Z value for metallicity, not log([Z/H])
        metallicity = {-1.35: 0.001,
                       -0.33: 0.01,
                       0.0: 0.02,
                       0.35: 0.04}[metallicity]
Yannick Roehlly's avatar
Yannick Roehlly committed
320

321
322
        base.add_m2005(M2005(imf, metallicity, time_grid, ssp_wave,
                             mass_table, ssp_lumin))
Yannick Roehlly's avatar
Yannick Roehlly committed
323
324


325
326
def build_bc2003(base, res):
    bc03_dir = os.path.join(os.path.dirname(__file__), 'bc03/')
327

328
329
    # Time grid (1 Myr to 14 Gyr with 1 Myr step)
    time_grid = np.arange(1, 14000)
330
    fine_time_grid = np.linspace(0.1, 13999, 139990)
331
332
333
334
335
336
337
338
339
340
341
342

    # Metallicities associated to each key
    metallicity = {
        "m22": 0.0001,
        "m32": 0.0004,
        "m42": 0.004,
        "m52": 0.008,
        "m62": 0.02,
        "m72": 0.05
    }

    for key, imf in itertools.product(metallicity, ["salp", "chab"]):
343
344
345
346
347
348
        ssp_filename = "{}bc2003_{}_{}_{}_ssp.ised_ASCII".format(bc03_dir, res,
                                                                 key, imf)
        color3_filename = "{}bc2003_lr_{}_{}_ssp.3color".format(bc03_dir, key,
                                                                imf)
        color4_filename = "{}bc2003_lr_{}_{}_ssp.4color".format(bc03_dir, key,
                                                                imf)
349

350
        print("Importing {}...".format(ssp_filename))
351
352
353
354
355

        # Read the desired information from the color files
        color_table = []
        color3_table = np.genfromtxt(color3_filename).transpose()
        color4_table = np.genfromtxt(color4_filename).transpose()
356
357
358
        color_table.append(color4_table[6])        # Mstar
        color_table.append(color4_table[7])        # Mgas
        color_table.append(10 ** color3_table[5])  # NLy
359
360
361
362
363

        color_table = np.array(color_table)

        ssp_time, ssp_wave, ssp_lumin = read_bc03_ssp(ssp_filename)

364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
        # 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.
        fill_value = (color_table[:, 0], color_table[:, -1])
        color_table = interpolate.interp1d(ssp_time, color_table,
                                           fill_value=fill_value,
                                           bounds_error=False,
                                           assume_sorted=True)(fine_time_grid)
        color_table = np.mean(color_table.reshape(3, -1, 10), axis=-1)

        # 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)
387

388
389
390
391
392
393
394
395
        # 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
        ssp_wave_resamp = np.around(np.logspace(np.log10(10000),
                                                np.log10(160000), 50))
        argmin = np.argmin(10000.-ssp_wave > 0)-1
        ssp_lumin_resamp = 10.**interpolate.interp1d(
                                    np.log10(ssp_wave[argmin:]),
396
                                    np.log10(ssp_lumin_interp[argmin:, :]),
397
398
399
400
                                    assume_sorted=True,
                                    axis=0)(np.log10(ssp_wave_resamp))

        ssp_wave = np.hstack([ssp_wave[:argmin+1], ssp_wave_resamp])
401
402
        ssp_lumin = np.vstack([ssp_lumin_interp[:argmin+1, :],
                               ssp_lumin_resamp])
403

404
        base.add_bc03(BC03(
405
406
407
408
409
410
411
412
            imf,
            metallicity[key],
            time_grid,
            ssp_wave,
            color_table,
            ssp_lumin
        ))

413

414
def build_dale2014(base):
415
    models = []
416
417
418
    dale2014_dir = os.path.join(os.path.dirname(__file__), 'dale2014/')

    # Getting the alpha grid for the templates
419
    d14cal = np.genfromtxt(dale2014_dir + 'dhcal.dat')
420
421
422
    alpha_grid = d14cal[:, 1]

    # Getting the lambda grid for the templates and convert from microns to nm.
423
    first_template = np.genfromtxt(dale2014_dir + 'spectra.0.00AGN.dat')
424
425
    wave = first_template[:, 0] * 1E3

Médéric Boquien's avatar
Médéric Boquien committed
426
427
428
429
    # Getting the stellar emission and interpolate it at the same wavelength
    # grid
    stell_emission_file = np.genfromtxt(dale2014_dir +
                                        'stellar_SED_age13Gyr_tau10Gyr.spec')
430
    # A -> to nm
Médéric Boquien's avatar
Médéric Boquien committed
431
    wave_stell = stell_emission_file[:, 0] * 0.1
432
    # W/A -> W/nm
Médéric Boquien's avatar
Médéric Boquien committed
433
434
    stell_emission = stell_emission_file[:, 1] * 10
    stell_emission_interp = np.interp(wave, wave_stell, stell_emission)
435
436
437
438
439
440
441
442
443
444
445
446

    # The models are in nuFnu and contain stellar emission.
    # We convert this to W/nm and remove the stellar emission.

    # Emission from dust heated by SB
    fraction = 0.0
    filename = dale2014_dir + "spectra.0.00AGN.dat"
    print("Importing {}...".format(filename))
    datafile = open(filename)
    data = "".join(datafile.readlines())
    datafile.close()

447
    for al in range(1, len(alpha_grid)+1, 1):
Médéric Boquien's avatar
Médéric Boquien committed
448
449
450
        lumin_with_stell = np.genfromtxt(io.BytesIO(data.encode()),
                                         usecols=(al))
        lumin_with_stell = pow(10, lumin_with_stell) / wave
451
452
        constant = lumin_with_stell[7] / stell_emission_interp[7]
        lumin = lumin_with_stell - stell_emission_interp * constant
Médéric Boquien's avatar
Médéric Boquien committed
453
454
455
        lumin[lumin < 0] = 0
        lumin[wave < 2E3] = 0
        norm = np.trapz(lumin, x=wave)
456
        lumin /= norm
457

458
        models.append(Dale2014(fraction, alpha_grid[al-1], wave, lumin))
459
    # Emission from dust heated by AGN - Quasar template
460
    filename = dale2014_dir + "shi_agn.regridded.extended.dat"
461
462
    print("Importing {}...".format(filename))

463
464
465
466
    wave, lumin_quasar = np.genfromtxt(filename, unpack=True)
    wave *= 1e3
    lumin_quasar = 10**lumin_quasar / wave
    norm = np.trapz(lumin_quasar, x=wave)
467
    lumin_quasar /= norm
468

469
470
471
    models.append(Dale2014(1.0, 0.0, wave, lumin_quasar))

    base.add_dale2014(models)
472

473

474
def build_dl2007(base):
475
    models = []
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
    dl2007_dir = os.path.join(os.path.dirname(__file__), 'dl2007/')

    qpah = {
        "00": 0.47,
        "10": 1.12,
        "20": 1.77,
        "30": 2.50,
        "40": 3.19,
        "50": 3.90,
        "60": 4.58
    }

    umaximum = ["1e3", "1e4", "1e5", "1e6"]
    uminimum = ["0.10", "0.15", "0.20", "0.30", "0.40", "0.50", "0.70",
                "0.80", "1.00", "1.20", "1.50", "2.00", "2.50", "3.00",
                "4.00", "5.00", "7.00", "8.00", "10.0", "12.0", "15.0",
                "20.0", "25.0"]

494
    # Mdust/MH used to retrieve the dust mass as models as given per atom of H
Médéric Boquien's avatar
Médéric Boquien committed
495
496
    MdMH = {"00": 0.0100, "10": 0.0100, "20": 0.0101, "30": 0.0102,
            "40": 0.0102, "50": 0.0103, "60": 0.0104}
497

498
499
500
501
502
503
504
505
506
507
508
509
510
511
    # Here we obtain the wavelength beforehand to avoid reading it each time.
    datafile = open(dl2007_dir + "U{}/U{}_{}_MW3.1_{}.txt".format(umaximum[0],
                                                                  umaximum[0],
                                                                  umaximum[0],
                                                                  "00"))
    data = "".join(datafile.readlines()[-1001:])
    datafile.close()

    wave = np.genfromtxt(io.BytesIO(data.encode()), usecols=(0))
    # For some reason wavelengths are decreasing in the model files
    wave = wave[::-1]
    # We convert wavelengths from μm to nm
    wave *= 1000.

512
513
    # Conversion factor from Jy cm² sr¯¹ H¯¹ to W nm¯¹ (kg of H)¯¹
    conv = 4. * np.pi * 1e-30 / (cst.m_p+cst.m_e) * cst.c / (wave*wave) * 1e9
514
515
516
517
518
519
520
521
522
523
524
525
526
527

    for model in sorted(qpah.keys()):
        for umin in uminimum:
            filename = dl2007_dir + "U{}/U{}_{}_MW3.1_{}.txt".format(umin,
                                                                     umin,
                                                                     umin,
                                                                     model)
            print("Importing {}...".format(filename))
            datafile = open(filename)
            data = "".join(datafile.readlines()[-1001:])
            datafile.close()
            lumin = np.genfromtxt(io.BytesIO(data.encode()), usecols=(2))
            # For some reason fluxes are decreasing in the model files
            lumin = lumin[::-1]
528
529
            # Conversion from Jy cm² sr¯¹ H¯¹to W nm¯¹ (kg of dust)¯¹
            lumin *= conv/MdMH[model]
530

531
            models.append(DL2007(qpah[model], umin, umin, wave, lumin))
532
533
534
535
536
537
538
539
540
541
542
543
544
            for umax in umaximum:
                filename = dl2007_dir + "U{}/U{}_{}_MW3.1_{}.txt".format(umin,
                                                                         umin,
                                                                         umax,
                                                                         model)
                print("Importing {}...".format(filename))
                datafile = open(filename)
                data = "".join(datafile.readlines()[-1001:])
                datafile.close()
                lumin = np.genfromtxt(io.BytesIO(data.encode()), usecols=(2))
                # For some reason fluxes are decreasing in the model files
                lumin = lumin[::-1]

545
546
                # Conversion from Jy cm² sr¯¹ H¯¹to W nm¯¹ (kg of dust)¯¹
                lumin *= conv/MdMH[model]
547

548
549
                models.append(DL2007(qpah[model], umin, umax, wave, lumin))
    base.add_dl2007(models)
550
551


552
def build_dl2014(base):
553
    models = []
554
555
    dl2014_dir = os.path.join(os.path.dirname(__file__), 'dl2014/')

Médéric Boquien's avatar
Médéric Boquien committed
556
557
558
    qpah = {"000": 0.47, "010": 1.12, "020": 1.77, "030": 2.50, "040": 3.19,
            "050": 3.90, "060": 4.58, "070": 5.26, "080": 5.95, "090": 6.63,
            "100": 7.32}
559
560
561
562
563
564
565

    uminimum = ["0.100", "0.120", "0.150", "0.170", "0.200", "0.250", "0.300",
                "0.350", "0.400", "0.500", "0.600", "0.700", "0.800", "1.000",
                "1.200", "1.500", "1.700", "2.000", "2.500", "3.000", "3.500",
                "4.000", "5.000", "6.000", "7.000", "8.000", "10.00", "12.00",
                "15.00", "17.00", "20.00", "25.00", "30.00", "35.00", "40.00",
                "50.00"]
566

567
568
569
570
    alpha = ["1.0", "1.1", "1.2", "1.3", "1.4", "1.5", "1.6", "1.7", "1.8",
             "1.9", "2.0", "2.1", "2.2", "2.3", "2.4", "2.5", "2.6", "2.7",
             "2.8", "2.9", "3.0"]

571
    # Mdust/MH used to retrieve the dust mass as models as given per atom of H
Médéric Boquien's avatar
Médéric Boquien committed
572
573
574
    MdMH = {"000": 0.0100, "010": 0.0100, "020": 0.0101, "030": 0.0102,
            "040": 0.0102, "050": 0.0103, "060": 0.0104, "070": 0.0105,
            "080": 0.0106, "090": 0.0107, "100": 0.0108}
575

576
577
578
579
580
581
582
583
584
585
586
587
588
    # Here we obtain the wavelength beforehand to avoid reading it each time.
    datafile = open(dl2014_dir + "U{}_{}_MW3.1_{}/spec_1.0.dat"
                    .format(uminimum[0], uminimum[0], "000"))

    data = "".join(datafile.readlines()[-1001:])
    datafile.close()

    wave = np.genfromtxt(io.BytesIO(data.encode()), usecols=(0))
    # For some reason wavelengths are decreasing in the model files
    wave = wave[::-1]
    # We convert wavelengths from μm to nm
    wave *= 1000.

589
590
    # Conversion factor from Jy cm² sr¯¹ H¯¹ to W nm¯¹ (kg of H)¯¹
    conv = 4. * np.pi * 1e-30 / (cst.m_p+cst.m_e) * cst.c / (wave*wave) * 1e9
591
592
593
594
595
596
597
598
599
600
601
602

    for model in sorted(qpah.keys()):
        for umin in uminimum:
            filename = (dl2014_dir + "U{}_{}_MW3.1_{}/spec_1.0.dat"
                        .format(umin, umin, model))
            print("Importing {}...".format(filename))
            with open(filename) as datafile:
                data = "".join(datafile.readlines()[-1001:])
            lumin = np.genfromtxt(io.BytesIO(data.encode()), usecols=(2))
            # For some reason fluxes are decreasing in the model files
            lumin = lumin[::-1]

603
604
            # Conversion from Jy cm² sr¯¹ H¯¹to W nm¯¹ (kg of dust)¯¹
            lumin *= conv/MdMH[model]
605

606
            models.append(DL2014(qpah[model], umin, umin, 1.0, wave, lumin))
607
608
609
610
611
612
613
614
615
616
            for al in alpha:
                filename = (dl2014_dir + "U{}_1e7_MW3.1_{}/spec_{}.dat"
                            .format(umin, model, al))
                print("Importing {}...".format(filename))
                with open(filename) as datafile:
                    data = "".join(datafile.readlines()[-1001:])
                lumin = np.genfromtxt(io.BytesIO(data.encode()), usecols=(2))
                # For some reason fluxes are decreasing in the model files
                lumin = lumin[::-1]

617
618
                # Conversion from Jy cm² sr¯¹ H¯¹to W nm¯¹ (kg of dust)¯¹
                lumin *= conv/MdMH[model]
619

620
                models.append(DL2014(qpah[model], umin, 1e7, al, wave, lumin))
621

622
    base.add_dl2014(models)
623

624
def build_fritz2006(base):
625
    models = []
626
    fritz2006_dir = os.path.join(os.path.dirname(__file__), 'fritz2006/')
627

628
629
    # Parameters of Fritz+2006
    psy = [0.001, 10.100, 20.100, 30.100, 40.100, 50.100, 60.100, 70.100,
630
631
           80.100, 89.990]  # Viewing angle in degrees
    opening_angle = ["20", "40", "60"]  # Theta = 2*(90 - opening_angle)
632
633
634
    gamma = ["0.0", "2.0", "4.0", "6.0"]
    beta = ["-1.00", "-0.75", "-0.50", "-0.25", "0.00"]
    tau = ["0.1", "0.3", "0.6", "1.0", "2.0", "3.0", "6.0", "10.0"]
635
    r_ratio = ["10", "30", "60", "100", "150"]
636
637

    # Read and convert the wavelength
638
639
640
    datafile = open(fritz2006_dir + "ct{}al{}be{}ta{}rm{}.tot"
                    .format(opening_angle[0], gamma[0], beta[0], tau[0],
                            r_ratio[0]))
641
642
643
644
    data = "".join(datafile.readlines()[-178:])
    datafile.close()
    wave = np.genfromtxt(io.BytesIO(data.encode()), usecols=(0))
    wave *= 1e3
Médéric Boquien's avatar
Médéric Boquien committed
645
    # Number of wavelengths: 178; Number of comments lines: 28
646
647
648
    nskip = 28
    blocksize = 178

649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
    iter_params = ((oa, gam, be, ta, rm)
                   for oa in opening_angle
                   for gam in gamma
                   for be in beta
                   for ta in tau
                   for rm in r_ratio)

    for params in iter_params:
        filename = fritz2006_dir + "ct{}al{}be{}ta{}rm{}.tot".format(*params)
        print("Importing {}...".format(filename))
        try:
            datafile = open(filename)
        except IOError:
            continue
        data = datafile.readlines()
        datafile.close()

666
        # Get intrinsic (de-reddened) AGN luminosity, i.e.,
667
668
669
670
671
672
673
674
675
676
677
678
        # central AGN luminosity at psy=89.990
        n = len(psy)-1
        block = data[nskip + blocksize * n + 4 * (n + 1) - 1:
                     nskip + blocksize * (n+1) + 4 * (n + 1) - 1]
        foo1, foo2, lumin_intrin_agn_unnormed = np.genfromtxt(
                io.BytesIO("".join(block).encode()), usecols=(2, 3, 4),
                unpack=True)
        # Remove NaN for AGN luminosity
        lumin_intrin_agn_unnormed = np.nan_to_num(lumin_intrin_agn_unnormed)
        # Conversion from erg/s/microns to W/nm
        lumin_intrin_agn_unnormed *= 1e-4

679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
        for n in range(len(psy)):
            block = data[nskip + blocksize * n + 4 * (n + 1) - 1:
                         nskip + blocksize * (n+1) + 4 * (n + 1) - 1]
            lumin_therm, lumin_scatt, lumin_agn = np.genfromtxt(
                io.BytesIO("".join(block).encode()), usecols=(2, 3, 4),
                unpack=True)
            # Remove NaN
            lumin_therm = np.nan_to_num(lumin_therm)
            lumin_scatt = np.nan_to_num(lumin_scatt)
            lumin_agn = np.nan_to_num(lumin_agn)
            # Conversion from erg/s/microns to W/nm
            lumin_therm *= 1e-4
            lumin_scatt *= 1e-4
            lumin_agn *= 1e-4
            # Normalization of the lumin_therm to 1W
            norm = np.trapz(lumin_therm, x=wave)
695
696
697
            lumin_therm /= norm
            lumin_scatt /= norm
            lumin_agn /= norm
698
699
            # The intrinsic (de-reddened) AGN luminosity
            lumin_intrin_agn = lumin_intrin_agn_unnormed/norm
700

701
            models.append(Fritz2006(params[4], params[3], params[2],
702
                                         params[1], params[0], psy[n], wave,
703
704
                                         lumin_therm, lumin_scatt, lumin_agn,
                                         lumin_intrin_agn))
705

706
    base.add_fritz2006(models)
707

708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
def build_skirtor2016(base):
    models = []
    skirtor2016_dir = os.path.join(os.path.dirname(__file__), 'skirtor2016/')

    files = glob.glob(skirtor2016_dir + '/*')
    files = [file.split('/')[-1] for file in files]
    params = [f.split('_')[:-1] for f in files]

    # Parameters of SKIRTOR 2016
    t = list({param[0][1:] for param in params})
    p = list({param[1][1:] for param in params})
    q = list({param[2][1:] for param in params})
    oa = list({param[3][2:] for param in params})
    R = list({param[4][1:] for param in params})
    Mcl = list({param[5][3:] for param in params})
    i = list({param[6][1:] for param in params})

    iter_params = ((p1, p2, p3, p4, p5, p6, p7)
                   for p1 in t
                   for p2 in p
                   for p3 in q
                   for p4 in oa
                   for p5 in R
                   for p6 in Mcl
                   for p7 in i)

734
735
    # Extract the intrinsic AGN disk luminosity
    filename = skirtor2016_dir+"t9_p1_q1_oa80_R30_Mcl0.97_i0_sed.dat"
736
737
738
739
740
741
742
743
744
745
746
747
748
749
    wl_uncut, intrin_disk_unnormed = np.genfromtxt(filename, unpack=True, usecols=(0, 6))
    # Convert to nm units
    wl_uncut *= 1e3
    # Convert from F_lam*lam to F_lam
    intrin_disk_unnormed /= wl_uncut

    # Define the new disk SED shape
    # Modify the disk SED to Feltre et al. (2012), following the method in
    # https://sites.google.com/site/skirtorus/sed-library
    # Set the broken power law shape
    # 1st part: < 50 nm
    disk_new_unnormed = wl_uncut.copy()
    # 2nd part: 50-125 nm
    wave_idxs = (wl_uncut>=50) & (wl_uncut<125)
Guang's avatar
Guang committed
750
751
752
    # Here, norm_fac is to re-nomalize for continuity
    norm_fac = 50**1.2
    disk_new_unnormed[wave_idxs] = disk_new_unnormed[wave_idxs]**-0.2 * norm_fac
753
    # 3rd part: 125-1e4 nm
Guang's avatar
Guang committed
754
    norm_fac *= 125**1.3
755
756
757
758
759
760
761
762
763
764
765
    wave_idxs = (wl_uncut>=125) & (wl_uncut<1e4)
    disk_new_unnormed[wave_idxs] = disk_new_unnormed[wave_idxs]**-1.5 * norm_fac
    # 4th part: >1e4 nm
    norm_fac *= 1e4**2.5
    wave_idxs = wl_uncut>=1e4
    disk_new_unnormed[wave_idxs] = disk_new_unnormed[wave_idxs]**-4 * norm_fac
    # Normalize the new disk SED, so that its total energy is 1
    disk_new_unnormed /= np.trapz(disk_new_unnormed, x=wl_uncut)

    # Change the SED shape of intrinsic disk luminosity
    intrin_disk_unnormed = disk_new_unnormed * np.trapz(intrin_disk_unnormed, x=wl_uncut)
766
    # Convert anisotropic to the isotropic disk luminosity
767
    intrin_disk_unnormed *= 7/18
768

769
770
771
772
    for params in iter_params:
        filename = skirtor2016_dir + \
                "t{}_p{}_q{}_oa{}_R{}_Mcl{}_i{}_sed.dat".format(*params)
        print("Importing {}...".format(filename))
773
774
        # Convert some useful parameters to float
        i_float, oa_float = float(params[6]), float(params[3])
775

776
777
778
779
        disk, scatt, dust = np.genfromtxt(filename, unpack=True, usecols=(2, 3, 4))
        disk /= wl_uncut
        scatt /= wl_uncut
        dust /= wl_uncut
780
781

        # Normalization of the lumin_therm to 1W
782
        norm = np.trapz(dust, x=wl_uncut)
783
        disk /= norm
784
        scatt /= norm
785
        dust /= norm
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
        intrin_disk = intrin_disk_unnormed / norm

        if i_float<=(90-oa_float):
            # For type-1 AGN, update the disk emission
            # Re-normalize the entire new disk SED to keep energy conservation
            disk = np.trapz(disk, x=wl_uncut) * disk_new_unnormed

        # Add the scatter component to disk for simplicity
        disk += scatt

        # Apply wavelength cut to avoid X-ray wavelength
        lam_cut = 5.
        lam_idxs = wl_uncut>=lam_cut
        # Calculate the re-normalization factor to keep energy conservation
        norm_fac = np.trapz(intrin_disk, x=wl_uncut) / \
                   np.trapz(intrin_disk[lam_idxs], x=wl_uncut[lam_idxs])
        # Perform the cut
        wl = wl_uncut[lam_idxs]
        disk = disk[lam_idxs]*norm_fac
        dust = dust[lam_idxs]
        intrin_disk = intrin_disk[lam_idxs]*norm_fac

        # Modify the model so that the central source is isotropic
        # Calculate the coverting factors disk emission
        # (dust emission is normalized to unity, so it remains same)
        disk_fac = 7/(12*np.cos(i_float*np.pi/180)**2 + 6*np.cos(i_float*np.pi/180))
        dust_fac = 7/(4*np.sin(oa_float*np.pi/180)**2 + 3*np.sin(oa_float*np.pi/180))
        # Only apply disk_fac to type 1 (for type 2 it is zero)
        if i_float<=(90-oa_float): disk*= disk_fac/dust_fac
        intrin_disk /= dust_fac
816
817
818

        models.append(SKIRTOR2016(params[0], params[1], params[2], params[3],
                                  params[4], params[5], params[6], wl, disk,
819
                                  dust, intrin_disk))
820
821
822

    base.add_skirtor2016(models)

823
def build_nebular(base):
824
825
    models_lines = []
    models_cont = []
826

827
828
829
    nebular_dir = os.path.join(os.path.dirname(__file__), 'nebular/')
    print("Importing {}...".format(nebular_dir + 'lines.dat'))
    lines = np.genfromtxt(nebular_dir + 'lines.dat')
830

831
832
833
834
    tmp = Table.read(nebular_dir + 'line_wavelengths.dat', format='ascii')
    wave_lines = tmp['col1'].data
    name_lines = tmp['col2'].data

835
836
    print("Importing {}...".format(nebular_dir + 'continuum.dat'))
    cont = np.genfromtxt(nebular_dir + 'continuum.dat')
837

838
839
840
    # Convert wavelength from Å to nm
    wave_lines *= 0.1
    wave_cont = cont[:1600, 0] * 0.1
841

842
843
    # Get the list of metallicities
    metallicities = np.unique(lines[:, 1])
844

845
846
847
    # Keep only the fluxes
    lines = lines[:, 2:]
    cont = cont[:, 1:]
848

849
850
851
    # We select only models with ne=100. Other values could be included later
    lines = lines[:, 1::3]
    cont = cont[:, 1::3]
852

853
854
    # Convert lines to W and to a linear scale
    lines = 10**(lines-7)
855

856
857
858
    # Convert continuum to W/nm
    cont *= np.tile(1e-7 * cst.c * 1e9 / wave_cont**2,
                    metallicities.size)[:, np.newaxis]
859

860
861
862
863
864
    # 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):
865
866
            models_lines.append(NebularLines(metallicity, logU, name_lines,
                                             wave_lines, spectrum))
867

868
869
870
871
872
873
874
    # 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))
875

876
    base.add_nebular_lines(models_lines)
877
    base.add_nebular_continuum(models_cont)
878

879

880
881
def build_schreiber2016(base):
    models = []
882
883
884
885
    schreiber2016_dir = os.path.join(os.path.dirname(__file__),
                                     'schreiber2016/')

    print("Importing {}...".format(schreiber2016_dir + 'g15_pah.fits'))
886
    pah = Table.read(schreiber2016_dir + 'g15_pah.fits')
887
    print("Importing {}...".format(schreiber2016_dir + 'g15_dust.fits'))
888
889
    dust = Table.read(schreiber2016_dir + 'g15_dust.fits')

890
    # Getting the lambda grid for the templates and convert from μm to nm.
891
    wave = dust['LAM'][0, 0, :].data * 1e3
892
893

    for td in np.arange(15., 100.):
894
        # Find the closest temperature in the model list of tdust
895
        tsed = np.argmin(np.absolute(dust['TDUST'][0].data-td))
896
897

        # The models are in νFν.  We convert this to W/nm.
898
899
        lumin_dust = dust['SED'][0, tsed, :].data / wave
        lumin_pah = pah['SED'][0, tsed, :].data / wave
900
901
902
903

        models.append(Schreiber2016(0, td, wave, lumin_dust))
        models.append(Schreiber2016(1, td, wave, lumin_pah))

904
905
    base.add_schreiber2016(models)

906

907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
def build_themis(base):
    models = []
    themis_dir = os.path.join(os.path.dirname(__file__), 'themis/')

    # Mass fraction of hydrocarbon solids i.e., a-C(:H) smaller than 1.5 nm,
    # also known as HAC
    qhac = {"000": 0.02, "010": 0.06, "020": 0.10, "030": 0.14, "040": 0.17,
            "050": 0.20, "060": 0.24, "070": 0.28, "080": 0.32, "090": 0.36,
            "100": 0.40}

    uminimum = ["0.100", "0.120", "0.150", "0.170", "0.200", "0.250", "0.300",
                "0.350", "0.400", "0.500", "0.600", "0.700", "0.800", "1.000",
                "1.200", "1.500", "1.700", "2.000", "2.500", "3.000", "3.500",
                "4.000", "5.000", "6.000", "7.000", "8.000", "10.00", "12.00",
                "15.00", "17.00", "20.00", "25.00", "30.00", "35.00", "40.00",
                "50.00", "80.00"]

    alpha = ["1.0", "1.1", "1.2", "1.3", "1.4", "1.5", "1.6", "1.7", "1.8",
             "1.9", "2.0", "2.1", "2.2", "2.3", "2.4", "2.5", "2.6", "2.7",
             "2.8", "2.9", "3.0"]

    # Mdust/MH used to retrieve the dust mass as models as given per atom of H
    MdMH = {"000": 7.4e-3, "010": 7.4e-3, "020": 7.4e-3, "030": 7.4e-3,
            "040": 7.4e-3, "050": 7.4e-3, "060": 7.4e-3, "070": 7.4e-3,
            "080": 7.4e-3, "090": 7.4e-3, "100": 7.4e-3}

    # Here we obtain the wavelength beforehand to avoid reading it each time.
    datafile = open(themis_dir + "U{}_{}_MW3.1_{}/spec_1.0.dat"
                    .format(uminimum[0], uminimum[0], "000"))

    data = "".join(datafile.readlines()[-576:])
    datafile.close()

    wave = np.genfromtxt(io.BytesIO(data.encode()), usecols=(0))

    # We convert wavelengths from μm to nm
    wave *= 1000.

    # Conversion factor from Jy cm² sr¯¹ H¯¹ to W nm¯¹ (kg of H)¯¹
    conv = 4. * np.pi * 1e-30 / (cst.m_p+cst.m_e) * cst.c / (wave*wave) * 1e9

    for model in sorted(qhac.keys()):
        for umin in uminimum:
            filename = (themis_dir + "U{}_{}_MW3.1_{}/spec_1.0.dat"
                        .format(umin, umin, model))
            print("Importing {}...".format(filename))
            with open(filename) as datafile:
                data = "".join(datafile.readlines()[-576:])
            lumin = np.genfromtxt(io.BytesIO(data.encode()), usecols=(2))

            # Conversion from Jy cm² sr¯¹ H¯¹to W nm¯¹ (kg of dust)¯¹
            lumin *= conv / MdMH[model]

            models.append(THEMIS(qhac[model], umin, umin, 1.0, wave, lumin))
            for al in alpha:
                filename = (themis_dir + "U{}_1e7_MW3.1_{}/spec_{}.dat"
                            .format(umin, model, al))
                print("Importing {}...".format(filename))
                with open(filename) as datafile:
                    data = "".join(datafile.readlines()[-576:])
                lumin = np.genfromtxt(io.BytesIO(data.encode()), usecols=(2))

                # Conversion from Jy cm² sr¯¹ H¯¹to W nm¯¹ (kg of dust)¯¹
                lumin *= conv/MdMH[model]

                models.append(THEMIS(qhac[model], umin, 1e7, al, wave, lumin))

    base.add_themis(models)


977
def build_base(bc03res='lr'):
978
979
980
981
982
983
    base = Database(writable=True)
    base.upgrade_base()

    print('#' * 78)
    print("1- Importing filters...\n")
    build_filters(base)
Médéric Boquien's avatar
Médéric Boquien committed
984
    build_filters_gazpar(base)
985
986
987
988
989
990
991
992
993
    print("\nDONE\n")
    print('#' * 78)

    print("2- Importing Maraston 2005 SSP\n")
    build_m2005(base)
    print("\nDONE\n")
    print('#' * 78)

    print("3- Importing Bruzual and Charlot 2003 SSP\n")
994
    build_bc2003(base, bc03res)
995
996
997
    print("\nDONE\n")
    print('#' * 78)

998
    print("4- Importing Draine and Li (2007) models\n")
999
1000
    build_dl2007(base)
    print("\nDONE\n")