%ddddddddddddddddddddd      基于選票分簇算法的能量消耗--非均勻分簇       ddddddddddddddddddddddddddd
%開始每個節點廣播ID號和當前能量；接收計算，每個節點廣播ID號和當前選票和適宜性（或節點度）
%接收，若本節點選票最多，作為簇頭，并廣播本節點ID號結束；簇成員節點接受簇頭信息，加入并廣播本節點ID號結束。

%能量消耗包括兩部分：簇形成階段的能耗和數據傳輸階段的能耗。
%dddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd

N=50;                    %節點個數 
f=1;                     %最后結果數組下標，50-500, 每隔50取一數，共有10

while N<=500

g=5;                     %相同條件下的重復次數    
for t=1:g
%dddddddddddddddddddddddddddddd     節點配置    dddddddddddddddddddddddddddddddddddddd
%   radius=25;                              %廣播距離
X_axis=100;                    %x軸長度 代表節點分布范圍中的一維 
Y_axis=100;                    %y軸長度 代表節點分布范圍中的另一維 
Loca=[X_axis*rand(1,N);Y_axis*rand(1,N)]';    %產生坐標矩陣 

for i=1:N
    e(i)=10;                  %節點的能量 
    Eprog(i)=0;
end  
er=0.533*10^(-6);
et=1.066*10^(-6);
eta=3.076*10^(-9);
efuse=1.8*10^(-7);
es=1.75*10^(-7);
D=-175;
c=0.2;
dopt=sqrt((er+et)/eta);
radius=dopt/2;                              %廣播距離
dT=210;                             %門限距離，小于直接傳輸到網關，大于多跳到網關

%ddddddddddddddddddddddddddddd         節點位置示意圖    dddddddddddddddddddddddddddddddd
 %plot(Loca(:,1),Loca(:,2),'o')

%ddddddddddddddddddddddddddd        計算鄰居節點矩陣   ddddddddddddddddddddddddddddddddddddddd
                                 
for i=1:N
    d(i)=0;                             %節點的度 
end    
for i=1:N
    for j=i:N
        dis(i,j)=sqrt((Loca(i,1)-Loca(j,1))*(Loca(i,1)-Loca(j,1))+(Loca(i,2)-Loca(j,2))*(Loca(i,2)-Loca(j,2)));          %存放節點間距離 
        dis(j,i)=dis(i,j);
    end
end
dmax=0;
dmin=250;
for i=1:N
    dis(i,N+1)=sqrt((Loca(i,1))*(Loca(i,1))+(Loca(i,2)-D)*(Loca(i,2)-D));          %存放節點到網關的距離 
    if dis(i,N+1)>dmax
       dmax=dis(i,N+1);
    end
    if dis(i,N+1)<dmin
       dmin=dis(i,N+1);
    end
end
  
for i=1:N   
    r(i)=radius*(1-c*(dmax-dis(i,N+1))/(dmax-dmin));
    for j=1:N
        if dis(i,j)<=r(i)&&i~=j
           d(i)=d(i)+1;
           neibor(i,d(i))=j;
       end 
    end
   if d(i)==0
     neibor(i,1)=0;
   end  
end
for i=1:N
    n1=ceil(dis(i,N+1)/dopt);
    syms k;
    n2=double(symsum(1/k,1,n1));
    dn(i)=dis(i,N+1)/n2;
end

%ddddddddddddddddddddddddddd        節點布置后，直到所有節點能量耗盡，所經歷的輪數   ddddddddddddddddddddddddddddddddddddddd
lunum=0;
flag2=1;
Ecluz=0;
Edataz=0;
flag3=0;
Nper=10;
dienote=0;

while flag2==1

%ddddddddddddddddddddddddddd        計算節點初始階段的能量消耗   ddddddddddddddddddddddddddddddddddddddd
for i=1:N
    Edata(i)=0;
    Einit(i)=0;
    Eprog(i)=0;
end    

for i=1:N
    if e(i)~=0
       Einit(i)=(6*8+8*8+8*8)*(eta*r(i)^2+et)+d(i)*(6*8+8*8+8*8)*er;
    end   
end    


      
%ddddddddddddddddddddddddddd        計算節點初始選票   ddddddddddddddddddddddddddddddddddddddd
for i=1:N
    if e(i)~=0
       E(i)=e(i);
       if d(i)~=0
          for k=1:d(i)
              E(i)=E(i)+e(neibor(i,k));
          end
       end  
       v(i,i)=e(i)/E(i);
       if d(i)~=0
          for k=1:d(i)
              v(i,neibor(i,k))=e(neibor(i,k))/E(i);
          end
       end   
    end
end    
for i=1:N
    if e(i)~=0
       if d(i)==0
          fitness(i)=e(i);
       else   
          fitness(i)=e(i)/d(i);           %節點的適應性
       end    
       curvote(i)=v(i,i);
       if d(i)~=0
          for k=1:d(i)
              curvote(i)=curvote(i)+v(neibor(i,k),i);          %節點的總選票
          end
       end
    end   
 end

