real L=5;
real h=2;

real omega=10.;
real c=1.;
real V=.3;

real k=omega/c;
real M=V/c;

int m=0;

real dn=15;

int n0=(k*h/pi)/sqrt(1-M^2);
cout<<"Nombre de modes propagatifs="<<n0<<endl;

// on tronque la serie pour T+ et T_ a n0
int nlimite=n0; 
real gamman=sqrt(2/h);

// Maillage sans la plaque
border Sigmaplus(t=0,h){x=L;y=t;label=2;};
border Dplus(t=L,-L){x=t;y=h;label=3;};
border Sigmamoins(t=h,0){x=-L;y=t;label=4;};
border Dmoins(t=-L,L){x=t;y=0;label=1;};



mesh Sh = buildmesh(Dmoins(2*L*dn) + Sigmaplus(h*dn)
+  Dplus(2*L*dn) + Sigmamoins(h*dn));

plot(Sh,wait=1);

// Construction d'EF P1

fespace Vh(Sh,P1);

// construction des matrices Tplus et Tmoins
////////////////////////////////////////////////////
// Traitement du bord x=L, construction de Tplus //
////////////////////////////////////////////////////

//////////////////////////////////////////////////
// On recupere la liste des noeuds sur le bord  //
//	Sigma moins (x=L)			//
//////////////////////////////////////////////////

// On compte le nombre d'aretes sur le bord Sigmaplus
int nbaretesbordplus=0;
for(int i=0;i<Sh.nt;i++){
	int nplus=0;
	for (int k=0;k<3;k++){
		if(Sh[i][k].x==L) {nplus++; }}
	if (nplus==2) {cout<<nbaretesbordplus<<endl; nbaretesbordplus++; }
};

cout<<"Nombre d'aretes sur Sigma+="<<nbaretesbordplus<<endl;

// On cree des tableau contenant la liste des aretes du bord (en fait deux tableaux, pour les deux extermites des aretes)
// remarque: on peut simplifier legerement cette partie avec les versions recentes de Freefem++.

int[int] aplus0(nbaretesbordplus);
int[int] aplus1(nbaretesbordplus);
int[int] eplus(nbaretesbordplus);

nbaretesbordplus=0;
for(int i=0;i<Sh.nt;i++){
	int nplus=0; int passurlebord=0;
	for (int k=0;k<3;k++){
		if(Sh[i][k].x==L) {nplus++;} else passurlebord=k;}
	if (nplus==2) {
		cout<<nbaretesbordplus<<endl;
		eplus[nbaretesbordplus]=i;
		if (passurlebord==0) {aplus0[nbaretesbordplus]=1; aplus1[nbaretesbordplus]=2;};
		if (passurlebord==1) {aplus0[nbaretesbordplus]=2; aplus1[nbaretesbordplus]=0;};
		if (passurlebord==2) {aplus0[nbaretesbordplus]=0; aplus1[nbaretesbordplus]=1;};
		nbaretesbordplus++;
		}
	}


// On affiche de resultat
cout<<"Liste des aretes sur le bord Sigma plus="<<endl;
// on en profite pour determiner le triangle inferieur
real miny=Sh[eplus[0]][aplus0[0]].y;
int Nymin=0;
for (int i=0;i<nbaretesbordplus;i++) {
	if (Sh[eplus[i]][aplus0[i]].y<miny) {miny=Sh[eplus[i]][aplus0[i]].y; Nymin=i;}
	cout<<"triangle "<<eplus[i]<<"="<<Sh[eplus[i]][aplus0[i]].x<<","<<Sh[eplus[i]][aplus0[i]].y<<" - "
					<<Sh[eplus[i]][aplus1[i]].x<<","<<Sh[eplus[i]][aplus1[i]].y<<endl;
}

// On construit la liste ordonnees des triangles du bord (oriente dans le sens croissant en y) //
int[int] tri(nbaretesbordplus);

