/* Integration using simpson's rule.
   Maybe you can make this code more elegant.

   Last edit: 30 jan 08
*/



// include the includes. you just need to know what to include.
// you learn this in C++ class.

#include <iostream>
#include <vector>
#include <string>
#include <stdexcept>
#include <cmath>   // you need this for math functions.

using namespace std;    // a useful trick

int main()
{

  // declare some variables

  int N;
  double a,b,h;
  double sum = 0.0;
  double weight = 0.0;
  double f(double); // declare function that will be integrated
                    // note function is not defined.
                    // you don't know what is does. that comes later.


  cout << " enter odd N: ";    // get the number of steps.
    cin >> N;

    // check that N is odd:

    if (N%2 == 0){     //the % operator is 'remainder division"
      cout << " \t --> we need N odd." << endl;
      return 1;
    }


    cout << " enter lower bound a: ";
    cin >> a;

    cout << " enter upper bound b: ";
    cin >> b;

    if (a>b){
      cout << " \t ====> we need b > a" << endl;
      return 1;
    }


    h = (b-a)/(double(N)); //get the step size
    
    for (int i = 1; i < N + 1; ++i){
      if ((i ==1) || (i==N)) weight = h/3.; // get the weights
      if ((i % 2) ==0) weight = 4*h/3.;
      if (((i % 2) == 1) && (i != N)) weight = 2*h/3.;
	  
      sum = sum + f(a +((i-1)*h))*weight;  
    }

    cout << "integral is : " << sum << endl;


  return 0;
}

// the function to be integrated.

double f (double x){
  return exp(-x);
}