using System;

namespace szfx
{
	//二重積分
	public class integral
	{
		public delegate double Calculation(double x,double y);
		public integral()
		{
		}
		public static void run()
		{
			double pi = Math.PI;
			double s;
			s = ff(0,pi/2,0,pi/4,100);
			Console.WriteLine(s);
		}
		public static double ff(double a,double b,double c,double d,int div)
		{
			double I,h;
			int n;
			I = 0.0;
			n = 2 * div;
			h = (b-a)/n;			
			double[] X = new double[n+1];
			double[] F = new double[n+1];
			//等分x
			for(int i = 0 ; i < n ; i++)
			{
				X[i] = a + i*h;
			}
			X[n] = b;
			//求得F[k] = F([Xi])
			Calculation f = new Calculation(mysin);

			for(int k = 0 ; k <= n ; k++)
			{
				F[k] = _integral1(f,X[k],c,d,100);
			}
//			F[n] = _integral1(f,b,c,d,100);
			for(int i=0;i < div ;i++)
			{
				I += F[2*i] + 4*F[2*i+1] + F[2*i+2];
			}
			I = h / 3.0 * I;
			return I;
		}
		public static double _integral1(Calculation f,double xk,double start,double end,int div)
		{
			double s,h,y;
			int i;
			if((div <= 0) ||(end - start <=0.0))
				Console.WriteLine("Start is larger than , or div <= 0 in _integral!");
			s = (f(xk,start) + f(xk,end))/2.00;
			h = (end - start)/div;
			for(i = 1; i < div ; i++)
			{
				s = s+f(xk,start+i*h);
			}
			y = s*h;
			return y;
		}
		public static double mysin(double x,double y)
		{
			double r_f;
			r_f = Math.Sin(x+y);
			return r_f;
		}
	}
}