sharp_ylmgen.cc 6.24 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke 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
/*
 *  This file is part of libsharp2.
 *
 *  libsharp2 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 2 of the License, or
 *  (at your option) any later version.
 *
 *  libsharp2 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 libsharp2; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

/* libsharp2 is being developed at the Max-Planck-Institut fuer Astrophysik */

/*
 *  Helper code for efficient calculation of Y_lm(theta,phi=0)
 *
 *  Copyright (C) 2005-2019 Max-Planck-Society
 *  Author: Martin Reinecke
 */

Martin Reinecke's avatar
Martin Reinecke committed
28
29
30
#include <cmath>
#include <cstdlib>
#include <algorithm>
Martin Reinecke's avatar
Martin Reinecke committed
31
#include "libsharp2/sharp_internal.h"
Martin Reinecke's avatar
Martin Reinecke committed
32
33
#include "mr_util/error_handling.h"

Martin Reinecke's avatar
Martin Reinecke committed
34
35
using namespace std;

Martin Reinecke's avatar
Martin Reinecke committed
36
37
#pragma GCC visibility push(hidden)

Martin Reinecke's avatar
Martin Reinecke committed
38
static inline void normalize (double &val, int &scale, double xfmax)
Martin Reinecke's avatar
Martin Reinecke committed
39
  {
Martin Reinecke's avatar
Martin Reinecke committed
40
41
42
  while (abs(val)>xfmax) { val*=sharp_fsmall; ++scale; }
  if (val!=0.)
    while (abs(val)<xfmax*sharp_fsmall) { val*=sharp_fbig; --scale; }
Martin Reinecke's avatar
Martin Reinecke committed
43
44
  }

