clear
%1.初始參數(shù)設(shè)定模塊
%.傳感器節(jié)點區(qū)域界限(單位 M)
xm=100;
ym=100;
%(1)匯聚節(jié)坐標(biāo)給定
sink.x=0.5*xm;
sink.y=0.5*ym;
%區(qū)域內(nèi)傳器節(jié)數(shù)
n=100
%簇頭優(yōu)化比例（當(dāng)選簇頭的概率）
p=0.1;
%能量模型（單位 焦）
%初始化能量模型
Eo=0.5;
%Eelec=Etx=Erx
ETX=50*0.000000001;
ERX=50*0.000000001;
%Transmit Amplifier types
Efs=10*0.000000000001;
Emp=0.0013*0.000000000001;
%Data Aggregation Energy
EDA=5*0.000000001;
%高能量節(jié)點超出一節(jié)點能量的百分比
a=1;
%最大循環(huán)次數(shù)
rmax=5000
%算出參數(shù) do
do=sqrt(Efs/Emp);
%2.無線傳感器網(wǎng)絡(luò)模型產(chǎn)生模塊
%構(gòu)建無線傳感器網(wǎng)絡(luò),在區(qū)域內(nèi)均勻投放100個節(jié)點,并畫出圖形
for i=1:1:n
    S1(i).xd=rand(1,1)*xm;
    S2(i).xd=S1(i).xd;
    XR2(i)=S2(i).xd;
    XR1(i)=S1(i).xd;
    S1(i).yd=rand(1,1)*ym;
    S2(i).yd=S1(i).yd;
    YR2(i)=S2(i).yd;
    YR1(i)=S1(i).yd;
    S1(i).G=0;
    S2(i).G=0;
    S1(i).E=Eo*(1+rand*a);
    S2(i).E=S1(i).E;
    %initially there are no cluster heads only nodes
    S1(i).type='N';
    S2(i).type='N';
end

S1(n+1).xd=sink.x;
S1(n+1).yd=sink.y;
S2(n+1).xd=sink.x;
S2(n+1).yd=sink.y;
%3.網(wǎng)絡(luò)運行模塊
%簇頭節(jié)點數(shù)
countCHs1=0;
cluster1=1;%此定義的目的僅僅是給定一個1開始的下標(biāo)參數(shù)，真正的簇頭數(shù)應(yīng)該還減去1
flag_first_dead1=0;
flag_teenth_dead1=0;
flag_all_dead1=0;
%死亡節(jié)點數(shù)
dead1=0;
first_dead1=0;
teenth_dead1=0;
all_dead1=0;
%活動節(jié)點數(shù)
allive1=n;
%counter for bit transmitted to Bases Station and to Cluster Heads
packets_TO_BS1=0;
packets_TO_CH1=0;
%(1)循環(huán)模式設(shè)定
for r=0:1:rmax     %該 for 循環(huán)將下面的所有程序包括在內(nèi)，直到最后一 end 才結(jié)束循環(huán)
    r
  %每過一個輪轉(zhuǎn)周期(本程序為10次)使各節(jié)點的S（i）.G參數(shù)（該參數(shù)用于后面的簇選舉，在該輪轉(zhuǎn)周期內(nèi)已當(dāng)選過簇頭的節(jié)點不能再當(dāng)選）恢復(fù)為零
  if(mod(r, round(1/p) )==0)
    for i=1:1:n
        S1(i).G=0;
        S1(i).cl=0;
    end
  end
%(2)死亡節(jié)點檢查模塊
dead1=0;
for i=1:1:n
    %檢查有無死亡節(jié)點
    if (S1(i).E<=0)
        dead1=dead1+1; 
        %(3)第一個死亡節(jié)點的產(chǎn)生時間(用輪次表示)
        %第一個節(jié)點死亡時間
        if (dead1==1)
           if(flag_first_dead1==0)
              first_dead1=r;
              flag_first_dead1=1;
           end
        end
        %10%的節(jié)點死亡時間
        if(dead1==0.1*n)
           if(flag_teenth_dead1==0)
              teenth_dead1=r;
              flag_teenth_dead1=1;
           end
        end
        if(dead1==n)
           if(flag_all_dead1==0)
              all_dead1=r;
              flag_all_dead1=1;
           end
        end
    end
    if S1(i).E>0
        S1(i).type='N';
    end
