// Mr. Minich
// CMPSC 201
// Ch. 13 Demo Program #7
// December 7, 2002 
// Purpose - to illustrate the use of the Multiple-Application Trapezoidal
//				and Simpsons 1/3 rules for approximating integrals

#include <iostream.h>
#include <math.h>

double f(double x, double A, double B, double C);
double trapezoid(double lowerLimit, double upperLimit, int n, double A, double B, double C);
double simpson1_3(double lowerLimit, double upperLimit, int n, double A, double B, double C);

int main()
{
		double sumT = 0.0;		// integral according to Trapezoidal Rule
		double sumS1_3 = 0.0;	// integral according to Simpson's 1/3 Rule
		double a = 0.0;			// lower limit of region to integrate
		double b = 0.0;			// upper limit of region to integrate
        int n = 0;				// number of segments
		double A = 0.0;			// coefficient for x squared term of user-inputed quadratic eq'n
		double B = 0.0;			// coefficient for x term of user-inputed quadratic eq'n
		double C = 0.0;			// constant term for user-inputed quadratic eq'n

		cout << "Enter the coefficients for the quadratic equation Ax^2 + Bx + C = 0\n";
		cout << "Enter coefficient A: ";
		cin >> A;
		cout << "Enter coefficient B: ";
		cin >> B;
		cout << "Enter coefficient C: ";
		cin >> C;
		cout << "\nEnter lower limit (a) of the region to integrate: ";
		cin >> a;
		cout << "Enter upper limit (b) of the region to integrate: ";
		cin >> b;
		cout << "Enter the number of segments (n; must be even): ";
		cin >> n;

        sumT = trapezoid(a, b, n, A, B, C);
        sumS1_3 = simpson1_3(a, b, n, A, B, C);

        cout << "Trapezoidal Rule: " << sumT << endl;
		cout << "Simpson's 1/3 Rule: " << sumS1_3 << endl << endl;

		return 0;
}// end of main

double trapezoid(double lowerLimit, double upperLimit, int n, double A, double B, double C)
{			
        double sum = 0.0;		// integral approximatoin
		double h = 0.0;			// average height of each segment
        int i = 0;				// loop variable

        h = (upperLimit - lowerLimit) / n;
        sum = f(lowerLimit, A, B, C) + f(upperLimit, A, B, C);

        for (i = 1; i < n; i++)
		{
			sum += 2 * f(lowerLimit + i * h, A, B, C);
        }

        return h / 2.0 * sum;
}// end of trapezoid

double simpson1_3(double lowerLimit, double upperLimit, int n, double A, double B, double C)
{
        // precondition: n is even

		double sum = 0.0;		// area of current "slice"
		double area = 0.0;		// integral approximation
		double a = 0.0;			// left endpoint of current "slice"
		double b = 0.0;			// right endpoint of current "slice"
		double h = 0.0;			// average height of each segment
        int i = 0;				// loop variable

        h = (upperLimit - lowerLimit) / n;
        a = lowerLimit;
        b = lowerLimit + 2 * h;

        for (i = 1; i <= n; i += 2)
		{
			area = (b - a)/6 * (f(a, A, B, C) + 4 * f((b + a) / 2, A, B, C) + f(b, A, B, C));
			sum = sum + area;
			a = b;
			b = a + 2 * h;
        }

        return sum;
}// end of simpson1_3

double f(double x, double A, double B, double C)
{
	return A * pow(x, 2) + B * x + C;
}// end of f