%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                                                      %% SEP: A Stable Election Protocol for clustered                        %%      heterogeneous wireless sensor networks                          %%                                                                      %% (c) Georgios Smaragdakis                                             %% WING group, Computer Science Department, Boston University           %%                                                                      %% You can find full documentation and related information at:          %% http://csr.bu.edu/sep                                                %%                                                                      %  % To report your comment or any bug please send e-mail to:             %% gsmaragd@cs.bu.edu                                                   %%                                                                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                                                      %% This is the SEP [1] code we have used.                               %%                                                                      %% [1]  Georgios Smaragdakis, Ibrahim Matta and Azer bestavros,         %      %      "SEP: A Stable Election Protocol for clustered                  %%      heterogeneous wireless sensor networks",                        %%      Second International Workshop on Sensor and Actor Network       %%      Protocols and Applications (SANPA 2004),Boston MA, August       %               					%      2004.                                                           %%                                                                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%Field Dimensions - x and y maximum (in meters)xm=100;ym=100;%x and y Coordinates of the Sinksink.x=0.5*xm;sink.y=0.5*ym;%Number of Nodes in the fieldn=100%Optimal Election Probability of a node%to become cluster headp=0.1;%Energy Model (all values in Joules)%Initial Energy Eo=0.5;%Eelec=Etx=ErxETX=50*0.000000001;ERX=50*0.000000001;%Transmit Amplifier typesEfs=10*0.000000000001;Emp=0.0013*0.000000000001;%Data Aggregation EnergyEDA=5*0.000000001;%Values for Hetereogeneity%Percentage of nodes than are advancedm=0.1;%\alphaa=1;%maximum number of roundsrmax=4999%%%%%%%%%%%%%%%%%%%%%%%%% END OF PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%Computation of dodo=sqrt(Efs/Emp); %Creation of the random Sensor Networkfigure(1);for i=1:1:n    S(i).xd=rand(1,1)*xm;    XR(i)=S(i).xd;    S(i).yd=rand(1,1)*ym;    YR(i)=S(i).yd;    S(i).G=0;    %initially there are no cluster heads only nodes    S(i).type='N';       temp_rnd0=i;    %Random Election of Normal Nodes    if (temp_rnd0>=m*n+1)         S(i).E=Eo;        S(i).ENERGY=0;        %%%%plot(S(i).xd,S(i).yd,'o');        hold on;    end    %Random Election of Advanced Nodes    if (temp_rnd0<m*n+1)          S(i).E=Eo*(1+a)        S(i).ENERGY=1;        %%%%plot(S(i).xd,S(i).yd,'+');         hold on;    endendS(n+1).xd=sink.x;S(n+1).yd=sink.y;%%%%plot(S(n+1).xd,S(n+1).yd,'x');            %First Iterationfigure(1);%counter for CHscountCHs=0;%counter for CHs per roundrcountCHs=0;cluster=1;countCHs;rcountCHs=rcountCHs+countCHs;flag_first_dead=0;for r=0:1:rmax    r  %Election Probability for Normal Nodes  pnrm=( p/ (1+a*m) );  %Election Probability for Advanced Nodes  padv= ( p*(1+a)/(1+a*m) );      %Operation for heterogeneous epoch  if(mod(r, round(1/pnrm) )==0)    for i=1:1:n        S(i).G=0;        S(i).cl=0;    end  end %Operations for sub-epochs if(mod(r, round(1/padv) )==0)    for i=1:1:n        if(S(i).ENERGY==1)            S(i).G=0;            S(i).cl=0;        end    end  end hold off;%Number of dead nodesdead=0;%Number of dead Advanced Nodesdead_a=0;%Number of dead Normal Nodesdead_n=0;%counter for bit transmitted to Bases Station and to Cluster Headspackets_TO_BS=0;packets_TO_CH=0;%counter for bit transmitted to Bases Station and to Cluster Heads %per roundPACKETS_TO_CH(r+1)=0;PACKETS_TO_BS(r+1)=0;figure(1);for i=1:1:n    %checking if there is a dead node    if (S(i).E<=0)        plot(S(i).xd,S(i).yd,'red .');        dead=dead+1;        if(S(i).ENERGY==1)            dead_a=dead_a+1;        end        if(S(i).ENERGY==0)            dead_n=dead_n+1;        end        hold on;        end    if S(i).E>0        S(i).type='N';        if (S(i).ENERGY==0)          plot(S(i).xd,S(i).yd,'o');        end        if (S(i).ENERGY==1)          plot(S(i).xd,S(i).yd,'+');        end        hold on;    endendplot(S(n+1).xd,S(n+1).yd,'x');STATISTICS(r+1).DEAD=dead;DEAD(r+1)=dead;DEAD_N(r+1)=dead_n;DEAD_A(r+1)=dead_a;%When the first node diesif (dead==1)    if(flag_first_dead==0)        first_dead=r        flag_first_dead=1;    endendcountCHs=0;cluster=1;for i=1:1:n   if(S(i).E>0)   temp_rand=rand;        if ( (S(i).G)<=0) %Election of Cluster Heads for normal nodes if( ( S(i).ENERGY==0 && ( temp_rand <= ( pnrm / ( 1 - pnrm * mod(r,round(1/pnrm)) )) ) )  )            countCHs=countCHs+1;            packets_TO_BS=packets_TO_BS+1;            PACKETS_TO_BS(r+1)=packets_TO_BS;                        S(i).type='C';            S(i).G=100;            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;                        %Calculation of Energy dissipated            distance;            if (distance>do)                S(i).E=S(i).E- ( (ETX+EDA)*(4000) + Emp*4000*( distance*distance*distance*distance ));             end            if (distance<=do)                S(i).E=S(i).E- ( (ETX+EDA)*(4000)  + Efs*4000*( distance * distance ));             end        end          %Election of Cluster Heads for Advanced nodes if( ( S(i).ENERGY==1 && ( temp_rand <= ( padv / ( 1 - padv * mod(r,round(1/padv)) )) ) )  )                    countCHs=countCHs+1;            packets_TO_BS=packets_TO_BS+1;            PACKETS_TO_BS(r+1)=packets_TO_BS;                        S(i).type='C';            S(i).G=100;            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;                        %Calculation of Energy dissipated            distance;            if (distance>do)                S(i).E=S(i).E- ( (ETX+EDA)*(4000) + Emp*4000*( distance*distance*distance*distance ));             end            if (distance<=do)                S(i).E=S(i).E- ( (ETX+EDA)*(4000)  + Efs*4000*( distance * distance ));             end        end             end  end endSTATISTICS(r+1).CLUSTERHEADS=cluster-1;CLUSTERHS(r+1)=cluster-1;%Election of Associated Cluster Head for Normal Nodesfor i=1:1:n   if ( S(i).type=='N' && S(i).E>0 )     if(cluster-1>=1)       min_dis=sqrt( (S(i).xd-S(n+1).xd)^2 + (S(i).yd-S(n+1).yd)^2 );       min_dis_cluster=1;       for c=1:1:cluster-1           temp=min(min_dis,sqrt( (S(i).xd-C(c).xd)^2 + (S(i).yd-C(c).yd)^2 ) );           if ( temp<min_dis )               min_dis=temp;               min_dis_cluster=c;           end       end              %Energy dissipated by associated Cluster Head            min_dis;            if (min_dis>do)                S(i).E=S(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis));             end            if (min_dis<=do)                S(i).E=S(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis));             end        %Energy dissipated        if(min_dis>0)            S(C(min_dis_cluster).id).E = S(C(min_dis_cluster).id).E- ( (ERX + EDA)*4000 );          PACKETS_TO_CH(r+1)=n-dead-cluster+1;         end       S(i).min_dis=min_dis;       S(i).min_dis_cluster=min_dis_cluster;              end endendhold on;countCHs;rcountCHs=rcountCHs+countCHs;%Code for Voronoi Cells%Unfortynately if there is a small%number of cells, Matlab's voronoi%procedure has some problems%[vx,vy]=voronoi(X,Y);%plot(X,Y,'r*',vx,vy,'b-');% hold on;% voronoi(X,Y);% axis([0 xm 0 ym]);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   STATISTICS    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                                                                     %%  DEAD  : a rmax x 1 array of number of dead nodes/round %  DEAD_A : a rmax x 1 array of number of dead Advanced nodes/round%  DEAD_N : a rmax x 1 array of number of dead Normal nodes/round%  CLUSTERHS : a rmax x 1 array of number of Cluster Heads/round%  PACKETS_TO_BS : a rmax x 1 array of number packets send to Base Station/round%  PACKETS_TO_CH : a rmax x 1 array of number of packets send to ClusterHeads/round%  first_dead: the round where the first node died                   %                                                                                     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%