aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim/extragalactic/clusters.py
blob: 8249e3a48dcb9bbbea78d33e465fee57512093d1 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
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
74
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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
# Copyright (c) 2016 Weitian LI <liweitianux@live.com>
# MIT license

"""
Simulation of the radio emissions from clusters of galaxies.

NOTE
----
Currently, only radio *halos* are considered with many simplifications.
Radio *relics* simulations need more investigations ...
"""

import os
import logging
from datetime import datetime, timezone

import numpy as np
import astropy.units as au
from astropy.io import fits
from astropy.cosmology import FlatLambdaCDM
import healpy as hp
import pandas as pd

from ..utils import write_fits_healpix
from ..utils.random import spherical_uniform
from ..utils.convert import Fnu_to_Tb
from ..utils.grid import make_grid_ellipse, map_grid_to_healpix


logger = logging.getLogger(__name__)


class GalaxyClusters:
    """
    Simulate the radio emissions from the clusters of galaxies, which
    host radio halos and relics (currently not considered).

    The simulation follows the method adopted by [Jelic2008]_, which uses
    the *ΛCDM deep wedge cluster catalog* derived from the *Hubble Volume
    Project* [HVP]_, [Evard2002]_.

    Every radio cluster is simulated as an *ellipse* of *uniform brightness*
    on a local coordinate grid with relatively higher resolution compared
    to the output HEALPix map, which is then mapped to the output HEALPix
    map by down-sampling, i.e., in a similar way as the simulations of SNRs.

    TODO: ???

    Parameters
    ----------
    configs : `ConfigManager`
        A `ConfigManager` instance containing default and user configurations.
        For more details, see the example configuration specifications.

    Attributes
    ----------
    TODO: ???

    NOTE
    ----
    Currently, only radio *halos* are considered with many simplifications.
    Radio *relics* simulations need more investigations ...

    References
    ----------
    .. [Jelic2008]
       Jelić, V. et al.,
       "Foreground simulations for the LOFAR-epoch of reionization experiment",
       2008, MNRAS, 389, 1319-1335,
       http://adsabs.harvard.edu/abs/2008MNRAS.389.1319J

    .. [Evard2002]
       Evard, A. E. et al.,
       "Galaxy Clusters in Hubble Volume Simulations: Cosmological Constraints
       from Sky Survey Populations",
       2002, ApJ, 573, 7-36,
       http://adsabs.harvard.edu/abs/2002ApJ...573....7E

    .. [EnBlin2002]
       Enßlin, T. A. & Röttgering, H.,
       "The radio luminosity function of cluster radio halos",
       2002, A&A, 396, 83-89,
       http://adsabs.harvard.edu/abs/2002A%26A...396...83E

    .. [Reiprich2002]
       Reiprich, Thomas H. & Böhringer, Hans,
       "The Mass Function of an X-Ray Flux-limited Sample of Galaxy Clusters",
       2002, ApJ, 567, 716-740,
       http://adsabs.harvard.edu/abs/2002ApJ...567..716R
    """
    # Component name
    name = "clusters of galaxies"

    def __init__(self, configs):
        self.configs = configs
        self._set_configs()

    def _set_configs(self):
        """Load the configs and set the corresponding class attributes."""
        comp = "extragalactic/clusters"
        self.catalog_path = self.configs.get_path(comp+"/catalog")
        self.catalog_outfile = self.configs.get_path(comp+"/catalog_outfile")
        self.halo_fraction = self.configs.getn(comp+"/halo_fraction")
        self.resolution = self.configs.getn(comp+"/resolution") * au.arcmin
        self.prefix = self.configs.getn(comp+"/prefix")
        self.save = self.configs.getn(comp+"/save")
        self.output_dir = self.configs.get_path(comp+"/output_dir")
        #
        self.filename_pattern = self.configs.getn("output/filename_pattern")
        self.use_float = self.configs.getn("output/use_float")
        self.clobber = self.configs.getn("output/clobber")
        self.nside = self.configs.getn("common/nside")
        self.freq_unit = au.Unit(self.configs.getn("frequency/unit"))
        # Cosmology model
        self.H0 = self.configs.getn("cosmology/H0")
        self.OmegaM0 = self.configs.getn("cosmology/OmegaM0")
        self.cosmo = FlatLambdaCDM(H0=self.H0, Om0=self.OmegaM0)
        #
        logger.info("Loaded and set up configurations")

    def _load_catalog(self):
        """Load the cluster catalog data and set up its properties."""
        self.catalog = pd.read_csv(self.catalog_path)
        nrow, ncol = self.catalog.shape
        logger.info("Loaded clusters catalog data from: {0}".format(
            self.catalog_path))
        logger.info("Clusters catalog data: {0} objects, {1} columns".format(
            nrow, ncol))
        # Set the properties for this catalog
        self.catalog_prop = {
            "omega_m": 0.3,
            "omega_lambda": 0.7,
            # Dimensionless Hubble constant
            "h": 0.7,
            "sigma8": 0.9,
            # Number of particles
            "n_particle": 1e9,
            # Particle mass of the simulation [ h^-1 ]
            "m_particle": 2.25e12 * au.solMass,
            # Cube side length [ h^-1 ]
            "l_side": 3000.0 * au.Mpc,
            # Overdensity adopted to derive the clusters
            "overdensity": 200,
            # Sky coverage
            "coverage": 10*au.deg * 10*au.deg
        }
        # Units for the catalog columns (also be populated by other methods)
        self.units = {}

    def _save_catalog_inuse(self):
        """Save the effective/inuse clusters catalog data to a CSV file."""
        if self.catalog_outfile is None:
            logger.warning("Catalog output file not set, so do NOT save.")
            return
        # Create directory if necessary
        dirname = os.path.dirname(self.catalog_outfile)
        if not os.path.exists(dirname):
            os.mkdir(dirname)
            logger.info("Created directory: {0}".format(dirname))
        # Save catalog data
        if os.path.exists(self.catalog_outfile):
            if self.clobber:
                os.remove(self.catalog_outfile)
            else:
                raise OSError("Output file already exists: {0}".format(
                    self.catalog_outfile))
        self.catalog.to_csv(self.catalog_outfile, header=True, index=False)
        logger.info("Save clusters catalog in use to: {0}".format(
            self.catalog_outfile))

    def _process_catalog(self):
        """Process the catalog to prepare for the simulation."""
        # Dimensionless Hubble parameter adopted in THIS simulation
        h = self.H0 / 100
        logger.info("Adopted dimensionless Hubble parameter: {0}".format(h))
        # Cluster masses, unit: solMass (NOTE: h dependence)
        self.catalog["mass"] = (self.catalog["m"] *
                                self.catalog_prop["m_particle"].value / h)
        self.units["mass"] = au.solMass
        logger.info("Catalog: calculated cluster masses")
        # Cluster distances from the observer, unit: Mpc
        dist = ((self.catalog["x"]**2 + self.catalog["y"]**2 +
                 self.catalog["z"]**2) ** 0.5 *
                self.catalog_prop["l_side"].value / h)
        self.catalog["distance"] = dist
        self.units["distance"] = au.Mpc
        logger.info("Catalog: calculated cluster distances")
        # Drop unnecessary columns to save memory
        columns_drop = ["m", "sigma", "ip", "x", "y", "z", "vx", "vy", "vz"]
        self.catalog.drop(columns_drop, axis=1, inplace=True)
        logger.info("Catalog: dropped unnecessary columns: {0}".format(
            ", ".join(columns_drop)))

    def _expand_catalog_fullsky(self):
        """Expand the catalog to be full sky, by assuming that clusters are
        uniformly distributed.  Also, the radio halo fraction is also
        considered to determine the final number of clusters on the full sky.
        """
        fullsky = 4*np.pi * au.sr
        factor = float(fullsky / self.catalog_prop["coverage"])
        n0_cluster = len(self.catalog)
        logger.info("Radio halo fraction in clusters: {0}".format(
            self.halo_fraction))
        # Total number of clusters on the full sky
        N_cluster = int(n0_cluster * factor * self.halo_fraction)
        logger.info("Total number of clusters on the full sky: {0:,}".format(
            N_cluster))
        logger.info("Expanding the catalog to be full sky ...")
        idx = np.round(np.random.uniform(low=0, high=n0_cluster-1,
                                         size=N_cluster)).astype(np.int)
        self.catalog = self.catalog.iloc[idx, :]
        self.catalog.reset_index(inplace=True)
        logger.info("Done expand the catalog to be full sky")

    def _add_random_position(self):
        """Add random positions for each cluster as columns "glon" and
        "glat" to the catalog data.

        Column "glon" is the Galactic longitudes, [0, 360) (degree).
        Column "glat" is the Galactic latitudes, [-90, 90] (degree).

        The positions are uniformly distributed on the spherical surface.
        """
        logger.info("Randomly generating positions for each cluster ...")
        num = len(self.catalog)
        theta, phi = spherical_uniform(num)
        glon = np.degrees(phi)
        glat = 90.0 - np.degrees(theta)
        self.catalog["glon"] = glon
        self.catalog["glat"] = glat
        self.units["coord"] = au.deg
        logger.info("Done add random positions for each cluster")

    def _add_random_eccentricity(self):
        """Add random eccentricities for each cluster as column
        "eccentricity" to the catalog data.

        The eccentricity of a ellipse is defined as:
            e = sqrt((a^2 - b^2) / a^2) = f / a
        where f is the distance from the center to either focus:
            f = sqrt(a^2 - b^2)

        NOTE
        ----
        The eccentricities are randomly generated from a *squared*
        standard normalization distribution, and with an upper limit
        at 0.9, i.e., the eccentricities are [0, 0.9].
        """
        logger.info("Adding random eccentricities for each cluster ...")
        num = len(self.catalog)
        eccentricity = np.random.normal(size=num) ** 2
        # Replace values beyond the upper limit by sampling from valid values
        ulimit = 0.9
        idx_invalid = (eccentricity > ulimit)
        num_invalid = idx_invalid.sum()
        eccentricity[idx_invalid] = np.random.choice(
            eccentricity[~idx_invalid], size=num_invalid)
        self.catalog["eccentricity"] = eccentricity
        logger.info("Done add random eccentricities to catalog")

    def _add_random_rotation(self):
        """Add random rotation angles for each cluster as column "rotation"
        to the catalog data.

        The rotation angles are uniformly distributed within [0, 360).

        The rotation happens on the spherical surface, i.e., not with respect
        to the line of sight, but to the Galactic frame coordinate axes.
        """
        logger.info("Adding random rotation angles for each cluster ...")
        num = len(self.catalog)
        rotation = np.random.uniform(low=0.0, high=360.0, size=num)
        self.catalog["rotation"] = rotation
        self.units["rotation"] = au.deg
        logger.info("Done add random rotation angles to catalog")

    def _calc_sizes(self):
        """Calculate the virial radii for each cluster from the masses,
        and then calculate the elliptical angular sizes by considering
        the added random eccentricities.

        Attributes
        ----------
        catalog["r_vir"] : 1D `~numpy.ndarray`
            The virial radii (unit: Mpc) calculated from the cluster masses
        catalog["size_major"], catalog["size_minor] : 1D `~numpy.ndarray`
            The major and minor axes (unit: degree) of the clusters calculated
            from the above virial radii and the random eccentricities.

        NOTE
        ----
        The elliptical major and minor axes are calculated by assuming
        the equal area between the ellipse and corresponding circle.
            theta2 = r_vir / distance
            pi * a * b = pi * (theta2)^2
            e = sqrt((a^2 - b^2) / a^2)
        thus,
            a = theta2 / (1-e^2)^(1/4)
            b = theta2 * (1-e^2)^(1/4)
        """
        logger.info("Calculating the virial radii ...")
        overdensity = self.catalog_prop["overdensity"]
        rho_crit = self.cosmo.critical_density(self.catalog["redshift"])
        mass = self.catalog["mass"].data * self.units["mass"]
        r_vir = (3 * mass / (4*np.pi * overdensity * rho_crit)) ** (1.0/3.0)
        self.catalog["r_vir"] = r_vir.to(au.Mpc).value
        self.units["r_vir"] = au.Mpc
        logger.info("Done calculate the virial radii")
        # Calculate (elliptical) angular sizes, i.e., major and minor axes
        logger.info("Calculating the elliptical angular sizes ...")
        distance = self.catalog["distance"].data * self.units["distance"]
        theta2 = (r_vir / distance).decompose().value  # [ rad ]
        size_major = theta2 / (1 - self.catalog["eccentricity"]**2) ** 0.25
        size_minor = theta2 * (1 - self.catalog["eccentricity"]**2) ** 0.25
        self.catalog["size_major"] = size_major * au.rad.to(au.deg)
        self.catalog["size_minor"] = size_minor * au.rad.to(au.deg)
        self.units["size"] = au.deg
        logger.info("Done calculate the elliptical angular sizes")

    def _calc_luminosity(self):
        """Calculate the radio luminosity (at 1.4 GHz) using empirical
        scaling relations.

        First, calculate the X-ray luminosity L_X using the empirical
        scaling relation between mass and X-ray luminosity.
        Then, derive the radio luminosity by employing the scaling
        relation between X-ray and radio luminosity.

        Attributes
        ----------
        catalog["luminosity"] : 1D `~numpy.ndarray`
            The luminosity density (at 1.4 GHz) of each cluster.
        catalog_prop["luminosity_freq"] : `~astropy.units.Quantity`
            The frequency (as an ``astropy`` quantity) where the above
            luminosity derived.
        units["luminosity"] : `~astropy.units.Unit`
            The unit used by the above luminosity.

        XXX/TODO
        --------
        The scaling relations used here may be outdated, and some of the
        quantities need trick conversions, which cause much confusion.

        Investigate for *more up-to-date scaling relations*, derived with
        new observation constraints.

        NOTE
        ----
        - The mass in the mass-X-ray luminosity scaling relation is NOT the
          cluster real mass, since [Reiprich2002]_ refer to the
          *critical density* ρ_c, while the scaling relation from
          Jenkins et al. (2001) requires mass refer to the
          *cosmic mean mass density ρ_m = Ω_m * ρ_c,
          therefore, the mass needs following conversion (which is an
          approximation):
              M_{R&B} ~= M * sqrt(OmegaM0)
        - The derived X-ray luminosity is for the 0.1-2.4 keV energy band.
        - The X-ray-radio luminosity scaling relation adopted here is
          derived at 1.4 GHz.
        - [EnBlin2002]_ assumes H0 = 50 h50 km/s/Mpc, so h50 = 1.

        References
        ----------
        - [Jelic2008], Eq.(13,14)
        - [EnBlin2002], Eq.(1,3)
        """
        # Dimensionless Hubble parameter adopted here and in the literature
        h_our = self.H0 / 100
        h_others = 50.0 / 100
        #
        logger.info("Calculating the radio luminosity (at 1.4 GHz) ...")
        # Calculate the X-ray luminosity from mass
        # XXX: why the mass conversion ??
        mass_RB = (self.catalog["mass"].data * self.units["mass"] *
                   self.catalog_prop["omega_m"]**0.5)
        a_X = 0.449
        b_X = 1.9
        # Hubble parameter conversion factor
        h_conv1 = (h_our / h_others) ** (b_X-2)
        # X-ray luminosity (0.1-2.4 keV) [ erg/s ]
        L_X = ((a_X * 1e45 *
                (mass_RB / (1e15*au.solMass)).decompose().value ** b_X) *
               h_conv1)
        # Calculate the radio luminosity from X-ray luminosity
        a_r = 2.78
        b_r = 1.94
        # Hubble parameter conversion factor
        h_conv2 = (h_our / h_others) ** (2*b_r-2)
        # Radio luminosity density (at 1.4 GHz) [ W/Hz ]
        L_r = (a_r * 1e24 * (L_X / 1e45)**b_r) * h_conv2
        self.catalog["luminosity"] = L_r
        self.catalog_prop["luminosity_freq"] = 1.4 * au.GHz
        self.units["luminosity"] = au.W / au.Hz
        logger.info("Done Calculate the radio luminosity")

    def _calc_specindex(self):
        """Calculate the radio spectral indexes for each cluster.

        Attributes
        ----------
        catalog["specindex"] : 1D `~numpy.ndarray`
            The radio spectral index of each cluster.

        XXX/TODO
        --------
        Currently, a common/uniform spectral index (1.2) is assumed for all
        clusters, which may be improved by investigating more recent results.
        """
        specindex = 1.2
        logger.info("Use common spectral index for all clusters: "
                    "{0}".format(specindex))
        self.catalog["specindex"] = specindex

    def _calc_Tb(self, luminosity, distance, specindex, frequency, size):
        """Calculate the brightness temperature at requested frequency
        by assuming a power-law spectral shape.

        Parameters
        ----------
        luminosity : float
            The luminosity density (unit: `self.units["luminosity"]`) at
            the reference frequency (i.e.,
            `self.catalog_prop["luminosity_freq"]`).
        distance : float
            The luminosity distance (unit: `self.units["distance"]`) to the
            object
        specindex : float
            The spectral index of the power-law spectrum.
            Note the *negative* sign in the formula.
        frequency : float
            The frequency (unit: `self.freq_unit`) where the brightness
            temperature requested.
        size : 2-float tuple
            The (major, minor) axes (unit: `self.units["size"]`).
            The order of major and minor can be arbitrary.

        Returns
        -------
        Tb : float
            Brightness temperature at the requested frequency, unit [ K ]

        NOTE
        ----
        The power-law spectral shape is assumed for *flux density* other
        than the *brightness temperature*.
        Therefore, the flux density at the requested frequency should first
        be calculated by extrapolating the spectrum, then convert the flux
        density to derive the brightness temperature.
        """
        freq = frequency * self.freq_unit
        freq_ref = self.catalog_prop["luminosity_freq"]
        luminosity = luminosity * self.units["luminosity"]
        Lnu = luminosity * float(freq / freq_ref) ** (-specindex)
        Fnu = Lnu / (4*np.pi * (distance*self.units["distance"])**2)
        omega = size[0]*self.units["size"] * size[1]*self.units["size"]
        Tb = Fnu_to_Tb(Fnu, omega, freq)
        return Tb.value

    def _simulate_templates(self):
        """Simulate the template (HEALPix) images for each cluster, and
        cache these templates within the class.

        The template images/maps have values of (or approximate) ones for
        these effective pixels, excluding the pixels corresponding to the
        edges of original rotated ellipse, which may have values of
        significantly less than 1 due to the rotation.

        Therefore, simulating the HEALPix map of one cluster at a requested
        frequency is simply multiplying the cached template image by the
        calculated brightness temperature (Tb) at that frequency.

        Furthermore, the total HEALPix map of all clusters are straightforward
        additions of all the maps of each cluster.

        Attributes
        ----------
        templates : list
            A List containing the simulated templates for each cluster.
            Each element is a `(hpidx, hpval)` tuple with `hpidx` the
            indexes of effective HEALPix pixels (RING ordering) and `hpval`
            the values of the corresponding pixels.
            e.g., ``[ (hpidx1, hpval1), (hpidx2, hpval2), ... ]``
        """
        logger.info("Simulating HEALPix templates for each cluster ...")
        templates = []
        resolution = self.resolution.to(au.deg).value
        # Make sure the index is reset, therefore, the *row indexes* can be
        # simply used to identify the corresponding template image.
        self.catalog.reset_index(inplace=True)
        # XXX/TODO: be parallel
        for row in self.catalog.itertuples():
            # TODO: progress bar
            center = (row.glon, row.glat)
            size = (row.size_major, row.size_minor)  # already [ degree ]
            rotation = row.rotation  # already [ degree ]
            grid = make_grid_ellipse(center, size, resolution, rotation)
            hpidx, hpval = map_grid_to_healpix(grid, self.nside)
            templates.append((hpidx, hpval))
        logger.info("Done simulate %d cluster templates" % len(templates))
        self.templates = templates

    def _simulate_single(self, data, frequency):
        """Simulate one single cluster at the specified frequency, based
        on the cached template image.

        Parameters
        ----------
        data : namedtuple
            The data of the SNR to be simulated, given in a ``namedtuple``
            object, from which can get the required properties by
            ``data.key``.
            e.g., elements of `self.catalog.itertuples()`
        frequency : float
            The simulation frequency (unit: `self.freq_unit`).

        Returns
        -------
        hpidx : 1D `~numpy.ndarray`
            The indexes (in RING ordering) of the effective HEALPix
            pixels for the SNR.
        hpval : 1D `~numpy.ndarray`
            The values (i.e., brightness temperature) of each HEALPix
            pixel with respect the above indexes.

        See Also
        --------
        `self._simulate_template()` for more detailed description.
        """
        index = data.Index
        hpidx, hpval = self.templates[index]
        # Calculate the brightness temperature
        luminosity = data.luminosity
        distance = data.distance
        specindex = data.specindex
        size = (data.size_major, data.size_minor)
        Tb = self._calc_Tb(luminosity, distance, specindex, frequency, size)
        hpval = hpval * Tb
        return (hpidx, hpval)

    def _make_filepath(self, **kwargs):
        """Make the path of output file according to the filename pattern
        and output directory loaded from configurations.
        """
        data = {
            "prefix": self.prefix,
        }
        data.update(kwargs)
        filename = self.filename_pattern.format(**data)
        filetype = self.configs.getn("output/filetype")
        if filetype == "fits":
            filename += ".fits"
        else:
            raise NotImplementedError("unsupported filetype: %s" % filetype)
        filepath = os.path.join(self.output_dir, filename)
        return filepath

    def _make_header(self):
        """Make the header with detail information (e.g., parameters and
        history) for the simulated products.
        """
        header = fits.Header()
        header["COMP"] = ("Extragalactic clusters of galaxies",
                          "Emission component")
        header["UNIT"] = ("Kelvin", "Map unit")
        header["CREATOR"] = (__name__, "File creator")
        # TODO:
        history = []
        comments = []
        for hist in history:
            header.add_history(hist)
        for cmt in comments:
            header.add_comment(cmt)
        self.header = header
        logger.info("Created FITS header")

    def output(self, hpmap, frequency):
        """Write the simulated free-free map to disk with proper header
        keywords and history.
        """
        if not os.path.exists(self.output_dir):
            os.mkdir(self.output_dir)
            logger.info("Created output dir: {0}".format(self.output_dir))
        #
        filepath = self._make_filepath(frequency=frequency)
        if not hasattr(self, "header"):
            self._make_header()
        header = self.header.copy()
        header["FREQ"] = (frequency, "Frequency [ MHz ]")
        header["DATE"] = (
            datetime.now(timezone.utc).astimezone().isoformat(),
            "File creation date"
        )
        if self.use_float:
            hpmap = hpmap.astype(np.float32)
        write_fits_healpix(filepath, hpmap, header=header,
                           clobber=self.clobber)
        logger.info("Write simulated map to file: {0}".format(filepath))

    def preprocess(self):
        """Perform the preparation procedures for the final simulations.

        Attributes
        ----------
        _preprocessed : bool
            This attribute presents and is ``True`` after the preparation
            procedures are performed, which indicates that it is ready to
            do the final simulations.
        """
        if hasattr(self, "_preprocessed") and self._preprocessed:
            return
        #
        logger.info("{name}: preprocessing ...".format(name=self.name))
        self._load_catalog()
        self._process_catalog()
        #
        self._expand_catalog_fullsky()
        self._add_random_position()
        self._add_random_eccentricity()
        self._add_random_rotation()
        #
        self._calc_sizes()
        self._calc_luminosity()
        self._calc_specindex()
        #
        self._simulate_templates()
        #
        self._preprocessed = True

    def simulate_frequency(self, frequency):
        """Simulate the emission (HEALPix) map of all Galactic SNRs at
        the specified frequency.

        Parameters
        ----------
        frequency : float
            The simulation frequency (unit: `self.freq_unit`).

        Returns
        -------
        hpmap_f : 1D `~numpy.ndarray`
            HEALPix map data in RING ordering

        See Also
        --------
        `self._simulate_template()` for more detailed description.
        """
        self.preprocess()
        #
        logger.info("Simulating {name} map at {freq} ({unit}) ...".format(
            name=self.name, freq=frequency, unit=self.freq_unit))
        hpmap_f = np.zeros(hp.nside2npix(self.nside))
        for row in self.catalog.itertuples():
            hpidx, hpval = self._simulate_single(row, frequency)
            hpmap_f[hpidx] += hpval
        #
        if self.save:
            self.output(hpmap_f, frequency)
        return hpmap_f

    def simulate(self, frequencies):
        """Simulate the emission (HEALPix) maps of all Galactic SNRs for
        every specified frequency.

        Parameters
        ----------
        frequency : list[float]
            List of frequencies (unit: `self.freq_unit`) where the
            simulation performed.

        Returns
        -------
        hpmaps : list[1D `~numpy.ndarray`]
            List of HEALPix maps (in RING ordering) at each frequency.
        """
        hpmaps = []
        for f in np.array(frequencies, ndmin=1):
            hpmap_f = self.simulate_frequency(f)
            hpmaps.append(hpmap_f)
        return hpmaps

    def postprocess(self):
        """Perform the post-simulation operations before the end."""
        logger.info("{name}: postprocessing ...".format(name=self.name))
        # Save the catalog actually used in the simulation
        self._save_catalog_inuse()