//計算單個矩形孔的菲涅耳衍射
#include "math.h"
#include "conio.h"
#include "iostream.h"
#include "fstream.h"
//全局變量
#define	pi 3.1415927
#define lambda	0.0006328
#define A	1.0
#define B	100
#define C	0.0
#define D	1.0
//菲涅耳積分
double Fun_cos(double);
double Fun_sin(double);
double Fun_sinc(double);
//振幅函數
double Fun_ER(double,double,double,double,double,double);
double Fun_EV(double,double,double,double,double,double);
//物面光強
double Fun_E0(double,double);
//坐標變換
double Fun_xi_1(double,double,double);
double Fun_xi_2(double,double,double);
double Fun_eta_1(double,double,double);
double Fun_eta_2(double,double,double);
int main()
{
	double x0,y0,Lx,Ly,lx,ly;
	double ER[200][200];
	double EV[200][200];
	double I;
	int i,j;
	ofstream fout("data.txt");
	for(i=0;i<200;i++)
	{
		for(j=0;j<200;j++)
		{
			ER[i][j]=0;
			EV[i][j]=0;
		}
	}
//	cout<<"輸入中心點（x0,y0），矩形的邊長一半（Lx，Ly），觀察屏的半邊長lx，ly"<<endl;
//	cin>>x0>>y0>>Lx>>Ly>>lx>>ly;
//	cout<<endl;
//	cout<<"你輸入的數據是："<<endl;
//	cout<<"中心點坐標："<<"("<<x0<<","<<y0<<")"<<endl;
//	cout<<"矩形半邊長："<<"Lx="<<Lx<<",Ly="<<Ly<<endl;
//	cout<<"觀察屏的半邊長:"<<"lx="<<lx<<","<<"ly="<<ly<<endl;
	lx=10;ly=10;
	x0=0;
	y0=-5;
	Lx=5;
	Ly=0.005;
	while(y0>=-5&&y0<5)
	{
		y0=y0+0.005;
	//	Lx=sqrt(25-y0*y0);
		for(i=0;i<200;i++)
		{
			for(j=0;j<200;j++)
			{
				if((i<2*lx/0.1)&&(j<(2*ly/0.1)))
				{
					double tempR=Fun_ER(x0,y0,Lx,Ly,lx-i*0.1,ly-j*0.1);
					double tempV=Fun_EV(x0,y0,Lx,Ly,lx-i*0.1,ly-j*0.1);
					ER[i][j]=ER[i][j]+tempR;
					EV[i][j]=EV[i][j]+tempV;
					//cout<<I[i][j]<<endl;
				}

			}


		}
		y0=y0+0.005;
	}

	//cout<<Fun_I(x0,y0,Lx,Ly,0,0)<<endl<<Fun_I(x0,y0,Lx,Ly,-5,0)<<endl<<Fun_I(x0,y0,Lx,Ly,5,0)<<endl;
	if(!fout)
	{
		cout<<"Can't open out put file\n";
		return 1;
	}
	
	for(i=0;i<200;i++)
	{
		for(j=0;j<200;j++)
		{
			if((i<2*lx/0.1)&&(j<(2*ly/0.1)))
			{
				I=ER[i][j]*ER[i][j]+EV[i][j]*EV[i][j];

				fout<<I<<" ";
			}

		}
		fout<<endl;
		
	}

	fout.close();
	return 0;
}


