spin_functions.py 6.04 KB
Newer Older
1
"""
Syl's avatar
Syl committed
2 3
Set of routines to compute the pixel covariance matrix using
Legendre polynomials
4 5 6 7 8 9 10 11 12 13

Author: Vanneste
"""
from __future__ import division

import math

import numpy as np
from scipy import special

Syl's avatar
Syl committed
14

15 16
def dlss(z, s1, s2, lmax):
    """
Syl's avatar
Syl committed
17
    Computes the reduced Wigner D-function d^l_ss'
18 19 20

    Parameters
    ----------
Syl's avatar
Syl committed
21 22
    z : float
        Cosine of the angle between two pixels
23
    s1 : int
Syl's avatar
Syl committed
24
        Spin number 1
25
    s2 : int
Syl's avatar
Syl committed
26
        Spin number 2
27 28 29 30 31 32 33 34
    lmax : int
        Maximum multipole

    Returns
    ----------
    d : 1D array of floats
        ???

Syl's avatar
Syl committed
35 36 37 38 39
    Example
    ----------
    >>> d = dlss(0.1, 2,  2, 5)
    >>> print(round(sum(d),5))
    0.24351
40 41 42 43 44 45
    """
    d = np.zeros((lmax + 1))
    if s1 < abs(s2):
        print("error spins, s1<|s2|")
        return

Julien's avatar
Julien committed
46
    # Conv: sign = -1 if (s1 + s2) and 1 else 1
47 48 49 50
    sign = (-1)**(s1 - s2)
    fs1 = math.factorial(2.0 * s1)
    fs1ps2 = math.factorial(1.0 * s1 + s2)
    fs1ms2 = math.factorial(1.0 * s1 - s2)
Julien's avatar
Julien committed
51
    num = (1.0 + z)**(0.5 * (s1 + s2)) * (1.0 - z)**(0.5 * (s1 - s2))
52 53 54 55 56 57 58 59 60

    # Initialise the recursion (l = s1 + 1)
    d[s1] = sign / 2.0**s1 * np.sqrt(fs1 / fs1ps2 / fs1ms2) * num

    l1 = s1 + 1.0
    rhoSSL1 = np.sqrt((l1 * l1 - s1 * s1) * (l1 * l1 - s2 * s2)) / l1
    d[s1+1] = (2 * s1 + 1.0)*(z - s2 / (s1 + 1.0)) * d[s1] / rhoSSL1

    # Build the recursion for l > s1 + 1
Julien's avatar
Julien committed
61
    for l in np.arange(s1 + 1, lmax, 1):
62 63 64 65 66 67 68 69 70 71 72 73 74
        l1 = l + 1.0
        numSSL = (l * l * 1.0 - s1 * s1) * (l * l * 1.0 - s2 * s2)
        rhoSSL = np.sqrt(numSSL) / (l * 1.0)
        numSSL1 = (l1 * l1 - s1 * s1) * (l1 * l1 - s2 * s2)
        rhoSSL1 = np.sqrt(numSSL1) / l1

        numd = (2.0 * l + 1.0) * (z - s1 * s2 / (l * 1.0) / l1) * d[l]
        d[l+1] = (numd - rhoSSL * d[l-1]) / rhoSSL1
    return d


def pl0(z, lmax):
    """
Syl's avatar
Syl committed
75 76
    Computes the associated Legendre function of the first kind of order 0
    Pn(z) from 0 to lmax (inclusive).
77 78 79

    Parameters
    ----------
Syl's avatar
Syl committed
80 81
    z : float
        Cosine of the angle between two pixels
82 83 84 85 86
    lmax : int
        Maximum multipole

    Returns
    ----------
Syl's avatar
Syl committed
87 88
    Pn : 1D array of floats
    Legendre function
89

Syl's avatar
Syl committed
90 91 92 93 94
    Example
    ----------
    >>> thepl0 = pl0(0.1, 5)
    >>> print(round(sum(thepl0),5))
    0.98427
95
    """
Syl's avatar
Syl committed
96 97 98
    Pn = special.lpn(lmax, z)[0]
    return Pn

99 100 101

def pl2(z, lmax):
    """
Syl's avatar
Syl committed
102 103
    Computes the associated Legendre function of the first kind of order 2
    from 0 to lmax (inclusive)
104 105 106

    Parameters
    ----------
Syl's avatar
Syl committed
107 108
    z : float
        Cosine of the angle between two pixels
109 110 111 112 113
    lmax : int
        Maximum multipole

    Returns
    ----------
Syl's avatar
Syl committed
114 115 116 117 118 119 120 121
    Pn2 : 1D array of floats
    Legendre function

    Example
    ----------
    >>> thepl2 = pl2(0.1, 5)
    >>> print(round(sum(thepl2),5))
    -7.49183
122
    """
Syl's avatar
Syl committed
123 124
    Pn2 = special.lpmn(2, lmax, z)[0][2]
    return Pn2
125 126


