/////////////////////////////////////////////////////////////////////////////////////////
//
//  projet 1: computation of the best reflectitivity spectrum with PSO optimization
//
//
///////////////////////////////////////////////////////////////////////////////////
//
//////////////direct problem computation  ///////////////
//
////////////////////////////////////////////////////////
//
function rhodot=f(z,rho,delta,lamb,zf,dn)         // right hand side of the ODE
d=splin(zf,dn);                                        
dnz=interp(z,zf,dn,d);
rhodot(1)=-2*delta*rho(2)-%pi*dnz*2*rho(1)*rho(2)/(lamb);
rhodot(2)=2*delta*rho(1)+%pi*dnz*(rho(1)^2-rho(2)^2)/(lamb)+%pi/lamb*dnz;
endfunction
//
////////////////////////////////////////////////////////////////////////////////////////
//
function err=fctbragg(x)
//
Nc=length(x);  // number of spline interpolation points of the apodisation function
//
n0=1.45;
L=0.01;
lambda0=1.55E-6;
Biglambda0=lambda0/(2*n0);
lambdamin=1.5497E-6;
lambdamax=1.5503E-6;
////////////////////////////////////////////////////////////
z0=L;rho0=[0;0];
z=L:(-L/(Ndis-1)):0;
zf=0:(L/(Nc-1)):L;
dn=-dnmax+2*dnmax*x;
//
lambda=lambdamin:(lambdamax-lambdamin)/(Np-1):lambdamax;
rtarget=zeros(1,Np);
ref=[];
////////////////////////////////////////////////////
for i=1:Np
  delta0=2*n0*%pi*1/lambda(i)-%pi/Biglambda0;     // delta=== beta
  if abs(lambda(i)-lambda0)<0.1E-9 then rtarget(i)=1;end
  rho=ode(rho0,z0,z,list(f,delta0,lambda0,zf,dn));
  R=rho(1,$)^2+rho(2,$)^2;
  ref=[ref,R];
end
err=(sum((abs(ref-rtarget)).^2))^(1/2);
endfunction;
////////////////////////////////////////////////////////
//
//////////////MAIN PROGRAM  //////////////////////
//
////////////////////////////////////////////////////////
//
//
n=evstr(x_dialog('parameter number (of the spline function)','4'));   // parameters number
//
// PSO parameters///////////////////

Vmax=evstr(x_dialog('maximal velocity of the swarm','1')); 
Npop=evstr(x_dialog('population number','14'));
Ngen=evstr(x_dialog('generation number','30'));
//
w=0.7;c1=0.7;c2=0.7;
xmin=0;xmax=1;
//
// Bragg parameters//////////////////
//
L=evstr(x_dialog('length of the FBG','0.01'));
Np=50;         // number of computed points of the spectrum
dnmax=2E-4;    // maximal value of the apodisation function
Ndis=200; // number of discretization points on [0,L] for the resolution of the ode
n0=1.45;
lambda0=1.55E-6;
Biglambda0=lambda0/(2*n0);
lambdamin=1.5497E-6;
lambdamax=1.5503E-6;
//
//
//initialisation//////////////////////
//
PSO=zeros(Npop,3*n+1);       
PSO(:,1:n)=xmin+(xmax-xmin)*rand(Npop,n);
PSO(:,(n+1:2*n))=-Vmax+2*Vmax*rand(Npop,n);
PSO(:,(2*n+1:3*n))=PSO(:,1:n);
PSO(:,3*n+1)=1E10;
pg=zeros(1,n+1);pg(n+1)=1E10;pgval=[];
for i=1:Ngen
  val=[];
  for j=1:Npop
    valj=fctbragg(PSO(j,1:n))
    val=[val;valj] // the evaluation of the swarm
    // computation of the new pi
    if (valj<PSO(j,3*n+1)) then
      PSO(j,2*n+1:3*n)=PSO(j,1:n);
      PSO(j,3*n+1)=valj;
    end 
  end
  // computation of pg
  [a,b]=gsort(PSO(:,3*n+1));ind=b($);
  pg=[PSO(ind,2*n+1:3*n),PSO(ind,3*n+1)];
  pgval=[pgval,pg(n+1)];
  // new velocitiy
PSO(:,(n+1:2*n))=w*PSO(:,(n+1:2*n))+c1*rand(Npop,n).*(PSO(:,2*n+1:3*n)-PSO(:,1:n))+c2*rand(Npop,n).*(ones(Npop,1)*pg(1:n)-PSO(:,1:n));
  // new value of x
  PSO(:,1:n)=PSO(:,1:n)+PSO(:,(n+1:2*n)); 
end
//
//
////////////////////////////////////////////////////////
//
//// Results display  ////////////////////////////////
//
////////////////////////////////////////////////////////
//
//
disp('best elemet found:');disp(pg(1:n))
xset('window',0);
clf()
plot2d(Npop*(1:Ngen),pgval)
xtitle('convergence history (best value vs evaluation number)')
/////////////////////
//
//
////////////////////////////////////////////////////////////
z0=L;rho0=[0;0];
z=L:(-L/(Ndis-1)):0;
zf=0:(L/(n-1)):L;
dn=-dnmax+2*dnmax*pg(1:n);     // optimal value
lambda=lambdamin:(lambdamax-lambdamin)/(Np-1):lambdamax;
rtarget=zeros(1,Np);
ref=[];
////////////////////////////////////////////////////
for i=1:Np
  delta0=2*n0*%pi*1/lambda(i)-%pi/Biglambda0;     // delta=== beta
  if abs(lambda(i)-lambda0)<0.1E-9 then rtarget(i)=1;end
  rho=ode(rho0,z0,z,list(f,delta0,lambda0,zf,dn));
  R=rho(1,$)^2+rho(2,$)^2;
  ref=[ref,R];
end
////////////////////////////////////////////////////
xset('window',1);
clf();
d=splin(zf,dn);
dnz=interp(0:L/Ndis:L,zf,dn,d);
plot2d(0:L/Ndis:L,dnz)
plot2d(zf,dn,-1);
xtitle('best apodisation function','z','dn(z)')
//////////////////////////////////////////////////////
xset('window',2);
clf();
plot2d(lambda,ref)
plot2d(lambda,rtarget,2)
xtitle('best reflectivity spectrum (black) and ideal spectrum (blue)','lambda','R(lambda)')
////////////////////////////////////////////////////////
error=(sum((abs(ref-rtarget)).^2))^(1/2);
disp('error function value:');disp(error)


