#include<math.h>
|
#include<stdlib.h>
|
#include "general.h"
|
#include "sh.h"
|
|
/* Gives you legendre polynomials. Taken from NR, p. 254 */
|
ts_double plgndr(ts_int l, ts_int m, ts_float x){
|
ts_double fact, pll, pmm, pmmp1, somx2;
|
ts_int i,ll;
|
|
#ifdef TS_DOUBLE_DOUBLE
|
if(m<0 || m>l || fabs(x)>1.0)
|
fatal("Bad arguments in routine plgndr",1);
|
#endif
|
#ifdef TS_DOUBLE_FLOAT
|
if(m<0 || m>l || fabsf(x)>1.0)
|
fatal("Bad arguments in routine plgndr",1);
|
#endif
|
#ifdef TS_DOUBLE_LONGDOUBLE
|
if(m<0 || m>l || fabsl(x)>1.0)
|
fatal("Bad arguments in routine plgndr",1);
|
#endif
|
pmm=1.0;
|
if (m>0) {
|
#ifdef TS_DOUBLE_DOUBLE
|
somx2=sqrt((1.0-x)*(1.0+x));
|
#endif
|
#ifdef TS_DOUBLE_FLOAT
|
somx2=sqrtf((1.0-x)*(1.0+x));
|
#endif
|
#ifdef TS_DOUBLE_LONGDOUBLE
|
somx2=sqrtl((1.0-x)*(1.0+x));
|
#endif
|
fact=1.0;
|
for (i=1; i<=m;i++){
|
pmm *= -fact*somx2;
|
fact +=2.0;
|
}
|
}
|
|
if (l == m) return pmm;
|
else {
|
pmmp1=x*(2*m+1)*pmm;
|
if(l==(m+1)) return(pmmp1);
|
else {
|
pll=0; /* so it can not be uninitialized */
|
for(ll=m+2;ll<=l;ll++){
|
pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);
|
pmm=pmmp1;
|
pmmp1=pll;
|
}
|
return(pll);
|
}
|
}
|
}
|
|
|
/*Computes Y(l,m,theta,fi) (Miha's definition that is different from common definition for factor srqt(1/(2*pi)) */
|
ts_double shY(ts_int l,ts_int m,ts_double theta,ts_double fi){
|
ts_double fac1, fac2, K;
|
int i;
|
|
if(l<0 || m>l || m<-l)
|
fatal("Error using shY function!",1);
|
|
fac1=1.0;
|
for(i=1; i<=l-m;i++){
|
fac1 *= i;
|
}
|
fac2=1.0;
|
for(i=1; i<=l+m;i++){
|
fac2 *= i;
|
}
|
|
if(m==0){
|
K=sqrt(1.0/(2.0*M_PI));
|
}
|
else if (m>0) {
|
K=sqrt(1.0/(M_PI))*cos(m*fi);
|
}
|
else {
|
//K=pow(-1.0,abs(m))*sqrt(1.0/(2.0*M_PI))*cos(m*fi);
|
if(abs(m)%2==0)
|
K=sqrt(1.0/(M_PI))*sin(m*fi);
|
else
|
K=-sqrt(1.0/(M_PI))*sin(m*fi);
|
}
|
|
return K*sqrt((2.0*l+1.0)/2.0*fac1/fac2)*plgndr(l,abs(m),cos(theta));
|
}
|