%ddddddddddddddddddddddddddd  查找本節點是否是鄰居集合中最多選票的節點，如是發簇頭廣播 ddddddddddddddddddddddddddddddddddddddd
for i=1:N
    CH(i)=0;
    cluster_head(i)=0;
end    
for i=1:N
    if e(i)~=0
       maxvote=curvote(i);
       maxvoteID(i)=i;
       flag=0;
       if d(i)~=0
          for k=1:d(i)
             if curvote(neibor(i,k))>maxvote
                flag=1;
                break
             end
          end
       end   
       if flag==0
          cluster_head(i)=i;
          CH(i)=1;
          Eprog(i)=Eprog(i)+4*8*(eta*r(i)^2+et);         %發簇頭信息
       end
    end   
end

    
%ddddddddddddddddddddddddddd       處理收到的簇頭信息的節點,按適宜性選擇簇頭  ddddddddddddddddddddddddddddddddddddddd
for i=1:N
    schnum(i)=0;
end    
for i=1:N
    if CH(i)==1
       if d(i)~=0 
          for k=1:d(i)
              schnum(neibor(i,k))=schnum(neibor(i,k))+1;
              sch(neibor(i,k),schnum(neibor(i,k)))=i;
              CH(neibor(i,k))=2;
              Eprog(neibor(i,k))=Eprog(neibor(i,k))+4*8*er+4*8*(eta*r(neibor(i,k))^2+et);       %接收簇頭信息并發撤消包
              if d(neibor(i,k))~=0
                 for b=1:d(neibor(i,k))
                     Eprog(neibor(neibor(i,k),b))=Eprog(neibor(neibor(i,k),b))+4*8*er;       %接收撤消包
                 end    
              end
          end
       end  
    end  
end   
     
 for i=1:N
     if CH(i)==2
        fitness1=0;
        if schnum(i)~=0
           for k=1:schnum(i)           %收到的簇頭集合中適宜性最高的簇頭
               if fitness(sch(i,k))>fitness1
                  fitness1=fitness(sch(i,k));
                  head=sch(i,k); 
              end
            end
        end 
        flag1=0; 
        if d(i)~=0
           for k=1:d(i)                     %簇頭的適宜性是否在未被覆蓋的節點中最高
               if CH(neibor(i,k))==0
                  if fitness(neibor(i,k))>fitness1
                     flag1=1;
                     break
                 end
              end
           end 
        end
        if flag1==0
           cluster_head(i)=head;
           CH(i)=1;
           Eprog(i)=Eprog(i)+6*8*(eta*r(i)^2+et);         %發加入信息
        end
    end   
  end
  
%ddddddddddddddddddddddddddd       處理未被覆蓋的節點  ddddddddddddddddddddddddddddddddddddddd 
flag=1;
while flag==1
      flag=0;
      chval=3;    
      for i=1:N            %重新計算未覆蓋節點的選票    
          if e(i)~=0 
             if CH(i)==0
                flag=1;
                E(i)=e(i);
                if d(i)~=0
                   for k=1:d(i)
                       if CH(neibor(i,k))==0
                          E(i)=E(i)+e(neibor(i,k));
                       end  
                   end   
                end  
                v(i,i)=e(i)/E(i);
                if d(i)~=0
                   for k=1:d(i)
                       if CH(neibor(i,k))==0
                          v(i,neibor(i,k))=e(neibor(i,k))/E(i);
                       else   
                          v(i,neibor(i,k))=0;
                       end 
                   end   
                end
            end
          end
      end   
         
        for i=1:N
            if e(i)~=0
               if CH(i)==0
                  Eprog(i)=Eprog(i)+6*8*(eta*r(i)^2+et);         %發新選票信息 
                  curvote(i)=v(i,i);
                  if d(i)~=0
                     for k=1:d(i)
                         if CH(neibor(i,k))==0
                            curvote(i)=curvote(i)+v(neibor(i,k),i);
                         end  
                     end  
                  end
               end
            end   
        end
        
        for i=1:N               %查找本節點是否是鄰居集合中最多選票的節點，如是發簇頭廣播
            if e(i)~=0
               if CH(i)==0
                  maxvote=curvote(i);
                  maxvoteID(i)=i;
                  flag1=0;
                  if d(i)~=0
                     for k=1:d(i)
                         if CH(neibor(i,k))==0
                            if curvote(neibor(i,k))>maxvote
                               flag1=1;
                               break
                            end
                         end  
                     end
                  end    
                  if flag1==0
                     cluster_head(i)=i;
                     CH(i)=chval;
                     Eprog(i)=Eprog(i)+4*8*(eta*r(i)^2+et);         %發簇頭信息
                  end     
               end
            end  
        end
        
       for i=1:N            %  處理收到的簇頭信息的節點