end
STATISTICS.DEAD1(r+1)=dead1;
STATISTICS.ALLIVE1(r+1)=allive1-dead1;
%(4)簇頭選舉模塊
countCHs1=0;
cluster1=1;
for i=1:1:n
   if(S1(i).E>0)
   temp_rand=rand;     
   if ( (S1(i).G)<=0)  
       %簇頭的選舉，當(dāng)選的簇頭會把各種相關(guān)信存入下面程序所給定的變量中
        if(temp_rand<= (p/(1-p*mod(r,round(1/p)))))
            countCHs1=countCHs1+1;
            packets_TO_BS1=packets_TO_BS1+1;
            PACKETS_TO_BS1(r+1)=packets_TO_BS1;
             S1(i).type='C';
            S1(i).G=round(1/p)-1;
            C1(cluster1).xd=S1(i).xd;
            C1(cluster1).yd=S1(i).yd;
           distance=sqrt( (S1(i).xd-(S1(n+1).xd) )^2 + (S1(i).yd-(S1(n+1).yd) )^2 );
            C1(cluster1).distance=distance;
            C1(cluster1).id=i;
            X1(cluster1)=S1(i).xd;
            Y1(cluster1)=S1(i).yd;
            cluster1=cluster1+1;
           %計算簇頭發(fā)送4000bit數(shù)據(jù)到基站的能量消耗（這里應(yīng)是所有節(jié)點包括簇頭每一輪發(fā)送4000bit數(shù)據(jù)）
           distance;
            if (distance>do)
                S1(i).E=S1(i).E- ( (ETX+EDA)*(4000) + Emp*4000*( distance*distance*distance*distance )); 
            end
            if (distance<=do)
                S1(i).E=S1(i).E- ( (ETX+EDA)*(4000)  + Efs*4000*( distance * distance )); 
            end
        end     
    
    end
    % S(i).G=S(i).G-1;  
   
 end 
