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
if countCHs2~=0
   u2=(n-y2)/countCHs2;
for c=1:1:cluster2-1
    z2=(x2(c)-u2)*(x2(c)-u2)+z2;
end
LBF2(r+1)=z2/countCHs2;
else  LBF2(r+1)=0;
end
STATISTICS.PACKETS_TO_CH2(r+1)=packets_TO_CH2;
STATISTICS.PACKETS_TO_BS2(r+1)=packets_TO_BS2;
end


d1=0.765*xm/2;
K=sqrt(0.5*n*do/pi)*xm/d1^2;
d2=xm/sqrt(2*pi*K);
Er=4000*(2*n*ETX+n*EDA+K*Emp*d1^4+n*Efs*d2^2);
S3(n+1).xd=sink.x;
S3(n+1).yd=sink.y;
S4(n+1).xd=sink.x;
S4(n+1).yd=sink.y;
%3.網絡運行模塊
%簇頭節點數
countCHs3=0;
cluster3=1;%此定義的目的僅僅是給定一個1開始的下標參數，真正的簇頭數應該還減去1
flag_first_dead3=0;
flag_teenth_dead3=0;
flag_all_dead3=0;
%死亡節點數
dead3=0;
first_dead3=0;
teenth_dead3=0;
all_dead3=0;
%活動節點數
allive3=n;
%counter for bit transmitted to Bases Station and to Cluster Heads
packets_TO_BS3=0;
packets_TO_CH3=0;
%(1)循環模式設定
for r=0:1:rmax     %該 for 循環將下面的所有程序包括在內，直到最后一 end 才結束循環
    r
  %每過一個輪轉周期(本程序為10次)使各節點的S（i）.G參數（該參數用于后面的簇選舉，在該輪轉周期內已當選過簇頭的節點不能再當選）恢復為零
  if(mod(r, round(1/P) )==0)
    for i=1:1:n
        S3(i).G=0;
        S3(i).cl=0;
    end
  end
Ea=Et*(1-r/rmax)/n;
%(2)死亡節點檢查模塊
dead3=0;
for i=1:1:n
    %檢查有無死亡節點
    if (S3(i).E<=0)
        dead3=dead3+1; 
        %(3)第一個死亡節點的產生時間(用輪次表示)
        %第一個節點死亡時間
        if (dead3==1)
           if(flag_first_dead3==0)
              first_dead3=r;
              flag_first_dead3=1;
           end
        end
        %10%的節點死亡時間
        if(dead3==0.1*n)
           if(flag_teenth_dead3==0)
              teenth_dead3=r;
              flag_teenth_dead3=1;
           end
        end
        if(dead3==n)
           if(flag_all_dead3==0)
              all_dead3=r;
              flag_all_dead3=1;
           end
        end
    end
    if S3(i).E>0
        S3(i).type='N';
    end
end
STATISTICS.DEAD3(r+1)=dead3;
STATISTICS.ALLIVE3(r+1)=allive3-dead3;
%(4)簇頭選舉模塊
countCHs3=0;
cluster3=1;
for i=1:1:n
 if Ea>0
 p(i)=P*n*S3(i).E*E3(i)/(Et*Ea);
 if(S3(i).E>0)
   temp_rand=rand;     
   if ( (S3(i).G)<=0)  
       %簇頭的選舉，當選的簇頭會把各種相關信存入下面程序所給定的變量中
        if(temp_rand<= (p(i)/(1-p(i)*mod(r,round(1/p(i))))))
            countCHs3=countCHs3+1;
            packets_TO_BS3=packets_TO_BS3+1;
            PACKETS_TO_BS3(r+1)=packets_TO_BS3;
             S3(i).type='C';
            S3(i).G=round(1/p(i))-1;
            C3(cluster3).xd=S3(i).xd;
            C3(cluster3).yd=S3(i).yd;
           distance=sqrt( (S3(i).xd-(S3(n+1).xd) )^2 + (S3(i).yd-(S3(n+1).yd) )^2 );
            C3(cluster3).distance=distance;
            C3(cluster3).id=i;
            X3(cluster3)=S3(i).xd;
            Y3(cluster3)=S3(i).yd;
            cluster3=cluster3+1;
           %計算簇頭發送4000bit數據到基站的能量消耗（這里應是所有節點包括簇頭每一輪發送4000bit數據）
           distance;
            if (distance>do)
                S3(i).E=S3(i).E- ( (ETX+EDA)*(4000) + Emp*4000*( distance*distance*distance*distance )); 
            end
            if (distance<=do)
                S3(i).E=S3(i).E- ( (ETX+EDA)*(4000)  + Efs*4000*( distance * distance )); 
            end
        end     
    
    end
    % S3(i).G=S3(i).G-1;  
   
 end 
 end
