Polylib.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: Polylib.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description:
//
///////////////////////////////////////////////////////////////////////////////
 
#include "Polylib.h"
#include <cfloat>
#include <cmath>
#include <complex>
#include <cstdio>
 
#include <boost/core/ignore_unused.hpp>
 
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
/// Maximum number of iterations in polynomial defalation routine Jacobz
#define STOP 30
/// Precision tolerance for two points to be similar
#define EPS 100 * DBL_EPSILON
/// return the sign(b)*a
#define sign(a, b) ((b) < 0 ? -fabs(a) : fabs(a))
 
namespace Polylib
{
 
/// The following function is used to circumvent/reduce "Subtractive
/// Cancellation" The expression 1/dz  is replaced by optinvsub(.,.) Added on 26
/// April 2017
double optdiff(double xl, double xr)
{
    double m_xln, m_xrn;
    int m_expn;
    int m_digits = static_cast<int>(fabs(floor(log10(DBL_EPSILON))) - 1);
 
    if (fabs(xl - xr) < 1.e-4)
    {
 
        m_expn = static_cast<int>(floor(log10(fabs(xl - xr))));
        m_xln  = xl * powl(10.0L, -m_expn) -
                floor(xl * powl(10.0L,
                                -m_expn)); // substract the digits overlap part
        m_xrn =
            xr * powl(10.0L, -m_expn) -
            floor(xl *
                  powl(10.0L,
                       -m_expn)); // substract the common digits overlap part
        m_xln =
            round(m_xln * powl(10.0L, m_digits + m_expn)); // git rid of rubbish
        m_xrn = round(m_xrn * powl(10.0L, m_digits + m_expn));
 
        return powl(10.0L, -m_digits) * (m_xln - m_xrn);
    }
    else
    {
        return (xl - xr);
    }
}
 
double laginterp(double z, int j, const double *zj, int np)
{
    double temp = 1.0;
    for (int i = 0; i < np; i++)
    {
        if (j != i)
        {
            temp *= optdiff(z, zj[i]) / (zj[j] - zj[i]);
        }
    }
    return temp;
}
/// Define whether to use polynomial deflation (1)  or tridiagonal solver (0).
#define POLYNOMIAL_DEFLATION 0
 
#ifdef POLYNOMIAL_DEFLATION
/// zero determination using Newton iteration with polynomial deflation
#define jacobz(n, z, alpha, beta) Jacobz(n, z, alpha, beta)
#else
/// zero determination using eigenvalues of tridiagaonl matrix
#define jacobz(n, z, alpha, beta) JacZeros(n, z, alpha, beta)
#endif
 
/* local functions */
static void Jacobz(const int n, double *z, const double alpha,
                   const double beta);
// static void   JacZeros (const int n, double *a, const double alpha,
//    const double beta);
// static void   TriQL    (const int n, double *d, double *e);
static void TriQL(const int, double *, double *, double **);
double gammaF(const double);
double gammaFracGammaF(const int, const double, const int, const double);
static void RecCoeff(const int, double *, double *, const double, const double);
void JKMatrix(int, double *, double *);
void chri1(int, double *, double *, double *, double *, double);
 
/**
\brief  Gauss-Jacobi zeros and weights.
 
\li Generate \a np Gauss Jacobi zeros, \a z, and weights,\a w,
associated with the Jacobi polynomial \f$ P^{\alpha,\beta}_{np}(z)
\f$,
 
\li Exact for polynomials of order \a 2np-1 or less
*/
 
void zwgj(double *z, double *w, const int np, const double alpha,
          const double beta)
{
    int i;
    double fac, one = 1.0, two = 2.0, apb = alpha + beta;
 
    jacobz(np, z, alpha, beta);
    jacobd(np, z, w, np, alpha, beta);
 
    fac = pow(two, apb + one) * gammaFracGammaF(np + 1, alpha, np + 1, 0.0) *
          gammaFracGammaF(np + 1, beta, np + 1, apb);
 
    for (i = 0; i < np; ++i)
        w[i] = fac / (w[i] * w[i] * (one - z[i] * z[i]));
 
    return;
}
 
/**
\brief  Gauss-Radau-Jacobi zeros and weights with end point at \a z=-1.
 
\li Generate \a np Gauss-Radau-Jacobi zeros, \a z, and weights,\a w,
associated with the  polynomial \f$(1+z) P^{\alpha,\beta+1}_{np-1}(z)
\f$.
 
\li  Exact for polynomials of order \a 2np-2 or less
*/
 
void zwgrjm(double *z, double *w, const int np, const double alpha,
            const double beta)
{
 
    if (np == 1)
    {
        z[0] = 0.0;
        w[0] = 2.0;
    }
    else
    {
        int i;
        double fac, one = 1.0, two = 2.0, apb = alpha + beta;
 
        z[0] = -one;
        jacobz(np - 1, z + 1, alpha, beta + 1);
        jacobfd(np, z, w, NULL, np - 1, alpha, beta);
 
        fac = pow(two, apb) * gammaFracGammaF(np, alpha, np, 0.0) *
              gammaFracGammaF(np, beta, np + 1, apb);
        fac /= (beta + np);
 
        for (i = 0; i < np; ++i)
            w[i] = fac * (1 - z[i]) / (w[i] * w[i]);
        w[0] *= (beta + one);
    }
 
    return;
}
 
/**
\brief  Gauss-Radau-Jacobi zeros and weights with end point at \a z=1
 
 
\li Generate \a np Gauss-Radau-Jacobi zeros, \a z, and weights,\a w,
associated with the  polynomial \f$(1-z) P^{\alpha+1,\beta}_{np-1}(z)
\f$.
 
\li Exact for polynomials of order \a 2np-2 or less
*/
 
void zwgrjp(double *z, double *w, const int np, const double alpha,
            const double beta)
{
 
    if (np == 1)
    {
        z[0] = 0.0;
        w[0] = 2.0;
    }
    else
    {
        int i;
        double fac, one = 1.0, two = 2.0, apb = alpha + beta;
 
        jacobz(np - 1, z, alpha + 1, beta);
        z[np - 1] = one;
        jacobfd(np, z, w, NULL, np - 1, alpha, beta);
 
        fac = pow(two, apb) * gammaFracGammaF(np, alpha, np, 0.0) *
              gammaFracGammaF(np, beta, np + 1, apb);
        fac /= (alpha + np);
 
        for (i = 0; i < np; ++i)
            w[i] = fac * (1 + z[i]) / (w[i] * w[i]);
        w[np - 1] *= (alpha + one);
    }
 
    return;
}
 
/**
\brief  Gauss-Lobatto-Jacobi zeros and weights with end point at \a z=-1,\a 1
 
 
\li Generate \a np Gauss-Lobatto-Jacobi points, \a z, and weights, \a w,
associated with polynomial \f$ (1-z)(1+z) P^{\alpha+1,\beta+1}_{np-2}(z) \f$
\li Exact for polynomials of order \a 2np-3 or less
*/
 
void zwglj(double *z, double *w, const int np, const double alpha,
           const double beta)
{
 
    if (np == 1)
    {
        z[0] = 0.0;
        w[0] = 2.0;
    }
    else if (np == 2)
    {
        z[0] = -1.0;
        z[1] = 1.0;
 
        w[0] = 1.0;
        w[1] = 1.0;
    }
    else
    {
        int i;
        double fac, one = 1.0, apb = alpha + beta, two = 2.0;
 
        z[0]      = -one;
        z[np - 1] = one;
        jacobz(np - 2, z + 1, alpha + one, beta + one);
        jacobfd(np, z, w, NULL, np - 1, alpha, beta);
 
        fac = pow(two, apb + 1) * gammaFracGammaF(np, alpha, np, 0.0) *
              gammaFracGammaF(np, beta, np + 1, apb);
        fac /= (np - 1);
 
        for (i = 0; i < np; ++i)
            w[i] = fac / (w[i] * w[i]);
        w[0] *= (beta + one);
        w[np - 1] *= (alpha + one);
    }
 
    return;
}
 
/**
\brief  Gauss-Kronrod-Jacobi zeros and weights.
 
\li Generate \a npt=2*np+1 Gauss-Kronrod Jacobi zeros, \a z, and weights,\a w,
associated with the Jacobi polynomial \f$ P^{\alpha,\beta}_{np}(z)
\f$,
 
\li Exact for polynomials of order \a 3np+1 or less
*/
void zwgk(double *z, double *w, const int npt, const double alpha,
          const double beta)
{
 
    int np = (npt - 1) / 2;
 
    int i, j;
 
    // number of kronrod points associated with the np gauss rule
    int kpoints = 2 * np + 1;
 
    // Define the number of required recurrence coefficents
    int ncoeffs = (int)floor(3.0 * (np + 1) / 2);
 
    // Define arrays  for the recurrence coefficients
    // We will use these arrays for the Kronrod results too, hence the
    // reason for the size of the arrays
    double *a = new double[kpoints];
    double *b = new double[kpoints];
 
    // Initialize a and b to zero
    for (i = 0; i < kpoints; i++)
    {
        a[i] = 0.0;
        b[i] = 0.0;
    }
 
    // Call the routine to calculate the recurrence coefficients
    RecCoeff(ncoeffs, a, b, alpha, beta);
 
    // Call the routine to caluclate the jacobi-Kronrod matrix
    JKMatrix(np, a, b);
 
    // Set up the identity matrix
    double **zmatrix = new double *[kpoints];
    for (i = 0; i < kpoints; i++)
    {
        zmatrix[i] = new double[kpoints];
        for (j = 0; j < kpoints; j++)
        {
            zmatrix[i][j] = 0.0;
        }
    }
    for (i = 0; i < kpoints; i++)
    {
        zmatrix[i][i] = 1.0;
    }
 
    // Calculte the points and weights
    TriQL(kpoints, a, b, zmatrix);
 
    for (i = 0; i < kpoints; i++)
    {
        z[i] = a[i];
        w[i] = b[i];
    }
    delete[] a;
    delete[] b;
    for (i = 0; i < kpoints; i++)
    {
        delete[] zmatrix[i];
    }
    delete[] zmatrix;
}
 
/**
\brief  Gauss-Radau-Kronrod-Jacobi zeros and weights.
 
\li Generate \a npt=2*np Radau-Kronrod Jacobi zeros, \a z, and weights,\a w,
associated with the Jacobi polynomial \f$ P^{\alpha,\beta}_{np}(z)
\f$,
*/
void zwrk(double *z, double *w, const int npt, const double alpha,
          const double beta)
{
 
    int np = npt / 2;
 
    if (np < 2)
    {
        fprintf(stderr, "too few points in formula\n");
        return;
    }
 
    double end0 = -1;
 
    int i, j;
 
    // number of kronrod points associated with the np gauss rule
    int kpoints = 2 * np;
 
    // Define the number of required recurrence coefficents
    int ncoeffs = (int)ceil(3.0 * np / 2);
 
    // Define arrays  for the recurrence coefficients
    double *a = new double[ncoeffs + 1];
    double *b = new double[ncoeffs + 1];
 
    // Initialize a and b to zero
    for (i = 0; i < ncoeffs + 1; i++)
    {
        a[i] = 0.0;
        b[i] = 0.0;
    }
 
    // Call the routine to calculate the recurrence coefficients
    RecCoeff(ncoeffs, a, b, alpha, beta);
 
    double *a0 = new double[ncoeffs];
    double *b0 = new double[ncoeffs];
 
    chri1(ncoeffs, a, b, a0, b0, end0);
 
    double s = b0[0] / fabs(b0[0]);
    b0[0]    = s * b0[0];
 
    // Finding the 2*np-1 gauss-kronrod points
    double *z1 = new double[2 * np - 1];
    double *w1 = new double[2 * np - 1];
    for (i = 0; i < ncoeffs; i++)
    {
        z1[i] = a0[i];
        w1[i] = b0[i];
    }
    JKMatrix(np - 1, z1, w1);
    // Set up the identity matrix
    double **zmatrix = new double *[2 * np - 1];
    for (i = 0; i < 2 * np - 1; i++)
    {
        zmatrix[i] = new double[2 * np - 1];
        for (j = 0; j < 2 * np - 1; j++)
        {
            zmatrix[i][j] = 0.0;
        }
    }
    for (i = 0; i < 2 * np - 1; i++)
    {
        zmatrix[i][i] = 1.0;
    }
 
    // Calculate the points and weights
    TriQL(2 * np - 1, z1, w1, zmatrix);
 
    double sumW1 = 0.0;
    for (i = 0; i < 2 * np - 1; i++)
    {
        w1[i] = s * w1[i] / (z1[i] - end0);
        sumW1 += w1[i];
    }
 
    z[0] = end0;
    w[0] = b[0] - sumW1;
    for (i = 1; i < kpoints; i++)
    {
        z[i] = z1[i - 1];
        w[i] = w1[i - 1];
    }
 
    delete[] a;
    delete[] b;
    delete[] a0;
    delete[] b0;
    delete[] z1;
    delete[] w1;
    for (i = 0; i < 2 * np - 1; i++)
    {
        delete[] zmatrix[i];
    }
    delete[] zmatrix;
}
 
/**
\brief  Gauss-Lobatto-Kronrod-Jacobi zeros and weights.
 
\li Generate \a npt=2*np-1 Lobatto-Kronrod Jacobi zeros, \a z, and weights,\a w,
associated with the Jacobi polynomial \f$ P^{\alpha,\beta}_{np}(z)
\f$,
*/
void zwlk(double *z, double *w, const int npt, const double alpha,
          const double beta)
{
 
    int np = (npt + 1) / 2;
 
    if (np < 4)
    {
        fprintf(stderr, "too few points in formula\n");
        return;
    }
 
    double endl = -1;
    double endr = 1;
    int i, j;
 
    // number of kronrod points associated with the np gauss rule
    int kpoints = 2 * np - 1;
 
    // Define the number of required recurrence coefficents
    int ncoeffs = (int)ceil(3.0 * np / 2) - 1;
 
    // Define arrays  for the recurrence coefficients
    double *a = new double[ncoeffs + 1];
    double *b = new double[ncoeffs + 1];
 
    // Initialize a and b to zero
    for (i = 0; i < ncoeffs + 1; i++)
    {
        a[i] = 0.0;
        b[i] = 0.0;
    }
 
    // Call the routine to calculate the recurrence coefficients
    RecCoeff(ncoeffs, a, b, alpha, beta);
 
    double *a0 = new double[ncoeffs];
    double *b0 = new double[ncoeffs];
 
    chri1(ncoeffs, a, b, a0, b0, endl);
 
    double *a1 = new double[ncoeffs - 1];
    double *b1 = new double[ncoeffs - 1];
 
    chri1(ncoeffs - 1, a0, b0, a1, b1, endr);
 
    double s = b1[0] / fabs(b1[0]);
    b1[0]    = s * b1[0];
 
    // Finding the 2*np-1 gauss-kronrod points
    double *z1 = new double[2 * np - 3];
    double *w1 = new double[2 * np - 3];
    for (i = 0; i < ncoeffs; i++)
    {
        z1[i] = a1[i];
        w1[i] = b1[i];
    }
    JKMatrix(np - 2, z1, w1);
    // Set up the identity matrix
    double **zmatrix = new double *[2 * np - 3];
    for (i = 0; i < 2 * np - 3; i++)
    {
        zmatrix[i] = new double[2 * np - 3];
        for (j = 0; j < 2 * np - 3; j++)
        {
            zmatrix[i][j] = 0.0;
        }
    }
    for (i = 0; i < 2 * np - 3; i++)
    {
        zmatrix[i][i] = 1.0;
    }
 
    // Calculate the points and weights
    TriQL(2 * np - 3, z1, w1, zmatrix);
 
    double sumW1   = 0.0;
    double sumW1Z1 = 0.0;
    for (i = 0; i < 2 * np - 3; i++)
    {
        w1[i] = s * w1[i] / (z1[i] - endl) / (z1[i] - endr);
        sumW1 += w1[i];
        sumW1Z1 += z1[i] * w1[i];
    }
 
    double c0 = b[0] - sumW1;
    double c1 = a[0] * b[0] - sumW1Z1;
 
    z[0]          = endl;
    z[2 * np - 2] = endr;
    w[0]          = (c0 * endr - c1) / (endr - endl);
    w[2 * np - 2] = (c1 - c0 * endl) / (endr - endl);
 
    for (i = 1; i < kpoints - 1; i++)
    {
        z[i] = z1[i - 1];
        w[i] = w1[i - 1];
    }
    delete[] a;
    delete[] b;
    delete[] a0;
    delete[] b0;
    delete[] a1;
    delete[] b1;
    delete[] z1;
    delete[] w1;
    for (i = 0; i < 2 * np - 3; i++)
    {
        delete[] zmatrix[i];
    }
    delete[] zmatrix;
}
 
/**
\brief Compute the Integration Matrix.
*/
void Qg(double *Q, const double *z, const int np, const int offset)
{
 
    if (np <= 0)
    {
        Q[0] = 0.0;
    }
    else
    {
        int i, j, k;
        double *pd;
 
        pd = (double *)malloc(np * sizeof(double));
 
        for (i = 0; i < np; i++)
        {
            polycoeffs(i, z, pd, np);
            for (j = 0; j < np; j++)
            {
                Q[j * (np + offset) + i + offset] = 0.0;
                for (k = 0; k < np; k++)
                {
                    Q[j * (np + offset) + i + offset] +=
                        pd[k] / (k + 1) * pow(z[j], k + 1);
                }
            }
        }
        free(pd);
    }
    return;
}
 
/**
\brief Compute the Derivative Matrix and its transpose associated
with the Gauss-Jacobi zeros.
 
\li Compute the derivative matrix, \a d, and its transpose, \a dt,
associated with the n_th order Lagrangian interpolants through the
\a np Gauss-Jacobi points \a z such that \n
\f$  \frac{du}{dz}(z[i]) =  \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$
 
*/
 
void Dgj(double *D, const double *z, const int np, const double alpha,
         const double beta)
{
 
    double one = 1.0, two = 2.0;
 
    if (np <= 0)
    {
        D[0] = 0.0;
    }
    else
    {
        int i, j;
        double *pd;
 
        pd = (double *)malloc(np * sizeof(double));
        jacobd(np, z, pd, np, alpha, beta);
 
        for (i = 0; i < np; i++)
        {
            for (j = 0; j < np; j++)
            {
 
                if (i != j)
                    D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
                else
                    D[i * np + j] =
                        (alpha - beta + (alpha + beta + two) * z[j]) /
                        (two * (one - z[j] * z[j]));
            }
        }
        free(pd);
    }
    return;
}
 
/**
\brief Compute the Derivative Matrix and its transpose associated
with the Gauss-Radau-Jacobi zeros with a zero at \a z=-1.
 
\li Compute the derivative matrix, \a d, associated with the n_th
order Lagrangian interpolants through the \a np Gauss-Radau-Jacobi
points \a z such that \n \f$ \frac{du}{dz}(z[i]) =
\sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$
 
*/
 
void Dgrjm(double *D, const double *z, const int np, const double alpha,
           const double beta)
{
 
    if (np <= 0)
    {
        D[0] = 0.0;
    }
    else
    {
        int i, j;
        double one = 1.0, two = 2.0;
        double *pd;
 
        pd = (double *)malloc(np * sizeof(double));
 
        pd[0] = pow(-one, np - 1) * gammaFracGammaF(np + 1, beta, np, 0.0);
        pd[0] /= gammaF(beta + two);
        jacobd(np - 1, z + 1, pd + 1, np - 1, alpha, beta + 1);
        for (i = 1; i < np; ++i)
            pd[i] *= (1 + z[i]);
 
        for (i = 0; i < np; i++)
        {
            for (j = 0; j < np; j++)
            {
                if (i != j)
                    D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
                else
                {
                    if (j == 0)
                        D[i * np + j] = -(np + alpha + beta + one) *
                                        (np - one) / (two * (beta + two));
                    else
                        D[i * np + j] =
                            (alpha - beta + one + (alpha + beta + one) * z[j]) /
                            (two * (one - z[j] * z[j]));
                }
            }
        }
        free(pd);
    }
 
    return;
}
 
/**
\brief Compute the Derivative Matrix  associated with the
Gauss-Radau-Jacobi zeros with a zero at \a z=1.
 
\li Compute the derivative matrix, \a d, associated with the n_th
order Lagrangian interpolants through the \a np Gauss-Radau-Jacobi
points \a z such that \n \f$ \frac{du}{dz}(z[i]) =
\sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$
*/
 
void Dgrjp(double *D, const double *z, const int np, const double alpha,
           const double beta)
{
 
    if (np <= 0)
    {
        D[0] = 0.0;
    }
    else
    {
        int i, j;
        double one = 1.0, two = 2.0;
        double *pd;
 
        pd = (double *)malloc(np * sizeof(double));
 
        jacobd(np - 1, z, pd, np - 1, alpha + 1, beta);
        for (i = 0; i < np - 1; ++i)
            pd[i] *= (1 - z[i]);
        pd[np - 1] = -gammaFracGammaF(np + 1, alpha, np, 0.0);
        pd[np - 1] /= gammaF(alpha + two);
 
        for (i = 0; i < np; i++)
        {
            for (j = 0; j < np; j++)
            {
                if (i != j)
                    D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
                else
                {
                    if (j == np - 1)
                        D[i * np + j] = (np + alpha + beta + one) * (np - one) /
                                        (two * (alpha + two));
                    else
                        D[i * np + j] =
                            (alpha - beta - one + (alpha + beta + one) * z[j]) /
                            (two * (one - z[j] * z[j]));
                }
            }
        }
        free(pd);
    }
 
    return;
}
 
/**
\brief Compute the Derivative Matrix associated with the
Gauss-Lobatto-Jacobi zeros.
 
\li Compute the derivative matrix, \a d, associated with the n_th
order Lagrange interpolants through the \a np
Gauss-Lobatto-Jacobi points \a z such that \n \f$
\frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$
 
*/
 
void Dglj(double *D, const double *z, const int np, const double alpha,
          const double beta)
{
    if (np <= 1)
    {
        D[0] = 0.0;
    }
    else
    {
        int i, j;
        double one = 1.0, two = 2.0;
        double *pd;
 
        pd = (double *)malloc(np * sizeof(double));
 
        pd[0] = two * pow(-one, np) * gammaFracGammaF(np, beta, np - 1, 0.0);
        pd[0] /= gammaF(beta + two);
        jacobd(np - 2, z + 1, pd + 1, np - 2, alpha + 1, beta + 1);
        for (i = 1; i < np - 1; ++i)
            pd[i] *= (one - z[i] * z[i]);
        pd[np - 1] = -two * gammaFracGammaF(np, alpha, np - 1, 0.0);
        pd[np - 1] /= gammaF(alpha + two);
 
        for (i = 0; i < np; i++)
        {
            for (j = 0; j < np; j++)
            {
                if (i != j)
                    D[i * np + j] = pd[j] / (pd[i] * (z[j] - z[i]));
                else
                {
                    if (j == 0)
                        D[i * np + j] =
                            (alpha - (np - 1) * (np + alpha + beta)) /
                            (two * (beta + two));
                    else if (j == np - 1)
                        D[i * np + j] =
                            -(beta - (np - 1) * (np + alpha + beta)) /
                            (two * (alpha + two));
                    else
                        D[i * np + j] = (alpha - beta + (alpha + beta) * z[j]) /
                                        (two * (one - z[j] * z[j]));
                }
            }
        }
        free(pd);
    }
 
    return;
}
 
/**
\brief Compute the value of the \a i th Lagrangian interpolant through
the \a np Gauss-Jacobi points \a zgj at the arbitrary location \a z.
 
\li \f$ -1 \leq z \leq 1 \f$
 
\li Uses the defintion of the Lagrangian interpolant:\n
%
\f$ \begin{array}{rcl}
h_j(z) =  \left\{ \begin{array}{ll}
\displaystyle \frac{P_{np}^{\alpha,\beta}(z)}
{[P_{np}^{\alpha,\beta}(z_j)]^\prime
(z-z_j)} & \mbox{if $z \ne z_j$}\\
& \\
1 & \mbox{if $z=z_j$}
\end{array}
\right.
\end{array}   \f$
*/
 
double hgj(const int i, const double z, const double *zgj, const int np,
           const double alpha, const double beta)
{
    boost::ignore_unused(alpha, beta);
    double zi, dz;
 
    zi = *(zgj + i);
    dz = z - zi;
    if (fabs(dz) < EPS)
        return 1.0;
 
    return laginterp(z, i, zgj, np);
}
 
/**
\brief Compute the value of the \a i th Lagrangian interpolant through the
\a np Gauss-Radau-Jacobi points \a zgrj at the arbitrary location
\a z. This routine assumes \a zgrj includes the point \a -1.
 
\li \f$ -1 \leq z \leq 1 \f$
 
\li Uses the defintion of the Lagrangian interpolant:\n
%
\f$ \begin{array}{rcl}
h_j(z) = \left\{ \begin{array}{ll}
\displaystyle \frac{(1+z) P_{np-1}^{\alpha,\beta+1}(z)}
{((1+z_j) [P_{np-1}^{\alpha,\beta+1}(z_j)]^\prime +
P_{np-1}^{\alpha,\beta+1}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\
& \\
1 & \mbox{if $z=z_j$}
\end{array}
\right.
\end{array}   \f$
*/
 
double hgrjm(const int i, const double z, const double *zgrj, const int np,
             const double alpha, const double beta)
{
    boost::ignore_unused(alpha, beta);
 
    double zi, dz;
 
    zi = *(zgrj + i);
    dz = z - zi;
    if (fabs(dz) < EPS)
        return 1.0;
 
    return laginterp(z, i, zgrj, np);
}
 
/**
\brief Compute the value of the \a i th Lagrangian interpolant through the
\a np Gauss-Radau-Jacobi points \a zgrj at the arbitrary location
\a z. This routine assumes \a zgrj includes the point \a +1.
 
\li \f$ -1 \leq z \leq 1 \f$
 
\li Uses the defintion of the Lagrangian interpolant:\n
%
\f$ \begin{array}{rcl}
h_j(z) = \left\{ \begin{array}{ll}
\displaystyle \frac{(1-z) P_{np-1}^{\alpha+1,\beta}(z)}
{((1-z_j) [P_{np-1}^{\alpha+1,\beta}(z_j)]^\prime -
P_{np-1}^{\alpha+1,\beta}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\
& \\
1 & \mbox{if $z=z_j$}
\end{array}
\right.
\end{array}   \f$
*/
 
double hgrjp(const int i, const double z, const double *zgrj, const int np,
             const double alpha, const double beta)
{
    boost::ignore_unused(alpha, beta);
    double zi, dz;
 
    zi = *(zgrj + i);
    dz = z - zi;
    if (fabs(dz) < EPS)
        return 1.0;
 
    return laginterp(z, i, zgrj, np);
}
 
/**
\brief Compute the value of the \a i th Lagrangian interpolant through the
\a np Gauss-Lobatto-Jacobi points \a zgrj at the arbitrary location
\a z.
 
\li \f$ -1 \leq z \leq 1 \f$
 
\li Uses the defintion of the Lagrangian interpolant:\n
%
\f$ \begin{array}{rcl}
h_j(z) = \left\{ \begin{array}{ll}
\displaystyle \frac{(1-z^2) P_{np-2}^{\alpha+1,\beta+1}(z)}
{((1-z^2_j) [P_{np-2}^{\alpha+1,\beta+1}(z_j)]^\prime -
2 z_j P_{np-2}^{\alpha+1,\beta+1}(z_j) ) (z-z_j)}&\mbox{if $z \ne z_j$}\\
& \\
1 & \mbox{if $z=z_j$}
\end{array}
\right.
\end{array}   \f$
*/
 
double hglj(const int i, const double z, const double *zglj, const int np,
            const double alpha, const double beta)
{
    boost::ignore_unused(alpha, beta);
 
    double zi, dz;
 
    zi = *(zglj + i);
    dz = z - zi;
    if (fabs(dz) < EPS)
        return 1.0;
 
    return laginterp(z, i, zglj, np);
}
 
/**
\brief Interpolation Operator from Gauss-Jacobi points to an
arbitrary distribution at points \a zm
 
\li Computes the one-dimensional interpolation matrix, \a im, to
interpolate a function from at Gauss-Jacobi distribution of \a nz
zeros \a zgrj to an arbitrary distribution of \a mz points \a zm, i.e.\n
\f$
u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j])
\f$
 
*/
 
void Imgj(double *im, const double *zgj, const double *zm, const int nz,
          const int mz, const double alpha, const double beta)
{
    double zp;
    int i, j;
 
    for (i = 0; i < nz; ++i)
    {
        for (j = 0; j < mz; ++j)
        {
            zp             = zm[j];
            im[i * mz + j] = hgj(i, zp, zgj, nz, alpha, beta);
        }
    }
 
    return;
}
 
/**
\brief Interpolation Operator from Gauss-Radau-Jacobi points
(including \a z=-1) to an arbitrary distrubtion at points \a zm
 
\li Computes the one-dimensional interpolation matrix, \a im, to
interpolate a function from at Gauss-Radau-Jacobi distribution of
\a nz zeros \a zgrj (where \a zgrj[0]=-1) to an arbitrary
distribution of \a mz points \a zm, i.e.
\n
\f$ u(zm[i]) =    \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$
 
*/
 
void Imgrjm(double *im, const double *zgrj, const double *zm, const int nz,
            const int mz, const double alpha, const double beta)
{
    double zp;
    int i, j;
 
    for (i = 0; i < nz; i++)
    {
        for (j = 0; j < mz; j++)
        {
            zp             = zm[j];
            im[i * mz + j] = hgrjm(i, zp, zgrj, nz, alpha, beta);
        }
    }
 
    return;
}
 
/**
\brief Interpolation Operator from Gauss-Radau-Jacobi points
(including \a z=1) to an arbitrary distrubtion at points \a zm
 
\li Computes the one-dimensional interpolation matrix, \a im, to
interpolate a function from at Gauss-Radau-Jacobi distribution of
\a nz zeros \a zgrj (where \a zgrj[nz-1]=1) to an arbitrary
distribution of \a mz points \a zm, i.e.
\n
\f$ u(zm[i]) =    \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$
 
*/
 
void Imgrjp(double *im, const double *zgrj, const double *zm, const int nz,
            const int mz, const double alpha, const double beta)
{
    double zp;
    int i, j;
 
    for (i = 0; i < nz; i++)
    {
        for (j = 0; j < mz; j++)
        {
            zp             = zm[j];
            im[i * mz + j] = hgrjp(i, zp, zgrj, nz, alpha, beta);
        }
    }
 
    return;
}
 
/**
\brief Interpolation Operator from Gauss-Lobatto-Jacobi points
to an arbitrary distrubtion at points \a zm
 
\li Computes the one-dimensional interpolation matrix, \a im, to
interpolate a function from at Gauss-Lobatto-Jacobi distribution of
\a nz zeros \a zgrj (where \a zgrj[0]=-1) to an arbitrary
distribution of \a mz points \a zm, i.e.
\n
\f$ u(zm[i]) =    \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$
 
*/
 
void Imglj(double *im, const double *zglj, const double *zm, const int nz,
           const int mz, const double alpha, const double beta)
{
    double zp;
    int i, j;
 
    for (i = 0; i < nz; i++)
    {
        for (j = 0; j < mz; j++)
        {
            zp             = zm[j];
            im[i * mz + j] = hglj(i, zp, zglj, nz, alpha, beta);
        }
    }
 
    return;
}
 
void polycoeffs(const int i, const double *z, double *c, const int np)
{
    int j, k, m;
 
    // Compute denominator
    double d = 1.0;
    for (j = 0; j < np; j++)
    {
        if (i != j)
        {
            d *= z[i] - z[j];
        }
        c[j] = 0.0;
    }
 
    // Compute coefficient
    c[0] = 1.0;
    m    = 0;
    for (j = 0; j < np; j++)
    {
        if (i != j)
        {
            m += 1;
            c[m] = c[m - 1];
            for (k = m - 1; k > 0; k--)
            {
                c[k] *= -1.0 * z[j];
                c[k] += c[k - 1];
            }
            c[0] *= -1.0 * z[j];
        }
    }
    for (j = 0; j < np; j++)
    {
        c[j] /= d;
    }
}
 
/**
\brief Routine to calculate Jacobi polynomials, \f$
P^{\alpha,\beta}_n(z) \f$, and their first derivative, \f$
\frac{d}{dz} P^{\alpha,\beta}_n(z) \f$.
 
\li This function returns the vectors \a poly_in and \a poly_d
containing the value of the \f$ n^th \f$ order Jacobi polynomial
\f$ P^{\alpha,\beta}_n(z) \alpha > -1, \beta > -1 \f$ and its
derivative at the \a np points in \a z[i]
 
- If \a poly_in = NULL then only calculate derivatice
 
- If \a polyd   = NULL then only calculate polynomial
 
- To calculate the polynomial this routine uses the recursion
relationship (see appendix A ref [4]) :
\f$ \begin{array}{rcl}
P^{\alpha,\beta}_0(z) &=& 1 \\
P^{\alpha,\beta}_1(z) &=& \frac{1}{2} [ \alpha-\beta+(\alpha+\beta+2)z] \\
a^1_n P^{\alpha,\beta}_{n+1}(z) &=& (a^2_n + a^3_n z)
P^{\alpha,\beta}_n(z) - a^4_n P^{\alpha,\beta}_{n-1}(z) \\
a^1_n &=& 2(n+1)(n+\alpha + \beta + 1)(2n + \alpha + \beta) \\
a^2_n &=& (2n + \alpha + \beta + 1)(\alpha^2 - \beta^2)  \\
a^3_n &=& (2n + \alpha + \beta)(2n + \alpha + \beta + 1)
(2n + \alpha + \beta + 2)  \\
a^4_n &=& 2(n+\alpha)(n+\beta)(2n + \alpha + \beta + 2)
\end{array} \f$
 
- To calculate the derivative of the polynomial this routine uses
the relationship (see appendix A ref [4]) :
\f$ \begin{array}{rcl}
b^1_n(z)\frac{d}{dz} P^{\alpha,\beta}_n(z)&=&b^2_n(z)P^{\alpha,\beta}_n(z)
+ b^3_n(z) P^{\alpha,\beta}_{n-1}(z) \hspace{2.2cm} \\
b^1_n(z) &=& (2n+\alpha + \beta)(1-z^2) \\
b^2_n(z) &=& n[\alpha - \beta - (2n+\alpha + \beta)z]\\
b^3_n(z) &=& 2(n+\alpha)(n+\beta)
\end{array} \f$
 
- Note the derivative from this routine is only valid for -1 < \a z < 1.
*/
void jacobfd(const int np, const double *z, double *poly_in, double *polyd,
             const int n, const double alpha, const double beta)
{
    int i;
    double zero = 0.0, one = 1.0, two = 2.0;
 
    if (!np)
        return;
 
    if (n == 0)
    {
        if (poly_in)
            for (i = 0; i < np; ++i)
                poly_in[i] = one;
        if (polyd)
            for (i = 0; i < np; ++i)
                polyd[i] = zero;
    }
    else if (n == 1)
    {
        if (poly_in)
            for (i = 0; i < np; ++i)
                poly_in[i] = 0.5 * (alpha - beta + (alpha + beta + two) * z[i]);
        if (polyd)
            for (i = 0; i < np; ++i)
                polyd[i] = 0.5 * (alpha + beta + two);
    }
    else
    {
        int k;
        double a1, a2, a3, a4;
        double two = 2.0, apb = alpha + beta;
        double *poly, *polyn1, *polyn2;
 
        if (poly_in)
        { // switch for case of no poynomial function return
            polyn1 = (double *)malloc(2 * np * sizeof(double));
            polyn2 = polyn1 + np;
            poly   = poly_in;
        }
        else
        {
            polyn1 = (double *)malloc(3 * np * sizeof(double));
            polyn2 = polyn1 + np;
            poly   = polyn2 + np;
        }
 
        for (i = 0; i < np; ++i)
        {
            polyn2[i] = one;
            polyn1[i] = 0.5 * (alpha - beta + (alpha + beta + two) * z[i]);
        }
 
        for (k = 2; k <= n; ++k)
        {
            a1 = two * k * (k + apb) * (two * k + apb - two);
            a2 = (two * k + apb - one) * (alpha * alpha - beta * beta);
            a3 =
                (two * k + apb - two) * (two * k + apb - one) * (two * k + apb);
            a4 = two * (k + alpha - one) * (k + beta - one) * (two * k + apb);
 
            a2 /= a1;
            a3 /= a1;
            a4 /= a1;
 
            for (i = 0; i < np; ++i)
            {
                poly[i]   = (a2 + a3 * z[i]) * polyn1[i] - a4 * polyn2[i];
                polyn2[i] = polyn1[i];
                polyn1[i] = poly[i];
            }
        }
 
        if (polyd)
        {
            a1 = n * (alpha - beta);
            a2 = n * (two * n + alpha + beta);
            a3 = two * (n + alpha) * (n + beta);
            a4 = (two * n + alpha + beta);
            a1 /= a4;
            a2 /= a4;
            a3 /= a4;
 
            // note polyn2 points to polyn1 at end of poly iterations
            for (i = 0; i < np; ++i)
            {
                polyd[i] = (a1 - a2 * z[i]) * poly[i] + a3 * polyn2[i];
                polyd[i] /= (one - z[i] * z[i]);
            }
        }
 
        free(polyn1);
    }
 
    return;
}
 
/**
\brief Calculate the  derivative of Jacobi polynomials
 
\li Generates a vector \a poly of values of the derivative of the
\a n th order Jacobi polynomial \f$ P^(\alpha,\beta)_n(z)\f$ at the
\a np points \a z.
 
\li To do this we have used the relation
\n
\f$ \frac{d}{dz} P^{\alpha,\beta}_n(z)
= \frac{1}{2} (\alpha + \beta + n + 1)  P^{\alpha,\beta}_n(z) \f$
 
\li This formulation is valid for \f$ -1 \leq z \leq 1 \f$
 
*/
void jacobd(const int np, const double *z, double *polyd, const int n,
            const double alpha, const double beta)
{
    int i;
    double one = 1.0;
    if (n == 0)
        for (i = 0; i < np; ++i)
            polyd[i] = 0.0;
    else
    {
        // jacobf(np,z,polyd,n-1,alpha+one,beta+one);
        jacobfd(np, z, polyd, NULL, n - 1, alpha + one, beta + one);
        for (i = 0; i < np; ++i)
            polyd[i] *= 0.5 * (alpha + beta + (double)n + one);
    }
    return;
}
 
/**
\brief Calculate the Gamma function , \f$ \Gamma(n)\f$, for integer
values and halves.
 
Determine the value of \f$\Gamma(n)\f$ using:
 
\f$ \Gamma(n) = (n-1)!  \mbox{ or  }  \Gamma(n+1/2) = (n-1/2)\Gamma(n-1/2)\f$
 
where \f$ \Gamma(1/2) = \sqrt(\pi)\f$
*/
 
double gammaF(const double x)
{
    double gamma = 1.0;
 
    if (x == -0.5)
        gamma = -2.0 * sqrt(M_PI);
    else if (!x)
        return gamma;
    else if ((x - (int)x) == 0.5)
    {
        int n      = (int)x;
        double tmp = x;
 
        gamma = sqrt(M_PI);
        while (n--)
        {
            tmp -= 1.0;
            gamma *= tmp;
        }
    }
    else if ((x - (int)x) == 0.0)
    {
        int n      = (int)x;
        double tmp = x;
 
        while (--n)
        {
            tmp -= 1.0;
            gamma *= tmp;
        }
    }
    else
        fprintf(stderr, "%lf is not of integer or half order\n", x);
    return gamma;
}
 
/**
\brief Calculate fraction of two Gamma functions, \f$
\Gamma(x+\alpha)/\Gamma(y+\beta) \f$, for integer values and halves.
 
Determine the value of \f$\Gamma(n)\f$ using:
 
\f$ \Gamma(n) = (n-1)!  \mbox{ or  }  \Gamma(n+1/2) = (n-1/2)\Gamma(n-1/2)\f$
 
where \f$ \Gamma(1/2) = \sqrt(\pi)\f$
 
Attempts simplification in cases like \f$ \Gamma(x+1)/\Gamma(x-1) = x*(x-1) \f$
*/
 
double gammaFracGammaF(const int x, const double alpha, const int y,
                       const double beta)
{
    double gamma = 1.0;
    double halfa = fabs(alpha - int(alpha));
    double halfb = fabs(beta - int(beta));
    if (halfa == 0.0 && halfb == 0.0) // integer value
    {
        int X = x + alpha;
        int Y = y + beta;
        if (X > Y)
        {
            for (int tmp = X - 1; tmp > Y - 1; tmp -= 1)
                gamma *= tmp;
        }
        else if (Y > X)
        {
            for (int tmp = Y - 1; tmp > X - 1; tmp -= 1)
                gamma *= tmp;
            gamma = 1. / gamma;
        }
    }
    else if (halfa == 0.5 && halfb == 0.5) // both are halves
    {
        double X = x + alpha;
        double Y = y + beta;
        if (X > Y)
        {
            for (int tmp = int(X); tmp > int(Y); tmp -= 1)
                gamma *= tmp - 0.5;
        }
        else if (Y > X)
        {
            for (int tmp = int(Y); tmp > int(X); tmp -= 1)
                gamma *= tmp - 0.5;
            gamma = 1. / gamma;
        }
    }
    else
    {
        double X = x + alpha;
        double Y = y + beta;
        while (X > 1 || Y > 1)
        {
            if (X > 1)
            {
                gamma *= X - 1.;
                X -= 1.;
            }
            if (Y > 1)
            {
                gamma /= Y - 1.;
                Y -= 1.;
            }
        }
        if (X == 0.5)
        {
            gamma *= sqrt(M_PI);
        }
        else if (Y == 0.5)
        {
            gamma /= sqrt(M_PI);
        }
        else
        {
            fprintf(stderr, "%lf or %lf is not of integer or half order\n", X,
                    Y);
        }
    }
 
    return gamma;
}
 
/**
\brief  Calculate the \a n zeros, \a z, of the Jacobi polynomial, i.e.
\f$ P_n^{\alpha,\beta}(z) = 0 \f$
 
This routine is only value for \f$( \alpha > -1, \beta > -1)\f$
and uses polynomial deflation in a Newton iteration
*/
 
static void Jacobz(const int n, double *z, const double alpha,
                   const double beta)
{
    int i, j, k;
    double dth = M_PI / (2.0 * (double)n);
    double poly, pder, rlast = 0.0;
    double sum, delr, r;
    double one = 1.0, two = 2.0;
 
    if (!n)
        return;
 
    for (k = 0; k < n; ++k)
    {
        r = -cos((two * (double)k + one) * dth);
        if (k)
            r = 0.5 * (r + rlast);
 
        for (j = 1; j < STOP; ++j)
        {
            jacobfd(1, &r, &poly, &pder, n, alpha, beta);
 
            for (i = 0, sum = 0.0; i < k; ++i)
                sum += one / (r - z[i]);
 
            delr = -poly / (pder - sum * poly);
            r += delr;
            if (fabs(delr) < EPS)
                break;
        }
        z[k]  = r;
        rlast = r;
    }
    return;
}
 
/**
\brief Zero and Weight determination through the eigenvalues and eigenvectors of
a tridiagonal matrix from the three term recurrence relationship.
 
Set up a symmetric tridiagonal matrix
 
\f$ \left [  \begin{array}{ccccc}
a[0] & b[0]   &        &        & \\
b[0] & a[1]   & b[1]   &        & \\
0   & \ddots & \ddots & \ddots &  \\
&        & \ddots & \ddots & b[n-2] \\
&        &        & b[n-2] & a[n-1] \end{array} \right ] \f$
 
Where the coefficients a[n], b[n] come from the  recurrence relation
 
\f$  b_j p_j(z) = (z - a_j ) p_{j-1}(z) - b_{j-1}   p_{j-2}(z) \f$
 
where \f$ j=n+1\f$ and \f$p_j(z)\f$ are the Jacobi (normalized)
orthogonal polynomials \f$ \alpha,\beta > -1\f$( integer values and
halves). Since the polynomials are orthonormalized, the tridiagonal
matrix is guaranteed to be symmetric. The eigenvalues of this
matrix are the zeros of the Jacobi polynomial.
*/
 
void JacZeros(const int n, double *a, double *b, const double alpha,
              const double beta)
{
 
    int i, j;
    RecCoeff(n, a, b, alpha, beta);
 
    double **z = new double *[n];
    for (i = 0; i < n; i++)
    {
        z[i] = new double[n];
        for (j = 0; j < n; j++)
        {
            z[i][j] = 0.0;
        }
    }
    for (i = 0; i < n; i++)
    {
        z[i][i] = 1.0;
    }
 
    // find eigenvalues and eigenvectors
    TriQL(n, a, b, z);
 
    delete[] z;
    return;
}
 
/**
\brief  The routine finds the recurrence coefficients \a a and
\a b of the orthogonal polynomials
*/
static void RecCoeff(const int n, double *a, double *b, const double alpha,
                     const double beta)
{
 
    int i;
    double apb, apbi, a2b2;
 
    if (!n)
        return;
 
    // generate normalised terms
    apb  = alpha + beta;
    apbi = 2.0 + apb;
 
    b[0] = pow(2.0, apb + 1.0) * gammaF(alpha + 1.0) * gammaF(beta + 1.0) /
           gammaF(apbi); // MuZero
    a[0] = (beta - alpha) / apbi;
    b[1] = (4.0 * (1.0 + alpha) * (1.0 + beta) / ((apbi + 1.0) * apbi * apbi));
 
    a2b2 = beta * beta - alpha * alpha;
 
    for (i = 1; i < n - 1; i++)
    {
        apbi     = 2.0 * (i + 1) + apb;
        a[i]     = a2b2 / ((apbi - 2.0) * apbi);
        b[i + 1] = (4.0 * (i + 1) * (i + 1 + alpha) * (i + 1 + beta) *
                    (i + 1 + apb) / ((apbi * apbi - 1) * apbi * apbi));
    }
 
    apbi     = 2.0 * n + apb;
    a[n - 1] = a2b2 / ((apbi - 2.0) * apbi);
}
 
/** \brief QL algorithm for symmetric tridiagonal matrix
 
This subroutine is a translation of an algol procedure,
num. math. \b 12, 377-383(1968) by martin and wilkinson, as modified
in num. math. \b 15, 450(1970) by dubrulle.  Handbook for
auto. comp., vol.ii-linear algebra, 241-248(1971).  This is a
modified version from numerical recipes.
 
This subroutine finds the eigenvalues and first components of the
eigenvectors of a symmetric tridiagonal matrix by the implicit QL
method.
 
on input:
- n is the order of the matrix;
- d contains the diagonal elements of the input matrix;
- e contains the subdiagonal elements of the input matrix
in its first n-2 positions.
    - z is the n by n identity matrix
 
on output:
 
- d contains the eigenvalues in ascending order.
- e contains the weight values - modifications of the first component
      of normalised eigenvectors
*/
 
static void TriQL(const int n, double *d, double *e, double **z)
{
    int m, l, iter, i, k;
    double s, r, p, g, f, dd, c, b;
 
    double MuZero = e[0];
 
    // Renumber the elements of e
    for (i = 0; i < n - 1; i++)
    {
        e[i] = sqrt(e[i + 1]);
    }
    e[n - 1] = 0.0;
 
    for (l = 0; l < n; l++)
    {
        iter = 0;
        do
        {
            for (m = l; m < n - 1; m++)
            {
                dd = fabs(d[m]) + fabs(d[m + 1]);
                if (fabs(e[m]) + dd == dd)
                    break;
            }
            if (m != l)
            {
                if (iter++ == STOP)
                {
                    fprintf(stderr, "triQL: Too many iterations in TQLI");
                    exit(1);
                }
                g = (d[l + 1] - d[l]) / (2.0 * e[l]);
                r = sqrt((g * g) + 1.0);
                g = d[m] - d[l] + e[l] / (g + sign(r, g));
                s = c = 1.0;
                p     = 0.0;
                for (i = m - 1; i >= l; i--)
                {
                    f = s * e[i];
                    b = c * e[i];
                    if (fabs(f) >= fabs(g))
                    {
                        c        = g / f;
                        r        = sqrt((c * c) + 1.0);
                        e[i + 1] = f * r;
                        c *= (s = 1.0 / r);
                    }
                    else
                    {
                        s        = f / g;
                        r        = sqrt((s * s) + 1.0);
                        e[i + 1] = g * r;
                        s *= (c = 1.0 / r);
                    }
                    g        = d[i + 1] - p;
                    r        = (d[i] - g) * s + 2.0 * c * b;
                    p        = s * r;
                    d[i + 1] = g + p;
                    g        = c * r - b;
 
                    // Calculate the eigenvectors
                    for (k = 0; k < n; k++)
                    {
                        f           = z[k][i + 1];
                        z[k][i + 1] = s * z[k][i] + c * f;
                        z[k][i]     = c * z[k][i] - s * f;
                    }
                }
                d[l] = d[l] - p;
                e[l] = g;
                e[m] = 0.0;
            }
        } while (m != l);
    }
 
    // order eigenvalues
    // Since we only need the first component of the eigenvectors
    // to calcualte the weight, we only swap the first components
    for (i = 0; i < n - 1; ++i)
    {
        k = i;
        p = d[i];
        for (l = i + 1; l < n; ++l)
            if (d[l] < p)
            {
                k = l;
                p = d[l];
            }
        d[k] = d[i];
        d[i] = p;
 
        double temp = z[0][k];
        z[0][k]     = z[0][i];
        z[0][i]     = temp;
    }
 
    // Calculate the weights
    for (i = 0; i < n; i++)
    {
        e[i] = MuZero * z[0][i] * z[0][i];
    }
}
 
/**
\brief Calcualtes the Jacobi-kronrod matrix by determining the
\a a and \b coefficients.
 
The first \a 3n+1 coefficients are already known
 
For more information refer to:
"Dirk P. Laurie, Calcualtion of Gauss-Kronrod quadrature rules"
*/
void JKMatrix(int n, double *a, double *b)
{
    int i, j, k, m;
    // Working storage
    int size  = (int)floor(n / 2.0) + 2;
    double *s = new double[size];
    double *t = new double[size];
 
    // Initialize s and t to zero
    for (i = 0; i < size; i++)
    {
        s[i] = 0.0;
        t[i] = 0.0;
    }
 
    t[1] = b[n + 1];
    for (m = 0; m <= n - 2; m++)
    {
        double u = 0.0;
        for (k = (int)floor((m + 1) / 2.0); k >= 0; k--)
        {
            int l = m - k;
            u     = u + (a[k + n + 1] - a[l]) * t[k + 1] + b[k + n + 1] * s[k] -
                b[l] * s[k + 1];
            s[k + 1] = u;
        }
 
        // Swap the contents of s and t
        double *hold = s;
        s            = t;
        t            = hold;
    }
 
    for (j = (int)floor(n / 2.0); j >= 0; j--)
    {
        s[j + 1] = s[j];
    }
 
    for (m = n - 1; m <= 2 * n - 3; m++)
    {
        double u = 0;
        for (k = m + 1 - n; k <= floor((m - 1) / 2.0); k++)
        {
            int l = m - k;
            j     = n - 1 - l;
            u = u - (a[k + n + 1] - a[l]) * t[j + 1] - b[k + n + 1] * s[j + 1] +
                b[l] * s[j + 2];
            s[j + 1] = u;
        }
 
        if (m % 2 == 0)
        {
            k = m / 2;
            a[k + n + 1] =
                a[k] + (s[j + 1] - b[k + n + 1] * s[j + 2]) / t[j + 2];
        }
        else
        {
            k            = (m + 1) / 2;
            b[k + n + 1] = s[j + 1] / s[j + 2];
        }
 
        // Swap the contents of s and t
        double *hold = s;
        s            = t;
        t            = hold;
    }
 
    a[2 * n] = a[n - 1] - b[2 * n] * s[1] / t[1];
}
 
/**
\brief
Given a weight function \f$w(t)\f$ through the first \a n+1
coefficients \a a and \a b of its orthogonal polynomials
this routine generates the first \a n recurrence coefficients for the orthogonal
polynomials relative to the modified weight function \f$(t-z)w(t)\f$.
 
The result will be placed in the array \a a0 and \a b0.
*/
 
void chri1(int n, double *a, double *b, double *a0, double *b0, double z)
{
 
    double q = ceil(3.0 * n / 2);
    int size = (int)q + 1;
    if (size < n + 1)
    {
        fprintf(stderr, "input arrays a and b are too short\n");
    }
    double *r = new double[n + 1];
    r[0]      = z - a[0];
    r[1]      = z - a[1] - b[1] / r[0];
    a0[0]     = a[1] + r[1] - r[0];
    b0[0]     = -r[0] * b[0];
 
    if (n == 1)
    {
        delete[] r;
        return;
    }
    int k = 0;
    for (k = 1; k < n; k++)
    {
        r[k + 1] = z - a[k + 1] - b[k + 1] / r[k];
        a0[k]    = a[k + 1] + r[k + 1] - r[k];
        b0[k]    = b[k] * r[k] / r[k - 1];
    }
    delete[] r;
}
 
/**
 
    \brief
 
Calcualte the bessel function of the first kind with complex double input y.
Taken from Numerical Recipies in C
 
Returns a complex double
*/
 
std::complex<Nektar::NekDouble> ImagBesselComp(
    int n, std::complex<Nektar::NekDouble> y)
{
    std::complex<Nektar::NekDouble> z(1.0, 0.0);
    std::complex<Nektar::NekDouble> zbes(1.0, 0.0);
    std::complex<Nektar::NekDouble> zarg;
    Nektar::NekDouble tol = 1e-15;
    int maxit             = 10000;
    int i                 = 1;
 
    zarg = -0.25 * y * y;
 
    while (abs(z) > tol && i <= maxit)
    {
        z = z * (1.0 / i / (i + n) * zarg);
        if (abs(z) <= tol)
            break;
        zbes = zbes + z;
        i++;
    }
    zarg = 0.5 * y;
    for (i = 1; i <= n; i++)
    {
        zbes = zbes * zarg;
    }
    return zbes;
}
} // namespace Polylib