Syl's avatar
Syl committed
127
# ####### F1 and F2 functions from Tegmark & De Oliveira-Costa, 2000  #########
128 129
def F1l0(z, lmax):
    """
Syl's avatar
Syl committed
130
    Compute the F1l0 function from Tegmark & De Oliveira-Costa, 2000
131 132 133

    Parameters
    ----------
Syl's avatar
Syl committed
134 135
    z : float
        Cosine of the angle between two pixels
136 137 138 139 140
    lmax : int
        Maximum multipole

    Returns
    ----------
Syl's avatar
Syl committed
141 142 143 144 145 146 147 148
    bla : 1D array of float
    F1l0 function from Tegmark & De Oliveira-Costa, 2000

    Example
    ----------
    >>> theF1l0= F1l0(0.1, 5)
    >>> print(round(sum(theF1l0),5))
    0.20392
149 150 151 152
    """
    if abs(z) == 1.0:
        return(np.zeros(lmax + 1))
    else:
Syl's avatar
Syl committed
153
        ell = np.arange(2, lmax + 1)
154
        thepl = pl0(z, lmax)
Syl's avatar
Syl committed
155 156 157
        theplm1 = np.append(0, thepl[:-1])
        thepl = thepl[2:]
        theplm1 = theplm1[2:]
158 159 160
        a0 = 2.0 / np.sqrt((ell - 1) * ell * (ell + 1) * (ell + 2))
        a1 = ell * z * theplm1 / (1 - z**2)
        a2 = (ell / (1 - z**2) + ell * (ell - 1) / 2) * thepl
Syl's avatar
Syl committed
161
        bla = np.append([0, 0], a0 * (a1 - a2))
162 163
        return bla

Syl's avatar
Syl committed
164

165 166
def F1l2(z, lmax):
    """
Syl's avatar
Syl committed
167
    Compute the F1l2 function from Tegmark & De Oliveira-Costa, 2000
168 169 170

    Parameters
    ----------
Syl's avatar
Syl committed
171 172
    z : float
        Cosine of the angle between two pixels
173 174 175 176 177
    lmax : int
        Maximum multipole

    Returns
    ----------
Syl's avatar
Syl committed
178 179 180 181 182 183 184 185
    bla : 1D array of float
    F1l2 function from Tegmark & De Oliveira-Costa, 2000

    Example
    ----------
    >>> theF1l2= F1l2(0.1, 5)
    >>> print(round(sum(theF1l2),5))
    0.58396
186 187 188 189 190 191 192
    """
    if z == 1.0:
        return np.append(np.zeros(2), np.ones(lmax - 1) * 0.5)
    elif z == -1.0:
        ell = np.arange(lmax + 1)
        return np.append(np.zeros(2), 0.5 * (-1)**ell[2:])
    else:
Syl's avatar
Syl committed
193
        ell = np.arange(2, lmax + 1)
194
        thepl2 = pl2(z, lmax)
Syl's avatar
Syl committed
195 196 197
        theplm1_2 = np.append(0, thepl2[:-1])
        thepl2 = thepl2[2:]
        theplm1_2 = theplm1_2[2:]
198 199 200
        a0 = 2.0 / ((ell - 1) * ell * (ell + 1) * (ell + 2))
        a1 = (ell + 2) * z * theplm1_2 / (1 - z**2)
        a2 = ((ell - 4) / (1 - z**2) + ell * (ell - 1) / 2) * thepl2
Syl's avatar
Syl committed
201
        bla = np.append([0, 0], a0 * (a1 - a2))
202 203
        return bla

Syl's avatar
Syl committed
204 205

def F2l2(z, lmax):
206
    """
Syl's avatar
Syl committed
207
    Compute the F2l2 function from Tegmark & De Oliveira-Costa, 2000
208 209 210

    Parameters
    ----------
Syl's avatar
Syl committed
211 212
    z : float
        Cosine of the angle between two pixels
213 214 215 216 217 218
    lmax : int
        Maximum multipole

    Returns
    ----------
    ???
Syl's avatar
Syl committed
219 220 221 222 223 224

    Example
    ----------
    >>> theF2l2= F2l2(0.1, 5)
    >>> print(round(sum(theF2l2),5))
    0.34045
225 226 227 228 229 230 231
    """
    if z == 1.0:
        return np.append(np.zeros(2), -0.5 * np.ones(lmax - 1))
    elif z == -1.0:
        ell = np.arange(lmax + 1)
        return np.append(np.zeros(2), 0.5 * (-1)**ell[2:])
    else:
Syl's avatar
Syl committed
232
        ell = np.arange(2, lmax + 1)
233
        thepl2 = pl2(z, lmax)
Syl's avatar
Syl committed
234 235 236
        theplm1_2 = np.append(0, thepl2[:-1])
        thepl2 = thepl2[2:]
        theplm1_2 = theplm1_2[2:]
237 238 239
        a0 = 4.0 / ((ell - 1) * ell * (ell + 1) * (ell + 2) * (1 - z**2))
        a1 = (ell + 2) * theplm1_2
        a2 = (ell - 1) * z * thepl2
Syl's avatar
Syl committed
240
        bla = np.append([0, 0], a0 * (a1 - a2))
241
        return bla
Syl's avatar
Syl committed
242 243 244 245 246 247 248 249 250 251 252 253 254 255

if __name__ == "__main__":
    """
    Run the doctest using

    python simulation.py

    If the tests are OK, the script should exit gracefuly, otherwise the
    failure(s) will be printed out.
    """
    import doctest
    if np.__version__ >= "1.14.0":
        np.set_printoptions(legacy="1.13")
    doctest.testmod()