// Simpsons Method for double integrals

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

#define true    1
#define false   0
#define PI  3.1415926535

// Y의 아래쪽 경계를 나타내는 함수
double lower(double x)
{
   return (x);
}

// Y의 위쪽 경계를 나타내는 함수
double upper(double x)
{
   return (2*x);
}

// 구하려는 함수
double func(double x, double y)
{
   return (pow((x),2.0)+pow((y),3.0));
}

void main()
{
   double h_out, sum_out=0.0, even_out=0.0, odd_out=0.0,
    	  x, h_in, sum_in, yLeft, yRight,
	  even_in, odd_in, y, temp_in, temp_out;
   double left=2.0;   // 바깥쪽 적분의 좌측경계
   double right=2.2;  // 바깥쪽 적분의 우측경계
   int inner_num=100; // 안쪽 적분의 구간 수   - M
   int outer_num=100; // 바깥쪽 적분의 구간 수 - N
   
   h_out = (right - left) / outer_num;         

   for (int i=0; i<=outer_num; i++) { 

      // x 에 대해서 Composite Simpson's Method 적용 
      x = left + i * h_out;
      yLeft = lower(x);
      yRight = upper(x);
      h_in = (yRight - yLeft) / inner_num;

      sum_in = func(x, yLeft) + func(x, yRight);
      even_in = 0.0;  
      odd_in = 0.0;   

      for (int j=1; j<=inner_num-1; j++) {
         y = yLeft + j * h_in;
	 temp_in = func(x, y);
	 if ((j%2) == 0) even_in = even_in + temp_in;
	 else odd_in = odd_in + temp_in; 
      }
      temp_out = (sum_in + 2.0 * even_in + 4.0 * odd_in) * h_in / 3.0;

      if ((i == 0) || (i == outer_num)) sum_out = sum_out + temp_out; 
      else {
         if ((i%2) == 0) even_out = even_out + temp_out;
	 else odd_out = odd_out + temp_out;
      }
   } 
   sum_out = (sum_out + 2.0 * even_out + 4.0 * odd_out) * h_out /3.0;

   // 결과 출력
   printf("Double integral of Function with Simpson's Method");
   printf("\nOuter Number : %3d", outer_num);
   printf("\nInner Number : %3d", inner_num);
   printf("\nLeft  End : %12.8f", left);
   printf("\nRight End : %12.8f ", right);
   printf("\nApproximation : %12.8f", sum_out);
}