%生成網絡部分
%第一步：隨即分布1000個點，并用三維立體圖繪出各點位置。
M=0.0000001;
for k=1:1000
    %if(k==1)
    %    P(1).x=5000;
    %    P(1).y=5000;
    %    P(1).z=5000;
    %else
        P(k).x=fix(10000*rand);
        P(k).y=fix(10000*rand);
        P(k).z=fix(10000*rand);
        hold on
        plot3(P(k).x,P(k).y,P(k).z,'.');
        if(k==1000)
            title('1000個傳感器隨機分布圖');
        end
        %end
P(k).energy=100;%初始的節點能量為100J
P(k).cid=0;%初始化節點網絡號
P(k).id=k;
P(k).next=0;%初始化下一點為0
P(k).level=1;%初始化節點的層次號為最低層號1，最高層層號為3
P(k).out=0;%初始化節點的初始發射功率
P(k).num=0;%初始化節點的鄰節點數

%新增的用于頭節點路由得參數
P(k).ort=0;%為簇頭跳躍的下一點地址(假設簇頭之間的跳躍最多3次)，一跳的平均距離為3000
P(k).sur=0;%為頭節點處在根節點的第幾層的外圍標示，因為最多3跳，所以范圍在1~3之間
%

end
fprintf('根節點的坐標為:');
disp([P(1).x P(1).y P(1).z]);%檢驗是否各點分布正常。

sur(1)=0;%紀錄外圍1的節點數
sur(2)=0;%紀錄外圍2的節點數
shot=50000;%用于紀錄最短的路徑,使用一個不可能存在的數值
shotnum=0;%紀錄最短路徑點的序號
time=0;%初始化時間記錄變量----（由于各網絡間是平行的所以其時間按含有最大子節點的網絡時間計算，按輪計數）
account=0;%初始化傳輸數據變量----(因為每次通信的信息量固定，所以以通信次數統計通信量)
cidnum=0;%網絡數

%第二步：構建網絡。
    %因為范圍在10000*10000*10000內，所以可以計算出每個網絡簇的大概范圍，每個節點的平均范圍為1000*1000*1000兩點間的平均距離為1417。設一個簇有N個節點，又要求兩跳內必須傳至根節點，所以
    %只能有兩層網絡。每個簇的平均大小應該為1+N+N*N=1000，所以N=31。所以認為每簇內節點數在27~35之間是合理的。最底層的網絡范圍
    %為2771
 %廣播生成最低層網絡 
h=1;%網絡頭鏈表的指針 
    for w=1:1000
       if(( P(w).cid==0)&(P(w).level==1))%若n點未處于任何一個網絡且位于最底層
           P(w).cid=w;%修改自己的網絡號
           head(h).id=w;%寫入網絡頭鏈表
           h=h+1;%向后爬行
           P(w).level= P(w).level+1;%修改自己的層號，以為作為簇頭所以層號上升1
           x=w+1;%用于向后續的點廣播該點
           next=w;
           cidnum=cidnum+1;%網絡數增加1
           fprintf('網絡%d的節點有：\n',P(w).id);
             for k=x:1000%向后續點廣播該點
                N1=fix(sqrt((P(w).x-P(k).x)^2+(P(w).y-P(k).y)^2+(P(w).z-P(k).z)^2));%距離公式
                if((N1<=2771)&(P(k).cid==0))%在廣播范圍內且未處于任何一個網絡中
                    P(next).next=k;
                    next=k;%該節點位于第P(n).level層的下一個鄰節點
                    P(next).cid=w;%修改該點的網絡號
                    P(w).num=P(w).num+1;
                    fprintf('%d ',P(next).id);
                    P(k).out=M*N1^2;%暫時假設最后的廣播功率
                end
             end
             fprintf('\n');
        end
    end
    fprintf('共有%d個網絡。\n',cidnum);
    fprintf('頭節點分別是:');
    fprintf('\n');
    for j=1:cidnum
        fprintf('%d  ',head(j).id);