//-------------------------------------------------------------------------------------
//計算像場的振幅函數
//輸入變量為：
//矩形物象的中心點（x_0，y_0）
//矩形邊長的一半L_x，L_y
//像點的坐標（x，y）
//輸出變量為像點的振幅實部
double Fun_ER(double x_0,double y_0,double L_x,double L_y,double x,double y)
{
	double E_R,C1,C2,S1,S2;
	C1=Fun_cos(Fun_xi_2(x_0,L_x,x))-Fun_cos(Fun_xi_1(x_0,L_x,x));
	S2=Fun_sin(Fun_eta_2(y_0,L_y,y))-Fun_sin(Fun_eta_1(y_0,L_y,y));
	C2=Fun_cos(Fun_eta_2(y_0,L_y,y))-Fun_cos(Fun_eta_1(y_0,L_y,y));
	S1=Fun_sin(Fun_xi_2(x_0,L_x,x))-Fun_sin(Fun_xi_1(x_0,L_x,x));
	E_R=C1*C2-S1*S2;
	return E_R;
}
//輸入變量為：
//矩形物象的中心點（x_0，y_0）
//矩形邊長的一半L_x，L_y
//像點的坐標（x，y）
//輸出變量為像點的振幅虛部
double Fun_EV(double x_0,double y_0,double L_x,double L_y,double x,double y)
{
	double E_V;
	double C1,C2,S1,S2;
	C1=Fun_cos(Fun_xi_2(x_0,L_x,x))-Fun_cos(Fun_xi_1(x_0,L_x,x));
	S2=Fun_sin(Fun_eta_2(y_0,L_y,y))-Fun_sin(Fun_eta_1(y_0,L_y,y));
	C2=Fun_cos(Fun_eta_2(y_0,L_y,y))-Fun_cos(Fun_eta_1(y_0,L_y,y));
	S1=Fun_sin(Fun_xi_2(x_0,L_x,x))-Fun_sin(Fun_xi_1(x_0,L_x,x));
	E_V=C1*S2+S1*C2;
	return E_V;
}
//-------------------------------------------------------------------------
//計算坐標為像點坐標是的輸入面的復振幅U0（x/A,y/A）
//輸入變量像點坐標（x，y）
//輸出變量為輸入面在該點的復振幅
//--------------------------------------------------------------
//計算正弦菲涅耳函數積分
double Fun_sin(double x)
{
	double return_sin;
	double mx;
	mx=fabs(x);
	if(mx<=sqrt(2))
	{
		return_sin=mx*sin(0.5567*exp(-(1.5545*mx-1.9941)*(1.5545*mx-1.9941)));
	}
	else
	{
		return_sin=0.5-(1-0.049*exp(-2*(mx-sqrt(2))))*cos(pi/2*mx*mx)/(pi*mx);
	}

	if(x<0)
	{
		return -return_sin;
	}
	else
	{
		return return_sin;
	}

}
//-----------------------------------------------------------------
//計算余弦菲涅耳函數積分
double Fun_cos(double x)
{
	double return_cos;
	double mx;
	mx=fabs(x);
	if(mx<=1)
	{
		return_cos=mx*cos(0.6855*mx*mx);
	}
	else
	{
		return_cos=0.5+(1-0.121*exp(-2*(mx-1)))/(pi*mx)*sin(pi/2*mx*mx);
	}
	if(x<0)
	{
		return -return_cos;
	}
	else
	{
		return return_cos;

	}
}
//----------------------------------------------------------------------
//計算sinc函數積分
double Fun_sinc(double x)
{
	double return_sinc;
	double mx;
	mx=fabs(x);


	if(mx<1)
	{
		return_sinc=(1-cos(pi*mx))/(3.39-3.39*pow(mx-1,1.6));
	}
	else
	{
		return_sinc=0.5+(1-0.1103*exp(-(mx-1)))/(pi*pi*mx)*cos(pi*(mx-1));
	}


	if(x<0)
	{
		return -return_sinc;
	}
	else
	{
		return return_sinc;
	}
}
//------------------------------------------------------------------------
//坐標變換
double Fun_xi_1(double x_0,double L_x,double x)
{
	return sqrt(2.0/fabs(lambda*B*A))*(A*(x_0-L_x)-x);
}
double Fun_xi_2(double x_0,double L_x,double x)
{
	return sqrt(2/fabs(lambda*B*A))*(A*(x_0+L_x)-x);
}
double Fun_eta_1(double y_0,double L_y,double y)
{
	return sqrt(2/fabs(lambda*B*A))*(A*(y_0-L_y)-y);
}
double Fun_eta_2(double y_0,double L_y,double y)
{
	return sqrt(2/fabs(lambda*B*A))*(A*(y_0+L_y)-y);
}