lm_space.py 4.75 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
# 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/>.
Theo Steininger's avatar
Theo Steininger committed
13
14
15
16
17
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
18

csongor's avatar
csongor committed
19
20
from __future__ import division
import numpy as np
Martin Reinecke's avatar
Martin Reinecke committed
21
from .space import Space
22
23
from .. import Field
from ..basic_arithmetics import exp
Martin Reinecke's avatar
Martin Reinecke committed
24
from ..dobj import from_ndarray as from_np
theos's avatar
theos committed
25

csongor's avatar
csongor committed
26

Theo Steininger's avatar
Theo Steininger committed
27
class LMSpace(Space):
Martin Reinecke's avatar
Martin Reinecke committed
28
29
30
31
32
33
34
35
    """NIFTY subclass for spherical harmonics components, for representations
    of fields on the two-sphere.

    Parameters
    ----------
    lmax : int
        The maximum :math:`l` value of any spherical harmonics
        :math:`Y_{lm}` that is represented in this Space.
36
37
38
39
40
41
42
        Must be >=0.

    mmax : int *optional*
        The maximum :math:`m` value of any spherical harmonics
        :math:`Y_{lm}` that is represented in this Space.
        If not supplied, it is set to lmax.
        Must be >=0 and <=lmax.
Martin Reinecke's avatar
Martin Reinecke committed
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

    See Also
    --------
    HPSpace : A class for the HEALPix discretization of the sphere [#]_.
    GLSpace : A class for the Gauss-Legendre discretization of the
        sphere [#]_.

    References
    ----------
    .. [#] K.M. Gorski et al., 2005, "HEALPix: A Framework for
           High-Resolution Discretization and Fast Analysis of Data
           Distributed on the Sphere", *ApJ* 622..759G.
    .. [#] M. Reinecke and D. Sverre Seljebotn, 2013, "Libsharp - spherical
           harmonic transforms revisited";
           `arXiv:1303.4945 <http://www.arxiv.org/abs/1303.4945>`_
csongor's avatar
csongor committed
58
59
    """

60
    def __init__(self, lmax, mmax=None):
Martin Reinecke's avatar
Martin Reinecke committed
61
        super(LMSpace, self).__init__()
62
63
64
65
66
67
68
69
70
        self._needed_for_hash += ["_lmax", "_mmax"]
        self._lmax = np.int(lmax)
        if self._lmax < 0:
            raise ValueError("lmax must be >=0.")
        if mmax is None:
            mmax = self._lmax
        self._mmax = np.int(mmax)
        if self._mmax < 0 or self._mmax > self._lmax:
            raise ValueError("mmax must be >=0 and <=lmax.")
csongor's avatar
csongor committed
71

72
    def __repr__(self):
73
        return ("LMSpace(lmax=%r, mmax=%r)" % (self.lmax, self.mmax))
74

csongor's avatar
csongor committed
75
76
77
    @property
    def harmonic(self):
        return True
csongor's avatar
csongor committed
78
79

    @property
80
81
    def shape(self):
        return (self.dim, )
csongor's avatar
csongor committed
82
83

    @property
84
    def dim(self):
85
86
        l = self._lmax
        m = self._mmax
87
        # the LMSpace consists of the full triangle (including -m's!),
theos's avatar
theos committed
88
        # minus two little triangles if mmax < lmax
89
        return (l+1)**2 - (l-m)*(l-m+1)
csongor's avatar
csongor committed
90

91
    def scalar_dvol(self):
92
        return 1.
csongor's avatar
csongor committed
93

94
    def get_k_length_array(self):
95
96
        lmax = self._lmax
        mmax = self._mmax
Martin Reinecke's avatar
stage1    
Martin Reinecke committed
97
98
99
100
101
102
        ldist = np.empty((self.dim,), dtype=np.float64)
        ldist[0:lmax+1] = np.arange(lmax+1, dtype=np.float64)
        tmp = np.empty((2*lmax+2), dtype=np.float64)
        tmp[0::2] = np.arange(lmax+1)
        tmp[1::2] = np.arange(lmax+1)
        idx = lmax+1
103
104
105
        for m in range(1, mmax+1):
            ldist[idx:idx+2*(lmax+1-m)] = tmp[2*m:]
            idx += 2*(lmax+1-m)
Martin Reinecke's avatar
Martin Reinecke committed
106
        return Field((self,), from_np(ldist))
107

108
    def get_unique_k_lengths(self):
Martin Reinecke's avatar
Martin Reinecke committed
109
110
        return np.arange(self.lmax+1, dtype=np.float64)

Martin Reinecke's avatar
Martin Reinecke committed
111
112
    @staticmethod
    def _kernel(x, sigma):
113
114
115
        # cf. "All-sky convolution for polarimetry experiments"
        # by Challinor et al.
        # http://arxiv.org/abs/astro-ph/0008228
Martin Reinecke's avatar
Martin Reinecke committed
116
117
118
        res = x+1.
        res *= x
        res *= -0.5*sigma*sigma
119
        exp(res, out=res)
Martin Reinecke's avatar
Martin Reinecke committed
120
121
        return res

122
    def get_fft_smoothing_kernel_function(self, sigma):
Martin Reinecke's avatar
Martin Reinecke committed
123
        return lambda x: self._kernel(x, sigma)
124

csongor's avatar
csongor committed
125
126
    @property
    def lmax(self):
Martin Reinecke's avatar
Martin Reinecke committed
127
        """ Returns the maximum :math:`l` value of any spherical harmonic
Theo Steininger's avatar
Theo Steininger committed
128
        :math:`Y_{lm}` that is represented in this Space.
129
        """
csongor's avatar
csongor committed
130
131
132
133
        return self._lmax

    @property
    def mmax(self):
Martin Reinecke's avatar
Martin Reinecke committed
134
135
        """ Returns the maximum :math:`m` value of any spherical harmonic
        :math:`Y_{lm}` that is represented in this Space.
136
        """
137
        return self._mmax
Martin Reinecke's avatar
Martin Reinecke committed
138
139
140
141
142
143
144
145

    def get_default_codomain(self):
        from .. import GLSpace
        return GLSpace(self.lmax+1, self.mmax*2+1)

    def check_codomain(self, codomain):
        from .. import GLSpace, HPSpace
        if not isinstance(codomain, (GLSpace, HPSpace)):
146
            raise TypeError("codomain must be a GLSpace or HPSpace.")