/**********************************************************************/
/*  Jeff Balsley                                                      */
/*                                                                    */
/*  This was Scott's code.                                            */
/*  It is an adaptive integration routine.  The interval will change  */
/*  based on the accuracy you want.                                   */
/*  gcc -lm adaptive_quad.c -o adaptive_quad                          */
/**********************************************************************/

#include<math.h>
#include<stdio.h>
#include"const.h"

/* adaptive_quad.c: Using adaptive quadrature and the trapazoidal rule,
   evaluate the integral the function f_x from x_min to x_max.  Here
   we choose to use recursion for easy of programming.  Although this is 
   not the most efficient method, it is sufficent to solve the problem
   at hand and easier to implement. */

/* Desired accuracy of each iteration. */
#define accuracy 1e-9

/* Function to integrate. */
double f_x(double x);
double adaptive_quad(double x_min, double x_max); 

int main(void){
     
     double xmin, xmax;              /*  limits of integration  */

     printf("Enter limits of integration > ");
     scanf("%lf %lf", &xmin, &xmax);

     printf("Integral = %e\n", adaptive_quad(xmin, xmax));

     return(0);
}

double adaptive_quad(double x_min,   /* I: Lower bound of integration. */
		     double x_max)   /* I: Upper bound of integration. */
{
     double I;      /* Approximate integral over given interval. */
     double f_max;  /* Function evaluated at x_max. */
     double f_min;  /* Function evaluated ay x_min. */
     double f_mid;  /* Function evaluated at x_mid. */
     double mid;    /* Midpoint of interval. */
     /* static double count = 0; /* Keep track of # of function calls. */

     /* count++;               */
     /* printf("%f\n", count); */

     f_max = f_x(x_max);
     f_min = f_x(x_min);

     /* Approximate the total integral using the trapezoidal rule. */
     I = 0.5*(x_max - x_min)*(f_max + f_min);
     
     /* Break the interval into to pieces and find the 
	integral of each piece. */
     mid = x_min + (x_max - x_min)/2.0;
     f_mid = f_x(mid);

     /* If we've reached our accuracy, then return our value.  Otherwise
	recursively act on interval. */
     if(I - (0.5*(mid - x_min)*(f_min + f_mid)
	     + 0.5*(x_max - mid)*(f_max + f_mid)) < accuracy)
	  return I;
     else
	  return (adaptive_quad(x_min, mid) + adaptive_quad(mid, x_max));
}

double f_x(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 );
}