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