/* Root finding and numerical integration by Joung-Mo Kang (mokang@mit.edu)
   1/7/00, submitted for pset question #4.1
   This program solves for x in the equation integral(from t=0 to 10) of
   x * exp(-xt^2)dt = 1/3 using the trapezoid rule to evaluate the integral
   and bisection to find the root. */

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

double integral(double,double);           /*function to be integrated*/
double func(double,double,int);               /*function to be solved*/
double traprule(double (*f)(double,double),double x,double t_init,double t_final,int N);
                                    /*trapezoid rule on a generic function*/
double bisect_root(double (*func)(double,double,int),double t,int N,double left_init,double right_init,double ATOL,double RTOL);
                                    /*bisect rule on a func with 1 param (x)*/
int sign(double num);              /*returns sign of a double*/

double integral(double t,double x) {
   return x * exp(-1*x*t*t);
}

double func(double t,double x,int N) {
   return traprule(integral,x,0,t,N)-1.0/3.0;
}

double traprule(double (*func)(double,double),double x,double t_init,double t_final,int N) {
   /* the area of a single trapezoid is h/2 * (f(x1) + f(x2)) */
   double h = (t_final - t_init) / N;                 /*h is the width of a single trapezoid*/
   double area = func(t_init/2.0,x) + func(t_final/2.0,x);    /*initialize sum with endpoint vals*/
   double t_current;                                  /*t-index for running sum*/
   for (t_current = t_init + h;t_current < t_final;t_current += h) {
      area += func(t_current,x);
   }
   area *= h;
   return area;
}

/*sign takes in any number up to a double and returns an int: 1 if positive, -1 if negative, 0 if neither*/
int sign(double num) {
   if (num > 0)
      return 1;
   if (num < 0)
      return -1;
   else return 0;
}

double bisect_root(double (*func)(double,double,int),double t,int N,double left_init,double right_init,double ATOL,double RTOL) {
   double left = left_init;
   double right = right_init;                   /*initialize iteration boundary params*/
   double midpoint = (left + right)/2.0;        /*stores midpoint between each pair of boundaries*/
   int left_sign = sign(func(t,left_init,N));
   int right_sign = sign((func(t,right_init,N)));  
   int midpoint_sign;                           /*store signs*/
   while (((right - left) > (ATOL + fabs(midpoint)*RTOL)) || ((func(t,midpoint,N)) > ATOL)) {
                               /*both x-variance and y-variance must satisfy tolerances*/
      midpoint_sign = sign(func(t,midpoint,N));
      if (midpoint_sign == left_sign)
         left = midpoint;                       /*If f(midpoint) has same sign as the left endpoint,*/
         else if (midpoint_sign == right_sign)  /*make midpoint new left endpoint and repeat.       */
              right = midpoint;                 /*Similar for right. If it happens to be 0 (very    */
              else left = right = midpoint;     /*unlikely) it sets both endpoints to the mid and   */
                                                /*loop should terminate on next iteration.          */

      midpoint = (left + right)/2.0;
   }
   return midpoint;
}


int main(void) {
   double ATOL,RTOL;                            /*User will define tolerances*/
   double t_final = 10;                         /*Limit of integral*/
   int N = 100;                                 /*iterations of traprule*/
   double trap_left = 0;
   double trap_right = 1;                       /*limits for traprule*/
   double answer;                               /*Is that your final answer?*/

   /*Prompt user for tolerances*/
   printf("Input ATOL\n");
   scanf("%lf",&ATOL);
   printf("Input RTOL\n");
   scanf("%lf",&RTOL);

   /*crunch crunch crunch*/
   answer = bisect_root(func,t_final,N,trap_left,trap_right,ATOL,RTOL);

   /*print final answer*/
   printf("\nx = %f\n",answer);

return 0;
}

/* sample input and output :
athena% hw4_1.exe
Input ATOL
1e-9
Input RTOL
1e-7

x = 0.135402 */