end
STATISTICS.COUNTCHS1(r+1)=countCHs1;
%(5)簇內(nèi)成員選擇簇頭模塊(即簇的形成模塊)
%簇內(nèi)成員對簇頭的選擇（即簇的形成）算法
for i=1:1:n
   if ( S1(i).type=='N' && S1(i).E>0 )
     if(cluster1-1>=1)
       min_dis=Inf;
       min_dis_cluster=0;
       for c=1:1:cluster1-1
           temp=min(min_dis,sqrt( (S1(i).xd-C1(c).xd)^2 + (S1(i).yd-C1(c).yd)^2 ) );
           if ( temp<min_dis )
               min_dis=temp;
               min_dis_cluster=c;
           end
       end
       %簇內(nèi)節(jié)點（發(fā)送4000bit數(shù)據(jù)）能量消耗
       
            min_dis;
            if (min_dis>do)
                S1(i).E=S1(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S1(i).E=S1(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
        %簇頭（接受和融合這一簇內(nèi)節(jié)點4000bit數(shù)據(jù)）的能量消耗
            S1(C1(min_dis_cluster).id).E = S1(C1(min_dis_cluster).id).E- ( (ERX + EDA)*4000 ); 
            packets_TO_CH1=packets_TO_CH1+1;    
 
        S1(i).min_dis=min_dis;
        S1(i).min_dis_cluster=min_dis_cluster;
    else
        min_dis=sqrt( (S1(i).xd-S1(n+1).xd)^2 + (S1(i).yd-S1(n+1).yd)^2 );
            if (min_dis>do)
                S1(i).E=S1(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S1(i).E=S1(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
            packets_TO_BS1=packets_TO_BS1+1;
     end
  end
end
STATISTICS.PACKETS_TO_CH1(r+1)=packets_TO_CH1;
STATISTICS.PACKETS_TO_BS1(r+1)=packets_TO_BS1;
end



%3.網(wǎng)絡(luò)運行模塊
%簇頭節(jié)點數(shù)
countCHs2=0;
cluster=1;%此定義的目的僅僅是給定一個1開始的下標(biāo)參數(shù)，真正的簇頭數(shù)應(yīng)該還減去1
flag_first_dead2=0;
flag_teenth_dead2=0;
flag_all_dead2=0;
%死亡節(jié)點數(shù)
dead2=0;
first_dead2=0;
teenth_dead2=0;
all_dead2=0;
%活動節(jié)點數(shù)
allive2=n;
%counter for bit transmitted to Bases Station and to Cluster Heads
packets_TO_BS2=0;
packets_TO_CH2=0;
%(1)循環(huán)模式設(shè)定
for r=0:1:rmax     %該 for 循環(huán)將下面的所有程序包括在內(nèi)，直到最后一 end 才結(jié)束循環(huán)
    r
  %每過一個輪轉(zhuǎn)周期(本程序為10次)使各節(jié)點的S（i）.G參數(shù)（該參數(shù)用于后面的簇選舉，在該輪轉(zhuǎn)周期內(nèi)已當(dāng)選過簇頭的節(jié)點不能再當(dāng)選）恢復(fù)為零
  if(mod(r, round(1/p) )==0)
    for i=1:1:n
        S2(i).G=0;
        S2(i).cl=0;
    end
  end
%(2)死亡節(jié)點檢查模塊
dead2=0;
for i=1:1:n
    %檢查有無死亡節(jié)點
    if (S2(i).E<=0)
        dead2=dead2+1; 
        %(3)第一個死亡節(jié)點的產(chǎn)生時間(用輪次表示)
        %第一個節(jié)點死亡時間
        if (dead2==1)
           if(flag_first_dead2==0)
              first_dead2=r;
              flag_first_dead2=1;
           end
        end
        %10%的節(jié)點死亡時間
        if(dead2==0.1*n)
           if(flag_teenth_dead2==0)
              teenth_dead2=r;
              flag_teenth_dead2=1;
           end
        end
        if(dead2==n)
           if(flag_all_dead2==0)
              all_dead2=r;
              flag_all_dead2=1;
           end
        end
    end
    if S2(i).E>0
        S2(i).type='N';
    end
end
STATISTICS.DEAD2(r+1)=dead2;
STATISTICS.ALLIVE2(r+1)=allive2-dead2;
%(4)簇頭選舉模塊
countCHs2=0;
cluster2=1;
for i=1:1:n
   if(S2(i).E>0)
   temp_rand=rand;     
   if ( (S2(i).G)<=0)  
       %簇頭的選舉，當(dāng)選的簇頭會把各種相關(guān)信存入下面程序所給定的變量中
        if(temp_rand<= (p/(1-p*mod(r,round(1/p)))))
            countCHs2=countCHs2+1;
            packets_TO_BS2=packets_TO_BS2+1;
            PACKETS_TO_BS2(r+1)=packets_TO_BS2;
             S2(i).type='C';
            S2(i).G=round(1/p)-1;
            C2(cluster2).xd=S2(i).xd;
            C2(cluster2).yd=S2(i).yd;
           distance=sqrt( (S2(i).xd-(S2(n+1).xd) )^2 + (S2(i).yd-(S2(n+1).yd) )^2 );
            C2(cluster2).distance=distance;
            C2(cluster2).id=i;
            X2(cluster2)=S2(i).xd;
            Y2(cluster2)=S2(i).yd;
            cluster2=cluster2+1;
           %計算簇頭發(fā)送4000bit數(shù)據(jù)到基站的能量消耗（這里應(yīng)是所有節(jié)點包括簇頭每一輪發(fā)送4000bit數(shù)據(jù)）
           distance;
            if (distance>do)
                S2(i).E=S2(i).E- ( (ETX+EDA)*(4000) + Emp*4000*( distance*distance*distance*distance )); 
            end
            if (distance<=do)
                S2(i).E=S2(i).E- ( (ETX+EDA)*(4000)  + Efs*4000*( distance * distance )); 
            end
        end     
    
    end
    % S2(i).G=S2(i).G-1;  
   
 end 
end
STATISTICS.COUNTCHS2(r+1)=countCHs2;
%(5)簇內(nèi)成員選擇簇頭模塊(即簇的形成模塊)
%簇內(nèi)成員對簇頭的選擇（即簇的形成）算法
for i=1:1:n
   if ( S2(i).type=='N' && S2(i).E>0 )
     if(cluster2-1>=1)
       min_dis=sqrt( (S2(i).xd-S2(n+1).xd)^2 + (S2(i).yd-S2(n+1).yd)^2 );
       min_dis_cluster=0;
       for c=1:1:cluster2-1
           temp=min(min_dis,sqrt( (S2(i).xd-C2(c).xd)^2 + (S2(i).yd-C2(c).yd)^2 ) );
           if ( temp<min_dis )
               min_dis=temp;
               min_dis_cluster=c;
           end
       end
       %簇內(nèi)節(jié)點（發(fā)送4000bit數(shù)據(jù)）能量消耗
       if(min_dis_cluster~=0)    
            min_dis;
            if (min_dis>do)
                S2(i).E=S2(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S2(i).E=S2(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
        %簇頭（接受和融合這一簇內(nèi)節(jié)點4000bit數(shù)據(jù)）的能量消耗
            S2(C2(min_dis_cluster).id).E = S2(C2(min_dis_cluster).id).E- ( (ERX + EDA)*4000 ); 
            packets_TO_CH2=packets_TO_CH2+1;
       else 
            min_dis;
            if (min_dis>do)
                S2(i).E=S2(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S2(i).E=S2(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
            packets_TO_BS2=packets_TO_BS2+1;
            
       end
        S2(i).min_dis=min_dis;
        S2(i).min_dis_cluster=min_dis_cluster;
    else
        min_dis=sqrt( (S2(i).xd-S2(n+1).xd)^2 + (S2(i).yd-S2(n+1).yd)^2 );
            if (min_dis>do)
                S2(i).E=S2(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S2(i).E=S2(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
            packets_TO_BS2=packets_TO_BS2+1;
     end
  end
end
STATISTICS.PACKETS_TO_CH2(r+1)=packets_TO_CH2;
STATISTICS.PACKETS_TO_BS2(r+1)=packets_TO_BS2;
end
r=0:5000;
subplot(1,2,1);
plot(r,STATISTICS.DEAD1,'--k',r,STATISTICS.DEAD2,'-r');
legend('leachm1','leachm2');
xlabel('x(time)');
ylabel('y(dead)');
title('\bf leachm1和leachm2的生命周期對比');
subplot(1,2,2);
plot(r,STATISTICS.PACKETS_TO_BS1,'--k',r,STATISTICS.PACKETS_TO_BS2,'-r');
legend('leachm1','leachm2');
xlabel('x(time)');
ylabel('y(data)');
title('\bf leachsm和leachm2的數(shù)據(jù)傳輸對比');