//© Copyright - All right reserved - Guillaume BURGNIES & Quentin MORINO


clear
//SINGLEMODE   SMF--------------------------------------------------------------------------------------------------------------
MAX=[]
V0=2;
a0=4;                                                                            //Rayon fibre SMF (µm))
aMMF=60.5                                                                        //Rayon fibre MMF (µm)

nc=1.4440                                                                       //Indice du ceour
ngVar=1.33:1:1.44                                                               //Variation de l'indice de la gaine 
lambda=1.45:0.001:1.65                                                          //Variation de la longueur d'onde 
b=.001:1e-5:1
aMpl=1;

z0=0000;                                                                        //Variation (ou fixe avec nBz=1 Z0=0 et zMax=valeur voulue) de la distance z 
zMax=58945;                                                                     //
nBz=1+0000;                                                                     //
z=linspace(z0,zMax,nBz);                                                        //
		
U0=V0*sqrt(1-b)
W0=V0*sqrt(b)
g0=U0.*besselj(1,U0).*besselk(0,W0)-W0.*besselk(1,W0).*besselj(0,U0)
Res0=b(find(g0(2:length(b)).*g0(1:$-1)<0))';                                    //Calcul des constantes de propagation
//clear g0 W0 U0

dR=1e-2
R=0.00001:dR:1.5                                                                // rayon normalisé R=r/a, a rayon du coeur.de la MMF
Ures0=V0*sqrt(1-Res0);
Wres0=V0*sqrt(Res0);
	


// partie radiale des modes-------------------------------------------------------------------------------------------------------------------

rSc=R(R/a0*aMMF<1)*aMMF/a0;   //Rayon coeur SMF
rSg=R(R/a0*aMMF>=1)*aMMF/a0;  //Rayon gaine SMF
Psi0=aMpl*[besselj(0,Ures0*rSc)./(besselj(0,Ures0)*ones(rSc)) , besselk(0,Wres0*rSg)./(besselk(0,Wres0)*ones(rSg))];//Champ mode SMF
Int0=(Psi0.^2).*(ones(Res0)*R)*dR
SumInt0=sum(Int0,"c")



    for j=1:length(ngVar)
		ng=ngVar(j)
		dn=nc-ng;  //delta n
	
		for m=1:length(lambda) 
		    k0=(2*%pi)/lambda(m);
			V=k0*aMMF*sqrt(nc^2-ng^2);
			
	
//MULTIMODE   MMF---------------------------------------------------------------------------------------------------------------
			
			U=V*sqrt(1-b)
			W=V*sqrt(b)
			g=U.*besselj(1,U).*besselk(0,W)-W.*besselk(1,W).*besselj(0,U)
			
			Res=b(find(g(2:length(b)).*g(1:$-1)<0))';                           // constantes de propagation
			
//clear g U W
			
			Ures=V*sqrt(1-Res);
			Wres=V*sqrt(Res);
			
// partie radiale des modes---------------------------------------------------------------------------------------------------------------

			Psi=aMpl*[besselj(0,Ures*R(R<1))./(besselj(0,Ures)*ones(R(R<1))) , besselk(0,Wres*R(R>=1))./(besselk(0,Wres)*ones(R(R>=1)))];
			
			Int=(Psi.^2).*(ones(Res)*R)*dR
			SumInt=sum(Int,"c")
			
//calcul Cm--------------------------------------------------------------------------------------------------------------------------------
			    //Num_cm=sum(ones(Res)*(Psi0.*R).*Psi*dR,"c")			
			    //Den_cm=SumInt
			
			Cm=(sum(ones(Res)*(Psi0.*R).*Psi*dR,"c"))./SumInt
			
			
//Calcul E(r,0) et E(r,z)------------------------------------------------------------------------------------------------------------------
			
			
			Beta=k0*sqrt(Res*(nc^2-ng^2)+ng^2)		
		
			Er0=sum(Cm(:,$)*(ones(R)).*Psi.*(exp(%i*Beta*0)*ones(R)),"r") //Calcul intégrale E(r,0)
			Er02=(abs(Er0)).^2;	  //
			
			
//calcul pertes----------------------------------------------------------------------------------------------------------------------------
	
			for i=1:length(z)
			
				Erz=sum(((Cm(:,$).*exp(%i*Beta*z(i)))*ones(R)).*Psi,"r")  //Calcul E(r,z(i))
				Erz2=abs(Erz).^2; //Calcul intégrale E(r,z)(i)) au carré 
				num_i=(abs(sum(Erz.*Psi0.*R*dR,"c"))).^2    //calcul du numérateur
				den_i=sum(Erz2.*R*dR,"c").*sum((abs(Psi0).^2).*R*dR,"c")  //calcul du dénominateur //
				
			//	den_i=max(den_i, 1D-300)		// contre singularités
			//	rApp=max(num_i/den_i,1D-300)
			//	
				rApp=((abs(sum(Erz.*Psi0.*R*dR,"c"))).^2)/(sum(Erz2.*R*dR,"c").*sum((abs(Psi0).^2).*R*dR,"c"))	//num_i/den_i
				L(i)=10*log10(rApp)
			end
			iMaxSpec=find(L==max(L))
			zMaxIndice(j)=z(iMaxSpec)
			LMaxIndice(j)=L(iMaxSpec)
		
		
			MAX=[MAX;[lambda(m),L(iMaxSpec),z(iMaxSpec)]]
		
		end //for lambda
end //for ng
 
 
 //Calcul final pertes
//clear U W h

// SORTIES----------------------------------------------------------------------------------------------------------------------

//disp(Res0,"b0=")
//disp(SumInt0,"SumInt0=")
//disp(Res,"b=")
//disp(SumInt,"SumInt=")
//disp(Cm,"Cm=")
//disp(E)

scf(0);clf(0);

//tItreGraph='Valeur de V = '+string(V)


//subplot(511)
//plot(R',Psi')		                                                            // champ
//plot(R,zeros(R),'k-')
//xtitle(tItreGraph,'Rayon normalisé','Champ E')
//xtitle(tItreGraph,'Rayon normalisé','Champ E (Courbe bleue) et Psi 0 (Courbe Verte)')


//subplot(512)                                                                  // intensité
//plot(R',Psi'.^2)
//xtitle(tItreGraph,'rayon normalisé','intensité Psi^2')


//subplot(513)
//plot(b,log10(abs(g)))                                                         //Résolution b
//xtitle(tItreGraph,'b','log(abs(g))')
//xstring(0,-5,(string(Res')));


//subplot(514)
//plot(Res($:-1:1),Cm','*')                                                     //b en fontion de Cm
//plot(Cm($:-1:1).^2,'d--')                                                     //Cm²
//xtitle(tItreGraph,'b','Cm')


//subplot(515)
//plot(R', 10*log10([Er0;Psi0]'));
//xtitle(tItreGraph,'Rayon normalisé','Er0 (bleu) et Psi0 (vert)')
//plot(R', ((Er0-Psi0)./Psi0)');
//plot(R',Psi0,'r');


scf(0);clf(0);
//plot (z,L','k-');
plot(MAX(:,1),MAX(:,2),'b.-')