/*******************************************************************/
/*   Jeff Balsley   Physics 551                                    */
/*                                                                 */
/*   gcc -lm integrator.c -o integrator                            */
/*******************************************************************/

#include<stdio.h>
#include<math.h>
#include"const.h"    /*  My own header with constants in it  */
#define accuracy 1.0e-6

double f(double x);

int main(void)
{
   double x0, x1;             /*  limits of integration */
   double h;                  /*  interval size  */
   double ans, ans_old;       /*  value of integral  */
   int i;
   printf("Enter limits of integration > ");
   scanf("%lf %lf", &x0, &x1);
   
   /*printf("Enter interval size > ");*/
   /*scanf("%lf", &h);*/
   
   h = (x1-x0)/2;
   ans = 0;
   ans_old = 1;

   while( fabs(ans - ans_old)/ans > accuracy){
      ans_old = ans;
      ans = 0;
      printf("%d   ",(int)((x1-x0)/h) ); 
      for(i=0; i<(int)((x1-x0)/h); ++i){
	 ans = .5*(f(x0 + i*h) + f(x0 + (i+1)*h))*h + ans;
      }
      printf("Integral = %.10e\n", ans);
      h = h/2;
   }
   printf("\nIntegral = %e\n", ans);
   return(0);
}

double f(double x)
{
   /*  The function here is dB/dT  */
   double T = 15000;
   double const1 = ( 2*h*h*pow(x,4) )/( k*c*c*T*T );
   double const2 = ( exp( h*x/(k*T) ) / pow(exp(h*x/(k*T)) - 1, 2));
   
   return ( const1*const2 );   
}