23 #define EPS 100*DBL_EPSILON
27 #define sign(a,b) ((b)<0 ? -fabs(a) : fabs(a))
41 int m_digits =
static_cast<int>(fabs(floor(log10(DBL_EPSILON)))-1);
43 if (fabs(xl-xr)<1.e-4){
45 m_expn =
static_cast<int>(floor(log10(fabs(xl-xr))));
46 m_xln = xl*powl(10.0
L,-m_expn)-floor(xl*powl(10.0
L,-m_expn));
47 m_xrn = xr*powl(10.0
L,-m_expn)-floor(xl*powl(10.0
L,-m_expn));
48 m_xln = round(m_xln*powl(10.0
L,m_digits+m_expn));
49 m_xrn = round(m_xrn*powl(10.0
L,m_digits+m_expn));
51 return powl(10.0
L,-m_digits)*(m_xln-m_xrn);
57 double laginterp(
double z,
int j,
const double *zj,
int np)
60 for (
int i=0; i<np; i++)
64 temp *=
optdiff(z,zj[i])/(zj[j]-zj[i]);
72 #define POLYNOMIAL_DEFLATION 0
76 #ifdef POLYNOMIAL_DEFLATION
80 #define jacobz(n,z,alpha,beta) Jacobz(n,z,alpha,beta)
86 #define jacobz(n,z,alpha,beta) JacZeros(n,z,alpha,beta)
98 static void Jacobz (
const int n,
double *z,
const double alpha,
108 static void TriQL(
const int,
double *,
double *,
double **);
110 double gammaF (
const double);
112 static void RecCoeff(
const int,
double *,
double *,
const double,
116 void JKMatrix(
int,
double *,
double *);
118 void chri1(
int,
double*,
double*,
double*,
double*,
double);
142 void zwgj (
double *z,
double *w,
const int np,
const double alpha,
150 double fac, one = 1.0, two = 2.0, apb = alpha + beta;
156 jacobd (np,z,w,np,alpha,beta);
160 fac = pow(two,apb + one)*
gammaF(alpha + np + one)*
gammaF(beta + np + one);
166 for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));
198 void zwgrjm(
double *z,
double *w,
const int np,
const double alpha,
218 double fac, one = 1.0, two = 2.0, apb = alpha + beta;
224 jacobz (np-1,z+1,alpha,beta+1);
226 jacobfd (np,z,w,NULL,np-1,alpha,beta);
230 fac = pow(two,apb)*
gammaF(alpha + np)*
gammaF(beta + np);
236 for(i = 0; i < np; ++i) w[i] = fac*(1-z[i])/(w[i]*w[i]);
238 w[0] *= (beta + one);
274 void zwgrjp(
double *z,
double *w,
const int np,
const double alpha,
294 double fac, one = 1.0, two = 2.0, apb = alpha + beta;
298 jacobz (np-1,z,alpha+1,beta);
302 jacobfd (np,z,w,NULL,np-1,alpha,beta);
306 fac = pow(two,apb)*
gammaF(alpha + np)*
gammaF(beta + np);
312 for(i = 0; i < np; ++i) w[i] = fac*(1+z[i])/(w[i]*w[i]);
314 w[np-1] *= (alpha + one);
346 void zwglj(
double *z,
double *w,
const int np,
const double alpha,
366 double fac, one = 1.0, apb = alpha + beta, two = 2.0;
374 jacobz (np-2,z + 1,alpha + one,beta + one);
376 jacobfd (np,z,w,NULL,np-1,alpha,beta);
380 fac = pow(two,apb + 1)*
gammaF(alpha + np)*
gammaF(beta + np);
382 fac /= (np-1)*
gammaF(np)*
gammaF(alpha + beta + np + one);
386 for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]);
388 w[0] *= (beta + one);
390 w[np-1] *= (alpha + one);
420 void zwgk(
double *z,
double *w,
const int npt ,
const double alpha,
438 int kpoints = 2*np + 1;
444 int ncoeffs = (int)floor(3.0*(np+1)/2);
454 double *a =
new double[kpoints];
456 double *b =
new double[kpoints];
462 for(i = 0; i < kpoints; i++)
488 double** zmatrix =
new double*[kpoints];
490 for(i = 0; i < kpoints; i++)
494 zmatrix[i] =
new double[kpoints];
496 for(j = 0; j < kpoints; j++)
506 for(i = 0; i < kpoints; i++)
518 TriQL(kpoints, a, b, zmatrix);
522 for(i = 0; i < kpoints; i++)
554 void zwrk(
double* z,
double* w,
const int npt ,
const double alpha,
570 fprintf(stderr,
"too few points in formula\n");
594 int ncoeffs = (int)ceil(3.0*np/2);
600 double *a =
new double[ncoeffs+1];
602 double *b =
new double[ncoeffs+1];
608 for(i = 0; i < ncoeffs+1; i++)
626 double* a0 =
new double[ncoeffs];
628 double* b0 =
new double[ncoeffs];
632 chri1(ncoeffs,a,b,a0,b0,end0);
636 double s = b0[0]/fabs(b0[0]);
644 double* z1 =
new double[2*np-1];
646 double* w1 =
new double[2*np-1];
648 for(i = 0; i < ncoeffs; i++)
662 double** zmatrix =
new double*[2*np-1];
664 for(i = 0; i < 2*np-1; i++)
668 zmatrix[i] =
new double[2*np-1];
670 for(j = 0; j < 2*np-1; j++)
680 for(i = 0; i < 2*np-1; i++)
692 TriQL(2*np-1, z1, w1, zmatrix);
698 for(i = 0; i < 2*np-1; i++)
702 w1[i] = s*w1[i]/(z1[i]-end0);
714 for(i = 1; i < kpoints; i++)
746 void zwlk(
double* z,
double* w,
const int npt,
748 const double alpha,
const double beta)
762 fprintf (stderr,
"too few points in formula\n");
780 int kpoints = 2*np-1;
786 int ncoeffs = (int)ceil(3.0*np/2)-1;
792 double *a =
new double[ncoeffs+1];
794 double *b =
new double[ncoeffs+1];
800 for(i = 0; i < ncoeffs+1; i++)
820 double* a0 =
new double[ncoeffs];
822 double* b0 =
new double[ncoeffs];
826 chri1(ncoeffs,a,b,a0,b0,endl);
830 double* a1 =
new double[ncoeffs-1];
832 double* b1 =
new double[ncoeffs-1];
836 chri1(ncoeffs-1,a0,b0,a1,b1,endr);
842 double s = b1[0]/fabs(b1[0]);
850 double* z1 =
new double[2*np-3];
852 double* w1 =
new double[2*np-3];
854 for(i = 0; i < ncoeffs; i++)
868 double** zmatrix =
new double*[2*np-3];
870 for(i = 0; i < 2*np-3; i++)
874 zmatrix[i] =
new double[2*np-3];
876 for(j = 0; j < 2*np-3; j++)
886 for(i = 0; i < 2*np-3; i++)
898 TriQL(2*np-3, z1, w1, zmatrix);
904 double sumW1Z1 = 0.0;
906 for(i = 0; i < 2*np-3; i++)
910 w1[i] = s*w1[i]/(z1[i]-endl)/(z1[i]-endr);
914 sumW1Z1 += z1[i]*w1[i];
920 double c0 = b[0]-sumW1;
922 double c1 = a[0]*b[0]-sumW1Z1;
930 w[0] = (c0*endr-c1)/(endr-endl);
932 w[2*np-2] = (c1-c0*endl)/(endr-endl);
936 for(i = 1; i < kpoints-1; i++)
972 void Dgj(
double *D,
const double *z,
const int np,
const double alpha,
980 double one = 1.0, two = 2.0;
998 pd = (
double *)malloc(np*
sizeof(
double));
1000 jacobd(np,z,pd,np,alpha,beta);
1004 for (i = 0; i < np; i++){
1006 for (j = 0; j < np; j++){
1012 D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1016 D[i*np+j] = (alpha - beta + (alpha + beta + two)*z[j])/
1018 (two*(one - z[j]*z[j]));
1058 void Dgrjm(
double *D,
const double *z,
const int np,
const double alpha,
1076 double one = 1.0, two = 2.0;
1082 pd = (
double *)malloc(np*
sizeof(
double));
1086 pd[0] = pow(-one,np-1)*
gammaF(np+beta+one);
1090 jacobd(np-1,z+1,pd+1,np-1,alpha,beta+1);
1092 for(i = 1; i < np; ++i) pd[i] *= (1+z[i]);
1096 for (i = 0; i < np; i++) {
1098 for (j = 0; j < np; j++){
1102 D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1108 D[i*np+j] = -(np + alpha + beta + one)*(np - one)/
1114 D[i*np+j] = (alpha - beta + one + (alpha + beta + one)*z[j])/
1116 (two*(one - z[j]*z[j]));
1158 void Dgrjp(
double *D,
const double *z,
const int np,
const double alpha,
1176 double one = 1.0, two = 2.0;
1182 pd = (
double *)malloc(np*
sizeof(
double));
1188 jacobd(np-1,z,pd,np-1,alpha+1,beta);
1190 for(i = 0; i < np-1; ++i) pd[i] *= (1-z[i]);
1192 pd[np-1] = -
gammaF(np+alpha+one);
1198 for (i = 0; i < np; i++) {
1200 for (j = 0; j < np; j++){
1204 D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1210 D[i*np+j] = (np + alpha + beta + one)*(np - one)/
1212 (two*(alpha + two));
1216 D[i*np+j] = (alpha - beta - one + (alpha + beta + one)*z[j])/
1218 (two*(one - z[j]*z[j]));
1260 void Dglj(
double *D,
const double *z,
const int np,
const double alpha,
1278 double one = 1.0, two = 2.0;
1284 pd = (
double *)malloc(np*
sizeof(
double));
1288 pd[0] = two*pow(-one,np)*
gammaF(np + beta);
1292 jacobd(np-2,z+1,pd+1,np-2,alpha+1,beta+1);
1294 for(i = 1; i < np-1; ++i) pd[i] *= (one-z[i]*z[i]);
1296 pd[np-1] = -two*
gammaF(np + alpha);
1302 for (i = 0; i < np; i++) {
1304 for (j = 0; j < np; j++){
1308 D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
1314 D[i*np+j] = (alpha - (np-1)*(np + alpha + beta))/(two*(beta+ two));
1318 D[i*np+j] =-(beta - (np-1)*(np + alpha + beta))/(two*(alpha+ two));
1322 D[i*np+j] = (alpha - beta + (alpha + beta)*z[j])/
1324 (two*(one - z[j]*z[j]));
1386 double hgj (
const int i,
const double z,
const double *zgj,
1388 const int np,
const double alpha,
const double beta)
1398 if (fabs(dz) <
EPS)
return 1.0;
1448 double hgrjm (
const int i,
const double z,
const double *zgrj,
const int np,
1450 const double alpha,
const double beta)
1461 if (fabs(dz) <
EPS)
return 1.0;
1513 double hgrjp (
const int i,
const double z,
const double *zgrj,
const int np,
1515 const double alpha,
const double beta)
1526 if (fabs(dz) <
EPS)
return 1.0;
1578 double hglj (
const int i,
const double z,
const double *zglj,
const int np,
1580 const double alpha,
const double beta)
1590 if (fabs(dz) <
EPS)
return 1.0;
1626 void Imgj(
double *im,
const double *zgj,
const double *zm,
const int nz,
1628 const int mz,
const double alpha,
const double beta){
1636 for (i = 0; i < nz; ++i) {
1638 for (j = 0; j < mz; ++j)
1644 im [i*mz+j] =
hgj(i, zp, zgj, nz, alpha, beta);
1684 void Imgrjm(
double *im,
const double *zgrj,
const double *zm,
const int nz,
1686 const int mz,
const double alpha,
const double beta)
1696 for (i = 0; i < nz; i++) {
1698 for (j = 0; j < mz; j++)
1704 im [i*mz+j] =
hgrjm(i, zp, zgrj, nz, alpha, beta);
1744 void Imgrjp(
double *im,
const double *zgrj,
const double *zm,
const int nz,
1746 const int mz,
const double alpha,
const double beta)
1756 for (i = 0; i < nz; i++) {
1758 for (j = 0; j < mz; j++)
1764 im [i*mz+j] =
hgrjp(i, zp, zgrj, nz, alpha, beta);
1806 void Imglj(
double *im,
const double *zglj,
const double *zm,
const int nz,
1808 const int mz,
const double alpha,
const double beta)
1818 for (i = 0; i < nz; i++) {
1820 for (j = 0; j < mz; j++)
1826 im[i*mz+j] =
hglj(i, zp, zglj, nz, alpha, beta);
1920 void jacobfd(
const int np,
const double *z,
double *poly_in,
double *polyd,
1922 const int n,
const double alpha,
const double beta){
1926 double zero = 0.0, one = 1.0, two = 2.0;
1940 for(i = 0; i < np; ++i)
1946 for(i = 0; i < np; ++i)
1956 for(i = 0; i < np; ++i)
1958 poly_in[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
1962 for(i = 0; i < np; ++i)
1964 polyd[i] = 0.5*(alpha + beta + two);
1974 double two = 2.0, apb = alpha + beta;
1976 double *poly, *polyn1,*polyn2;
1982 polyn1 = (
double *)malloc(2*np*
sizeof(
double));
1992 polyn1 = (
double *)malloc(3*np*
sizeof(
double));
2002 for(i = 0; i < np; ++i){
2006 polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
2012 for(k = 2; k <= n; ++k){
2014 a1 = two*k*(k + apb)*(two*k + apb - two);
2016 a2 = (two*k + apb - one)*(alpha*alpha - beta*beta);
2018 a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb);
2020 a4 = two*(k + alpha - one)*(k + beta - one)*(two*k + apb);
2032 for(i = 0; i < np; ++i){
2034 poly [i] = (a2 + a3*z[i])*polyn1[i] - a4*polyn2[i];
2036 polyn2[i] = polyn1[i];
2038 polyn1[i] = poly [i];
2048 a1 = n*(alpha - beta);
2050 a2 = n*(two*n + alpha + beta);
2052 a3 = two*(n + alpha)*(n + beta);
2054 a4 = (two*n + alpha + beta);
2056 a1 /= a4; a2 /= a4; a3 /= a4;
2062 for(i = 0; i < np; ++i){
2064 polyd[i] = (a1- a2*z[i])*poly[i] + a3*polyn2[i];
2066 polyd[i] /= (one - z[i]*z[i]);
2120 void jacobd(
const int np,
const double *z,
double *polyd,
const int n,
2122 const double alpha,
const double beta)
2132 for(i = 0; i < np; ++i) polyd[i] = 0.0;
2138 jacobfd(np,z,polyd,NULL,n-1,alpha+one,beta+one);
2140 for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (
double)n + one);
2180 if (x == -0.5) gamma = -2.0*sqrt(M_PI);
2182 else if (!x)
return gamma;
2184 else if ((x-(
int)x) == 0.5){
2204 else if ((x-(
int)x) == 0.0){
2224 fprintf(stderr,
"%lf is not of integer or half order\n",x);
2248 static void Jacobz(
const int n,
double *z,
const double alpha,
2254 double dth = M_PI/(2.0*(double)n);
2256 double poly,pder,rlast=0.0;
2260 double one = 1.0, two = 2.0;
2270 for(k = 0; k < n; ++k){
2272 r = -cos((two*(
double)k + one) * dth);
2274 if(k) r = 0.5*(r + rlast);
2278 for(j = 1; j <
STOP; ++j){
2280 jacobfd(1,&r,&poly, &pder, n, alpha, beta);
2284 for(i = 0, sum = 0.0; i < k; ++i) sum += one/(r - z[i]);
2288 delr = -poly / (pder - sum * poly);
2292 if( fabs(delr) <
EPS )
break;
2358 void JacZeros(
const int n,
double *a,
double*b,
const double alpha,
2370 double **z =
new double*[n];
2372 for(i = 0; i < n; i++)
2376 z[i] =
new double[n];
2378 for(j = 0; j < n; j++)
2388 for(i = 0; i < n; i++)
2418 static void RecCoeff(
const int n,
double *a,
double *b,
const double alpha,
2426 double apb, apbi,a2b2;
2446 a[0] = (beta-alpha)/apbi;
2448 b[1] = (4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
2452 a2b2 = beta*beta-alpha*alpha;
2456 for(i = 1; i < n-1; i++){
2458 apbi = 2.0*(i+1) + apb;
2460 a[i] = a2b2/((apbi-2.0)*apbi);
2462 b[i+1] = (4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
2464 ((apbi*apbi-1)*apbi*apbi));
2472 a[n-1] = a2b2/((apbi-2.0)*apbi);
2534 static void TriQL(
const int n,
double *d,
double *e,
double **z){
2538 double s,r,
p,g,f,dd,c,b;
2542 double MuZero = e[0];
2548 for(i = 0; i < n-1; i++)
2552 e[i] = sqrt(e[i+1]);
2568 for (m=l;m<n-1;m++) {
2570 dd=fabs(d[m])+fabs(d[m+1]);
2572 if (fabs(e[m])+dd == dd)
break;
2578 if (iter++ ==
STOP){
2580 fprintf(stderr,
"triQL: Too many iterations in TQLI");
2586 g=(d[l+1]-d[l])/(2.0*e[l]);
2590 g=d[m]-d[l]+e[l]/(g+
sign(r,g));
2596 for (i=m-1;i>=l;i--) {
2602 if (fabs(f) >= fabs(g)) {
2626 r=(d[i]-g)*s+2.0*c*b;
2638 for(k = 0; k < n; k++)
2644 z[k][i+1] = s*z[k][i] + c*f;
2646 z[k][i] = c*z[k][i] - s*f;
2674 for(i = 0; i < n-1; ++i){
2680 for(l = i+1; l < n; ++l)
2696 double temp = z[0][k];
2708 for(i =0 ; i < n; i++)
2712 e[i] = MuZero*z[0][i]*z[0][i];
2746 int size = (int)floor(n/2.0)+2;
2748 double *s =
new double[size];
2750 double *t =
new double[size];
2756 for(i = 0; i < size; i++)
2770 for(m = 0; m <= n-2; m++)
2776 for(k = (
int)floor((m+1)/2.0); k >= 0; k--)
2782 u = u+(a[k+n+1]-a[l])*t[k+1] + b[k+n+1]*s[k] - b[l]*s[k+1];
2804 for(j = (
int)floor(n/2.0); j >= 0; j--)
2814 for(m = n-1; m <= 2*n-3; m++)
2820 for(k = m+1-n; k <= floor((m-1)/2.0); k++)
2828 u = u-(a[k+n+1]-a[l])*t[j+1] - b[k+n+1]*s[j+1] + b[l]*s[j+2];
2842 a[k+n+1] = a[k] + (s[j+1]-b[k+n+1]*s[j+2])/t[j+2];
2852 b[k+n+1] = s[j+1]/s[j+2];
2872 a[2*n ] = a[n-1]-b[2*n]*s[1]/t[1];
2900 void chri1(
int n,
double* a,
double* b,
double* a0,
2902 double* b0,
double z)
2908 double q = ceil(3.0*n/2);
2910 int size = (int)q+1;
2916 fprintf(stderr,
"input arrays a and b are too short\n");
2920 double* r =
new double[n+1];
2924 r[1] = z - a[1] - b[1]/r[0];
2926 a0[0] = a[1] + r[1] - r[0];
2942 for(k = 1; k < n; k++)
2946 r[k+1] = z - a[k+1] - b[k+1]/r[k];
2948 a0[k] = a[k+1] + r[k+1] - r[k];
2950 b0[k] = b[k] * r[k]/r[k-1];
2968 std::complex<Nektar::NekDouble>
ImagBesselComp(
int n,std::complex<Nektar::NekDouble> y)
2970 std::complex<Nektar::NekDouble> z (1.0,0.0);
2971 std::complex<Nektar::NekDouble> zbes (1.0,0.0);
2972 std::complex<Nektar::NekDouble> zarg;
2979 while (abs(z) > tol && i <= maxit){
2980 z = z*(1.0/i/(i+n)*zarg);
2981 if (abs(z) <= tol)
break;
void Dgj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated.
#define sign(a, b)
return the sign(b)*a
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and.
void zwgrjp(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=1.
double hgj(const int i, const double z, const double *zgj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through.
void Imgrjp(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Radau-Jacobi points.
void zwglj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.
double hgrjp(const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the.
void zwlk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Lobatto-Kronrod-Jacobi zeros and weights.
void JacZeros(const int n, double *a, double *b, const double alpha, const double beta)
Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal.
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1.
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
void Imglj(double *im, const double *zglj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Lobatto-Jacobi points.
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the.
void Dgrjm(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated.
void zwgk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Kronrod-Jacobi zeros and weights.
void Imgrjm(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Radau-Jacobi points.
void Dgrjp(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the.
void zwgrjm(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
The namespace associated with the the Polylib library (Polylib introduction)
double optdiff(double xl, double xr)
The following function is used to circumvent/reduce "Subtractive Cancellation" The expression 1/dz is...
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
static void Jacobz(const int n, double *z, const double alpha, const double beta)
Calculate the n zeros, z, of the Jacobi polynomial, i.e.
void zwrk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Radau-Kronrod-Jacobi zeros and weights.
double hglj(const int i, const double z, const double *zglj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the.
std::complex< Nektar::NekDouble > ImagBesselComp(int n, std::complex< Nektar::NekDouble > y)
Calcualte the bessel function of the first kind with complex double input y. Taken from Numerical Rec...
double hgrjm(const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the.
#define EPS
Precision tolerance for two points to be similar.
void Dglj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the.
void Imgj(double *im, const double *zgj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Jacobi points to an.
double laginterp(double z, int j, const double *zj, int np)
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
double gammaF(const double)
Calculate the Gamma function , , for integer.
void zwgj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Jacobi zeros and weights.