/* Legendre polynomial */

#include <math.h>
#include "../iolib/message.h"

float *assolegendre();
float Legpolyn();

void legpolyn_(f,n,x)
float *f,*x;
int *n;
{
  *f=Legpolyn(*n,*x);
}

float Legpolyn(n,x)
float x;
int n;
{
float theta,f;
	if(x<-1.0 || x> 1.0) STOP(Legpolyn: illegal argument);
	if(!n)        f=0.7071068;
	else if(n==1) f=1.2247449*x;
	else if(n==2) f=0.7905694*(3.0*x*x-1.0);
	else if(n==3) f=0.9354144*x*(5.*x*x-3.0);
	else if(n==4) f=0.2651650*(x*x*(35.*x*x-30.)+3.0);
	else if(n==5) f=0.2931510*x*(x*x*(63.*x*x-70.)+15.0);
	else
	{ theta=acos((double)x);
	  f=*assolegendre(theta,n,"f");
	}
	return(f);
}

float dLeg_dx(n,x)  /*derivative of Legpolyn*/
float x;
int n;
{
float y;
float Legpolyn();

  if(x<-1.0 || x> 1.0) STOP(dLeg_dx: illegal argument);
  if(!n) y=0.0;
  else if(n==1) y=1.2247449;
  else if(n==2) y=4.7434164*x;
  else if(n==3) y=2.8062432*(5.*x*x-1.);
  else if(n==4) y=5.3033000*(7.*x*x*x-3.*x);
  else if(n==5) y=4.397265*(21.*x*x*x*x-14.*x*x+1.);
  else 
  { double arg,nn;
    nn=n;
    arg=(2.*nn+1.0)/(2.*nn-1.0);
    y=sqrt(arg)*(nn*Legpolyn(n-1,x)+x*dLeg_dx(n-1,x));
  }
  
  return(y);
  
}


float d2Leg_dx2(n,x) /* 2nd derivative of Legpolyn*/
float x;
int n;
{
float y;
float dLeg_dx();

  if(x<-1.0 || x> 1.0) STOP(d2Leg_dx2: illegal argument);
  if(n<2)       y=0.0;
  else if(n==2) y=4.7434164;
  else if(n==3) y=28.062432*x;
  else if(n==4) y=15.909900*(7.*x*x-1.);
  else if(n==5) y=123.123420*(3.*x*x*x-x);
  else
  { double arg,nn;
    nn=n;
    arg=(2.*nn+1.0)/(2.*nn-1.0);
    y=sqrt(arg)*((nn+1.)*dLeg_dx(n-1,x)+x*d2Leg_dx2(n-1,x));
  }
  return(y);
}