

clear;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARAMETRES %%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Dimensions du terrain - Maximum Axe des abscisses (X) et Axe des ordonnées (Y) (en metre)
xm=50;
ym=50;
warning off

%coordonnées du puits (station base)
puits.x=0.5*xm;
puits.y=-2*ym;

%Nombre de noeuds (capteurs)
n=100

%Probabilité optimale de l'éléction d'un noeud comme étant CH
p=0.05;

%Modéle énergetique toutes les valeurs sont en JOULES
%Energie initial (de chaque noeud):
E0=0.5;
%longueur des paquets ou messages transmit
k=2000;
%energie de transmission (Eelec) = énergie nécessaire pour émettre (Etx) =  énergie nécessaire pour recevoir (ERx)
ETX=50*0.000000001;
ERX=50*0.000000001;
%Energie d'amplification
Eamp=100*0.000000000001;

%energie nécessaire pour l'agrégation des données
EAg=5*0.000000001;
%Le nombre maximum de tours
max_tours=1300

%%%%%%%%%%%%%%%%%%%%%%%%% Fin des parametres %%%%%%%%%%%%%%%%%%%%%%%%

%Creation du RCSF
figure(1);
for i=1:1:n
    S(i).xd=rand(1,1)*xm;    
    S(i).yd=rand(1,1)*ym;    
    S(i).G=0;
    %initialement, il n'y a aucun CH
    S(i).type='N';     
    S(i).E=E0;
    plot(S(i).xd,S(i).yd,'o');
    hold on;      
end

% Création du puits (Station de Base)
S(n+1).xd=puits.x;
S(n+1).yd=puits.y;
plot(S(n+1).xd,S(n+1).yd,'x');
    
        
%Premiere Iteration
figure(1);

%compteur de CHs
cpt_CHs=0;
%compteur de CHs par tour
cpt_tour_CHs=0;
cluster=1;

cpt_CHs;
cpt_tour_CHs=cpt_tour_CHs+cpt_CHs;
indic_premier_mort=0;

for tour=0:1:max_tours
    tour

  %Ajustement de la probabilite d'election d'un CH
  if(mod(tour, round(1/p) )==0)
    for i=1:1:n
       S(i).G=0;       
    end
  end  

hold off;

%Nombre de noeuds morts
mort=0;

%compteurs de bits transmis vers la Station de Base et les CHs
messages_pour_SB=0;
messages_pour_CH=0;
%compteurs de bits transmis vers la Station de Base et les CHs par tour
MESSAGES_POUR_SB(tour+1)=0;
MESSAGES_POUR_CH(tour+1)=0;

figure(1);

for i=1:1:n
    %verification d'existence de noeuds morts
    if (S(i).E<=0)
        plot(S(i).xd,S(i).yd,'red .');
        mort=mort+1;
        S(i).type='N'; % fin de mission        
        hold on;    
    end
    if S(i).E>0
        S(i).type='N';
        plot(S(i).xd,S(i).yd,'o');
        hold on;
    end
end
plot(S(n+1).xd,S(n+1).yd,'x');


STATISTIQUE(tour+1).MORTS=mort;
MORTS(tour+1)=mort;

%Mort du premier noeud
if (mort==1)
    if(indic_premier_mort==0)
        premier_mort=tour
        indic_premier_mort=1;
    end
end

cpt_CHs=0;
cluster=1;
for i=1:1:n
   if(S(i).E>0)
   nombre_aleatoire=rand;     
   if ((S(i).G)<=0)

 %Election des CHs
  if(nombre_aleatoire<= (p/(1-p*mod(tour,round(1/p)))))&&(cluster<=(n*p))
             cpt_CHs=cpt_CHs+1;
            messages_pour_SB=messages_pour_SB+1;
            MESSAGES_POUR_SB(tour+1)=messages_pour_SB;
            
            S(i).type='C';
            S(i).G=round(1/p)-1;
            C(cluster).xd=S(i).xd;
            C(cluster).yd=S(i).yd;
            plot(S(i).xd,S(i).yd,'k*');
            
            distance=sqrt( (S(i).xd-(S(n+1).xd) )^2 + (S(i).yd-(S(n+1).yd) )^2 );
            C(cluster).distance=distance;
            C(cluster).id=i;
            X(cluster)=S(i).xd;
            Y(cluster)=S(i).yd;
            cluster=cluster+1;
            
            %Calcul de la quantite d'energie dissipee
            distance;         
                S(i).E=S(i).E- ( (ETX+EAg)*k  + Eamp*k*( distance * distance ));            
        end     
    
    end
  end 
end

STATISTIQUE(tour+1).CLUSTERHEADS=cluster-1;
CLUSTERHS(tour+1)=cluster-1;

%Ralliement des CHs
for i=1:1:n
   if ( S(i).type=='N' && S(i).E>0 )
     if(cluster-1>=1)
       distance_minimale=sqrt( (S(i).xd-S(n+1).xd)^2 + (S(i).yd-S(n+1).yd)^2 );
       plus_proche_cluster=1;
       for c=1:1:cluster-1
           temp=min(distance_minimale,sqrt( (S(i).xd-C(c).xd)^2 + (S(i).yd-C(c).yd)^2 ) );
           if ( temp<distance_minimale )
               distance_minimale=temp;
               plus_proche_cluster=c;
           end
       end
       
       %Energie dissipee par le CH associe
            distance_minimale;     
                S(i).E=S(i).E- ( ETX*(k) + Eamp*k*( distance_minimale * distance_minimale));         
        %Energie dissipee
        if(distance_minimale>0)
          S(C(plus_proche_cluster).id).E = S(C(plus_proche_cluster).id).E- ( (ERX + EAg)*k ); 
         MESSAGES_POUR_CH(tour+1)=n-mort-cluster+1; 
        end

       S(i).distance_minimale=distance_minimale;
       S(i).plus_proche_cluster=plus_proche_cluster;
           
   end
 end
end
hold on;

cpt_CHs;
cpt_tour_CHs=cpt_tour_CHs+cpt_CHs;

%Code de division en cellules de Voronoi

[vx,vy]=voronoi(X,Y);
plot(X,Y,'r*',vx,vy,'b-');
 hold on;
 voronoi(X,Y);
 axis([0 xm 0 ym]);

end

%Creation du graphe representant le nombre de noeuds morts au fil des tours
figure(2);
plot(MORTS);
hold on;