/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  Gaussian Quadrature Routines                                           //
//                                                                         //
//  B. Militzer                                     Berkeley 04-18-17      //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#include "Quadrature.h"

double gammln(double xx) {
  double x,y,tmp,ser;
  static double cof[6]={76.18009172947146,-86.50532032941677,
			24.01409824083091,-1.231739572450155,
			0.1208650973866179e-2,-0.5395239384953e-5};
  y=x=xx;
  tmp=x+5.5;
  tmp -= (x+0.5)*log(tmp);
  ser=1.000000000190015;
  for (int j=0; j<=5; j++) 
    ser += cof[j]/++y;

  return -tmp+log(2.5066282746310005*ser/x);
}

void GaussJacobiQuadraturePoints(Array1 <double> & x, Array1 <double> & w, const int n, const double alf, const double bet) {
  double an,bn,r1,r2,r3;
  double a,b,c,p1,p2,pp,z;

  x.Resize(n);
  w.Resize(n);
  for (int i=1; i<=n; i++) {
    if (i == 1) {
      an=alf/n;
      bn=bet/n;
      r1=(1.0+alf)*(2.78/(4.0+n*n)+0.768*an/n);
      r2=1.0+1.48*an+0.96*bn+0.452*an*an+0.83*an*bn;
      z=1.0-r1/r2;
    } else if (i == 2) {
      r1=(4.1+alf)/((1.0+alf)*(1.0+0.156*alf));
      r2=1.0+0.06*(n-8.0)*(1.0+0.12*alf)/n;
      r3=1.0+0.012*bet*(1.0+0.25*fabs(alf))/n;
      z -= (1.0-z)*r1*r2*r3;
    } else if (i == 3) {
      r1=(1.67+0.28*alf)/(1.0+0.37*alf);
      r2=1.0+0.22*(n-8.0)/n;
      r3=1.0+8.0*bet/((6.28+bet)*n*n);
      z -= (x[1]-z)*r1*r2*r3;
    } else if (i == n-1) {
      r1=(1.0+0.235*bet)/(0.766+0.119*bet);
      r2=1.0/(1.0+0.639*(n-4.0)/(1.0+0.71*(n-4.0)));
      r3=1.0/(1.0+20.0*alf/((7.5+alf)*n*n));
      z += (z-x[n-3])*r1*r2*r3;
    } else if (i == n) {
      r1=(1.0+0.37*bet)/(1.67+0.28*bet);
      r2=1.0/(1.0+0.22*(n-8.0)/n);
      r3=1.0/(1.0+8.0*alf/((6.28+alf)*n*n));
      z += (z-x[n-2])*r1*r2*r3;
    } else {
      z=3.0*x[i-1]-3.0*x[i-2]+x[i-3];
    }

    double alfbet=alf+bet;
    const int maxIt = 10;
    int its=0;
    double temp;
    for (; its<=maxIt; its++) {
      temp=2.0+alfbet;
      p1=(alf-bet+temp*z)/2.0;
      p2=1.0;
      for (int j=2; j<=n; j++) {
	double p3=p2;
	p2=p1;
	temp=2*j+alfbet;
	a=2*j*(j+alfbet)*(temp-2.0);
	b=(temp-1.0)*(alf*alf-bet*bet+temp*(temp-2.0)*z);
	c=2.0*(j-1+alf)*(j-1+bet)*temp;
	p1=(b*p2-c*p3)/a;
      }
      pp=(n*(alf-bet-temp*z)*p1+2.0*(n+alf)*(n+bet)*p2)/(temp*(1.0-z*z));
      double z1=z;
      z=z1-p1/pp;
      const double eps = 3.0e-14;
      if (fabs(z-z1) <= eps) break;
    }
    if (its > maxIt) error("Too many iterations in GaussJacobiQuadraturePoints()",its);
    x[i]=z;
    w[i]=exp(gammln(alf+n)+gammln(bet+n)-gammln(n+1.0)-gammln(n+alfbet+1.0))*temp*pow(2.0,alfbet)/(pp*p2);
  }
}

void GaussLegendreQuadraturePoints(const double x1, const double x2, const int n, Array1 <double> & x, Array1 <double> & w) {
  x.Resize(n);
  w.Resize(n);
  double z1,pp,p3,p2,p1;

  int m=(n+1)/2;
  double xm=0.5*(x2+x1);
  double xl=0.5*(x2-x1);

  //  const double eps = 3.0e-11;
  const double eps = 3.0e-15; // needed to get good agreement for Sean's F90 code
  for (int i=0; i<m; i++) {
    double z=cos(pi*(i+0.75)/(n+0.5));
    do {
      p1=1.0;
      p2=0.0;
      for (int j=1; j<=n; j++) {
	p3=p2;
	p2=p1;
	p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
      }
      pp=n*(z*p1-p2)/(z*z-1.0);
      z1=z;
      z=z1-p1/pp;
    } while (fabs(z-z1) > eps);
    x[i]    =xm-xl*z;
    x[n-1-i]=xm+xl*z; // in fortran x[n+1-i]
    w[i]    =2.0*xl/((1.0-z*z)*pp*pp);
    w[n-1-i]=w[i];
  }
}