end
STATISTICS.COUNTCHS3(r+1)=countCHs3;
%(5)簇內成員選擇簇頭模塊(即簇的形成模塊)
%簇內成員對簇頭的選擇（即簇的形成）算法
for c=1:1:cluster3-1
    x3(c)=0;
end
y3=0;
z3=0;
for i=1:1:n
   if ( S3(i).type=='N' && S3(i).E>0 )
     if(cluster3-1>=1)
       min_dis=Inf;
       min_dis_cluster=0;
       for c=1:1:cluster3-1
           temp=min(min_dis,sqrt( (S3(i).xd-C3(c).xd)^2 + (S3(i).yd-C3(c).yd)^2 ) );
           if ( temp<min_dis )
               min_dis=temp;
               min_dis_cluster=c;
               x3(c)=x3(c)+1;
           end
       end
       %簇內節點（發送4000bit數據）能量消耗
 
            min_dis;
            if (min_dis>do)
                S3(i).E=S3(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S3(i).E=S3(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
        %簇頭（接受和融合這一簇內節點4000bit數據）的能量消耗
            S3(C3(min_dis_cluster).id).E = S3(C3(min_dis_cluster).id).E- ( (ERX + EDA)*4000 ); 
            packets_TO_CH3=packets_TO_CH3+1;

            
      
        S3(i).min_dis=min_dis;
        S3(i).min_dis_cluster=min_dis_cluster;
   else
       y3=y3+1;
            min_dis=sqrt( (S3(i).xd-S3(n+1).xd)^2 + (S3(i).yd-S3(n+1).yd)^2 );
            if (min_dis>do)
                S3(i).E=S3(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S3(i).E=S3(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
           packets_TO_BS3=packets_TO_BS3+1;
   end
  end
end
if countCHs3~=0
   u3=(n-y3)/countCHs3;
for c=1:1:cluster3-1
    z3=(x3(c)-u3)*(x3(c)-u3)+z3;
end
LBF3(r+1)=z3/countCHs3;
else  LBF3(r+1)=0;
end
 STATISTICS.PACKETS_TO_CH3(r+1)=packets_TO_CH3;
 STATISTICS.PACKETS_TO_BS3(r+1)=packets_TO_BS3;
end



countCHs4=0;
cluster=1;%此定義的目的僅僅是給定一個1開始的下標參數，真正的簇頭數應該還減去1
flag_first_dead4=0;
flag_teenth_dead4=0;
flag_all_dead4=0;
%死亡節點數
dead4=0;
first_dead4=0;
teenth_dead4=0;
all_dead4=0;
%活動節點數
allive4=n;
%counter for bit transmitted to Bases Station and to Cluster Heads
packets_TO_BS4=0;
packets_TO_CH4=0;
%(1)循環模式設定
for r=0:1:rmax     %該 for 循環將下面的所有程序包括在內，直到最后一 end 才結束循環
    r
  %每過一個輪轉周期(本程序為10次)使各節點的S（i）.G參數（該參數用于后面的簇選舉，在該輪轉周期內已當選過簇頭的節點不能再當選）恢復為零
  if(mod(r, round(1/P) )==0)
    for i=1:1:n
        S4(i).G=0;
        S4(i).cl=0;
    end
  end
Ea=Et*(1-r/rmax)/n;
%(2)死亡節點檢查模塊
dead4=0;
for i=1:1:n
    %檢查有無死亡節點
    if (S4(i).E<=0)
        dead4=dead4+1; 
        %(3)第一個死亡節點的產生時間(用輪次表示)
        %第一個節點死亡時間
        if (dead4==1)
           if(flag_first_dead4==0)
              first_dead4=r;
              flag_first_dead4=1;
           end
        end
        %10%的節點死亡時間
        if(dead4==0.1*n)
           if(flag_teenth_dead4==0)
              teenth_dead4=r;
              flag_teenth_dead4=1;
           end
        end
        if(dead4==n)
           if(flag_all_dead4==0)
              all_dead4=r;
              flag_all_dead4=1;
           end
        end
    end
    if S4(i).E>0
        S4(i).type='N';
    end
end
STATISTICS.DEAD4(r+1)=dead4;
STATISTICS.ALLIVE4(r+1)=allive4-dead4;
%(4)簇頭選舉模塊
countCHs4=0;
cluster4=1;
for i=1:1:n
 if Ea>0
 p(i)=P*n*S4(i).E*E4(i)/(Et*Ea);
 if(S4(i).E>0)
   temp_rand=rand;     
   if ( (S4(i).G)<=0)  
       %簇頭的選舉，當選的簇頭會把各種相關信存入下面程序所給定的變量中
        if(temp_rand<= (p(i)/(1-p(i)*mod(r,round(1/p(i))))))
            countCHs4=countCHs4+1;
            packets_TO_BS4=packets_TO_BS4+1;
            PACKETS_TO_BS4(r+1)=packets_TO_BS4;
             S4(i).type='C';
            S4(i).G=round(1/p(i))-1;
            C4(cluster4).xd=S4(i).xd;
            C4(cluster4).yd=S4(i).yd;
           distance=sqrt( (S4(i).xd-(S4(n+1).xd) )^2 + (S4(i).yd-(S4(n+1).yd) )^2 );
            C4(cluster4).distance=distance;
            C4(cluster4).id=i;
            X4(cluster4)=S4(i).xd;
            Y4(cluster4)=S4(i).yd;
            cluster4=cluster4+1;
           %計算簇頭發送4000bit數據到基站的能量消耗（這里應是所有節點包括簇頭每一輪發送4000bit數據）
           distance;
            if (distance>do)
                S4(i).E=S4(i).E- ( (ETX+EDA)*(4000) + Emp*4000*( distance*distance*distance*distance )); 
            end
            if (distance<=do)
                S4(i).E=S4(i).E- ( (ETX+EDA)*(4000)  + Efs*4000*( distance * distance )); 
            end
        end     
    
    end
    % S4(i).G=S4(i).G-1;  
   
 end 
 end
end
STATISTICS.COUNTCHS4(r+1)=countCHs4;
%(5)簇內成員選擇簇頭模塊(即簇的形成模塊)
%簇內成員對簇頭的選擇（即簇的形成）算法
for c=1:1:cluster4-1
    x4(c)=0;
end
y4=0;
z4=0;
for i=1:1:n
   if ( S4(i).type=='N' && S4(i).E>0 )
     if(cluster4-1>=1)
       min_dis=sqrt( (S4(i).xd-S4(n+1).xd)^2 + (S4(i).yd-S4(n+1).yd)^2 );
       min_dis_cluster=0;
       for c=1:1:cluster4-1
           temp=min(min_dis,sqrt( (S4(i).xd-C4(c).xd)^2 + (S4(i).yd-C4(c).yd)^2 ) );
           if ( temp<min_dis )
               min_dis=temp;
               min_dis_cluster=c;
               x4(c)=x4(c)+1;
           end
       end
       %簇內節點（發送4000bit數據）能量消耗
       if(min_dis_cluster~=0)    
            min_dis;
            if (min_dis>do)
                S4(i).E=S4(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S4(i).E=S4(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
        %簇頭（接受和融合這一簇內節點4000bit數據）的能量消耗
            S4(C4(min_dis_cluster).id).E = S4(C4(min_dis_cluster).id).E- ( (ERX + EDA)*4000 ); 
            packets_TO_CH4=packets_TO_CH4+1;
       else 
            min_dis;
            if (min_dis>do)
                S4(i).E=S4(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S4(i).E=S4(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
            packets_TO_BS4=packets_TO_BS4+1;
            
       end
        S4(i).min_dis=min_dis;
       S4(i).min_dis_cluster=min_dis_cluster;
   else
       y4=y4+1;
            min_dis=sqrt( (S4(i).xd-S4(n+1).xd)^2 + (S4(i).yd-S4(n+1).yd)^2 );
            if (min_dis>do)
                S4(i).E=S4(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S4(i).E=S4(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis)); 
            end
            packets_TO_BS4=packets_TO_BS4+1;
   end
  end
end
if countCHs4~=0
   u4=(n-y4)/countCHs4;
for c=1:1:cluster4-1
    z4=(x4(c)-u4)*(x4(c)-u4)+z4;
end
LBF4(r+1)=z4/countCHs4;
else  LBF4(r+1)=0;
end
STATISTICS.PACKETS_TO_CH4(r+1)=packets_TO_CH4;
STATISTICS.PACKETS_TO_BS4(r+1)=packets_TO_BS4;
end

r=0:5000;
subplot(2,1,1);
plot(r,STATISTICS.DEAD1,':k',r,STATISTICS.DEAD2,'-r',r,STATISTICS.DEAD3,'--b',r,STATISTICS.DEAD4,'-.m');
legend('leachm1','leachm2','DEEC1','DEEC2');
xlabel('x(time)');
ylabel('y(dead)');
title('\bf leach(1,2)和DEEC(1,2)的生命周期對比');
subplot(2,1,2);
plot(r,STATISTICS.PACKETS_TO_BS1,':k',r,STATISTICS.PACKETS_TO_BS2,'-r',r,STATISTICS.PACKETS_TO_BS3,'--b',r,STATISTICS.PACKETS_TO_BS4,'-.m');
legend('leachm1','leachm2','DEEC1','DEEC2');
xlabel('x(time)');
ylabel('y(data)');
title('\bf leach(1,2)和DEEC(1,2)的數據傳輸對比');