移動通信基站布局規(guī)劃分析

隨著社會的快速發(fā)展，我們的生活越來越離不開手機，手機已成為我們生活的必需品。大家有沒有遇到這樣的問題？就是在高樓大廈林立的地方，我們會發(fā)現(xiàn)手機信號特別差，通話斷斷續(xù)續(xù)并伴有雜音，視頻卡成狗，地下車庫更是盲區(qū)，為什么會出現(xiàn)這樣的問題呢？

其實手機信號的傳輸是一種無線電的發(fā)射與接收過程，傳播路徑主要是以直線傳播為主，若在傳播過程中遇到阻礙（即與高樓大廈墻壁發(fā)生多次反射、繞射等物理現(xiàn)象。其中，繞射是指在電磁波傳播路徑上，當電波被尺寸較大(與波長相比)的障礙物遮擋時，電磁波改變傳播方向的現(xiàn)象），就會大大削弱其信號強度。那么，該如何解決呢？

在現(xiàn)實生活中，運營商主要通過安裝基站的方式增強手機信號。那么問題來了，如何最優(yōu)選擇基站的位置將成為重點考慮的問題？

目前，比較有代表性的就是射線跟蹤模型。射線跟蹤是一種被廣泛用于移動通信中的預(yù)測無線電波傳播特性的技術(shù)，由于移動通信中使用的超高頻微波和光同屬電磁波，有一定近似性，按光學(xué)方法辨認出多路徑信道中收、發(fā)射機間所有主要的傳播路徑。一旦這些傳播路徑被辨認后，就可根據(jù)電波傳播理論來計算每條傳播路徑信號的幅度、相位、延遲和極化，然后結(jié)合天線方向圖和系統(tǒng)帶寬就可得到到達接收點的所有傳播路徑的相干合成結(jié)果。

對于城市基站的二維模型，建筑群可被劃分為一定的“塊”，建筑物(即下圖中帶有灰色陰影的多邊形)則被定義為“多邊形”，多邊形的“邊”代表建筑物的表面，多邊形的“頂點”則代表了建筑物的拐角。這種簡化了的市區(qū)平面圖大致反映出城市的主體結(jié)構(gòu)，利用它進行射線跟蹤，可以得到較為準確的路徑損耗。

在多邊形的頂點上僅能產(chǎn)生繞射，而在多邊形的邊上僅能產(chǎn)生反射，這些多次的反射、繞射及其組合便是收、發(fā)射機間的傳播路徑。二維射線跟蹤模型可以通過以下兩種規(guī)律分別確定反射傳播路徑和繞射傳播路徑：

(1)反射傳播路徑，如下圖(a)所示，產(chǎn)生反射時入射角 等于反射角 ;

(2)繞射傳播路徑，如下圖(b)所示，不論入射線以任意角度入射到建筑物頂點上，繞射射線都會以任意出射角向沒有建筑物覆蓋的區(qū)域傳播。

現(xiàn)在我們要求出，在發(fā)射機Tx 即坐標為(500, 200)、接收機Rx 即坐標為(250, 350)之間，通過多次反射與繞射的組合，從而形成的主要傳播路徑。

拿到這樣的題目，我們通過“枚舉法”，即發(fā)射機Tx，從0°逐步加上步長（比如0.1°，甚至可以更小，步長越小精度越高，但是程序運行時間更長），從這個角度發(fā)出一條射線，然后考慮是反射還是繞射情況，直到最后的射線與接收機Rx非常接近（可以利用點到直線的距離公式求得），此次循環(huán)才能結(jié)束，進入下一個循環(huán)，最終當Tx角度超過360°，程序才能結(jié)束。

下面，我們通過Matlab編程，運行主程序main.m文件，由于該程序涉及到大量矩陣求逆運算等過程，計算量大，在筆者電腦上面運行該程序大約需要15分鐘（取決于電腦配置），執(zhí)行過程中請大家耐心等待哈。比如，點擊運行程序后，可以先去洗個澡，回來后就能看到運行結(jié)果啦

最后，我們得到運行結(jié)果為：

詳細代碼如下：

（1）主程序main.m詳細代碼

clc;
clear;
close all;


