00001 #include <math.h>
00002 #include <stdio.h>
00003 #include "plgndr.h"
00004 #ifndef PI
00005 #define PI 3.14159265359
00006 #endif
00007
00008
00009 double plgndr(int l, int m, double x)
00010 {
00011 double fact,pll,pmm,pmmp1,somx2;
00012 int i,ll;
00013
00014 if (m < 0 || m > l || fabs(x) > 1.0) {
00015 fprintf(stderr, "Bad arguments in routine plgndr()\n");
00016 fprintf(stderr, "l: %d m: %d x: %g\n", l, m, x);
00017 exit(1);
00018 }
00019 pmm=1.0;
00020 if (m > 0) {
00021 somx2=sqrt((1.0-x)*(1.0+x));
00022 fact=1.0;
00023 for (i=1;i<=m;i++) {
00024 pmm *= -fact*somx2;
00025 fact += 2.0;
00026 }
00027 }
00028 if (l == m)
00029 return pmm;
00030 else {
00031 pmmp1=x*(2*m+1)*pmm;
00032 if (l == (m+1))
00033 return pmmp1;
00034 else {
00035 for (ll=m+2;ll<=l;ll++) {
00036 pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);
00037 pmm=pmmp1;
00038 pmmp1=pll;
00039 }
00040 return pll;
00041 }
00042 }
00043 }
00044
00045
00046 double sphericalY(int l, int m, double theta)
00047 {
00048 int i;
00049 double lmm = 1.0, lpm;
00050
00051 for (i = 2; i <= l - m; i++)
00052 lmm *= (double)i;
00053 lpm = lmm;
00054 for (i = l - m + 1; i <= l + m; i++)
00055 lpm *= (double)i;
00056 return sqrt((2.*l + 1.)*lmm/(4.*PI*lpm))*plgndr(l, m, cos(theta));
00057 }
00058
00059
00060