legendre.cpp

The Spectral Toolkit is a C++ spectral transform library written by Rodney James and Chuck Panaccion

//
// 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;
    }
}