%加入網絡外圍路由
        P(head(1).id).sur=100;%定義一個超范圍的數以區分根節點
        N2=fix(sqrt((P(head(j).id).x-P(1).x)^2+(P(head(j).id).y-P(1).y)^2+(P(head(j).id).z-P(1).z)^2));%距離公式
        if(N2<=6928)
            P(head(j).id).ort=head(1).id;
            P(head(j).id).sur=1;
            P(head(j).id).out=M*N2^2;%修改頭節點到根節點的傳輸功率      
        end
    end
    fprintf('\n');
    fprintf('外圍一層節點頭節點分別是:');
    for j=1:cidnum
        if(P(head(j).id).sur==1)
            sur(1)=sur(1)+1;
            fprintf('%d  ',head(j).id);
        end
    end
    fprintf('\n');
    
    %建立第二層外圍簇頭節點
    fprintf('外圍二層節點頭節點分別是:');
    for j=1:cidnum
        if(P(head(j).id).sur==0)
            sur(2)=sur(2)+1;  
            for v=1:cidnum
                if((P(head(v).id).sur==100)|(P(head(v).id).sur==1))
                    N4=fix(sqrt((P(head(j).id).x-P(head(v).id).x)^2+(P(head(j).id).y-P(head(v).id).y)^2+(P(head(j).id).z-P(head(v).id).z)^2));%距離公式
                    if(N4<shot)
                        shot=N4;
                        shotnum=head(v).id;
                    end
                end
            end
            if(shotnum==1)
                P(head(j).id).sur=1;
                sur(2)=sur(2)-1;
                sur(1)=sur(1)+1;
            else
                P(head(j).id).sur=2;
            end
            P(head(j).id).out=M*shot^2;
            P(head(j).id).ort=shotnum;
            fprintf('%d(%d)  ',head(j).id,shotnum);
        end
    end
    fprintf('\n');
    h=h-1;%最后多加的一點要減去
   
    %繪圖部分
    figure
    hold on
    plot3(P(1).x,P(1).y,P(1).z,'r:p');
    for i=1:1000
        hold on
        if(P(i).sur==1)
            plot3(P(i).x,P(i).y,P(i).z,'b:s');
        end
        if(P(i).sur==2)
            plot3(P(i).x,P(i).y,P(i).z,'k:o');
        end
        if(i==1000)
            title('各頭節點的層次圖');
        end
    end
    figure
    hold on
    plot3(P(1).x,P(1).y,P(1).z,'r:p');
    for i=1:1000
        hold on
        if(P(i).sur==1)
            plot3(P(i).x,P(i).y,P(i).z,'b:s');
        end
        if(P(i).sur==2)
            plot3(P(i).x,P(i).y,P(i).z,'k:o');
        end
        if(P(i).sur==0)
            plot3(P(i).x,P(i).y,P(i).z,'g:.');
        end
        if(i==1000)
            title('最終的網絡層次圖');
        end
    end
    
    while(1)
        n10=input('請輸入您想觀察的網絡號(退出請輸入"0")\n');
        if(n10==0)
            break
        else
            while(P(n10).level~=2)
                fprintf('輸入錯誤!請從新輸入\n');
                n10=input('請輸入您想觀察的網絡號\n');
            end
            z=n10;
            figure
            hold on
            plot3(P(z).x,P(z).y,P(z).z,'r:d');
            fprintf('該簇的節點有: %d ',n10);
            while(P(z).next~=0)
                hold on
                plot3(P(P(z).next).x,P(P(z).next).y,P(P(z).next).z,'b:x');
                fprintf('%d ',P(z).next);
                z=P(z).next;
                if(P(z).next==0)
                    title('該簇的結構圖');
                    fprintf('\n');
                end
            end
        end
    end 
    fprintf('網絡運行中,請稍等……\n');
   
   
