function?Power_all?=?get_rss_by_ray_tracing（room_x?room_y?room_z?source_x?source_y?source_z?grid_size?f）
%?射線跟蹤：?輸入房間大小，信號源的坐標，設置網格大小和信號頻率，輸出每個網格點上接收到的信號強度
%?輸入參數分別為：房間尺寸x，?y，?z，信號源（天線）坐標（x，?y，?z），網格大小（輸出數據的密度），單位（m）；信號的發射頻率，單位（MHz）

????if?nargin?==?0?%無輸入參數時的默認設置
????????room_x?=?20;
????????room_y?=?15;
????????room_z?=?4;
????????source_x?=?10;
????????source_y?=?7.5;
????????source_z?=?1;
????????grid_size?=?0.1;
????????f?=?2400;
????end

????room_x?=?1000?*?room_x;
????room_y?=?1000?*?room_y;
????room_z?=?1000?*?room_z;
????source_x?=?1000?*?source_x;
????source_y?=?1000?*?source_y;
????source_z?=?1000?*?source_z;
????grid_size?=?1000?*?grid_size;
????
????%介電常數和導電系數，參考射線跟蹤相關的論文
????epsilon_c=10-1.2j;
????epsilon_w=6.1-1.2j;
????%?epsilon_c=7.9-0.89j;
????%?epsilon_w=6.2-0.69j;

????if?room_x?*?room_y?/?grid_size^2?>?100000000
????????disp（‘提示：程序使用大矩陣同時計算所有位置點上的信號強度，位置點設置的過多可能導致內存不足，matlab可能會卡死。‘）;
????????input（‘繼續請回車；退出請ctrl+c‘）;
????end

????T?=?1?/?（f?*?10^6）;
????c?=?3.0e8;
????lambda=c?/?（f?*?10^6）;

????%得到空間中所有的網格點位置坐標把z固定，和信號源一個高度
????[X?Y]?=?meshgrid（grid_size:grid_size:（room_x-grid_size）?grid_size:grid_size:（room_y-grid_size））;
????L?=?[X（:）?Y（:）];
????L?=?[L?zeros（size（X（:）））?+?source_z];
????
????%%?直射路徑
????%將反射引起的相位變化加入電場E的計算里面。后面注釋所說的相位只是由于距離引起的相位。
????d_direct?=?sqrt（（L（:1）?-?source_x）.^2?+?（L（:2）-source_y）.^2?+?（L（:3）?-?source_z）.^2）;%每個網格點距發射源的歐式距離
????t_direct_0?=?d_direct./1000./c;%直射時延
????p_direct?=?mod（t_direct_0*2*pi/T2*pi）;%直射相位
????E_direct?=?（lambda./（4.*pi.*d_direct./1000））;%這里和下面的E不一定正好是場強大小，但和場強成正比
????E0=E_direct.*?exp（1i.*（-p_direct））;

????%%
????%上面計算的L為所有網格點的坐標。每行為一組坐標。
????%下面的Li相應的為所有網格點所對應的鏡像點
????%下面對六組反射路徑分別計算

????%%?前平面反射路徑（前后左右上下的意思是：人站才這個長方體中，面朝y軸，這時的六個面分別稱為前后左右上下）
????Li=[L（:1）??2.*room_y-L（:2）??L（:3）];?%計算鏡像點
????d_reflect?=?sqrt（（Li（:1）-source_x）.^2+（Li（:2）-source_y）.^2+（Li（:3）-source_z）.^2）;%反射路徑總長度
????t_reflect_1?=?d_reflect./1000./c;%時延
????p_reflect?=?mod（t_reflect_1*2*pi/T2*pi）;%相位
????thet?=?abs（atan（（Li（:2）-source_y）./（Li（:1）-source_x）））;%入射角
????reflect_coefficient?=?（sin（thet）-sqrt（epsilon_w-（cos（thet））.^2））./（sin（thet）+sqrt（epsilon_w-（cos（thet））.^2））;%反射系數也是矩陣
????E_reflect?=?（lambda./（4.*pi.*d_reflect./1000））?.*??reflect_coefficient;
????E1=E_reflect?.*?exp（1i.*（-p_reflect））;%將延遲的相位加進來，與反射造成的衰減與相位變化合在一起。

????%%?后平面反射路徑
????Li=[L（:1）??-L（:2）??L（:3）];?%計算鏡像點
????d_reflect?=?sqrt（（Li（:1）-source_x）.^2+（Li（:2）-source_y）.^2+（Li（:3）-source_z）.^2）;%反射路徑總長度
????t_reflect_2?=?d_reflect./1000./c;%時延
????p_reflect?=?mod（t_reflect_2*2*pi/T2*pi）;%相位
????thet?=?abs（atan（（Li（:2）-source_y）./（Li（:1）-source_x）））;%入射角
????reflect_coefficient?=?（sin（thet）-sqrt（epsilon_w-（cos（thet））.^2））./（sin（thet）+sqrt（epsilon_w-（cos（thet））.^2））;
????E_reflect?=?（lambda./（4.*pi.*d_reflect./1000））?.*??reflect_coefficient;
????E2=E_reflect?.*?exp（1i.*（-p_reflect））;%將延遲的相位加進來，與反射造成的衰減與相位變化合在一起。

????