There is a maintenance of MPCDF Gitlab on Thursday, April 22st 2020, 9:00 am CEST - Expect some service interruptions during this time

demo_faraday.py 4.65 KB
Newer Older
Marco Selig's avatar
Marco Selig committed
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
## NIFTY (Numerical Information Field Theory) has been developed at the
## Max-Planck-Institute for Astrophysics.
##
## Copyright (C) 2013 Max-Planck-Society
##
## Author: Marco Selig
## Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
## See the GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.

"""
    ..                  __   ____   __
    ..                /__/ /   _/ /  /_
    ..      __ ___    __  /  /_  /   _/  __   __
    ..    /   _   | /  / /   _/ /  /   /  / /  /
    ..   /  / /  / /  / /  /   /  /_  /  /_/  /
    ..  /__/ /__/ /__/ /__/    \___/  \___   /  demo
    ..                               /______/

    NIFTY demo applying transformations to the "Faraday Map" [#]_.

    References
    ----------
    .. [#] N. Opermann et. al., "An improved map of the Galactic Faraday sky",
        Astronomy & Astrophysics, Volume 542, id.A93, 06/2012;
        `arXiv:1111.6186 <http://www.arxiv.org/abs/1111.6186>`_

"""
from __future__ import division
from nifty import *



##-----------------------------------------------------------------------------

Marco Selig's avatar
Marco Selig committed
47
# (global) Faraday map
48
m = field(HPSpace(128), val=np.load(os.path.join(get_demo_dir(),
csongor's avatar
csongor committed
49
                                           "demo_faraday_map.npy")))
Marco Selig's avatar
Marco Selig committed
50 51 52 53 54 55 56 57 58 59 60 61

##-----------------------------------------------------------------------------

##=============================================================================

def run(projection=False, power=False):
    """
        Runs the demo.

        Parameters
        ----------
        projection : bool, *optional*
Marco Selig's avatar
Marco Selig committed
62
            Whether to additionaly show projections or not (default: False).
Marco Selig's avatar
Marco Selig committed
63
        power : bool, *optional*
Marco Selig's avatar
Marco Selig committed
64
            Whether to additionaly show power spectra or not (default: False).
Marco Selig's avatar
Marco Selig committed
65 66

    """
Marco Selig's avatar
Marco Selig committed
67 68 69
    nicely = {"vmin":-4, "vmax":4, "cmap":ncmap.fm()}

    # (a) representation on HEALPix grid
Marco Selig's avatar
Marco Selig committed
70
    m0 = m
Marco Selig's avatar
Marco Selig committed
71 72
    m0.plot(title=r"$m$ on a HEALPix grid", **nicely)
    nicely.update({"cmap":ncmap.fm()}) # healpy bug workaround
Marco Selig's avatar
Marco Selig committed
73

Marco Selig's avatar
Marco Selig committed
74 75 76 77 78 79 80
    # (b) representation in spherical harmonics
    k_space = m0.target # == lm_space(383, mmax=383)
    m1 = m0.transform(k_space) # == m.transform()
#    m1.plot(title=r"$m$ in spherical harmonics")

    if(power):
        m1.plot(title=r"angular power spectrum of $m$", vmin=1E-2, vmax=1E+1, mono=False)
Marco Selig's avatar
Marco Selig committed
81
    if(projection):
Marco Selig's avatar
Marco Selig committed
82 83 84
        # define projection operator
        Sk = projection_operator(m1.domain)
        # project quadrupole
Marco Selig's avatar
Marco Selig committed
85
        m2 = Sk(m0, band=2)
Marco Selig's avatar
Marco Selig committed
86
        m2.plot(title=r"angular quadrupole of $m$ on a HEALPix grid", **nicely)
Marco Selig's avatar
Marco Selig committed
87

Marco Selig's avatar
Marco Selig committed
88 89 90 91
    # (c) representation on Gauss-Legendre grid
    y_space = m.target.get_codomain(coname="gl") # == gl_space(384, nlon=767)
    m3 = m1.transform(y_space) # == m0.transform().transform(y_space)
    m3.plot(title=r"$m$ on a spherical Gauss-Legendre grid", **nicely)
Marco Selig's avatar
Marco Selig committed
92

Marco Selig's avatar
Marco Selig committed
93 94 95 96 97
    if(projection):
        m4 = Sk(m3, band=2)
        m4.plot(title=r"angular quadrupole of $m$ on a Gauss-Legendre grid", **nicely)

    # (d) representation on regular grid
98 99
    y_space = GLSpace(384, nlon=768) # auxiliary gl_space
    z_space = RGSpace([768, 384], dist=[1/360, 1/180])
Marco Selig's avatar
Marco Selig committed
100 101
    m5 = m1.transform(y_space)
    m5.cast_domain(z_space)
102
    m5.set_val(np.roll(m5.val[::-1, ::-1], y_space.paradict['nlon']//2, axis=1)) # rearrange value array
Marco Selig's avatar
Marco Selig committed
103
    m5.plot(title=r"$m$ on a regular 2D grid", **nicely)
Marco Selig's avatar
Marco Selig committed
104 105

    if(power):
Marco Selig's avatar
Marco Selig committed
106 107
        m5.target.set_power_indices(log=False)
        m5.plot(power=True, title=r"Fourier power spectrum of $m$", vmin=1E-3, vmax=1E+0, mono=False)
Marco Selig's avatar
Marco Selig committed
108
    if(projection):
Marco Selig's avatar
Marco Selig committed
109 110 111 112 113 114
        m5.target.set_power_indices(log=False)
        # define projection operator
        Sk = projection_operator(m5.target)
        # project quadrupole
        m6 = Sk(m5, band=2)
        m6.plot(title=r"Fourier quadrupole of $m$ on a regular 2D grid", **nicely)
Marco Selig's avatar
Marco Selig committed
115 116 117

##=============================================================================

Marco Selig's avatar
Marco Selig committed
118 119 120 121 122 123 124 125 126 127 128 129
##-----------------------------------------------------------------------------

if(__name__=="__main__"):
#    pl.close("all")

    # run demo
    run(projection=False, power=False)
    # define projection operator
    Sk = projection_operator(m.target)

##-----------------------------------------------------------------------------