Martin Reinecke's avatar
Martin Reinecke committed
45
sharp_Ylmgen::sharp_Ylmgen (size_t l_max, size_t m_max, size_t spin)
Martin Reinecke's avatar
Martin Reinecke committed
46
  {
Martin Reinecke's avatar
Martin Reinecke committed
47
  constexpr double inv_sqrt4pi = 0.2820947917738781434740397257803862929220;
Martin Reinecke's avatar
Martin Reinecke committed
48

Martin Reinecke's avatar
Martin Reinecke committed
49
50
  lmax = l_max;
  mmax = m_max;
Martin Reinecke's avatar
Martin Reinecke committed
51
52
  MR_assert(l_max>=spin,"incorrect l_max: must be >= spin");
  MR_assert(l_max>=m_max,"incorrect l_max: must be >= m_max");
Martin Reinecke's avatar
Martin Reinecke committed
53
  s = spin;
Martin Reinecke's avatar
Martin Reinecke committed
54
55
  MR_assert((sharp_minscale<=0)&&(sharp_maxscale>0),
    "bad value for min/maxscale");
Martin Reinecke's avatar
Martin Reinecke committed
56
57
  cf.resize(sharp_maxscale-sharp_minscale+1);
  cf[-sharp_minscale]=1.;
Martin Reinecke's avatar
Martin Reinecke committed
58
  for (int sc=-sharp_minscale-1; sc>=0; --sc)
Martin Reinecke's avatar
Martin Reinecke committed
59
    cf[size_t(sc)]=cf[size_t(sc+1)]*sharp_fsmall;
Martin Reinecke's avatar
Martin Reinecke committed
60
  for (int sc=-sharp_minscale+1; sc<(sharp_maxscale-sharp_minscale+1); ++sc)
Martin Reinecke's avatar
Martin Reinecke committed
61
    cf[size_t(sc)]=cf[size_t(sc-1)]*sharp_fbig;
Martin Reinecke's avatar
Martin Reinecke committed
62
63
64
65
  powlimit.resize(m_max+spin+1);
  powlimit[0]=0.;
  constexpr double ln2 = 0.6931471805599453094172321214581766;
  constexpr double expo=-400*ln2;
Martin Reinecke's avatar
Martin Reinecke committed
66
67
  for (size_t i=1; i<=m_max+spin; ++i)
    powlimit[i]=exp(expo/i);
Martin Reinecke's avatar
Martin Reinecke committed
68

Martin Reinecke's avatar
Martin Reinecke committed
69
  m = ~size_t(0);
Martin Reinecke's avatar
Martin Reinecke committed
70
71
  if (spin==0)
    {
Martin Reinecke's avatar
Martin Reinecke committed
72
73
    mfac.resize(mmax+1);
    mfac[0] = inv_sqrt4pi;
Martin Reinecke's avatar
Martin Reinecke committed
74
75
    for (size_t i=1; i<=mmax; ++i)
      mfac[i] = mfac[i-1]*sqrt((2*i+1.)/(2*i));
Martin Reinecke's avatar
Martin Reinecke committed
76
77
    root.resize(2*lmax+8);
    iroot.resize(2*lmax+8);
Martin Reinecke's avatar
Martin Reinecke committed
78
    for (size_t i=0; i<2*lmax+8; ++i)
Martin Reinecke's avatar
Martin Reinecke committed
79
      {
Martin Reinecke's avatar
Martin Reinecke committed
80
81
      root[i] = sqrt(i);
      iroot[i] = (i==0) ? 0. : 1./root[i];
Martin Reinecke's avatar
Martin Reinecke committed
82
      }
Martin Reinecke's avatar
Martin Reinecke committed
83
84
85
    eps.resize(lmax+4);
    alpha.resize(lmax/2+2);
    coef.resize(lmax/2+2);
Martin Reinecke's avatar
Martin Reinecke committed
86
87
88
    }
  else
    {
Martin Reinecke's avatar
Martin Reinecke committed
89
    m=mlo=mhi=~size_t(0);
Martin Reinecke's avatar
Martin Reinecke committed
90
    coef.resize(lmax+3);
Martin Reinecke's avatar
Martin Reinecke committed
91
92
    for (size_t i=0; i<lmax+3; ++i)
      coef[i].a=coef[i].b=0.;
Martin Reinecke's avatar
Martin Reinecke committed
93
94
95
    alpha.resize(lmax+3);
    inv.resize(lmax+2);
    inv[0]=0;
Martin Reinecke's avatar
Martin Reinecke committed
96
    for (size_t i=1; i<lmax+2; ++i) inv[i]=1./i;
Martin Reinecke's avatar
Martin Reinecke committed
97
98
    flm1.resize(2*lmax+3);
    flm2.resize(2*lmax+3);
Martin Reinecke's avatar
Martin Reinecke committed
99
    for (size_t i=0; i<2*lmax+3; ++i)
Martin Reinecke's avatar
Martin Reinecke committed
100
      {
Martin Reinecke's avatar
Martin Reinecke committed
101
102
      flm1[i] = sqrt(1./(i+1.));
      flm2[i] = sqrt(i/(i+1.));
Martin Reinecke's avatar
Martin Reinecke committed
103
      }
Martin Reinecke's avatar
Martin Reinecke committed
104
105
106
107
    prefac.resize(mmax+1);
    fscale.resize(mmax+1);
    vector<double> fac(2*lmax+1);
    vector<int> facscale(2*lmax+1);
Martin Reinecke's avatar
Martin Reinecke committed
108
    fac[0]=1; facscale[0]=0;
Martin Reinecke's avatar
Martin Reinecke committed
109
    for (size_t i=1; i<2*lmax+1; ++i)
Martin Reinecke's avatar
Martin Reinecke committed
110
      {
Martin Reinecke's avatar
Martin Reinecke committed
111
112
113
      fac[i]=fac[i-1]*sqrt(i);
      facscale[i]=facscale[i-1];
      normalize(fac[i],facscale[i],sharp_fbighalf);
Martin Reinecke's avatar
Martin Reinecke committed
114
      }
Martin Reinecke's avatar
Martin Reinecke committed
115
    for (size_t i=0; i<=mmax; ++i)
Martin Reinecke's avatar
Martin Reinecke committed
116
      {
Martin Reinecke's avatar
Martin Reinecke committed
117
      size_t mlo_=s, mhi_=i;
Martin Reinecke's avatar
Martin Reinecke committed
118
119
120
      if (mhi_<mlo_) swap(mhi_,mlo_);
      double tfac=fac[2*mhi_]/fac[mhi_+mlo_];
      int tscale=facscale[2*mhi_]-facscale[mhi_+mlo_];
Martin Reinecke's avatar
Martin Reinecke committed
121
      normalize(tfac,tscale,sharp_fbighalf);
Martin Reinecke's avatar
Martin Reinecke committed
122
123
      tfac/=fac[mhi_-mlo_];
      tscale-=facscale[mhi_-mlo_];
Martin Reinecke's avatar
Martin Reinecke committed
124
      normalize(tfac,tscale,sharp_fbighalf);
Martin Reinecke's avatar
Martin Reinecke committed
125
126
      prefac[i]=tfac;
      fscale[i]=tscale;
Martin Reinecke's avatar
Martin Reinecke committed
127
128
129
130
      }
    }
  }

