// CMPSC 201
// Integration Demo #1

#include <iostream>
#include <cmath>
using namespace std;

double f(double x);
double trapezoid1(double lowerLimit, double upperLimit, int n);
double trapezoid2(double lowerLimit, double upperLimit, int n);
double simpson1_3(double lowerLimit, double upperLimit, int n);
double simpson3_8(double lowerLimit, double upperLimit, int n);

int main()
{
   double a = 0.0;         // lower limit
   double b = 0.0;         // upper limit
   int n = 0;              // number of segments

   cout << "The function to integrate is f(x) = x * x" << endl;
   cout << "Enter lower limit (a): ";
   cin >> a;
   cout << "Enter upper limit (b): ";
   cin >> b;
   cout << "Enter even number of segments: ";
   cin >> n;
                
   cout << "Trapezoid #1" << "\t" << trapezoid1(a, b, n) << endl; 
   cout << "Trapezoid #2" << "\t" << trapezoid2(a, b, n) << endl; 
   cout << "Simpsons 1/3" << "\t" << simpson1_3(a, b, n) << endl;
   cout << "Simpsons 3/8" << "\t" << simpson3_8(a, b, n) << endl;

   return 0;
}// end of main

// It is more important that CMPSC 201 students understand the trapezoid1
// version of using the Multiple Application Trapezoidal Rule than the
// trapezoid2 version below.
double trapezoid1(double lowerLimit, double upperLimit, int n)
{
   double sum = 0.0;
   double area = 0.0;
   double a = 0.0;
   double b = 0.0;
   double h = 0.0;
   int i = 0; 

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

   for (i = 1; i <= n; i++)
   {
      area = (b - a)/2 * (f(a) + f(b));
      sum = sum + area;							// adding the areas of each small trapezoid together
      a = b;
      b = a + h;
   }

   return sum;
}// end of trapezoid1

// This trapezoid2 version of using the Multiple Application Trapezoidal Rule is less important
// than the version in trapezoid1 above.
double trapezoid2(double lowerLimit, double upperLimit, int n)
{
   double sum = 0.0;
   double h = 0.0;
   int i = 0; 

   h = (upperLimit - lowerLimit) / n;
   sum = f(lowerLimit) + f(upperLimit);

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

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

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

   double sum = 0.0;
   double area = 0.0;
   double a = 0.0;
   double b = 0.0;
   double h = 0.0;
   int i = 0; 

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

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

   return sum;
}// end of simpson1_3

// we may not be studying the simpson 3/8 method in CMPSC 201 this semester
// ask the instructor if you are responsible for undertanding this function
double simpson3_8(double lowerLimit, double upperLimit, int n)
{
   // precondition: n is odd

   double sum = 0.0;
   double area = 0.0;
   double a = 0.0;
   double b = 0.0;
   double h = 0.0;
   int i = 0; 

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

   for (i = 1; i <= n; i += 3)
   {
      area = (b - a)/8 * (f(a) + 3 * f(a + h) + 3 * f(a+2 * h) + f(b));
      sum = sum + area;
      a = b;
      b = a + 3 * h;
   }

   return sum;
}// end of simpson3_8

double f(double x)
{
   double y = 0.0;  
	
   y = x * x;
	
   return y;
}// end of f