clear
%1.初始參數設定模塊
%.傳感器節點區域界限(單位 M)
xm=100;
ym=100;
%(1)匯聚節坐標給定
sink.x=0.5*xm;
sink.y=0.5*ym;
%區域內傳器節數
n=100
%簇頭優化比例（當選簇頭的概率）
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;
%最大循環次數
rmax=5000
%算出參數 do
do=sqrt(Efs/Emp);
%2.無線傳感器網絡模型產生模塊
%構建無線傳感器網絡,在區域內均勻投放100個節點,并畫出圖形
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;
    S(i).E=Eo;
    %initially there are no cluster heads only nodes
    S(i).type='N';
end

S(n+1).xd=sink.x;
S(n+1).yd=sink.y;
%3.網絡運行模塊
%簇頭節點數
countCHs=0;
cluster=1;%此定義的目的僅僅是給定一個1開始的下標參數，真正的簇頭數應該還減去1
flag_first_dead=0;
flag_teenth_dead=0;
flag_all_dead=0;
%死亡節點數
dead=0;
first_dead=0;
teenth_dead=0;
all_dead=0;
%活動節點數
allive=n;
%counter for bit transmitted to Bases Station and to Cluster Heads
packets_TO_BS=0;
packets_TO_CH=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
        S(i).G=0;
        S(i).cl=0;
    end
  end
%(2)死亡節點檢查模塊
dead=0;
for i=1:1:n
    %檢查有無死亡節點
    if (S(i).E<=0)
        dead=dead+1; 
        %(3)第一個死亡節點的產生時間(用輪次表示)
        %第一個節點死亡時間
        if (dead==1)
           if(flag_first_dead==0)
              first_dead=r;
              flag_first_dead=1;
           end
        end
        %10%的節點死亡時間
        if(dead==0.1*n)
           if(flag_teenth_dead==0)
              teenth_dead=r;
              flag_teenth_dead=1;
           end
        end
        if(dead==n)
           if(flag_all_dead==0)
              all_dead=r;
              flag_all_dead=1;
           end
        end
    end
    if S(i).E>0
        S(i).type='N';
    end
end
STATISTICS.DEAD(r+1)=dead;
STATISTICS.ALLIVE(r+1)=allive-dead;
%(4)簇頭選舉模塊
countCHs=0;
cluster=1;
for i=1:1:n
 if(S(i).E>0)
   temp_rand=rand;     
   if ( (S(i).G)<=0) 
       %簇頭的選舉，當選的簇頭會把各種相關信存入下面程序所給定的變量中
        if(temp_rand<= (p/(1-p*mod(r,round(1/p)))))
            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=round(1/p)-1;
            C(cluster).xd=S(i).xd;
            C(cluster).yd=S(i).yd;
           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;
           %計算簇頭發送4000bit數據到基站的能量消耗（這里應是所有節點包括簇頭每一輪發送4000bit數據）
           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
    % S(i).G=S(i).G-1;  
   
 end 
end
STATISTICS.COUNTCHS(r+1)=countCHs;
%(5)簇內成員選擇簇頭模塊(即簇的形成模塊)
%簇內成員對簇頭的選擇（即簇的形成）算法
for c=1:1:cluster-1
    x(c)=0;
end
y=0;
z=0;
for i=1:1:n
   if ( S(i).type=='N' && S(i).E>0 )
     if(cluster-1>=1)
       min_dis=Inf;
       min_dis_cluster=0;
       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;
               x(c)=x(c)+1;
           end
       end
       %簇內節點（發送4000bit數據）能量消耗
            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
        %簇頭（接受和融合這一簇內節點4000bit數據）的能量消耗
            S(C(min_dis_cluster).id).E = S(C(min_dis_cluster).id).E- ( (ERX + EDA)*4000 ); 
            packets_TO_CH=packets_TO_CH+1;

       S(i).min_dis=min_dis;
       S(i).min_dis_cluster=min_dis_cluster;
   else
         y=y+1;
          min_dis=sqrt( (S(i).xd-S(n+1).xd)^2 + (S(i).yd-S(n+1).yd)^2 );
            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
            packets_TO_BS=packets_TO_BS+1;
    end
  end
end
if countCHs~=0
   u=(n-y)/countCHs;
for c=1:1:cluster-1
    z=(x(c)-u)*(x(c)-u)+z;
end
LBF(r+1)=countCHs/z;
else  LBF(r+1)=0;
end
STATISTICS.PACKETS_TO_CH(r+1)=packets_TO_CH;
STATISTICS.PACKETS_TO_BS(r+1)=packets_TO_BS;
end
first_dead
teenth_dead
all_dead
STATISTICS.DEAD(r+1)
STATISTICS.ALLIVE(r+1)
STATISTICS.PACKETS_TO_CH(r+1)
STATISTICS.PACKETS_TO_BS(r+1)
STATISTICS.COUNTCHS(r+1)
r=0:5000;
subplot(2,2,1);
plot(r,STATISTICS.DEAD);
subplot(2,2,2);
plot(r,STATISTICS.ALLIVE);
subplot(2,2,3);
plot(r,STATISTICS.PACKETS_TO_BS);
subplot(2,2,4);
plot(r,STATISTICS.COUNTCHS);
%STATISTICS，結構體數組，包括下面的5個變量；
%countCHs(r+1）,每一輪所選出的簇頭數目;
%packets_TO_BS(r+1),基站收到的數據包總數;
%PACKETS_TO_CH(r+1),簇頭收到的數據包總數;
%first_dead,第一個節點死亡的時間;
%teenth_dead=r,10%的節點死亡的時間；
%dead(r+1),每一輪的死亡節點數；
%allive(r+1),每一輪的活動節點數。