/* calculate associated Legendre functions of various conventions*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "const.h"
#include "message.h"

float **farray2();
void free_farray2();
float factor();
void rotmatrix();

float *assolegendre(theta,l,conventn)
float theta;
int l;
char *conventn; /*available convention: d(l,m), 0<=m<=l;
                 *                      p(l,m), 0<=m<=l;
                 *                      f(l)           ;
                 *                      X(l,m), 0<=m<=l;
                 *                      X^(l,m),0<=m<=l;
                 *                      P(l,m), 0<=m<=l;
                 */
{
static int current_dim=0;
static float *x,**d;
static int mem;
int dim,m;
void *test;
float fac,sign;

	dim=2*l+1;
	if(dim>current_dim)
	{ if(current_dim)
	  { free(x);
	    free_farray2(d,0,0,mem);
	  }
	  if( (x=(float *)malloc(sizeof(float)*dim)) ==NULL)
	    stop("assolegendre: malloc failed 1");
	  d=farray2(0,0,-l,l);
	  current_dim=dim;
	  mem=-l;
	}
	
	rotmatrix(0,l,theta,d);
	
	if(!strcmp(conventn,"d"))
	{ for(m=0;m<=l;m++) x[m]=d[0][m]; 
	}
	else if(!strcmp(conventn,"p"))
	{ for(sign=1.0,m=0;m<=l;sign=-sign,m++)
	  { fac=(float)(2*l+1);
	    if(m) fac*=2.0;
	    fac=sign*(float) sqrt( (double)fac );
	    x[m]=fac*d[0][m];
	  }
	}
	else if(!strcmp(conventn,"f"))
	{ fac=0.5*(float)(2*l+1);
	  fac=(float) sqrt( (double)fac );
	  *x=fac*d[0][0];
	}
	else if(!strcmp(conventn,"X"))
	{ for(m=0;m<=l;m++)
	  { fac=((float)(2*l+1))/(4.*PI);
	    fac=(float)sqrt((double)fac);
	    x[m]=fac*d[0][m];
	  }
	}
	else if(!strcmp(conventn,"X^"))
	{ for(m=0;m<=l;m++)
	  { fac=((float)(2*l+1))/(8.*PI);
	    if(m) fac*=2.0;
	    fac=(float)sqrt((double)fac);
	    x[m]=fac*d[0][m];
	  }
	}
	else if(!strcmp(conventn,"P"))
	{ for(sign=1.0,m=0;m<=l;sign=-sign,m++)
	  { fac=factor(l,m);
	    fac=sign*(float)sqrt((double)fac);
	    x[m]=fac*d[0][m];
	  }
	}
	else stop("assolegendre: illegal convention");
	
	return(x);
}

float factor(l,m) /* (l+m)!/(l-m)! | m>=0 */
int l,m;
{
float f;
int i;  if(m<0 || m>l) stop("assolegendre_factor: illegal m");
	f=1.0;
	for(i=m;i>=-m+1;i--) f*=(float)(l+i);
	return(f);
}