load "iovtk"; //to export the results in .vtk format
bool display=false;
int[int] fforder=[1];
int m=24; //parameter of the meshsize 
real L=2; 	//parameter of the geometry
real xRight=0; //position of the right slit
real xLeft=-2.5; //position of the left slit
real eps=0.05; //width of the slits
real w=0.8*pi; //wavenumber
real lMin=1.0;
real l2Min=1; // the expected critical length is 1.25=pi/w
real l2Max=1.34; 
real l1Min=1; // the expected critical length is 1.25=pi/w
real l1Max=1.34; 
real l1; //length of the first slit
real l2; //length of the second slit 
int start=1; //Loop on the lengths of the slits
int end1=10;
int end2=10;
int size1=end1-start+1;
int size2=end2-start+1;
real [int] Tablel1(size1*size2); // to store the results
real [int] Tablel2(size1*size2);
real [int] TableRreal(size1*size2);
real [int] TableRimag(size1*size2);
real [int] TableTpreal(size1*size2);
real [int] TableTpimag(size1*size2);
real [int] TableTmreal(size1*size2);
real [int] TableTmimag(size1*size2);
func  uifonc=exp(1i*w*x);
func  uiconjfonc=exp(-1i*w*x);
func  uiconjfoncH=exp(-1i*w*y);
/******************************************************
We make a double loop on the lengths of the two slits
/*****************************************************/
for (int parEps1=start;parEps1<=end1;parEps1++)
{
l1=l1Min+(parEps1-1)*(l1Max-l1Min)/(end1-start+1); //length of the first slit
for (int parEps2=start;parEps2<=end2;parEps2++)
{
l2=l2Min+(parEps2-1)*(l2Max-l2Min)/(end2-start+1); //length of the second slit
// Construction of the mesh
border a1(t=-L-2,eps/2){x =t;y=0; label = 0;}
border a2(t=0,1){x =eps/2;y=t; label = 0;}
border a3(t=0,l1){x =xRight+eps/2;y=1+t; label = 0;}
border a4(t=0,0.5-eps/2){x =xRight+eps/2+t;y=1+l1; label = 0;}
border a5(t=0,L-l1+lMin){x =xRight+0.5;y=1+l1+t; label = 0;}
border a6(t=0,1){x =xRight+0.5-t;y=1+lMin+L; label = 22;} // top right
border a7(t=0,L-l1+lMin){x =xRight-0.5;y=1+lMin+L-t; label = 0;}
border a8(t=0,0.5-eps/2){x =xRight-0.5+t;y=1+l1; label = 0;}
border a9(t=0,l1){x =xRight-eps/2;y=1+l1-t; label = 0;}
border a11(t=-xRight+eps/2,-xLeft-eps/2){x =-t;y=1; label = 0;} 
border a12(t=0,l2){x =xLeft+eps/2;y=1+t; label = 0;}
border a13(t=0,0.5-eps/2){x =xLeft+eps/2+t;y=1+l2; label = 0;}
border a14(t=0,L-l2+lMin){x =xLeft+0.5;y=1+l2+t; label = 0;}
border a15(t=0,1){x =xLeft+0.5-t;y=1+lMin+L; label = 23;} //top left
border a16(t=0,L-l2+lMin){x =xLeft-0.5;y=1+lMin+L-t; label = 0;}
border a17(t=0,0.5-eps/2){x =xLeft-0.5+t;y=1+l2; label = 0;}
border a18(t=0,l2){x =xLeft-eps/2;y=1+l2-t; label = 0;}
border a19(t=0,L-(-xLeft+eps/2)+2){x =xLeft-eps/2-t;y=1; label = 0;} 
border a20(t=0,1){x =-L-2;y=1-t; label = 21;} //entrance of the waveguide
mesh Th = buildmesh(a1(m*(L+2))+a2(m)+a3(m*l1)+a4(m*(0.5-eps/2))+a5(m*(L-l1+lMin))+a6(m)+a7(m*(L-l1+lMin))+a8(m*(0.5-eps/2))+a9(m*l1)+a11(m*(-xLeft+eps))+a12(m*l2)+a13(m*(0.5-eps/2))+a14(m*(L-l2+lMin))+a15(m*1)+a16(m*(L-l2+lMin))+a17(m*(0.5-eps/2))+a18(m*l2)+a19(m*(L-(-xLeft+eps)+2))+a20(m));
//plot(Th,wait=1,cmm="Press Enter to continue"); // uncomment to visualize the mesh
fespace Vh(Th,P2); // finite element space
/*************************************************
Construction of the Dirichlet-to-Neumann operators
*************************************************/
int nbfpro = 15; //we put nbfpro terms in the Dirichlet-to-Neumann
//Fourier basis
func complex expin(real x1,real x2, int n)
{
	if (n==0)
	return 1.;
	else
	return (sqrt(2.)*cos(n*pi*x2));	
}
/******************************************************
Left Dirichlet-to-Neumann 
/*****************************************************/
complex[int,int] vDtNG( Vh.ndof, nbfpro);
matrix<complex> DtNG;
for (int n=0;n<nbfpro;n++)
{
	func f= expin(x,y,n);
	varf FiniFourier(u,v) = int1d(Th,21)(v*f);
	complex[int] temp = FiniFourier(0,Vh);
	vDtNG(:,n)=temp;
}
DtNG=vDtNG;
// Construction of the diagonal matrix for the EFL
matrix<complex> DG;
complex[int] diagofDG(nbfpro);
for (int n =0;n<nbfpro;n++)
{
	if (n==0)
	diagofDG[n] = -1i*w; 
	else
	diagofDG[n] = sqrt((n*pi)^2-w^2); 
}
DG = [diagofDG];
matrix<complex> EFLG ;
EFLG = DtNG*DG;
EFLG= EFLG*DtNG';

/******************************************************
Top right Dirichlet-to-Neumann 
/*****************************************************/
//Fourier basis
func complex expinHD(real x1,real x2, int n)
{
	if (n==0)
	return 1.;
	else
	return (sqrt(2.)*cos(n*pi*(x1+0.5)));	
}
complex[int,int] vDtND( Vh.ndof, nbfpro);
matrix<complex> DtND;
for (int n=0;n<nbfpro;n++)
{
	func f= expinHD(x,y,n);
	varf FiniFourier(u,v) = int1d(Th,22)(v*f); 
	complex[int] temp = FiniFourier(0,Vh);
	vDtND(:,n)=temp;
}
DtND=vDtND;
matrix<complex> DD;
complex[int] diagofDD(nbfpro);
for (int n =0;n<nbfpro;n++)
{
	if (n==0)
	diagofDD[n] = -1i*w; 
	else
	diagofDD[n] = sqrt((n*pi)^2-w^2); 
}
DD = [diagofDD];
matrix<complex> EFLD ;
EFLD = DtND*DD;
EFLD= EFLD*DtND';


/******************************************************
Top left Dirichlet-to-Neumann 
/*****************************************************/
//Fourier basis
func complex expinH(real x1,real x2, int n)
{
	if (n==0)
	return 1.;
	else
	return (sqrt(2.)*cos(n*pi*(x1-(xLeft-0.5))));	
}

complex[int,int] vDtNH( Vh.ndof, nbfpro);
matrix<complex> DtNH;
for (int n=0;n<nbfpro;n++)
{
	func f= expinH(x,y,n);
	varf FiniFourier(u,v) = int1d(Th,23)(v*f); 
	complex[int] temp = FiniFourier(0,Vh);
	vDtNH(:,n)=temp;
}
DtNH=vDtNH;
matrix<complex> DH;
complex[int] diagofDH(nbfpro);
for (int n =0;n<nbfpro;n++)
{
	if (n==0)
	diagofDH[n] = -1i*w; 
	else
	diagofDH[n] = sqrt((n*pi)^2-w^2); 
}
DH = [diagofDH];
matrix<complex> EFLH ;
EFLH = DtNH*DH;
EFLH= EFLH*DtNH';
// End of the construction of the Dirichlet-to-Neumann operators

/***************************************
We solve the variational formulation 
The unknown is the total  field
****************************************/
Vh<complex> u,v,ui=uifonc;
Vh uiReal,uiImag,usReal,usImag,uReal,uImag;
matrix<complex> A,B,C;
complex[int] G(Vh.ndof);
varf aForme(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))-int2d(Th)(w^2*u*v);				
varf gForme(u,v) =  int1d(Th,21)(-2*1i*w*uifonc*v); //the source term coming from the incident field
A= aForme(Vh,Vh);
C= A+EFLD+EFLH+EFLG; //we add the terms involving the Dirichlet-to-Neumann
set(C,solver=UMFPACK);
G= gForme(0,Vh);
u[]=C^-1*G; // we solve the system
uReal=real(u);
plot(uReal,fill=1,dim=2,nbiso=40, value=30,wait=0,cmm="Real part of u");
Vh<complex> uconj=conj(u);
if (display==true)
{
uReal=real(u);
uImag=imag(u);
uiReal=real(ui);
uiImag=imag(ui);
//plot(usReal,fill=1,dim=2,nbiso=40, value=30,wait=1,cmm="Real part of us");
//plot(usImag,fill=1,dim=2,nbiso=40, value=10,wait=1,cmm="Imaginary part of us");
//plot(uReal,fill=1,dim=2,nbiso=40, value=20,wait=1,cmm="Real part of u");
//plot(uImag,fill=1,dim=2,nbiso=40, value=10,wait=1,cmm="Imaginary part of u");
/**
// If one wishes to store the results in vtk format (then to display with paraview)
savevtk("uScaReal.vtk",Th,usReal,order=fforder,dataname="solution");
savevtk("uScaImag.vtk",Th,usImag,order=fforder,dataname="solution");
savevtk("uTotReal.vtk",Th,uReal,order=fforder,dataname="solution");
savevtk("uTotImag.vtk",Th,uImag,order=fforder,dataname="solution");
***/
}

/***************************************
Computation of the scattering coefficients
****************************************/
Vh<complex> uiconj=conj(ui);
Vh<complex> us=u-ui;
complex Tp=int1d(Th,23)(u*uiconjfoncH); //left transmission coefficient
complex Tm=int1d(Th,22)(u*uiconjfoncH); //right transmission coefficient
complex R=int1d(Th,21)(us*uifonc); //reflection coefficient
cout<<"*****************************************************" <<endl;
cout<<"Progression = " << (parEps2-1+(parEps1-1)*size2)*100./(size1*size2) <<"%"<<endl;
cout<<"parEps1 = " << parEps1 <<endl;
cout<<"parEps2 = " << parEps2 <<endl;
cout<<"Coefficient R = " << R <<endl; //reflection coefficient
cout<<"Coefficient Tp = " << Tp <<endl; //first transmission coefficient
cout<<"Coefficient Tm = " << Tm <<endl; //second transmission coefficient
cout<<"Energy conservation |R|^2+|Tp|^2+|Tm|^2 = " << abs(R)*abs(R)+abs(Tp)*abs(Tp)+abs(Tm)*abs(Tm) <<endl;
cout<<"l1 = " << l1 <<endl; //length of the right slit 
cout<<"l2 = " << l2 <<endl; //length of the left slit 
cout<<"*****************************************************" <<endl;
TableRreal(parEps2-1+(parEps1-1)*size2)=real(R);
TableRimag(parEps2-1+(parEps1-1)*size2)=imag(R);
TableTpreal(parEps2-1+(parEps1-1)*size2)=real(Tp);
TableTpimag(parEps2-1+(parEps1-1)*size2)=imag(Tp);
TableTmreal(parEps2-1+(parEps1-1)*size2)=real(Tm);
TableTmimag(parEps2-1+(parEps1-1)*size2)=imag(Tm);
Tablel1(parEps2-1+(parEps1-1)*size2)=l1; 
Tablel2(parEps2-1+(parEps1-1)*size2)=l2;
}
}

/*********************************************************
We store the results in a file
********************************************************/
ofstream f1("Results.txt");
int nold=f1.precision(15);
for (int i=start;i<=end1;i++)
{	
for (int j=start;j<=end2;j++)	
{
f1 << Tablel1(j-1+(i-1)*size2) << " " << Tablel2(j-1+(i-1)*size2) << " " << TableRreal(j-1+(i-1)*size2)<< " " << TableRimag(j-1+(i-1)*size2) << " " << TableTpreal(j-1+(i-1)*size2)<< " " << TableTpimag(j-1+(i-1)*size2) << " " << TableTmreal(j-1+(i-1)*size2)<< " " << TableTmimag(j-1+(i-1)*size2)  << endl;
};
};