Martin Reinecke's avatar
Martin Reinecke committed
131
void sharp_Ylmgen::prepare (size_t m_)
Martin Reinecke's avatar
Martin Reinecke committed
132
  {
Martin Reinecke's avatar
Martin Reinecke committed
133
134
  if (m_==m) return;
  m = m_;
Martin Reinecke's avatar
Martin Reinecke committed
135

Martin Reinecke's avatar
Martin Reinecke committed
136
  if (s==0)
Martin Reinecke's avatar
Martin Reinecke committed
137
    {
Martin Reinecke's avatar
Martin Reinecke committed
138
    eps[m] = 0.;
Martin Reinecke's avatar
Martin Reinecke committed
139
    for (size_t l=m+1; l<lmax+4; ++l)
Martin Reinecke's avatar
Martin Reinecke committed
140
141
142
      eps[l] = root[l+m]*root[l-m]*iroot[2*l+1]*iroot[2*l-1];
    alpha[0] = 1./eps[m+1];
    alpha[1] = eps[m+1]/(eps[m+2]*eps[m+3]);
Martin Reinecke's avatar
Martin Reinecke committed
143
    for (size_t il=1, l=m+2; l<lmax+1; ++il, l+=2)
Martin Reinecke's avatar
Martin Reinecke committed
144
      alpha[il+1]= ((il&1) ? -1 : 1) / (eps[l+2]*eps[l+3]*alpha[il]);
Martin Reinecke's avatar
Martin Reinecke committed
145
    for (size_t il=0, l=m; l<lmax+2; ++il, l+=2)
Martin Reinecke's avatar
Martin Reinecke committed
146
      {
Martin Reinecke's avatar
Martin Reinecke committed
147
148
149
      coef[il].a = ((il&1) ? -1 : 1)*alpha[il]*alpha[il];
      double t1 = eps[l+2], t2 = eps[l+1];
      coef[il].b = -coef[il].a*(t1*t1+t2*t2);
Martin Reinecke's avatar
Martin Reinecke committed
150
151
152
153
      }
    }
  else
    {
Martin Reinecke's avatar
Martin Reinecke committed
154
    size_t mlo_=m, mhi_=s;
Martin Reinecke's avatar
Martin Reinecke committed
155
    if (mhi_<mlo_) swap(mhi_,mlo_);
Martin Reinecke's avatar
Martin Reinecke committed
156
    bool ms_similar = ((mhi==mhi_) && (mlo==mlo_));
Martin Reinecke's avatar
Martin Reinecke committed
157

Martin Reinecke's avatar
Martin Reinecke committed
158
    mlo = mlo_; mhi = mhi_;
Martin Reinecke's avatar
Martin Reinecke committed
159
160
161

    if (!ms_similar)
      {
Martin Reinecke's avatar
Martin Reinecke committed
162
163
      alpha[mhi] = 1.;
      coef[mhi].a = coef[mhi].b = 0.;
Martin Reinecke's avatar
Martin Reinecke committed
164
      for (size_t l=mhi; l<=lmax; ++l)
Martin Reinecke's avatar
Martin Reinecke committed
165
        {
Martin Reinecke's avatar
Martin Reinecke committed
166
        double t = flm1[l+m]*flm1[l-m]*flm1[l+s]*flm1[l-s];
Martin Reinecke's avatar
Martin Reinecke committed
167
168
169
        double lt = 2*l+1;
        double l1 = l+1;
        double flp10=l1*lt*t;
Martin Reinecke's avatar
Martin Reinecke committed
170
171
172
173
174
        double flp11=m*s*inv[l]*inv[l+1];
        t = flm2[l+m]*flm2[l-m]*flm2[l+s]*flm2[l-s];
        double flp12=t*l1*inv[l];
        if (l>mhi)
          alpha[l+1] = alpha[l-1]*flp12;
Martin Reinecke's avatar
Martin Reinecke committed
175
        else
Martin Reinecke's avatar
Martin Reinecke committed
176
177
178
          alpha[l+1] = 1.;
        coef[l+1].a = flp10*alpha[l]/alpha[l+1];
        coef[l+1].b = flp11*coef[l+1].a;
Martin Reinecke's avatar
Martin Reinecke committed
179
180
181
        }
      }

Martin Reinecke's avatar
Martin Reinecke committed
182
    preMinus_p = preMinus_m = false;
Martin Reinecke's avatar
Martin Reinecke committed
183
    if (mhi==m)
Martin Reinecke's avatar
Martin Reinecke committed
184
      {
Martin Reinecke's avatar
Martin Reinecke committed
185
186
      cosPow = mhi+s; sinPow = mhi-s;
      preMinus_p = preMinus_m = ((mhi-s)&1);
Martin Reinecke's avatar
Martin Reinecke committed
187
188
189
      }
    else
      {
Martin Reinecke's avatar
Martin Reinecke committed
190
191
      cosPow = mhi+m; sinPow = mhi-m;
      preMinus_m = ((mhi+m)&1);
Martin Reinecke's avatar
Martin Reinecke committed
192
193
194
195
      }
    }
  }

