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=1;
       S(i).child_leaf=0;
       S(i).child_middle=0;
       S(i).child=[];
    end
    sink.child=[];
    sink.CHO=[];
    sink.num_neighbor=0;

    for i=1:1:node_number
        dis=sqrt((sink.x-S(i).xd)^2+(sink.y-S(i).yd)^2);
        if dis<R
            communication_overhead=communication_overhead+beacon_length;
            energy_consumption=energy_consumption+beacon_length*(en+FS*dis^2);
            energy_consumption=energy_consumption+beacon_length*en;
            energy_consumption=energy_consumption+beacon_length*(en+FS*dis^2);
            energy_consumption=energy_consumption+beacon_length*en;
            S(i).en=S(i).en-beacon_length*en;
            S(i).father=0;
            sink.num_neighbor=sink.num_neighbor+1;
            sink.child(sink.num_neighbor)=i;
            S(i).intree=1;
        end
    end
    line=[];
    line_head=0;
    line_tail=0;
    Dis=[];
    for i=1:1:sink.num_neighbor
        Dis(i)=sqrt((sink.x-S(sink.child(i)).xd)^2+(sink.y-S(sink.child(i)).yd)^2);
    end

    for i=1:1:sink.num_neighbor
    [min_value, min_index] = min(Dis);
    line_tail=line_tail+1;
    line(line_tail)=sink.child(min_index);
    Dis(min_index)=9999;
    end

    line_head=line_head+1;
    while(line_head<=line_tail)
        Dis=[];
        for i=1:1:node_number
            if S(i).intree==0
                dis=sqrt((S(line(line_head)).xd-S(i).xd)^2+(S(line(line_head)).yd-S(i).yd)^2);
                if dis<R
                    communication_overhead=communication_overhead+beacon_length;
                    energy_consumption=energy_consumption+beacon_length*(en+FS*dis^2);
                    energy_consumption=energy_consumption+beacon_length*en;
                    energy_consumption=energy_consumption+beacon_length*(en+FS*dis^2);
                    energy_consumption=energy_consumption+beacon_length*en;
                    S(line(line_head)).en=S(line(line_head)).en-beacon_length*(en+FS*dis^2);
                    S(i).en=S(i).en-beacon_length*en;
                    S(line(line_head)).child_num=S(line(line_head)).child_num+1;
                    S(line(line_head)).child(S(line(line_head)).child_num)=i;
                    S(i).father=line(line_head);
                    S(i).intree=1;
                end
            end
        end
        if S(line(line_head)).child_num==0
            S(line(line_head)).ifleaf=1;
        else
            for j=1:1:S(line(line_head)).child_num
                Dis(j)=sqrt((S(line(line_head)).xd-S(S(line(line_head)).child(j)).xd)^2 ...
                    +(S(line(line_head)).yd-S(S(line(line_head)).child(j)).yd)^2);
            end
            for j=1:1:S(line(line_head)).child_num
                [min_value, min_index] = min(Dis);
                line_tail=line_tail+1;
                line(line_tail)=S(line(line_head)).child(min_index);
                Dis(min_index)=9999;
            end
        end
        line_head=line_head+1;
    end
    for i=1:1:node_number
        if S(i).child_num>0
            for j=1:1:S(i).child_num
                if S(S(i).child(j)).ifleaf==1
                    S(i).child_leaf= S(i).child_leaf+1;
                else
                    S(i).child_middle=S(i).child_middle+1;
                end
            end
        end
    end
    
    for i=1:1:node_number
        if S(i).ifleaf==1
            S(i).PIOC=0;
        end
    end
    count=0;
    sink.num_child=0;
    for i=1:1:sink.num_neighbor
        if S(sink.child(i)).ifleaf==1
               count=count+1;
        else
            sink.num_child=sink.num_child+1;
        end
    end
    over_flag_1=1;
    while 1
        over_flag_1=1;
        for i=1:1:sink.num_neighbor
            if S(sink.child(i)).PIOC==0
                list_head=sink.child(i);
                S(list_head).father=0;
                sink.num_child=sink.num_child+1;
                over_flag_1=0;
                break;
            end
        end
        if over_flag_1==1
            break;
        end
        over_flag_2=1;
        while 1
            S(list_head).PIOC=1;
            father=list_head;
            over_flag_2=1;
            for i=1:1:sink.num_neighbor
                if S(sink.child(i)).ifleaf==1&&S(sink.child(i)).PIOC==0
                    dis=sqrt((S(sink.child(i)).xd-S(list_head).xd)^2+ ...
                        (S(sink.child(i)).yd-S(list_head).yd)^2);
                    if dis<R
                        communication_overhead=communication_overhead+beacon_length;
                        energy_consumption=energy_consumption+beacon_length*(en+FS*dis^2);
                        energy_consumption=energy_consumption+beacon_length*en;
                        energy_consumption=energy_consumption+beacon_length*(en+FS*dis^2);
                        energy_consumption=energy_consumption+beacon_length*en;
                        S(list_head).en=S(list_head).en-beacon_length*(en+FS*dis^2);
                        S(sink.child(i)).en=S(sink.child(i)).en-beacon_length*en;
                        S(list_head).child_num=1;
                        S(list_head).child(1)=sink.child(i);
                        S(list_head).ifleaf=0;
                        list_head=sink.child(i);
                        S(list_head).father=father;
                        over_flag_2=0;
                        break;
                    end
                end
            end
            if over_flag_2==1
                break;
            end
        end
    end

    line=[];
    line_head=1;
    line_tail=0;
    for i=1:1:sink.num_neighbor
        if S(sink.child(i)).ifleaf==0
                line_tail=line_tail+1;
                line(line_tail)=sink.child(i);
        end
    end
    
    while line_head<=line_tail
            for i=1:1:S(line(line_head)).child_num
                if S(S(line(line_head)).child(i)).ifleaf==0
                    line_tail=line_tail+1;
                    line(line_tail)=S(line(line_head)).child(i);
                end
            end
            over_flag_3=1;
        while 1
            over_flag_3=1;
            for i=1:1:S(line(line_head)).child_num
                if S(S(line(line_head)).child(i)).PIOC==0
                    list_head=S(line(line_head)).child(i);
                    S(list_head).father=line(line_head);
                    over_flag_3=0;
                end
            end
            if over_flag_3==1
                break;
            end
           over_flag_4=0;
            while 1
                S(list_head).PIOC=1;
                father=list_head;
                over_flag_4=1;
                for i=1:1:S(line(line_head)).child_num
                    if S(S(line(line_head)).child(i)).ifleaf==1&&S(S(line(line_head)).child(i)).PIOC==0
                        dis=sqrt((S(S(line(line_head)).child(i)).xd-S(list_head).xd)^2 ...
                            +(S(S(line(line_head)).child(i)).yd-S(list_head).yd)^2);
                        if dis<R
                            communication_overhead=communication_overhead+beacon_length;
                            energy_consumption=energy_consumption+beacon_length*(en+FS*dis^2);
                            energy_consumption=energy_consumption+beacon_length*en;
                            S(S(line(line_head)).child(i)).en=S(S(line(line_head)).child(i)).en-beacon_length*en;
                            S(list_head).en=S(list_head).en-beacon_length*(en+FS*dis^2);
                            energy_consumption=energy_consumption+beacon_length*(en+FS*dis^2);
                            energy_consumption=energy_consumption+beacon_length*en;
                            S(list_head).en=S(list_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*dis^2);
                            S(list_head).child_num=1;
                            S(list_head).child(1)=S(line(line_head)).child(i);
                            S(S(line(line_head)).child(i)).father=list_head;
                            S(list_head).ifleaf=0;
                            list_head=S(line(line_head)).child(i);
                            over_flag_4=0;
                            break;
                        end
                    end
                end
                if over_flag_4==1
                    break;
                end
            end
        end
        line_head=line_head+1;
    end

    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(S(i).neighbor(rem_j(l))).en=S(S(i).neighbor(rem_j(l))).en-S(i).data*en;
                   S(i).en=S(i).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