/*     Problem 1                                   */
/*     Romberg Extrapolation of Simpson's Rule     */
/*     Aptil 19 1999                               */
/*     Calculats the gamma function                */
/*     gcc -lm  integrater.c -o integrater         */

#include <stdio.h>
#include <math.h>
# define SIZE 21         /*  defines the size of the array
                          *  that will hold the integrals
                          *  also equal to the number of 
                          *  different h's to be used     */

double get_value(double x, douple p);
double integrate(double x1, double x2, double h, double p);

int
main(void)
{
   double x1, x2, dist, h, x, integral, slope, y_int;
   int j, a, r, i;
   double p;                /*  gamma of p  */
   double stuff[SIZE][2];
   double line[100][2];
   FILE *out_p;
   FILE *line_p;

   printf("Enter the limits of integration> "); 
   scanf("%lf%lf", &x1, &x2);
   printf("This program calculates the gamma function.\n");
   printf("what is your argument for the gamma function >");
   scanf("%lf", &p);
   dist = x2-x1;
   h = dist/SIZE;

   /* initiate the first column of the array to all the h's */

   stuff[0][0] = dist/2;
   for(j=1; j < SIZE; ++j)
      stuff[j][0] = (stuff[j-1][0] / 2);
   /* set the second column of the array equal to the integral done 
    * with the corresponding h in the first column   */

   for(a=0; a < SIZE; ++a)
      stuff[a][1] = integrate(x1, x2, stuff[a][0], p);

   out_p = fopen("output.dat", "w");
   line_p = fopen("line.dat", "w");

   /*  print out the array with h's and integrals.  To be plotted on sm */

   for(r=0; r<SIZE; ++r)
      fprintf(out_p, "%f, %f\n", stuff[r][0], stuff[r][1]);
   fclose(out_p);

   /*  fit last two points to a straight line  */ 

   slope = ( (stuff[SIZE-1][1] - stuff[SIZE-2][1])/(stuff[SIZE-1][0] - stuff[SIZE-2][0]) );
   y_int = stuff[SIZE-1][1] - (slope * stuff[SIZE-1][0]);

   /*  now that we have the slope and y_int of line, make a 100x2 array
    *  of separate points on line to be plotted on sm       */

   line[0][0] = 0;
   line[0][1] = y_int;
   for (i=1; i<100; ++i)
     {
      line[i][0] = line[i-1][0] + dist/100;
      line[i][1] = slope*line[i][0] + y_int;
      fprintf(line_p, "%f %f\n", line[i][0], line[i][1]);  
     }  
    fclose(line_p); 

    /*  print value of integral  */

    printf("The integral is %f\n", stuff[SIZE-1][1]);
 return(0);
}

double 
integrate(double x1, double x2, double h, double p)
{
   int i;
   double a, integral, x;

   integral=(h/3)*(get_value(x1, p) + get_value(x2, p));

   i = 1;
   x = x1 + h;
   
   for (i=1; i*h < x2-x1; ++i )
     {
       if ( i%2 == 0 ){
          integral = integral + (4*h/3)*(get_value(x, p));
          x = x + h;
       }else{
          integral = integral + (2*h/3)*(get_value(x, p));
          x = x + h;
       }
     }
  
return(integral);

}

double 
get_value(double x, double p)
{
   double omega = .4;
   double lamda = .4;

   return( pow(x, p-1)*exp(-x) );

}