#include "iostream.h"
#include "math.h"

void qgausx(double, double, double&);
void qgausy(double, double, double&);
void qgausz(double, double, double&);

double g(double);
double h(double);
double f(double);
double xmax, x, y, z;

double func(double x, double y, double z)
{
    return x*x + y*y + z*z;
}

double z1(double x, double y)
{
    return -sqrt(fabs(xmax * xmax - x * x - y * y));
}

double z2(double x, double y)
{
    return sqrt(fabs(xmax * xmax - x * x - y * y));
}

double y1(double x)
{
    return -sqrt(fabs(xmax * xmax - x * x));
}

double y2(double x)
{
    return sqrt(fabs(xmax * xmax - x * x));
}

void qgausx(double a, double b, double& ss)
{
    double x1[6], w[6];
	double xm,xr,dx;
    x1[1] = 0.1488743389; x1[2] = 0.4333953941; x1[3] = 0.6794095682;
    x1[4] = 0.8650633666; x1[5] = 0.9739065285;
    w[1] = 0.2955242247; w[2] = 0.2692667193; w[3] = 0.2190863625;
    w[4] = 0.1494513491; w[5] = 0.0666713443;
    xm = 0.5 * (b + a);
    xr = 0.5 * (b - a);
    ss = 0.0;
    for (int j = 1; j<=5; j++)
	{
        dx = xr * x1[j];
        ss = ss + w[j] * (h(xm + dx) + h(xm - dx));
    }
    ss = xr * ss;
}

void qgausy(double a, double b, double& ss)
{
    double x1[6], w[6];
	double xm,xr,dx;
    x1[1] = 0.1488743389; x1[2] = 0.4333953941; x1[3] = 0.6794095682;
    x1[4] = 0.8650633666; x1[5] = 0.9739065285;
    w[1] = 0.2955242247; w[2] = 0.2692667193; w[3] = 0.2190863625;
    w[4] = 0.1494513491; w[5] = 0.0666713443;
    xm = 0.5 * (b + a);
    xr = 0.5 * (b - a);
    ss = 0.0;
    for (int j = 1; j<=5; j++)
	{
        dx = xr * x1[j];
        ss = ss + w[j] * (g(xm + dx) + g(xm - dx));
    }
    ss = xr * ss;
}

void qgausz(double a, double b, double& ss)
{
    double x1[6], w[6];
	double xm,xr,dx;
    x1[1] = 0.1488743389; x1[2] = 0.4333953941; x1[3] = 0.6794095682;
    x1[4] = 0.8650633666; x1[5] = 0.9739065285;
    w[1] = 0.2955242247; w[2] = 0.2692667193; w[3] = 0.2190863625;
    w[4] = 0.1494513491; w[5] = 0.0666713443;
    xm = 0.5 * (b + a);
    xr = 0.5 * (b - a);
    ss = 0.0;
    for (int j = 1; j<=5; j++)
	{
        dx = xr * x1[j];
        ss = ss + w[j] * (f(xm + dx) + f(xm - dx));
    }
    ss = xr * ss;
}

void main()
{
    //program d3r9
    //driver for routine quad3d
    const double pi = 3.1415926;
    int nval = 10;
	double xmin,s;
    cout<<endl;
    cout<<"integral of r^2 over a spherical volume"<<endl;
    cout<<endl;
    cout<<"      radius      quad3d      actual"<<endl;
	cout.setf(ios::fixed|ios::right);
	cout.precision(4);
    for (int i = 1; i<=nval; i++)
	{
        xmax = 0.1 * i;
        xmin = -xmax;
        quad3d(xmin, xmax, s);
		cout.width(2);
        cout<<"      "<<xmax;
		cout.width(12);
        cout<<s;
		cout.width(12);
        cout<<(4.0 * pi * pow(xmax , 5)) / 5.0<<endl;
    }
}