alf.cpp

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

//
// spectral toolkit 
// copyright 2005 university corporation for atmospheric research
// licensed under the gnu general public license
//
// $Id: alf.cpp,v 1.5 2005/04/13 17:06:25 rodney Exp $
//

#include "alf.h"

namespace spectral
{

  /// Creates an instance of the associated Legendre function \f$ P_n^m \f$.
  /// \param N order of Legendre polynomial
  /// \param M order of Fourier mode, \f$ 0 \leq M \leq N \f$
  /// \throws bad_parameter if m>n

  alf::alf(int N,int M) 
  {
    real fden,fnum,fnmh,fnmsq,fnnp1,pm1,fk,a1,b1,c1,t1,t2,cp2,twoton;
    int  nmms2,nex,i,l,ma;
    n=N;
    m=M;
    cp=new real[n];
    if(m>n)
      throw bad_parameter();
    if(m<0)
      ma=-m;
    else
      ma=m;
    cp[0]=0.0;
    if(n==0)
      {
	cp[0]=root2; 
      }
    else if(n==1)
      {
	if(ma==0) 
	  cp[0]=root3/root2;
	else
	  cp[0]=root3/2.0;
	if(m == -1) 
	  cp[0] = -cp[0];
      }
    else
      {
        
	real sc10=1024.0;
	real sc20=sc10*sc10;
	real sc40=sc20*sc20;
	if((n+ma)%2 == 0)
	  {
	    nmms2 = (n-ma)/2;
	    fnum = n+ma+1;
	    fnmh = n-ma+1;
	    pm1 = 1.0;
	  }
	else
	  {
	    nmms2 = (n-ma-1)/2;
	    fnum = n+ma+2;
	    fnmh = n-ma+2;
	    pm1 = -1.0;
	  }
	t1 = 1.0/sc20;
	nex = 20;
	fden = 2.0;
	if(nmms2>=1)
	  {
	    for(i=0;i<nmms2;i++)
	      {
		t1 = fnum*t1/fden;
		if(t1>sc20)
		  {
		    t1 = t1/sc40;
		    nex += 40;
		  }
		fnum += 2.0;
		fden += 2.0;
	      }
	  }
	twoton= std::pow(2.0,n-1-nex);
	t1 = t1/twoton;
	if((ma/2)%2 != 0) 
	  t1 = -t1;
	t2 = 1.0;
	if(ma != 0)
	  {
	    for(i=0;i<ma;i++)
	      {
		t2 = fnmh*t2/(fnmh+pm1);
		fnmh += 2.0;
	      }
	  }
	cp2 = t1*std::sqrt((n+0.5)*t2);
	fnnp1 = n*(n+1);
	fnmsq = fnnp1-2.0*ma*ma;
	l = (n+1)/2;
	if(n%2 == 0 && ma%2==0) 
	  l++;
	cp[l-1] = cp2;
	if( (m < 0) && (ma%2 != 0)) 
	  cp[l-1] = -cp[l-1];
	if(l<=1) 
	  return;
	fk = n;
	a1 = (fk-2.0)*(fk-1.0)-fnnp1;
	b1 = 2.0*(fk*fk-fnmsq);
	cp[l-2] = b1*cp[l-1]/a1;
	l--;
	while(l>1)
	  {
	    fk = fk-2.0;
	    a1 = (fk-2.0)*(fk-1.0)-fnnp1;
	    b1 = -2.0*(fk*fk-fnmsq);
	    c1 = (fk+1.0)*(fk+2.0)-fnnp1;
	    cp[l-2] = -(b1*cp[l-1]+c1*cp[l])/a1;
	    l--;
	  }
      }
  }

  /// Destroys associated Legendre function object.

  alf::~alf()
  {
    delete [] cp;
  }

  /// Evaluates the associated Legendre function at a point.
  /// The result is the normalized associated Legendre function 
  /// \f[ \overline{P_n^m}(x) = (-1)^m \sqrt{\frac{(2n+1)(n-m)!}{2(n+m)!}} P_n^m(x) \f]
  /// \param x real number in [-1,1]
  /// \returns \f$ \overline{P_n^m}(x) \f$

  real alf::operator() (real x)
  {
    real theta=std::acos(x);
    real cth,sth,chh,pmn;
    real cdt = std::cos(theta+theta);
    real sdt = std::sin(theta+theta);
    int k,kdo;
    
    if(n%2==0)
      {   
	if(m%2==0)
	  {
	    /* n even and m even */
	    kdo=n/2;
	    pmn = 0.5*cp[0];
	    if(n==0) return(pmn);
	    cth = cdt;
	    sth = sdt;
	    for(k=0;k<kdo;k++)
	      {
		pmn += cp[k+1]*cth;
		chh = cdt*cth-sdt*sth;
		sth = sdt*cth+cdt*sth;
		cth = chh;
	      }
	  }
	else
	  {  
	    /* n even and m odd */
	    kdo = n/2;
	    pmn = 0.0;
	    cth = cdt;
	    sth = sdt;
	    for(k=0;k<kdo;k++)
	      {
		pmn  += cp[k]*sth;
		chh = cdt*cth-sdt*sth;
		sth = sdt*cth+cdt*sth;
		cth = chh;
	      }
	  }
      }
    else
      {
	if(m%2==0)
	  {
	    /* n odd and m even */
	    kdo = (n+1)/2;
	    pmn = 0.0;
	    cth = std::cos(theta);
	    sth = std::sin(theta);
	    for(k=0;k<kdo;k++)
	      {
		pmn += cp[k]*cth;
		chh = cdt*cth-sdt*sth;
		sth = sdt*cth+cdt*sth;
		cth = chh;
	      }
	  }
	else
	  { 
	    /* n odd and m odd */
	    kdo = (n+1)/2;
	    pmn = 0.0;
	    cth = std::cos(theta);
	    sth = std::sin(theta);
	    for(k=0;k<kdo;k++)
	      {
		pmn += cp[k]*sth;
		chh = cdt*cth-sdt*sth;
		sth = sdt*cth+cdt*sth;
		cth = chh;
	      }
	  }
      }
    return(pmn);
  }

}