LevelS SHT library  3.50
sharp_ylmgen_c.c
1 /*
2  * This file is part of libsharp.
3  *
4  * libsharp is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation; either version 2 of the License, or
7  * (at your option) any later version.
8  *
9  * libsharp is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with libsharp; if not, write to the Free Software
16  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
17  */
18 
19 /*
20  * libsharp is being developed at the Max-Planck-Institut fuer Astrophysik
21  * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
22  * (DLR).
23  */
24 
25 /*
26  * Helper code for efficient calculation of Y_lm(theta,phi=0)
27  *
28  * Copyright (C) 2005-2016 Max-Planck-Society
29  * Author: Martin Reinecke
30  */
31 
32 #include <math.h>
33 #include <stdlib.h>
34 #include "sharp_ylmgen_c.h"
35 #include "c_utils.h"
36 
37 static inline void normalize (double *val, int *scale, double xfmax)
38  {
39  while (fabs(*val)>xfmax) { *val*=sharp_fsmall; ++*scale; }
40  if (*val!=0.)
41  while (fabs(*val)<xfmax*sharp_fsmall) { *val*=sharp_fbig; --*scale; }
42  }
43 
44 void sharp_Ylmgen_init (sharp_Ylmgen_C *gen, int l_max, int m_max, int spin)
45  {
46  const double inv_sqrt4pi = 0.2820947917738781434740397257803862929220;
47 
48  gen->lmax = l_max;
49  gen->mmax = m_max;
50  UTIL_ASSERT(spin>=0,"incorrect spin: must be nonnegative");
51  UTIL_ASSERT(l_max>=spin,"incorrect l_max: must be >= spin");
52  UTIL_ASSERT(l_max>=m_max,"incorrect l_max: must be >= m_max");
53  gen->s = spin;
54  UTIL_ASSERT((sharp_minscale<=0)&&(sharp_maxscale>0),
55  "bad value for min/maxscale");
56  gen->cf=RALLOC(double,sharp_maxscale-sharp_minscale+1);
57  gen->cf[-sharp_minscale]=1.;
58  for (int m=-sharp_minscale-1; m>=0; --m)
59  gen->cf[m]=gen->cf[m+1]*sharp_fsmall;
60  for (int m=-sharp_minscale+1; m<(sharp_maxscale-sharp_minscale+1); ++m)
61  gen->cf[m]=gen->cf[m-1]*sharp_fbig;
62  gen->powlimit=RALLOC(double,m_max+spin+1);
63  gen->powlimit[0]=0.;
64  const double ln2 = 0.6931471805599453094172321214581766;
65  const double expo=-400*ln2;
66  for (int m=1; m<=m_max+spin; ++m)
67  gen->powlimit[m]=exp(expo/m);
68 
69  gen->m = -1;
70  if (spin==0)
71  {
72  gen->rf = RALLOC(sharp_ylmgen_dbl2,gen->lmax+1);
73  gen->mfac = RALLOC(double,gen->mmax+1);
74  gen->mfac[0] = inv_sqrt4pi;
75  for (int m=1; m<=gen->mmax; ++m)
76  gen->mfac[m] = gen->mfac[m-1]*sqrt((2*m+1.)/(2*m));
77  gen->root = RALLOC(double,2*gen->lmax+5);
78  gen->iroot = RALLOC(double,2*gen->lmax+5);
79  for (int m=0; m<2*gen->lmax+5; ++m)
80  {
81  gen->root[m] = sqrt(m);
82  gen->iroot[m] = (m==0) ? 0. : 1./gen->root[m];
83  }
84  }
85  else
86  {
87  gen->m=gen->mlo=gen->mhi=-1234567890;
88  ALLOC(gen->fx,sharp_ylmgen_dbl3,gen->lmax+2);
89  for (int m=0; m<gen->lmax+2; ++m)
90  gen->fx[m].f[0]=gen->fx[m].f[1]=gen->fx[m].f[2]=0.;
91  ALLOC(gen->inv,double,gen->lmax+1);
92  gen->inv[0]=0;
93  for (int m=1; m<gen->lmax+1; ++m) gen->inv[m]=1./m;
94  ALLOC(gen->flm1,double,2*gen->lmax+1);
95  ALLOC(gen->flm2,double,2*gen->lmax+1);
96  for (int m=0; m<2*gen->lmax+1; ++m)
97  {
98  gen->flm1[m] = sqrt(1./(m+1.));
99  gen->flm2[m] = sqrt(m/(m+1.));
100  }
101  ALLOC(gen->prefac,double,gen->mmax+1);
102  ALLOC(gen->fscale,int,gen->mmax+1);
103  double *fac = RALLOC(double,2*gen->lmax+1);
104  int *facscale = RALLOC(int,2*gen->lmax+1);
105  fac[0]=1; facscale[0]=0;
106  for (int m=1; m<2*gen->lmax+1; ++m)
107  {
108  fac[m]=fac[m-1]*sqrt(m);
109  facscale[m]=facscale[m-1];
110  normalize(&fac[m],&facscale[m],sharp_fbighalf);
111  }
112  for (int m=0; m<=gen->mmax; ++m)
113  {
114  int mlo=gen->s, mhi=m;
115  if (mhi<mlo) SWAP(mhi,mlo,int);
116  double tfac=fac[2*mhi]/fac[mhi+mlo];
117  int tscale=facscale[2*mhi]-facscale[mhi+mlo];
118  normalize(&tfac,&tscale,sharp_fbighalf);
119  tfac/=fac[mhi-mlo];
120  tscale-=facscale[mhi-mlo];
121  normalize(&tfac,&tscale,sharp_fbighalf);
122  gen->prefac[m]=tfac;
123  gen->fscale[m]=tscale;
124  }
125  DEALLOC(fac);
126  DEALLOC(facscale);
127  }
128  }
129 
130 void sharp_Ylmgen_destroy (sharp_Ylmgen_C *gen)
131  {
132  DEALLOC(gen->cf);
133  DEALLOC(gen->powlimit);
134  if (gen->s==0)
135  {
136  DEALLOC(gen->rf);
137  DEALLOC(gen->mfac);
138  DEALLOC(gen->root);
139  DEALLOC(gen->iroot);
140  }
141  else
142  {
143  DEALLOC(gen->fx);
144  DEALLOC(gen->prefac);
145  DEALLOC(gen->fscale);
146  DEALLOC(gen->flm1);
147  DEALLOC(gen->flm2);
148  DEALLOC(gen->inv);
149  }
150  }
151 
152 void sharp_Ylmgen_prepare (sharp_Ylmgen_C *gen, int m)
153  {
154  if (m==gen->m) return;
155  UTIL_ASSERT(m>=0,"incorrect m");
156  gen->m = m;
157 
158  if (gen->s==0)
159  {
160  gen->rf[m].f[0] = gen->root[2*m+3];
161  gen->rf[m].f[1] = 0.;
162  for (int l=m+1; l<=gen->lmax; ++l)
163  {
164  double tmp=gen->root[2*l+3]*gen->iroot[l+1+m]*gen->iroot[l+1-m];
165  gen->rf[l].f[0] = tmp*gen->root[2*l+1];
166  gen->rf[l].f[1] = tmp*gen->root[l+m]*gen->root[l-m]*gen->iroot[2*l-1];
167  }
168  }
169  else
170  {
171  int mlo_=m, mhi_=gen->s;
172  if (mhi_<mlo_) SWAP(mhi_,mlo_,int);
173  int ms_similar = ((gen->mhi==mhi_) && (gen->mlo==mlo_));
174 
175  gen->mlo = mlo_; gen->mhi = mhi_;
176 
177  if (!ms_similar)
178  {
179  for (int l=gen->mhi; l<gen->lmax; ++l)
180  {
181  double t = gen->flm1[l+gen->m]*gen->flm1[l-gen->m]
182  *gen->flm1[l+gen->s]*gen->flm1[l-gen->s];
183  double lt = 2*l+1;
184  double l1 = l+1;
185  gen->fx[l+1].f[0]=l1*lt*t;
186  gen->fx[l+1].f[1]=gen->m*gen->s*gen->inv[l]*gen->inv[l+1];
187  t = gen->flm2[l+gen->m]*gen->flm2[l-gen->m]
188  *gen->flm2[l+gen->s]*gen->flm2[l-gen->s];
189  gen->fx[l+1].f[2]=t*l1*gen->inv[l];
190  }
191  }
192 
193  gen->preMinus_p = gen->preMinus_m = 0;
194  if (gen->mhi==gen->m)
195  {
196  gen->cosPow = gen->mhi+gen->s; gen->sinPow = gen->mhi-gen->s;
197  gen->preMinus_p = gen->preMinus_m = ((gen->mhi-gen->s)&1);
198  }
199  else
200  {
201  gen->cosPow = gen->mhi+gen->m; gen->sinPow = gen->mhi-gen->m;
202  gen->preMinus_m = ((gen->mhi+gen->m)&1);
203  }
204  }
205  }
206 
207 double *sharp_Ylmgen_get_norm (int lmax, int spin)
208  {
209  const double pi = 3.141592653589793238462643383279502884197;
210  double *res=RALLOC(double,lmax+1);
211  /* sign convention for H=1 (LensPix paper) */
212 #if 1
213  double spinsign = (spin>0) ? -1.0 : 1.0;
214 #else
215  double spinsign = 1.0;
216 #endif
217 
218  if (spin==0)
219  {
220  for (int l=0; l<=lmax; ++l)
221  res[l]=1.;
222  return res;
223  }
224 
225  spinsign = (spin&1) ? -spinsign : spinsign;
226  for (int l=0; l<=lmax; ++l)
227  res[l] = (l<spin) ? 0. : spinsign*0.5*sqrt((2*l+1)/(4*pi));
228  return res;
229  }
230 
231 double *sharp_Ylmgen_get_d1norm (int lmax)
232  {
233  const double pi = 3.141592653589793238462643383279502884197;
234  double *res=RALLOC(double,lmax+1);
235 
236  for (int l=0; l<=lmax; ++l)
237  res[l] = (l<1) ? 0. : 0.5*sqrt(l*(l+1.)*(2*l+1.)/(4*pi));
238  return res;
239  }
#define RALLOC(type, num)
void sharp_Ylmgen_destroy(sharp_Ylmgen_C *gen)
double * sharp_Ylmgen_get_norm(int lmax, int spin)
double * sharp_Ylmgen_get_d1norm(int lmax)
void sharp_Ylmgen_prepare(sharp_Ylmgen_C *gen, int m)
#define ALLOC(ptr, type, num)
#define DEALLOC(ptr)
#define UTIL_ASSERT(cond, msg)
void sharp_Ylmgen_init(sharp_Ylmgen_C *gen, int l_max, int m_max, int spin)

Generated on Mon Dec 10 2018 10:24:20 for LevelS SHT library