%運行部分
%(加入路由部分后，因為認為路由結構不變了，所以本層只需在簇頭能量用完時修改即可)
tr=2;
ti(1)=0;
ac(1)=0;
act(1)=0;
maxnum=0;%用于記錄還有最多子節點的簇的子節點數
num=1000;%num為記錄能量不為零的節點數
receive=0.0005;  %接收所消耗的能量
head=1;%紀錄根節點位置(因為開始時根節點必為1)為0時同樣可以運行
while(num>1)
    maxnum=0;
    %底層通信（網絡的路由結構不改變）
    for k=1:1000 %子節點與簇頭節點的通信，不會影響簇的數目（一次運行為一簇）
        
        if(head==0)%若頭節點無，則跳出本層循環。重新進行網絡重組。
            break;
        end
        if((P(k).level==2)&(P(k).energy>0 )&(P(k).next~=0))%若k為頭節點，開始這一簇得的廣播
            m=k;%m用于移動
            while((P(m).next~=0)&(P(k).energy>0))%簇內的各節點逐個的運行（負責一簇內的運行，一次運行為一簇內的一個子節點）
                P(P(m).next).energy=P(P(m).next).energy-P(P(m).next).out;%發射節點能量損耗
                if(P(P(m).next).energy<=0)%若最底層的子節點的能量耗盡，則刪除該節點并修改鄰接點鏈表
                    num=num-1;%總數減少了1
                    P(k).num=P(k).num-1;%該簇頭子節點減少1
                    P(P(m).next).level=0;%作為能量耗盡的標示
                    fprintf('子節點%d  ',P(m).next);%顯示出能量耗盡的點的序號
                    P(m).next=P(P(m).next).next;%從鏈表中刪除該點P(m).next
                else   %簇內的該子點的能量沒有用完
                    P(k).energy=P(k).energy-receive;%底層能量未耗盡，則簇頭減去接收所需的能量,receive為接收信號的固定損耗
                    if(P(k).energy<=0) %若頭節點的能量耗盡；同時認為本簇的本次循環終止，系統會按順序轉入新的簇運行<修改頭節點時不會影響到網絡結構>
                        num=num-1;%總數減少了1
                        fprintf('頭節點%d  ',P(k).id);%輸出該頭節點
                        P(k).level=0;%作為能量耗盡的標示
                        P(P(k).next).level=2;%改變頭節點。不用考慮P(k).next是否存在，因為如果不存在則不會出現頭節點能量為零的情況
                        z=k;%z用于對頭節點改造的計序
                        while(P(z).next~=0)%（內二層循環）用于對頭節點的改造
                            P(P(k).next).num=P(P(k).next).num+1;%鄰節點數加1
                            P(P(z).next).cid=P(k).next;%改變網絡中各節點的頭節點標識
                            P(P(z).next).out=M*(fix(sqrt((P(P(z).next).x-P(P(k).next).x)^2+(P(P(z).next).y-P(P(k).next).y)^2+(P(P(z).next).z-P(P(k).next).z)^2)))^2;%距離公式%改變至總簇頭的傳輸能量%修改傳輸功率
                            z=P(z).next;%向后爬動
                        end%完成對各子節點的改造
                        if(P(k).next~=0)%如果上面的循環運行了，減去紀錄自己的一次
                            P(P(k).next).num=P(P(k).next).num-1;%減去自己的一點
                        end
                        if(k==head)%如果是頭節點則結構重組，跳過本次運行
                            head=0;
                            break;
                        end
                        %修改路由
                        P(P(k).next).sur=P(k).sur;
                        P(P(k).next).ort=P(k).ort;
                        %該點代替頭節點鏈表中的位置,因為有后續節點代替所以簇數不會減少
                        %用后續節點取代原節點作中繼節點
                        if(P(k).sur==1)%第一層的節點才有可能為中繼節點
                            for i=1:1000
                                if(P(i).ort==k)
                                    P(i).ort=P(k).next;
                                    P(i).out=M*(fix(sqrt((P(P(k).next).x-P(i).x)^2+(P(P(k).next).y-P(i).y)^2+(P(P(k).next).z-P(i).z)^2)))^2;%距離公式，改變傳輸功率為至總簇頭的傳輸能量
                                end
                            end
                        end
                        Q=P(P(k).next).ort;
                        P(P(k).next).out=M*(fix(sqrt((P(P(k).next).x-P(Q).x)^2+(P(P(k).next).y-P(Q).y)^2+(P(P(k).next).z-P(Q).z)^2)))^2;%距離公式，改變傳輸功率為至總簇頭的傳輸能量
                        break;%因為簇頭節點的能量用完了，所以需要跳出本簇該層的循環，即 while((P(m).next~=0)&(P(k).energy>0))；提前結束本簇運行
                    else%若頭節點的能量也沒有用完,完成一次通信，并向后爬行
                        account=account+1;%完成一次通信  
                    end%完成簇內的一次通信(頭節點的能量沒有用完)
                end%子節點的能量沒有用完
                if(P(m).next~=0)
                    m=P(m).next;%向后爬行?????(存在問題:當改為在下個end之后就不對了)