figure
building={ [5.54892205638474e+002    1.53797468354430e+002;
            4.03648424543947e+002    1.53797468354430e+002;
            4.03648424543947e+002    1.00632911392405e+002;
            5.54892205638474e+002    1.00632911392405e+002;
            5.54892205638474e+002    1.53797468354430e+002];

            [5.54892205638474e+002    1.75738396624473e+002;
            5.54892205638474e+002    1.90928270042194e+002;
            5.29684908789386e+002    1.90928270042194e+002;
            5.29684908789386e+002    1.75738396624473e+002;
            5.54892205638474e+002    1.75738396624473e+002];

            [5.54892205638474e+002    1.95147679324895e+002;
            5.54892205638474e+002    2.59282700421941e+002;
            5.36981757877280e+002    2.59282700421941e+002;
            5.36981757877280e+002    1.95147679324895e+002;
            5.54892205638474e+002    1.95147679324895e+002];

            [4.42786069651741e+002    1.75738396624473e+002;
            4.42786069651741e+002    2.59282700421941e+002;
            4.03648424543947e+002    2.59282700421941e+002;
            4.03648424543947e+002    1.75738396624473e+002;
            4.42786069651741e+002    1.75738396624473e+002];

            [5.54892205638474e+002    2.81223628691983e+002;
            5.54892205638474e+002    3.00632911392405e+002;
            4.99170812603648e+002    3.00632911392405e+002;      %%%%%%淇敼
            4.99170812603648e+002    2.81223628691983e+002;
            5.54892205638474e+002    2.81223628691983e+002];

            [4.89220563847430e+002    2.81223628691983e+002;
            4.89220563847430e+002    2.95569620253165e+002;
            4.56053067993366e+002    2.95569620253165e+002;
            4.56053067993366e+002    2.81223628691983e+002;
            4.89220563847430e+002    2.81223628691983e+002];

            [4.32172470978441e+002    2.81223628691983e+002;
            4.32172470978441e+002    3.39451476793249e+002;
            4.03648424543947e+002    3.39451476793249e+002;
            4.03648424543947e+002    2.81223628691983e+002;
            4.32172470978441e+002    2.81223628691983e+002];

            [5.54892205638474e+002    3.09071729957806e+002;
            5.54892205638474e+002    3.39451476793249e+002;
            4.47429519071310e+002    3.39451476793249e+002;
            4.47429519071310e+002    3.09071729957806e+002;
            5.54892205638474e+002    3.09071729957806e+002];

            [4.75953565505804e+002    3.62236286919831e+002;
            4.75953565505804e+002    3.79113924050633e+002;
            4.71310116086236e+002    3.79113924050633e+002;
            4.71310116086236e+002    3.62236286919831e+002;
            4.75953565505804e+002    3.62236286919831e+002];

            [4.62686567164179e+002    3.62236286919831e+002;
            4.62686567164179e+002    4.18776371308017e+002;
            4.01658374792703e+002    4.18776371308017e+002;
            4.01658374792703e+002    3.62236286919831e+002;
            4.62686567164179e+002    3.62236286919831e+002];

            [4.99834162520730e+002    3.88396624472574e+002;
            4.99834162520730e+002    4.18776371308017e+002;
            4.80597014925373e+002    4.18776371308017e+002;
            4.80597014925373e+002    3.88396624472574e+002;
            4.99834162520730e+002    3.88396624472574e+002];

            [3.79104477611940e+002    4.18776371308017e+002;
            2.94195688225539e+002    4.18776371308017e+002;
            2.94195688225539e+002    3.62236286919831e+002;
            3.79104477611940e+002    3.62236286919831e+002;
            3.79104477611940e+002    4.18776371308017e+002];

            [3.79104477611940e+002    2.81223628691983e+002;
            3.79104477611940e+002    3.39451476793249e+002;
            2.12603648424544e+002    3.39451476793249e+002;
            2.12603648424544e+002    2.81223628691983e+002;
            3.79104477611940e+002    2.81223628691983e+002];

            [3.79104477611940e+002    1.53797468354430e+002;
            3.51243781094527e+002    1.53797468354430e+002;
            3.51243781094527e+002    1.00632911392405e+002;
            3.79104477611940e+002    1.00632911392405e+002;
            3.79104477611940e+002    1.53797468354430e+002];

            [3.79104477611940e+002    1.75738396624473e+002;
            3.79104477611940e+002    2.59282700421941e+002;
            3.32669983416252e+002    2.59282700421941e+002;
            3.32669983416252e+002    2.16244725738397e+002;
            3.57213930348259e+002    2.16244725738397e+002;
            3.57213930348259e+002    1.75738396624473e+002;
            3.79104477611940e+002    1.75738396624473e+002];

            [3.38640132669983e+002    1.53797468354430e+002;
            2.41791044776119e+002    1.53797468354430e+002;
            2.41791044776119e+002    1.25949367088608e+002;
            3.02819237147595e+002    1.25949367088608e+002;
            3.02819237147595e+002    1.00632911392405e+002;
            3.27363184079602e+002    1.00632911392405e+002;
            3.27363184079602e+002    1.25949367088608e+002;
            3.38640132669983e+002    1.25949367088608e+002;
            3.38640132669983e+002    1.53797468354430e+002];

            [2.90215588723051e+002    1.75738396624473e+002;
            2.90215588723051e+002    2.08649789029536e+002;
            2.65008291873964e+002    2.08649789029536e+002;
            2.65008291873964e+002    2.59282700421941e+002;
            2.12603648424544e+002    2.59282700421941e+002;
            2.12603648424544e+002    2.33966244725738e+002;
            2.43117744610282e+002    2.33966244725738e+002;
            2.43117744610282e+002    2.01054852320675e+002;
            2.12603648424544e+002    2.01054852320675e+002;
            2.12603648424544e+002    1.75738396624473e+002;
            2.90215588723051e+002    1.75738396624473e+002];

            };