Vh next=0;
for (int i=0;i<nbaretesbordplus;i++) {
	next[][Sh[eplus[i]][aplus0[i]]]=i; 
}

cout<<"next est defini"<<endl;

int j=Nymin;
for (int i=0;i<nbaretesbordplus;i++){
	tri[i]=j;
	j=next[][Sh[eplus[j]][aplus1[j]]];
}

cout<<"La liste triees des triangles sur le bord plus"<<endl;

for (int i=0;i<nbaretesbordplus;i++){
	cout<<Sh[eplus[tri[i]]]<<endl;
}

cout<<"On cree la liste des noeuds du bord sigma plus"<<endl;
int nbnoeudsbordplus=nbaretesbordplus+1;

int[int] Nnoeudbordplus(nbnoeudsbordplus);
Nnoeudbordplus[0]=Sh[eplus[tri[0]]][aplus0[tri[0]]];
for(int i=0;i<nbaretesbordplus;i++) Nnoeudbordplus[i+1]=Sh[eplus[tri[i]]][aplus1[tri[i]]];

cout<<"Noeuds du bord Sigma plus="<<endl;
for(int i=0;i<nbnoeudsbordplus;i++) cout<<Nnoeudbordplus[i]<<endl;

cout<<"test"<<endl;
Vh test=0;
for(int i=0;i<nbnoeudsbordplus;i++) test[][Nnoeudbordplus[i]]=i+1;
plot(test,wait=1);


///////////////////////////////////
//                               //
// Creation de la matrice T plus //
//                               //
///////////////////////////////////
complex[int,int] Tplus(nbnoeudsbordplus,nbnoeudsbordplus);
Tplus=0.;

real[int] pscnwi(nbnoeudsbordplus);
real y0,ym,yp;

// A completer 
// Debut


// Fin


////////////////////////////////////////////////////
// Traitement du bord y=-L, construction de Tmoins //
////////////////////////////////////////////////////


//////////////////////////////////////////////////
// On recupere la liste des noeuds sur le bord  //
//	Sigma moins (x=-L)			//
/////////////////////////////////////////////////

// On compte le nombre d'aretes sur le bord Sigmamoins
int nbaretesbordmoins=0;
for(int i=0;i<Sh.nt;i++){
	int nmoins=0;
	for (int k=0;k<3;k++){
		if(Sh[i][k].x==-L) {nmoins++; cout<<"x,y="<<Sh[i][k].x<<","<<Sh[i][k].y<<endl;}}
	if (nmoins==2) nbaretesbordmoins++;
};

cout<<"Nombre d'aretes sur Sigma-="<<nbaretesbordmoins<<endl;

// On cree des tableau contenant la liste des aretes du bord (en fait deux tableaux, pour les deux extermites de chaque arete)
// remarque: on peut simplifier legerement cette partie avec les versions recentes de Freefem++.

int[int] amoins0(nbaretesbordmoins);
int[int] amoins1(nbaretesbordmoins);
int[int] emoins(nbaretesbordmoins);

nbaretesbordmoins=0;
for(int i=0;i<Sh.nt;i++){
	int nmoins=0; int passurlebord=0;
	for (int k=0;k<3;k++){
		if(Sh[i][k].x==-L) {nmoins++;} else passurlebord=k;}
	if (nmoins==2) {
		emoins[nbaretesbordmoins]=i;
		if (passurlebord==0) {amoins0[nbaretesbordmoins]=1; amoins1[nbaretesbordmoins]=2;};
		if (passurlebord==1) {amoins0[nbaretesbordmoins]=2; amoins1[nbaretesbordmoins]=0;};
		if (passurlebord==2) {amoins0[nbaretesbordmoins]=0; amoins1[nbaretesbordmoins]=1;};
		nbaretesbordmoins++;
		}
	}


