demo_faraday.py 3.96 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 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
## 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 *
from nifty.nifty_cmaps import *


about.warnings.off()
about.infos.off()


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

## spaces
h  = hp_space(128)
l  = lm_space(383)
g  = gl_space(384) ## nlon == 767
g_ = gl_space(384, nlon=768)
r  = rg_space([768, 384], dist=[1/360, 1/180])
r_ = rg_space([256, 128], dist=[1/120, 1/60])

## map
m = field(h, val=np.load("demo_faraday_map.npy"))

## projection operator
Sk = None

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



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

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

        Parameters
        ----------
        projection : bool, *optional*
            Whether to additionaly include projections in the demo or not. If
            ``projection == True`` the projection operator `Sk` will be
            defined. (default: False)
        power : bool, *optional*
            Whether to additionaly show power spectra in the demo or not.
            (default: False)

    """
    global Sk
    ## start in hp_space
    m0 = m

    ## transform to lm_space
    m1 = m0.transform(l)
    if(projection):
        ## define projection operator
        powerindex = l.get_power_index()
        Sk = projection_operator(l, assign=powerindex)
        ## project quadrupole
        m2 = Sk(m0, band=2)

    ## transform to gl_space
    m3 = m1.transform(g)

    ## transform to rg_space
    m4 = m1.transform(g_) ## auxiliary gl_space
    m4.cast_domain(r) ## rg_space cast
    m4.set_val(np.roll(m4.val[::-1, ::-1], g.nlon()//2, axis=1)) ## rearrange
    if(power):
        ## restrict to central window
        m5 = field(r_, val=m4[128:256, 256:512]).transform()

    ## plots
    m0.plot(title=r"$m$ on a HEALPix grid", vmin=-4, vmax=4, cmap=ncmap.fm())
    if(power):
        m1.plot(title=r"angular power spectrum of $m$", vmin=1E-2, vmax=1E+1, mono=False)
    if(projection):
        m2.plot(title=r"quadrupole of $m$ on a HEALPix grid", vmin=-4, vmax=4, cmap=ncmap.fm())
    m3.plot(title=r"$m$ on a spherical Gauss-Legendre grid", vmin=-4, vmax=4, cmap=ncmap.fm())
    m4.plot(title=r"$m$ on a regular 2D grid", vmin=-4, vmax=4, cmap=ncmap.fm())
    if(power):
        m5.plot(title=r"(restricted) Fourier power spectrum of $m$", vmin=1E-3, vmax=1E+0, mono=False)

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