37 static inline void normalize (
double *val,
int *scale,
double xfmax)
39 while (fabs(*val)>xfmax) { *val*=sharp_fsmall; ++*scale; }
41 while (fabs(*val)<xfmax*sharp_fsmall) { *val*=sharp_fbig; --*scale; }
46 const double inv_sqrt4pi = 0.2820947917738781434740397257803862929220;
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");
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);
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);
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)
81 gen->root[m] = sqrt(m);
82 gen->iroot[m] = (m==0) ? 0. : 1./gen->root[m];
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);
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)
98 gen->flm1[m] = sqrt(1./(m+1.));
99 gen->flm2[m] = sqrt(m/(m+1.));
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)
108 fac[m]=fac[m-1]*sqrt(m);
109 facscale[m]=facscale[m-1];
110 normalize(&fac[m],&facscale[m],sharp_fbighalf);
112 for (
int m=0; m<=gen->mmax; ++m)
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);
120 tscale-=facscale[mhi-mlo];
121 normalize(&tfac,&tscale,sharp_fbighalf);
123 gen->fscale[m]=tscale;
154 if (m==gen->m)
return;
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)
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];
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_));
175 gen->mlo = mlo_; gen->mhi = mhi_;
179 for (
int l=gen->mhi; l<gen->lmax; ++l)
181 double t = gen->flm1[l+gen->m]*gen->flm1[l-gen->m]
182 *gen->flm1[l+gen->s]*gen->flm1[l-gen->s];
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];
193 gen->preMinus_p = gen->preMinus_m = 0;
194 if (gen->mhi==gen->m)
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);
201 gen->cosPow = gen->mhi+gen->m; gen->sinPow = gen->mhi-gen->m;
202 gen->preMinus_m = ((gen->mhi+gen->m)&1);
209 const double pi = 3.141592653589793238462643383279502884197;
210 double *res=
RALLOC(
double,lmax+1);
213 double spinsign = (spin>0) ? -1.0 : 1.0;
215 double spinsign = 1.0;
220 for (
int l=0; l<=lmax; ++l)
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));
233 const double pi = 3.141592653589793238462643383279502884197;
234 double *res=
RALLOC(
double,lmax+1);
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));
#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 UTIL_ASSERT(cond, msg)
void sharp_Ylmgen_init(sharp_Ylmgen_C *gen, int l_max, int m_max, int spin)