Martin Reinecke's avatar
Martin Reinecke committed
196
vector<double> sharp_Ylmgen::get_norm (size_t lmax, size_t spin)
Martin Reinecke's avatar
Martin Reinecke committed
197
198
199
200
201
202
203
204
205
206
  {
  const double pi = 3.141592653589793238462643383279502884197;
  /* sign convention for H=1 (LensPix paper) */
#if 1
   double spinsign = (spin>0) ? -1.0 : 1.0;
#else
   double spinsign = 1.0;
#endif

  if (spin==0)
Martin Reinecke's avatar
Martin Reinecke committed
207
    return vector<double>(lmax+1,1.);
Martin Reinecke's avatar
Martin Reinecke committed
208

Martin Reinecke's avatar
Martin Reinecke committed
209
  vector<double> res(lmax+1);
Martin Reinecke's avatar
Martin Reinecke committed
210
  spinsign = (spin&1) ? -spinsign : spinsign;
Martin Reinecke's avatar
Martin Reinecke committed
211
  for (size_t l=0; l<=lmax; ++l)
Martin Reinecke's avatar
Martin Reinecke committed
212
213
214
215
    res[l] = (l<spin) ? 0. : spinsign*0.5*sqrt((2*l+1)/(4*pi));
  return res;
  }

Martin Reinecke's avatar
Martin Reinecke committed
216
vector<double> sharp_Ylmgen::get_d1norm (size_t lmax)
Martin Reinecke's avatar
Martin Reinecke committed
217
218
  {
  const double pi = 3.141592653589793238462643383279502884197;
Martin Reinecke's avatar
Martin Reinecke committed
219
  vector<double> res(lmax+1);
Martin Reinecke's avatar
Martin Reinecke committed
220

Martin Reinecke's avatar
Martin Reinecke committed
221
  for (size_t l=0; l<=lmax; ++l)
Martin Reinecke's avatar
Martin Reinecke committed
222
223
224
225
226
    res[l] = (l<1) ? 0. : 0.5*sqrt(l*(l+1.)*(2*l+1.)/(4*pi));
  return res;
  }

#pragma GCC visibility pop