clc;
clear;
close all;
load('new_point');
round=1;
death_nodes=[];
communication=[];
energy=[];
while round<=100
    for i=1:1:node_number
       S(i).father=-1;
       S(i).child_num=0;
       S(i).intree=0;
       S(i).data=data_length;
       S(i).ifleaf=0;
       S(i).child_leaf=0;
       S(i).child_middle=0;
       S(i).child=[];
    end
    sink.child=[];
    sink.child_num=0;
    sink.CHO=[];
     sink.num_neighbor=0;
    line=[];
    line_head=1;
    line_tail=0;
    l=0;
    sink.child_num=0;
    for i=1:1:node_number
        l=sqrt((sink.x-S(i).xd)^2+(sink.y-S(i).yd));
        if l<R
            communication_overhead=communication_overhead+beacon_length;
            energy_consumption=energy_consumption+beacon_length*(en+FS*l^2);
            energy_consumption=energy_consumption+beacon_length*en;
            S(i).en=S(i).en-beacon_length*en;
            communication_overhead=communication_overhead+beacon_length;
            energy_consumption=energy_consumption+beacon_length*(en+FS*l^2);
            energy_consumption=energy_consumption+beacon_length*en;
            S(i).en=S(i).en-beacon_length*(en+FS*l^2);
            sink.child_num=sink.child_num+1;
            sink.child(sink.child_num)=i;
            S(i).father=0;
            S(i).intree=1;
        end
    end
    
    for i=1:1:sink.child_num
        sink.CHO(i)=S(sink.child(i)).en/sqrt((sink.x-S(sink.child(i)).xd)^2+(sink.y-S(sink.child(i)).yd)^2);
    end
        bs_child=sink.CHO;
    for i=1:1:sink.child_num
        [max_value, max_index] = max(bs_child);
        line_tail=line_tail+1;
        line(line_tail)=sink.child(max_index);
        bs_child(max_index)=-9999;
    end
    
    while line_head<=line_tail
            bs_child=[];
            cho=[];
            count=0;
            for i=1:1:S(line(line_head)).num_neighbor
                l1=sqrt((S(line(line_head)).xd-S(S(line(line_head)).neighbor(i)).xd)^2+ ...
                    (S(line(line_head)).yd-S(S(line(line_head)).neighbor(i)).yd)^2);
                communication_overhead=communication_overhead+beacon_length;
                energy_consumption=energy_consumption+beacon_length*(en+FS*l1^2);
                energy_consumption=energy_consumption+beacon_length*en;
                S(line(line_head)).en=S(line(line_head)).en-beacon_length*(en+FS*l1^2);
                S(S(line(line_head)).neighbor(i)).en=S(S(line(line_head)).neighbor(i)).en-beacon_length*en;
                if S(S(line(line_head)).neighbor(i)).intree==0
                    count=count+1;
                    bs_child(count)=S(line(line_head)).neighbor(i);
                    communication_overhead=communication_overhead+beacon_length;
                    energy_consumption=energy_consumption+beacon_length*(en+FS*l1^2);
                    energy_consumption=energy_consumption+beacon_length*en;
                    S(S(line(line_head)).neighbor(i)).en=S(S(line(line_head)).neighbor(i)).en-beacon_length*(en+FS*l1^2);
                    S(line(line_head)).en=S(line(line_head)).en-beacon_length*en;
                end
            end
            if count>0
                if count==1
                    S(line(line_head)).child_num=1;
                else
                    S(line(line_head)).child_num=randi(count);
                end
                for i=1:1:count
                    cho(i)=S(bs_child(i)).en/sqrt((S(line(line_head)).xd-S(bs_child(i)).xd)^2+(S(line(line_head)).yd-S(bs_child(i)).yd)^2);
                end
                for i=1:1:S(line(line_head)).child_num
                    [max_value, max_index] = max(cho);
                    S(line(line_head)).child(i)=bs_child(max_index);
                    S(S(line(line_head)).child(i)).intree=1;
                    S(S(line(line_head)).child(i)).father=line(line_head);
                    line_tail=line_tail+1;
                    line(line_tail)=S(line(line_head)).child(i);
                    cho(max_index)=-9999;
                    l2=sqrt((S(line(line_head)).xd-S(S(line(line_head)).child(i)).xd)^2 ...
                        +(S(line(line_head)).yd-S(S(line(line_head)).child(i)).yd)^2);
                    communication_overhead=communication_overhead+beacon_length;
                    energy_consumption=energy_consumption+beacon_length*(en+FS*l2^2);
                    S(line(line_head)).en=S(line(line_head)).en-beacon_length*en;
                    S(S(line(line_head)).child(i)).en=S(S(line(line_head)).child(i)).en-beacon_length*(en+FS*l2^2);
                    energy_consumption=energy_consumption+beacon_length*en;
                end
            else
                S(line(line_head)).ifleaf=1;
            end
        line_head=line_head+1;
    end
    count_intree=0;
 
        for i=1:1:node_number
               cho=[];
            if S(i).intree==1
                count_intree=count_intree+1;
            else
                for j=1:1:S(i).num_neighbor
                    cho(j)=sqrt((S(i).xd-S(S(i).neighbor(j)).xd)^2+(S(i).yd-S(S(i).neighbor(j)).yd)^2);
                    if S(S(i).neighbor(j)).intree==0
                        cho(j)=-9999;
                    end
                end
                [max_value, max_index] = max(cho);
                S(i).father=S(i).neighbor(max_index);
                S(i).intree=1;
                communication_overhead=communication_overhead+beacon_length;
                l3=sqrt((S(i).xd-S(S(i).father).xd)^2+(S(i).yd-S(S(i).father).yd));
                energy_consumption=energy_consumption+beacon_length*(en+FS*l3^2);
                energy_consumption=energy_consumption+beacon_length*en;
                S(S(i).neighbor(max_index)).en=S(S(i).neighbor(max_index)).en-beacon_length;
                S(i).en=S(i).en-beacon_length*(en+FS*l3^2);
                S(S(i).neighbor(max_index)).child_num=S(S(i).neighbor(max_index)).child_num+1;
                S(S(i).neighbor(max_index)).child(S(S(i).neighbor(max_index)).child_num)=i;
                count_intree=count_intree+1;
            end
        end
    
        for i=1:1:node_number
            if S(i).child_num>0
                S(i).ifleaf=0;
            end
        end