// On affiche de resultat
cout<<"Liste des aretes sur le bord Sigma moins="<<endl;
// on en profite pour determiner le triangle superieur
real maxy=Sh[emoins[0]][amoins0[0]].y;
int Nymax=0;
for (int i=0;i<nbaretesbordmoins;i++) {
	if (Sh[emoins[i]][amoins0[i]].y>maxy) {maxy=Sh[emoins[i]][amoins0[i]].y; Nymax=i;}
	cout<<"triangle "<<emoins[i]<<"="<<Sh[emoins[i]][amoins0[i]].x<<","<<Sh[emoins[i]][amoins0[i]].y<<" - "
					<<Sh[emoins[i]][amoins1[i]].x<<","<<Sh[emoins[i]][amoins1[i]].y<<endl;
}


// On construit la liste ordonnees des triangles du bord (oriente dans le sens croissant en y) //
int[int] trim(nbaretesbordmoins);

next=0;
for (int i=0;i<nbaretesbordmoins;i++) {
	next[][Sh[emoins[i]][amoins0[i]]]=i; 
}

cout<<"next est defini"<<endl;

 j=Nymax;
for (int i=0;i<nbaretesbordmoins;i++){
	trim[i]=j;
	j=next[][Sh[emoins[j]][amoins1[j]]];
}

cout<<"La liste trimees des trimangles sur le bord moins"<<endl;

for (int i=0;i<nbaretesbordmoins;i++){
	cout<<Sh[emoins[trim[i]]]<<endl;
}


cout<<"On cree la liste des noeuds du bord sigma moins"<<endl;
int nbnoeudsbordmoins=nbaretesbordmoins+1;

int[int] Nnoeudbordmoins(nbnoeudsbordmoins);
Nnoeudbordmoins[0]=Sh[emoins[trim[0]]][amoins0[trim[0]]];
for(int i=0;i<nbaretesbordmoins;i++) Nnoeudbordmoins[i+1]=Sh[emoins[trim[i]]][amoins1[trim[i]]];

cout<<"Noeuds du bord Sigma moins="<<endl;
for(int i=0;i<nbnoeudsbordmoins;i++) cout<<Nnoeudbordmoins[i]<<endl;

cout<<"test"<<endl;
test=0;
for(int i=0;i<nbnoeudsbordmoins;i++) test[][Nnoeudbordmoins[i]]=i+1;
plot(test,wait=1);


//////////////////////////////////////
//                                  //
// Creation de la matrice T moins  //
//                                  //
//////////////////////////////////////
complex[int,int] Tmoins(nbnoeudsbordmoins,nbnoeudsbordmoins);
Tmoins=0.;

real[int] pscnwim(nbnoeudsbordmoins);

for(int n=0;n<nlimite;n++){
// On cree un tableau contenant les valeurs                   
// des produits scalaires cn,w_I pour tout les noeuds du bord 
	// noeud inferieur
	y0=Sh[emoins[trim[0]]][amoins0[trim[0]]].y;
	yp=Sh[emoins[trim[0]]][amoins1[trim[0]]].y;

	if (n==0) {pscnwim[0]=1./sqrt(h)/2.*(yp-y0);} else
		pscnwim[0]=gamman*h^2/((pi*n)^2)*(cos(pi*n*yp/h)-cos(pi*n*y0/h))/(y0-yp);

	// noeuds intermediaires
	for(int i=1;i<nbnoeudsbordmoins-1;i++){
		ym=y0;
		y0=yp;
		yp=Sh[emoins[trim[i]]][amoins1[trim[i]]].y;
		if (n==0) {pscnwim[i]=1./sqrt(h)/2.*(yp-ym);} else
			pscnwim[i]=gamman*h^2/((pi*n)^2)*(cos(pi*n*yp/h)/(y0-yp)+cos(pi*n*y0/h)*(1./(y0-ym)-1./(y0-yp))-cos(pi*n*ym/h)/(y0-ym));
	};

	// noeud superieur
	ym=y0;
	y0=yp;
	if (n==0) {pscnwim[nbnoeudsbordplus-1]=1./sqrt(h)/2.*(y0-ym);} else
		pscnwim[nbnoeudsbordmoins-1]=gamman*h^2/((pi*n)^2)*(cos(pi*n*y0/h)-cos(pi*n*ym/h))/(y0-ym);

	real delta=k^2-(n*pi/h)^2*(1-M^2);
	complex sqrtdelta;
	if (delta<0) {sqrtdelta=sqrt(-delta)*(1.i);} else {sqrtdelta=sqrt(delta);}
	complex betamoinsn=(-k*M-sqrtdelta)/(1-M^2);
	for(int i=0;i<nbnoeudsbordmoins;i++)
		for(int j=0;j<nbnoeudsbordmoins;j++)
			Tmoins(i,j)+=betamoinsn*pscnwim[i]*pscnwim[j]*(1.i);
};

