| commit | author | age | ||
| 88f451 | 1 | #include<math.h> |
| SP | 2 | #include<stdlib.h> |
| 3 | #include "general.h" | |
| 4 | #include "sh.h" | |
| 5 | ||
| 6 | /* Gives you legendre polynomials. Taken from NR, p. 254 */ | |
| 7 | ts_double plgndr(ts_int l, ts_int m, ts_float x){ | |
| 8 | ts_double fact, pll, pmm, pmmp1, somx2; | |
| 9 | ts_int i,ll; | |
| 10 | ||
| 11 | #ifdef TS_DOUBLE_DOUBLE | |
| 12 | if(m<0 || m>l || fabs(x)>1.0) | |
| 13 | fatal("Bad arguments in routine plgndr",1); | |
| 14 | #endif | |
| 15 | #ifdef TS_DOUBLE_FLOAT | |
| 16 | if(m<0 || m>l || fabsf(x)>1.0) | |
| 17 | fatal("Bad arguments in routine plgndr",1); | |
| 18 | #endif | |
| 19 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 20 | if(m<0 || m>l || fabsl(x)>1.0) | |
| 21 | fatal("Bad arguments in routine plgndr",1); | |
| 22 | #endif | |
| 23 | pmm=1.0; | |
| 24 | if (m>0) { | |
| 25 | #ifdef TS_DOUBLE_DOUBLE | |
| 26 | somx2=sqrt((1.0-x)*(1.0+x)); | |
| 27 | #endif | |
| 28 | #ifdef TS_DOUBLE_FLOAT | |
| 29 | somx2=sqrtf((1.0-x)*(1.0+x)); | |
| 30 | #endif | |
| 31 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 32 | somx2=sqrtl((1.0-x)*(1.0+x)); | |
| 33 | #endif | |
| 34 | fact=1.0; | |
| 35 | for (i=1; i<=m;i++){ | |
| 36 | pmm *= -fact*somx2; | |
| 37 | fact +=2.0; | |
| 38 | } | |
| 39 | } | |
| 40 | ||
| 41 | if (l == m) return pmm; | |
| 42 | else { | |
| 43 | pmmp1=x*(2*m+1)*pmm; | |
| 44 | if(l==(m+1)) return(pmmp1); | |
| 45 | else { | |
| 46 | pll=0; /* so it can not be uninitialized */ | |
| 47 | for(ll=m+2;ll<=l;ll++){ | |
| 48 | pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m); | |
| 49 | pmm=pmmp1; | |
| 50 | pmmp1=pll; | |
| 51 | } | |
| 52 | return(pll); | |
| 53 | } | |
| 54 | } | |
| 55 | } | |
| 56 | ||
| 57 | ||
| 58 | /*Computes Y(l,m,theta,fi) (Miha's definition that is different from common definition for factor srqt(1/(2*pi)) */ | |
| 59 | ts_double shY(ts_int l,ts_int m,ts_double theta,ts_double fi){ | |
| 60 | ts_double fac1, fac2, K; | |
| 61 | int i; | |
| 62 | ||
| 63 | if(l<0 || m>l || m<-l) | |
| 64 | fatal("Error using shY function!",1); | |
| 65 | ||
| 66 | fac1=1.0; | |
| 67 | for(i=1; i<=l-m;i++){ | |
| 68 | fac1 *= i; | |
| 69 | } | |
| 70 | fac2=1.0; | |
| 71 | for(i=1; i<=l+m;i++){ | |
| 72 | fac2 *= i; | |
| 73 | } | |
| 74 | ||
| 75 | if(m==0){ | |
| 76 | K=sqrt(1.0/(2.0*M_PI)); | |
| 77 | } | |
| 78 | else if (m>0) { | |
| 79 | K=sqrt(1.0/(M_PI))*cos(m*fi); | |
| 80 | } | |
| 81 | else { | |
| 82 | //K=pow(-1.0,abs(m))*sqrt(1.0/(2.0*M_PI))*cos(m*fi); | |
| 83 | if(abs(m)%2==0) | |
| 84 | K=sqrt(1.0/(M_PI))*sin(m*fi); | |
| 85 | else | |
| 86 | K=-sqrt(1.0/(M_PI))*sin(m*fi); | |
| 87 | } | |
| 88 | ||
| 89 | return K*sqrt((2.0*l+1.0)/2.0*fac1/fac2)*plgndr(l,abs(m),cos(theta)); | |
| 90 | } | |
