/***************************************************
Scattering of an incident wave for the Dirichlet pb
Delta u+k^2u=0 in \Om
u=0 on \partial\Om
u-ui ougoing
****************************************************/
/***************************************************
Physical parameters
****************************************************/
real k=1.8*pi; // frequency
int N=int(k/pi); //nb of propagating modes
real beta=sqrt(k^2-pi^2);
func uiF=exp(1i*beta*x)*sin(pi*y); //incident field
func uiFconj=exp(-1i*beta*x)*sin(pi*y); 
/***************************************************
Mesh construction
****************************************************/
real L=4; 	//domain [-L;L]x[0;1]
int Nmesh=20;
border aMesh(t=-L,L){x =t; y = 0; label = 0;}
border bMesh(t=0,1){x =L; y = t; label = 3;} //reference 3 at x=L 
border cMesh(t=0,L-1){x =L-t; y = 1; label = 0;}
border dMesh(t=-1,1){x =t; y = 1+0.4*sin(t*pi); label = 0;}
border eMesh(t=0,L-1){x =-1-t; y = 1; label = 0;} 
border fMesh(t=0,1){x =-L; y = 1-t; label = 2;} //reference 2  at x=-L
border inclusion(t=0,2*pi){x=0.3+0.15*cos(t);y=0.7+0.2*sin(t);label=0;}
mesh Th=buildmesh(aMesh(2*L*Nmesh)+bMesh(Nmesh)+cMesh(int((L-1)*Nmesh))+dMesh(-int(2*Nmesh+10))+eMesh(int((L-1)*Nmesh))+fMesh(Nmesh)+inclusion(-int(0.2*2*pi*Nmesh)));
plot(Th,wait=1,cmm="Press Enter to continue");
fespace Vh(Th,P2); // Finite element space P2

/***************************************************
Construction of the Dirichlet-to-Neumann
****************************************************/
int nbfpro = 15; //Number of terms in  the Dirichlet-to-Neumann series
//Fourier basis in the section
func complex expin(real x1,real x2, int n)
{	
	return (sqrt(2.)*sin(n*pi*x2));	
}
/***************************************************
Dirichlet-to-Neumann at x=-L
****************************************************/
complex[int,int] vDtNG(Vh.ndof, nbfpro);
matrix<complex> DtNG;
for (int n=0;n<nbfpro;n++)
{
	func f= expin(x,y,n+1);
	varf FiniFourier(u,v) = int1d(Th,2)(v*f);  // integration at x=-L
	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<N)
	diagofDG[n] = 1i*sqrt(k^2-pi^2); 
	else
	diagofDG[n] = -sqrt(((n+1)*pi)^2-k^2); 
}
DG = [diagofDG];
// Final construction 
matrix<complex> EFLG ;
EFLG = DtNG*DG;
EFLG= EFLG*DtNG';

/***************************************************
Dirichlet-to-Neumann at x=L
****************************************************/
complex[int,int] vDtND( Vh.ndof, nbfpro);
matrix<complex> DtND;
for (int n=0;n<nbfpro;n++)
{
	func f= expin(x,y,n+1);
	varf FiniFourier(u,v) = int1d(Th,3)(v*f); // integration at x=L
	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<N)
	diagofDD[n] = 1i*sqrt(k^2-pi^2); 
	else
	diagofDD[n] = -sqrt(((n+1)*pi)^2-k^2); 
}
DD = [diagofDD];
matrix<complex> EFLD ;
EFLD = DtND*DD;
EFLD= EFLD*DtND';

/***************************************************
Variational formulation for the total field
****************************************************/
Vh<complex> u,v;
matrix<complex> A,C;
complex[int] F(Vh.ndof);
varf aForme(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))-int2d(Th)(k^2*u*v)+on(0,u=0);
varf fForme(u,v)=-int1d(Th,2)(2*1i*beta*uiF*v); //the rightgoing incident wave is not outgoing at x=-L (referece =2)
A= aForme(Vh,Vh);
F= fForme(0,Vh);
C= A-EFLG-EFLD; //here we add the DTN
set(C,solver=UMFPACK);
u[]=C^-1*F; //we solve the system
Vh uReel=real(u);
plot(uReel,fill=1,dim=2,nbiso=40, value=10,wait=0,cmm="Real part of u");

/***************************************************
Computation of the reflexion/transmission coefficients
****************************************************/
Vh<complex> ui=uiF,us=u-ui;
complex R=2*int1d(Th,2)(us*uiF);
complex T=2*int1d(Th,3)(u*uiFconj);

/****************************************************
We display the results in the terminal
*****************************************************/
cout<<"*****************************************************" <<endl;
cout<<"Dirichlet problem" <<endl;
cout<<"k = " << k <<endl;
cout<<"N = " << N <<endl;
cout<<"Coefficient R = " << R <<endl;
cout<<"Coefficient T = " << T <<endl;
cout<<"Energy conservation |R|^2+|T|^2 = " << abs(R)*abs(R)+abs(T)*abs(T) <<endl;  //should be equal to one in monomode regime (k\in(pi;2pi))
cout<<"*****************************************************" <<endl;