%%%%%%%%%%%%%%%%%%%cut%%%%%%%%%%%%%%%%%%    
    for i=1:1:node_number
        rem_j=[];
        if S(i).ifleaf==1
           if(S(i).num_neighbor>J)
               j=1;
               rem_j(j)=randi(S(i).num_neighbor);
                while j<J-1
                    rem=randi(S(i).num_neighbor);
                    flag=0;
                    for k=1:1:j
                        if rem_j(k)==rem||S(i).neighbor(rem)==S(i).father
                             flag=1; 
                        end
                    end
                    if flag==0
                        j=j+1;
                        rem_j(j)=rem;
                    end      
                end
                S(i).data=S(i).data/J;
                communication_overhead=communication_overhead+(J-1)*S(i).data;
                for l=1:1:J-1
                    S(S(i).neighbor(rem_j(l))).data=S(S(i).neighbor(rem_j(l))).data+S(i).data;
                    l4=sqrt((S(i).xd-S(S(i).neighbor(rem_j(l))).xd)^2 ...
                        +(S(i).yd-S(S(i).neighbor(rem_j(l))).yd)^2);
                    energy_consumption=energy_consumption+S(i).data*(en+FS*l4^2);
                    energy_consumption=energy_consumption+S(i).data*en;
                    S(i).en=S(i).en-S(i).data*(en+FS*l4^2);
                    S(S(i).neighbor(rem_j(l))).en=S(S(i).neighbor(rem_j(l))).en-S(i).data*en;
                end
           end
        end
    end
    
    sentinel=1;
    while sentinel==1
        sentinel=0;
        for i=1:1:node_number
            if S(i).data~=0
               sentinel=1;
               communication_overhead=communication_overhead+S(i).data;
               if S(i).father==0
                   l4=sqrt((S(i).xd-sink.x)^2+(S(i).yd-sink.y)^2);
                   energy_consumption=energy_consumption+S(i).data*(en+FS*l4^2);
                   energy_consumption=energy_consumption+S(i).data*en;
                   S(i).en=S(i).en-S(i).data*(en+FS*l4^2);
                    S(i).data=0;
               else 
                   S(S(i).father).data=S(S(i).father).data+S(i).data;
                   l4=sqrt((S(i).xd-S(S(i).father).xd)^2+(S(i).yd-S(S(i).father).yd)^2);
                   energy_consumption=energy_consumption+S(i).data*(en+FS*l4^2);
                    energy_consumption=energy_consumption+S(i).data*en;
                    S(i).en=S(i).en-S(i).data*(en+FS*l4^2);
                    S(S(i).father).en=S(S(i).father).en-S(i).data*en;
                   S(i).data=0;
               end
            end
        end
    end
    death_node=0;
    for i=1:1:node_number
        if S(i).en<0
            death_node=death_node+1;
        end
    end
    death_nodes(round)=death_node;
    communication(round)=communication_overhead;
    energy(round)=energy_consumption;
    round=round+1;
end