 	0.000000000000000000000000000000000
 };
 
-__device__ void qgauss(double a, double b, double d, double ep1, double ep2, double sigma, double k1, double k2, double tau, double *out){
-	// translated from Fortran from Quadpack's dqk21.f
+__device__ void qg10k21(double a, double b, double *out, double d, double ep1, double ep2, double sigma, double k1, double k2, double tau) {
+	// Do quadrature with the Gauss-10, Kronrod-21 rule
+	// translated from Fortran from Quadpack's dqk21.f...
 
 	double center = 0.5*(b+a);
 	double halflength = 0.5*(b-a);

			* (cyl_bessel_i0(y) + (k2*f+k1*(1.-f))*tau
				* cyl_bessel_i1(y)/y); // Full Expression
	}

}


/*
	Parameters ripped from Quadpack dqk21.f on 3/20/2017:

	c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
	c bell labs, nov. 1981.
*/

__constant__ static double wg[5] = {
	0.066671344308688137593568809893332,
