//// spectral toolkit // copyright (c) 2005 university corporation for atmospheric research// licensed under the gnu general public license//#include "legendre.h"using namespace spectral;/// Constructor for Gauss-Legendre quadrature./// \param M total length of stencil on \f$ (-1,1) \f$legendre::legendre(int M) : gauss(M){  int i,ns2,nhalf,nix;  real x,dtheta,cmax;  real zero,zprev,zhold,zlast;  real dthalf,cz,pb,dpb,dcor,sgnd,sum,temp;    const real pis2=0.5L*pi;    switch(n)    {    case (1):      points[0] = pis2;      weights[0] = 2.0;      break;    case (2):      x = 1.0/root3;      points[0] = x;      points[1] = -x;      weights[0] = 1.0;      weights[1] = 1.0;      break;    default:      ns2 = M/2;      nhalf =(M+1)/2;      cz=cpdp(points,weights);      dtheta = pi/(real)(M+1);      dthalf = dtheta*0.5;      cmax = 0.2*dtheta;      if(n%2)	{	  zero = pis2-dtheta;	  zprev = pis2;	  zlast = zprev;	  nix = nhalf-1;	}      else	{	  zero = pis2-dthalf;	  zprev = pis2;	  zlast = zprev;	  nix = nhalf;	}          do	{	  int iteration=0;	  while((std::abs(zero-zlast) > eps*std::abs(zero))&&iteration<10)	    {	      zlast = zero;	      tpdp(zero,cz,points,weights,pb,dpb);	      dcor = sgn(pb/dpb)*min(std::abs(pb/dpb),cmax);	      zero -= dcor;	      iteration++;	    }	  points[nix-1] = zero;	  zhold = zero;                    	  //	  // yakimiw's formula permits using old pb and dpb	  //                    	  temp = (dpb+pb*std::cos(zlast)/std::sin(zlast));	  weights[nix-1] = (n+n+1)/(temp*temp);	  nix--;	  if(nix==0)	    break;	  if(nix==nhalf-1)	    zero = 3.0*zero-pi;	  if(nix <nhalf-1)	    zero= zero+zero-zprev;	  zprev = zhold;	} while (nix!=0);      if(n%2)	{	  points[nhalf-1] = pis2;	  tpdp(pis2,cz,points,weights,pb,dpb);	  weights[nhalf-1] = (n+n+1)/(dpb*dpb);	}      for(i=1;i<ns2+1;i++)	{	  weights[n-i] = weights[i-1];	  points[n-i] = pi - points[i-1];	}      for(i=0;i<n;i++)	points[i]=std::cos(points[i]);                          sum = 0.0;      for(i=1;i<n+1;i++)	sum += weights[i-1];      for(i=1;i<n+1;i++)	weights[i-1] *= 2.0/sum;      if(n%2)	points[nhalf-1]=0;          break;    }}/// Internal method for computing intermediate values.real legendre::cpdp(real *cp, real *dcp){  real t1,t2,t3,t4;  real cz;  int ncp,j,i;  for(i=0;i<n;i++)    {      cp[i]=0.0;      dcp[i]=0.0;    }  ncp = (n+1)/2;  t1 = -1.0;  t2 = n+1.0;  t3 = 0.0;  t4 = n+n+1.0;  if(n%2)    {      cp[n-1] = 1.0;      for(j=n-1;j>ncp-1;j--)        {	  t1 +=2.0;	  t2 -=1.0;	  t3 +=1.0;	  t4 -=2.0;	  cp[j-1] = (t1*t2)/(t3*t4)*cp[j];        }      for(j=ncp;j<n+1;j++)	dcp[j-1] = ((j-ncp+1)*2-1)*cp[j-1];    }  else    {      cp[n-1] = 1.0;      for(j=n;j>ncp+1;j--)        {	  t1 +=2.0;	  t2 -=1.0;	  t3 +=1.0;	  t4 -=2.0;	  cp[j-2] = (t1*t2)/(t3*t4)*cp[j-1];        }      t1 +=2.0;      t2 -=1.0;      t3 +=1.0;      t4 -=2.0;      for(j=ncp;j<n+1;j++)	dcp[j-1] = (j-ncp)*2*cp[j-1];    }  cz = (t1*t2)/(t3*t4)*cp[ncp];  return(cz);}/// Internal routine that/// computes \f$ pn(\theta) \f$ and its derivative \f$ dpb(\theta) \f$/// with respect to \f$ \theta \f$.void legendre::tpdp(real theta,real cz,real *cp,real *dcp,real &pb,real &dpb){  real cdt,sdt,sth,cth,chh;  int kdo,k;  cdt=std::cos(theta+theta);  sdt=std::sin(theta+theta);  dpb = 0.0;      if(n%2)    {      kdo=(n+1)/2;      pb = 0.0;      cth = std::cos(theta);      sth = std::sin(theta);    }  else    {      kdo = n/2;      pb = 0.5 * cz;      cth = cdt;      sth = sdt;    }  for(k=1;k<kdo+1;k++)    {      pb += cp[kdo+k-1-n%2]*cth;      dpb -= dcp[kdo+k-1-n%2]*sth;      chh = cdt*cth-sdt*sth;      sth = sdt*cth+cdt*sth;      cth = chh;    }}