//// spectral toolkit // copyright (c) 2005 university corporation for atmospheric research// licensed under the gnu general public license//#include "lobatto.h"using namespace spectral;/// Constructor for Gauss-Lobatto quadrature.  /// \param N total length of stencil on \f$ [-1,1] \f$lobatto::lobatto(int N) : gauss(N){  int m=N-1;  jpx = new real[m+2];  djpx= new real[m+2];      real *xj  = new real[m+2];  real *glw = new real[m+2];  real *dxx = new real[m+2];  real *xjr = new real[m+2];  real *xpt = new real[m+2];  real *dpt = new real[m+2];  real *xpt1= new real[m+2];  real *dpt1= new real[m+2];      real alpha, beta;  real rr, sm,  djp, xjp;      real dth,dd, aa,bb, cd,cs,ss,sd,tmp;           const int itn = 25;         // Max iteration number  alpha = 0.0;            beta = 0.0;                xj[0] = 1.0;  xj[m] = -1.0;        int nh = (m+1)/2;      jacobi(m+1,alpha,beta, 1.0,xpt, dpt);  jacobi(m+1,alpha,beta,-1.0,xpt1,dpt1);      dd =  xpt[m]*xpt1[m-1] - xpt1[m]*xpt[m-1];  aa = (xpt1[m+1]*xpt[m-1] - xpt1[m-1]*xpt[m+1] )/dd;  bb = (xpt1[m]*xpt[m+1] - xpt1[m+1]*xpt[m] )/dd;      dth = pi/(real)(2*m+1);  cd = std::cos(2.0*dth);  sd = std::sin(2.0*dth);  cs = std::cos(dth);  ss = std::sin(dth);      for(int k=1;k<nh;k++)    {      rr = cs;          real dl = 1.0;       int iteration=0;      dl=2.0*eps;      while((std::abs(dl) > eps)&&(iteration<itn))        {	  jacobi(m+1,alpha,beta,rr,xpt,dpt);	  xjp =  xpt[m+1]+aa* xpt[m]+bb* xpt[m-1];	  djp =  dpt[m+1]+aa* dpt[m]+bb* dpt[m-1];	  sm = 0.0;  	  for(int j=0;j<k;j++)	    sm = sm + 1.0 /(rr - xj[j]);	  dl = -(xjp /(djp - xjp*sm));  	  rr = rr + dl; 	  iteration++;        }       xj[k] = rr;      tmp = cs*cd-ss*sd;      ss = cs*sd+ss*cd;      cs = tmp;    }  for(int k=1;k<=nh;k++)    xj[m-k] = -xj[k];  if(m%2==0)     xj[nh] = 0.0;     for(int k=0;k<m+1;k++)    {      rr = xj[k];      jacobi(m,alpha,beta,rr,xpt,dpt);      xjp = xpt[m];      glw[k] = 2.0 / ((real)(m*(m+1)) * xjp*xjp);      dxx[m-k] = xjp;    }  for(int j = 0;j<m+1;j++)    {      xjr[m-j] = xj[j];      points[m-j] = xj[j];      weights[m-j] = glw[j];    }  delete [] xj;  delete [] glw;  delete [] dxx;  delete [] xjr;  delete [] xpt;  delete [] dpt;  delete [] xpt1;  delete [] dpt1;}/// Destructor for Gauss-Lobatto quadrature class.lobatto::~lobatto(){  delete [] jpx;  delete [] djpx;}    /// Internal routine.void lobatto::jacobi(int nox,real alpha,real beta,real x,real *px,real *dpx){  real k1,dk, ab, ac1,ac2,ac3,acx;    for(int k = 0;k<nox+1;k++)    {      px[k] = 0.0;          dpx[k] = 0.0;       }  switch(nox)    {    case 0:      {	px[nox] = 1.0;	dpx[nox] = 0.0;	break;      }    case 1:      {	px[0] = 1.0;	px[1] =  (1.0 + alpha)*x;	dpx[0] = 0.0; 	dpx[1] = 1.0 + alpha;	break;      }    default:      {	ab = alpha + beta; 	jpx[0] = 1.0; 	jpx[1] = (1.0 + alpha)*x;	djpx[0] = 0.0; 	djpx[1] = (1.0 + alpha);	for(int k = 1;k< nox;k++)	  {	    dk = (real)k;	    k1 = (real)(k+1);	    ac1 =  2.0*k1* (dk+ ab +1.0) * (2.0*dk + ab);	    ac2 =  (2.0*k1+ ab) * (2.0*dk +ab + 1.0) * (2.0*dk + ab);	    acx =  (2.0*dk + ab + 1.0) * (alpha*alpha - beta*beta) + x *ac2; 	    ac3 =  2.0*(dk + alpha) * (dk + beta) * (2.0*k + ab + 2.0);	    jpx[k+1] = (acx*jpx[k] - ac3*jpx[k-1]) / ac1;	    djpx[k+1] = (acx*djpx[k] + ac2*jpx[k] - ac3*djpx[k-1]) / ac1;	  }      }    }  for(int k = 0;k< nox+1;k++)                  {      px[k] = jpx[k];      dpx[k] = djpx[k];    }}