N=17;


for i=1:N
    plot(building{i}(:,1),building{i}(:,2));
    fill(building{i}(:,1),building{i}(:,2),[0.5 0.5 0.5],'FaceAlpha',0.8);
    hold on
end
axis equal
xlabel('X(m)');
ylabel('Y(m)');
%grid on


TX=[500;200];
plot(TX(1),TX(2),'*r');
h_text1 = text(TX(1)-6,TX(2)-10,['Tx'],'color','r');
set(h_text1,'fontsize',14)    %設(shè)置字體大小


RX=[250;350];
plot(RX(1),RX(2),'or');     
h_text2 = text(RX(1)-5,RX(2)+10,['Rx'],'color','r');
set(h_text2,'fontsize',14)    %設(shè)置字體大小    


[aa]=calculate_point(TX(1),TX(2),building);
num1=size(aa);
num=num1(1);
for i=1:num
    [bb]=calculate_point(aa(i,1),aa(i,2),building);
    num2=size(bb);
    number=num2(1);
    for j=1:number
         for thea=0:0.1:2*pi
            x1=bb(j,1);y1=bb(j,2);x2=x1+cos(thea);y2=y1+sin(thea); 
            node=calculate_node(x1,y1,building);
            if x2<=max(building{node}(:,1))&&x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
               [insect1,insect2]=insection_point(x1,y1,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   continue
               end
               value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
               zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
               if zhi >=0.95&&zhi<=1&&value1==1
                   plot([TX(1),aa(i,1)],[TX(2),aa(i,2)],'k');
                   plot([aa(i,1),bb(j,1)],[aa(i,2),bb(j,2)],'k');
                   plot([bb(j,1),insect1(1)],[bb(j,2),insect1(2)],'k');
                   plot([insect1(1),RX(1)],[insect1(2),RX(2)],'k');
               end
        end
    end     
end


for thea=0:0.1:2*pi
    x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
    [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
    if insect1(1)==0&&insect1(2)==0
        continue
    end
    x1=insect1(1);y1=insect1(2);
    x2=insect2(1);y2=insect2(2);    
    [x,y]=line_to_point(x1,y1,x2,y2,building);
    if x==0&&y==0
        continue
    end
    [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
    if x1==0&&y1==0
        continue
    end
   [bb]=calculate_point(x,y,building);
    num2=size(bb);
    number=num2(1);
    for j=1:number     
        value1=calculate(bb(j,1),bb(j,2),RX(1),RX(2),building);
        if value1==1
            plot([TX(1),x1],[TX(2),y1],'r');
            plot([x1,x],[y1,y],'r');
            plot([x,bb(j,1)],[y,bb(j,2)],'r');
            plot([bb(j,1),RX(1)],[bb(j,2),RX(2)],'r');
        end
    end  
end        


[aa]=calculate_point(TX(1),TX(2),building);
num1=size(aa);
num=num1(1);
for i=1:num
    for thea=0:0.1:2*pi
        x1=aa(i,1);y1=aa(i,2);x2=x1+cos(thea);y2=y1+sin(thea);     
        node=calculate_node(x1,y1,building);
        if x2<=max(building{node}(:,1))&&x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
            continue
        end
        [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
        if insect1(1)==0&&insect1(2)==0
            continue
        end
        x1=insect1(1);y1=insect1(2);
        x2=insect2(1);y2=insect2(2);    
        [x,y]=line_to_point(x1,y1,x2,y2,building);
        if x==0&&y==0
            continue
        end
        [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
        if x1==0&&y1==0
            continue
        end  
        value1=calculate(x,y,RX(1),RX(2),building);
        if value1==1
            plot([TX(1),aa(i,1)],[TX(2),aa(i,2)],'r');
            plot([aa(i,1),x1],[aa(i,2),y1],'r');
            plot([x1,x],[y1,y],'r');
            plot([x,RX(1)],[y,RX(2)],'r');
        end
    end  
end

[aa]=calculate_point(TX(1),TX(2),building);
num1=size(aa);
num=num1(1);
for i=1:num
    [bb]=calculate_point(aa(i,1),aa(i,2),building); 
    num2=size(bb);
    number=num2(1);
    for j=1:number
        value1=calculate(bb(j,1),bb(j,2),RX(1),RX(2),building);
        if value1==1
            plot([TX(1),aa(i,1)],[TX(2),aa(i,2)],'g');
            plot([aa(i,1),bb(j,1)],[aa(i,2),bb(j,2)],'g');
            plot([bb(j,1),RX(1)],[bb(j,2),RX(2)],'g');
        end
    end
end      


[aa]=calculate_point(TX(1),TX(2),building);
num1=size(aa);
num=num1(1);
for i=1:num  
    for thea=0:0.1:2*pi
         for k=1:4
            x1=aa(i,1);y1=aa(i,2);
            x2=x1+cos(thea);y2=y1+sin(thea);   
            node=calculate_node(x1,y1,building);
            if x2<=max(building{node}(:,1)) && x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
            for j=1:k
               [insect1,insect2]=insection_point(x1,y1,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   break
               end
               xm{j}=insect1;ym{j}=insect2;
               x1=insect1(1);y1=insect1(2);
               x2=insect2(1);y2=insect2(2); 
            end
            value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
            zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
            if zhi >=0.95&&zhi<=1&&value1==1
               plot([TX(1),aa(i,1)],[TX(2),aa(i,2)],'k');
               xx=aa(i,1);yy=aa(i,2);
               for m=1:j
                  plot([xx,xm{m}(1)],[yy,xm{m}(2)],'k');
                  xx=xm{m}(1);yy=xm{m}(2);
               end
               plot([xx,RX(1)],[yy,RX(2)],'k');
           end            
         end
    end
end


 for thea=0:0.1:2*pi
    for thea1=0:0.1:2*pi
         for k=1:3

             x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);    
            [x,y]=line_to_point(x1,y1,x2,y2,building);
            if x==0&&y==0
                continue
            end
            [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
            if x1==0&&y1==0
                continue
            end
            aa1=x1;bb1=y1;
            aa2=x;bb2=y;
            x2=x+cos(thea1);y2=y+sin(thea1);   
            node=calculate_node(x,y,building);
            if x2<=max(building{node}(:,1)) && x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
            for j=1:k
               [insect1,insect2]=insection_point(x,y,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   break
               end
               xm{j}=insect1;ym{j}=insect2;
               x=insect1(1);y=insect1(2);
               x2=insect2(1);y2=insect2(2); 
            end
            value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
            zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
            if zhi >=0.95&&zhi<=1&&value1==1
               plot([TX(1),aa1],[TX(2),bb1],'k');
               plot([aa1,aa2],[bb1,bb2],'k');
               xx=aa2;yy=bb2;
               for m=1:j
                  plot([xx,xm{m}(1)],[yy,xm{m}(2)],'k');
                  xx=xm{m}(1);yy=xm{m}(2);
               end
               plot([xx,RX(1)],[yy,RX(2)],'k');
           end            
         end
    end    
 end           


 for thea=0:0.1:2*pi
    for thea1=0:0.1:2*pi
         for k=1:2

             x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);  
            aa3=x1;bb3=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 

            [x,y]=line_to_point(x1,y1,x2,y2,building);
            if x==0&&y==0
                continue
            end
            [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
            if x1==0&&y1==0
                continue
            end
            aa1=x1;bb1=y1;
            aa2=x;bb2=y;
            x2=x+cos(thea1);y2=y+sin(thea1);   
            node=calculate_node(x,y,building);
            if x2<=max(building{node}(:,1)) && x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
            for j=1:k
               [insect1,insect2]=insection_point(x,y,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   break
               end
               xm{j}=insect1;ym{j}=insect2;
               x=insect1(1);y=insect1(2);
               x2=insect2(1);y2=insect2(2); 
            end
            value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
            zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
            if zhi >=0.95&&zhi<=1&&value1==1
               plot([TX(1),aa3],[TX(2),bb3],'k');
               plot([aa3,aa1],[bb3,bb1],'k');
               plot([aa1,aa2],[bb1,bb2],'k');
               xx=aa2;yy=bb2;
               for m=1:j
                  plot([xx,xm{m}(1)],[yy,xm{m}(2)],'k');
                  xx=xm{m}(1);yy=xm{m}(2);
               end
               plot([xx,RX(1)],[yy,RX(2)],'k');
           end            
         end
    end    
 end 
 
 for thea=0:0.1:2*pi
    for thea1=0:0.1:2*pi
         for k=1:1

             x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);  
            aa4=x1;bb4=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 
            aa3=x1;bb3=y1;  

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 

            [x,y]=line_to_point(x1,y1,x2,y2,building);
            if x==0&&y==0
                continue
            end
            [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
            if x1==0&&y1==0
                continue
            end
            aa1=x1;bb1=y1;
            aa2=x;bb2=y;
            x2=x+cos(thea1);y2=y+sin(thea1);   
            node=calculate_node(x,y,building);
            if x2<=max(building{node}(:,1)) && x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
            for j=1:k
               [insect1,insect2]=insection_point(x,y,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   break
               end
               xm{j}=insect1;ym{j}=insect2;
               x=insect1(1);y=insect1(2);
               x2=insect2(1);y2=insect2(2); 
            end
            value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
            zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
            if zhi >=0.95&&zhi<=1&&value1==1
               plot([TX(1),aa4],[TX(2),bb4],'k');
               plot([aa4,aa3],[bb4,bb3],'k');
               plot([aa3,aa1],[bb3,bb1],'k');
               plot([aa1,aa2],[bb1,bb2],'k');
               xx=aa2;yy=bb2;
               for m=1:j
                  plot([xx,xm{m}(1)],[yy,xm{m}(2)],'k');
                  xx=xm{m}(1);yy=xm{m}(2);
               end
               plot([xx,RX(1)],[yy,RX(2)],'k');
           end            
         end
    end    
 end         


 for thea=0:0.1:2*pi
    for thea1=0:0.1:2*pi
             x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);  
            aa5=x1;bb5=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);  
            aa4=x1;bb4=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 
            aa3=x1;bb3=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 

            [x,y]=line_to_point(x1,y1,x2,y2,building);
            if x==0&&y==0
                continue
            end
            [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
            if x1==0&&y1==0
                continue
            end
            aa2=x1;bb2=y1;
            aa1=x;bb1=y;
            value1=calculate(x,y,RX(1),RX(2),building);
            if value1==1
                plot([TX(1),aa5],[TX(2),bb5],'r');
                plot([aa5,aa4],[bb5,bb4],'r');
                plot([aa4,aa3],[bb4,bb3],'r');
                plot([aa3,aa2],[bb3,bb2],'r');
                plot([aa2,aa1],[bb2,bb1],'r');
                plot([aa1,RX(1)],[bb1,RX(2)],'r');
            end
    end      
 end

（2）calculate.m詳細代碼

function value=calculate(x2,y2,x3,y3,build)
    N=17;
    k=1;mm=0;
    value=0;
    while k >=1
        row_column1=size(build{k});
        row1=row_column1(1);
        for p=1:row1-1
            X1=[x3,y3];Y1=[x2,y2];X2=build{k}(p,:);Y2=build{k}(p+1,:);


            if X1(1)==Y1(1)
                X=X1(1);
                k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                b2=X2(2)-k2*X2(1); 
                Y=k2*X+b2;
            end
            if X2(1)==Y2(1)
                X=X2(1);
                k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                b1=X1(2)-k1*X1(1);
                Y=k1*X+b1;
            end
            if X1(1)~=Y1(1)&&X2(1)~=Y2(1)
                k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                b1=X1(2)-k1*X1(1);
                b2=X2(2)-k2*X2(1);
                if k1==k2
                    continue
                else
                    X=(b2-b1)/(k1-k2);
                    Y=k1*X+b1;
                end
            end


            Q=([X,Y]-X1)*(Y1-[X,Y])';
            W=([X,Y]-X2)*(Y2-[X,Y])';
            if (Q >0.001)&&(W >0.001)
                mm=1;
            end
        end 
        if mm==1
            mm=0;
            k=1;
            break
        end    
        k=k+1;     
        if k >N
            break
        end
    end
    if k >N
        value=1;
        k=1;
    end  
   end

（3）calculate_node.m詳細代碼

function node=calculate_node(x1,y1,build)
    N=17;m=0;node=0;
    for k=1:N
        row_column1=size(build{k});
        row1=row_column1(1);
        for p=1:row1-1
            x2=build{k}(p,1);y2=build{k}(p,2);
         if abs(x1-x2)<=0.1&&abs(y1-y2)<=0.1
             node=k;
             m=1;
             break
         end
         if m==1
             break 
         end
        end  
    end
end

（4）calculate_point.m詳細代碼

function [x1]=calculate_point(x2,y2,build)
    N=17;
    k=1;m=1;mm=0;dd=0;
    node=calculate_node(x2,y2,build);
    for i=1:N
        if i==node
            continue
        end
        row_column=size(build{i});
        row=row_column(1);

        value=findway(x2,y2,build);
        if value==i
            continue
        end

        for j=1:row-1
            x=build{i}(j,1);
            y=build{i}(j,2);
            while k >=1
                row_column1=size(build{k});
                row1=row_column1(1);
                for p=1:row1-1
                    X1=[x,y];Y1=[x2,y2];X2=build{k}(p,:);Y2=build{k}(p+1,:);

                    if (X1(1)==Y1(1)&&X1(2)==Y1(2))  ||  (X2(1)==Y2(1)&&X2(2)==Y2(2))
                        continue
                    end
                    if X1(1)==Y1(1)
                        X=X1(1);
                        k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                        b2=X2(2)-k2*X2(1); 
                        Y=k2*X+b2;
                    end
                    if X2(1)==Y2(1)
                        X=X2(1);
                        k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                        b1=X1(2)-k1*X1(1);
                        Y=k1*X+b1;
                    end
                    if X1(1)~=Y1(1)&&X2(1)~=Y2(1)
                        k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                        k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                        b1=X1(2)-k1*X1(1);
                        b2=X2(2)-k2*X2(1);
                        if k1==k2
                            continue
                        else
                            X=(b2-b1)/(k1-k2);
                            Y=k1*X+b1;
                        end
                    end

                    Q=([X,Y]-X1)*(Y1-[X,Y])';
                    W=([X,Y]-X2)*(Y2-[X,Y])';
                    if (Q >0.001)&&(W >0.001)
                        mm=1;
                    end
                end 
                if mm==1
                    mm=0;
                    k=1;
                    break
                end    
                k=k+1;     
                if k >N
                    break
                end
            end
            if k >N
                for ii=1:N
                    row_column11=size(build{ii});
                    row11=row_column11(1);
                    for pp=1:row11-1
                        x22=build{ii}(pp,1);y22=build{ii}(pp,2);
                        QQ=([x22,y22]-Y1)*([x,y]-[x22,y22])'/norm([x22,y22]-Y1)/norm([x,y]-[x22,y22]);
                        if QQ==1
                            dd=1;
                            break
                        end
                    end
                    if dd==1
                            break
                    end
                end
                if dd==0
                    x1(m,1)=x;x1(m,2)=y;
                    m=m+1;
                end
                k=1;
                dd=0;
            end          
        end
    end
end

（5）findway.m詳細代碼

function value=findway(xx,yy,build)
    N=17;
    value=0;
    for i=1:N
        row_column=size(build{i});
        row=row_column(1);
        for j=1:row-1
            x=build{i}(j,1);
            y=build{i}(j,2);
            if x==xx&&y==yy
                value=i;
                break
            end
        end
        if value==i
            break
        end
    end
end

（6）insection_point.m詳細代碼

function [point1,point2,point3,a,b,c,d]=insection_point(x2,y2,x3,y3,build)
    N=17;min_num=10^10;
    point1(1)=0;point1(2)=0;point2(1)=0;point2(2)=0;point3(1)=0;point3(2)=0;a=0;b=0;c=0;d=0;
    for k=1:N
        row_column1=size(build{k});
        row1=row_column1(1);
        for p=1:row1-1
            X1=[x3,y3];Y1=[x2,y2];X2=build{k}(p,:);Y2=build{k}(p+1,:);
            if X1(1)==Y1(1)
                X=X1(1);
                k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                b2=X2(2)-k2*X2(1); 
                Y=k2*X+b2;
            end
            if X2(1)==Y2(1)
                X=X2(1);
                k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                b1=X1(2)-k1*X1(1);
                Y=k1*X+b1;
            end
            if X1(1)~=Y1(1)&&X2(1)~=Y2(1)
                k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                b1=X1(2)-k1*X1(1);
                b2=X2(2)-k2*X2(1);
                if k1==k2
                    continue
                else
                    X=(b2-b1)/(k1-k2);
                    Y=k1*X+b1;
                end
            end


            Q=(X1-Y1)*([X,Y]-Y1)';
            W=([X,Y]-X2)*(Y2-[X,Y])';
            if (Q >0.001)&&(W >0.001)
                num1=(x2-X)^2+(y2-Y)^2;
               if min_num >=num1
                   min_num=num1;
                   point1=[X,Y];
                   a=build{k}(p,1);b=build{k}(p,2);
                   c=build{k}(p+1,1);d=build{k}(p+1,2);
                       if a==c
                            point2=[2*a-x2,y2];
                            point3=point2;
                            point2=2.*point1-point2;
                       end
                       if b==d
                            point2=[x2,2*d-y2];
                            point3=point2;
                            point2=2.*point1-point2;
                       end
               end              
            end            
        end  
    end
end

（7）line_to_point.m詳細代碼

function [x,y]=line_to_point(x1,y1,x2,y2,build)
    x=0;y=0;
    N=17;
    for k=1:N
        row_column1=size(build{k});
        row1=row_column1(1);
        for p=1:row1-1
            x3=build{k}(p,1);y3=build{k}(p,2);
            value1=calculate(x1,y1,x3,y3,build);
            zhi=([x3,y3]-[x1,y1])*([x2,y2]-[x1,y1])'/norm([x3,y3]-[x1,y1])/norm([x2,y2]-[x1,y1]);
            if zhi >=0.99&&zhi<=1&&value1==1
                x=x3;y=y3;
            end
        end
    end
end

（8）point_in.m詳細代碼

function [value1,value2]=point_in(x2,y2,x3,y3,a,b,c,d)
    value1=0;value2=0;
    X1=[x3,y3];Y1=[x2,y2];X2=[a,b];Y2=[c,d];
    if X1(1)==Y1(1)
        X=X1(1);
        k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
        b2=X2(2)-k2*X2(1); 
        Y=k2*X+b2;
    end
    if X2(1)==Y2(1)
        X=X2(1);
        k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
        b1=X1(2)-k1*X1(1);
        Y=k1*X+b1;
    end
    if X1(1)~=Y1(1)&&X2(1)~=Y2(1)
        k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
        k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
        b1=X1(2)-k1*X1(1);
        b2=X2(2)-k2*X2(1);
        if k1==k2
            value1=0;value2=0;
        else
            X=(b2-b1)/(k1-k2);
            Y=k1*X+b1;
        end
    end
    Q=([X,Y]-X1)*(Y1-[X,Y])';
    W=([X,Y]-X2)*(Y2-[X,Y])';
    if (Q >0.001)&&(W >0.001)
       value1=X;value2=Y;
    end
end