///////////////////////////////////////
// construction des matrices creuses //
///////////////////////////////////////

// petit hack pour pouvoir faire des matrices creuses de la taille voulue
// Inutile si une version recente de Freefem++ est utilisee
	bool maxplus=false;
	for (int i=0;i<Nnoeudbordplus.n;i++) if (Nnoeudbordplus[i]==Sh.nv-1){maxplus=true;};
	if (!maxplus){
		Nnoeudbordplus.resize(Nnoeudbordplus.n+1);
		Nnoeudbordplus[Nnoeudbordplus.n-1]=Sh.nv-1;
		Tplus.resize(Tplus.n+1,Tplus.n+1);
		Tplus(Tplus.n-1,Tplus.n-1)=0.;
	}
	bool maxmoins=false;
	for (int i=0;i<Nnoeudbordmoins.n;i++) if (Nnoeudbordmoins[i]==Sh.nv-1){maxmoins=true;};
	if (!maxmoins){
		Nnoeudbordmoins.resize(Nnoeudbordmoins.n+1);
		Nnoeudbordmoins[Nnoeudbordmoins.n-1]=Sh.nv-1;
		Tmoins.resize(Tmoins.n+1,Tmoins.n+1);
		Tmoins(Tmoins.n-1,Tmoins.n-1)=0.;
	}

matrix<complex> TTplus;
TTplus=Tplus(Nnoeudbordplus^-1,Nnoeudbordplus^-1);
cout<<"N bord plus="<<Nnoeudbordplus<<endl;
cout<<"valeurs bords plus="<<Tplus<<endl;

matrix<complex> TTmoins;
TTmoins=Tmoins(Nnoeudbordmoins^-1,Nnoeudbordmoins^-1);

varf mat(u,v)=int2d(Sh)((1-M^2)*dx(u)*dx(v)+dy(u)*dy(v)-2i*k*M*dx(u)*v-k^2*u*v);
matrix<complex> A=mat(Vh,Vh);

varf matb(u,v)=int1d(Sh,4)(u*v);
matrix<complex> B=matb(Vh,Vh);

// Construction du second membre
real Discriminant=k^2-(pi*m/h)^2*(1.-M^2);
complex deltabeta;
if (Discriminant<0) {deltabeta=(-2.i)*sqrt(abs(Discriminant))/(1.-M^2);} else deltabeta=-2.*sqrt(Discriminant)/(1.-M^2);
Vh<complex> g=(1.-M^2)*(1.i)*deltabeta*gamman*cos(pi*m*y/h);

Vh<complex> b=0.; // le second membre
b[]=B'*g[];

cout<<"A="<<A.n<<endl;
cout<<"TTplus="<<TTplus.n<<endl;
cout<<"TTmoins="<<TTmoins.n<<endl;

matrix<complex> C=A+(1.-M^2)*TTplus+(1.-M^2)*TTmoins;//A+(1-M^2)*(TTplus+TTmoins);


set(C,solver=UMFPACK);

Vh<complex> p;
p[]=C^-1*b[];

Vh pr=real(p);
Vh pim=imag(p);
plot(pr,wait=1);
plot(pim,wait=1);