sharp_ylmgen.cc 6.28 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
cleanup    
Martin Reinecke committed
31
#include "libsharp2/sharp_ylmgen.h"
Martin Reinecke's avatar
Martin Reinecke committed
32
33
34
#include "libsharp2/sharp_utils.h"
#include "mr_util/error_handling.h"

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

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

Martin Reinecke's avatar
Martin Reinecke committed
39
static inline void normalize (double &val, int &scale, double xfmax)
Martin Reinecke's avatar
Martin Reinecke committed
40
  {
Martin Reinecke's avatar
Martin Reinecke committed
41
42
43
  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
44
45
  }

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

Martin Reinecke's avatar
Martin Reinecke committed
50
51
  lmax = l_max;
  mmax = m_max;
Martin Reinecke's avatar
Martin Reinecke committed
52
53
54
  MR_assert(spin>=0,"incorrect spin: must be nonnegative");
  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
55
  s = spin;
Martin Reinecke's avatar
Martin Reinecke committed
56
57
  MR_assert((sharp_minscale<=0)&&(sharp_maxscale>0),
    "bad value for min/maxscale");
Martin Reinecke's avatar
Martin Reinecke committed
58
59
  cf.resize(sharp_maxscale-sharp_minscale+1);
  cf[-sharp_minscale]=1.;
Martin Reinecke's avatar
Martin Reinecke committed
60
  for (int m=-sharp_minscale-1; m>=0; --m)
Martin Reinecke's avatar
Martin Reinecke committed
61
    cf[m]=cf[m+1]*sharp_fsmall;
Martin Reinecke's avatar
Martin Reinecke committed
62
  for (int m=-sharp_minscale+1; m<(sharp_maxscale-sharp_minscale+1); ++m)
Martin Reinecke's avatar
Martin Reinecke committed
63
64
65
66
67
    cf[m]=cf[m-1]*sharp_fbig;
  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
68
  for (int m=1; m<=m_max+spin; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
69
    powlimit[m]=exp(expo/m);
Martin Reinecke's avatar
Martin Reinecke committed
70

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

Martin Reinecke's avatar
Martin Reinecke committed
133
void sharp_Ylmgen::prepare (int m)
Martin Reinecke's avatar
Martin Reinecke committed
134
  {
Martin Reinecke's avatar
Martin Reinecke committed
135
  if (m==this->m) return;
Martin Reinecke's avatar
Martin Reinecke committed
136
  MR_assert(m>=0,"incorrect m");
Martin Reinecke's avatar
Martin Reinecke committed
137
  this->m = m;
Martin Reinecke's avatar
Martin Reinecke committed
138

Martin Reinecke's avatar
Martin Reinecke committed
139
  if (s==0)
Martin Reinecke's avatar
Martin Reinecke committed
140
    {
Martin Reinecke's avatar
Martin Reinecke committed
141
142
143
144
145
146
147
148
    eps[m] = 0.;
    for (int l=m+1; l<lmax+4; ++l)
      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]);
    for (int il=1, l=m+2; l<lmax+1; ++il, l+=2)
      alpha[il+1]= ((il&1) ? -1 : 1) / (eps[l+2]*eps[l+3]*alpha[il]);
    for (int il=0, l=m; l<lmax+2; ++il, l+=2)
Martin Reinecke's avatar
Martin Reinecke committed
149
      {
Martin Reinecke's avatar
Martin Reinecke committed
150
151
152
      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
153
154
155
156
      }
    }
  else
    {
Martin Reinecke's avatar
Martin Reinecke committed
157
158
159
    int mlo_=m, mhi_=s;
    if (mhi_<mlo_) swap(mhi_,mlo_);
    int ms_similar = ((mhi==mhi_) && (mlo==mlo_));
Martin Reinecke's avatar
Martin Reinecke committed
160

Martin Reinecke's avatar
Martin Reinecke committed
161
    mlo = mlo_; mhi = mhi_;
Martin Reinecke's avatar
Martin Reinecke committed
162
163
164

    if (!ms_similar)
      {
Martin Reinecke's avatar
Martin Reinecke committed
165
166
167
      alpha[mhi] = 1.;
      coef[mhi].a = coef[mhi].b = 0.;
      for (int l=mhi; l<=lmax; ++l)
Martin Reinecke's avatar
Martin Reinecke committed
168
        {
Martin Reinecke's avatar
Martin Reinecke committed
169
        double t = flm1[l+m]*flm1[l-m]*flm1[l+s]*flm1[l-s];
Martin Reinecke's avatar
Martin Reinecke committed
170
171
172
        double lt = 2*l+1;
        double l1 = l+1;
        double flp10=l1*lt*t;
Martin Reinecke's avatar
Martin Reinecke committed
173
174
175
176
177
        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
178
        else
Martin Reinecke's avatar
Martin Reinecke committed
179
180
181
          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
182
183
184
        }
      }

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

Martin Reinecke's avatar
Martin Reinecke committed
199
vector<double> sharp_Ylmgen::get_norm (int lmax, int spin)
Martin Reinecke's avatar
Martin Reinecke committed
200
201
  {
  const double pi = 3.141592653589793238462643383279502884197;
Martin Reinecke's avatar
Martin Reinecke committed
202
  vector<double> res(lmax+1);
Martin Reinecke's avatar
Martin Reinecke committed
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
  /* 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)
    {
    for (int l=0; l<=lmax; ++l)
      res[l]=1.;
    return res;
    }

  spinsign = (spin&1) ? -spinsign : spinsign;
  for (int l=0; l<=lmax; ++l)
    res[l] = (l<spin) ? 0. : spinsign*0.5*sqrt((2*l+1)/(4*pi));
  return res;
  }

Martin Reinecke's avatar
Martin Reinecke committed
223
vector<double> sharp_Ylmgen::get_d1norm (int lmax)
Martin Reinecke's avatar
Martin Reinecke committed
224
225
  {
  const double pi = 3.141592653589793238462643383279502884197;
Martin Reinecke's avatar
Martin Reinecke committed
226
  vector<double> res(lmax+1);
Martin Reinecke's avatar
Martin Reinecke committed
227
228
229
230
231
232
233

  for (int l=0; l<=lmax; ++l)
    res[l] = (l<1) ? 0. : 0.5*sqrt(l*(l+1.)*(2*l+1.)/(4*pi));
  return res;
  }

#pragma GCC visibility pop