【电磁波】基于matlab大气波导电磁波射频追踪仿真【含Matlab源码 3559期】.zip

%% 大气传播射线跟踪法 %clc;clear all; function xsitax = Main_ray(M0,hmax,ht,HH1,HH2,slopeElvated1,slopeElvated2,xmaxkm,sitamax) %% 参数说明 % hmax: 最大显示高度(m)；dh：高度步长(m)；ht：发射天线高度(m)； % HH1：波导基础层高度(m)；HH2：波导层+基础层高度/蒸发波导高度(m)； % slopeElvated1：基础层斜率(/km)；slopeElvated2：陷获层斜率(/km)； % xmaxkm：最大显示距离(km);ae:等效地球半径(m) % sita0：天线初始仰角范围(°) %% 参数输入 hmax = round(hmax);%%四舍五入返回 %hmax= 500; %M0=330;sitamax=4; dh = 1; h = 0:dh:hmax; %ht = 20; %HH1=0;%基础层高度(m) %HH2=50;%波导层+基础层高度(m) %slopeElvated1=118;%基础层斜率/km %slopeElvated2=-120;%陷获层斜率/km %xmaxkm = 200; %% 折射率剖面和折射率梯度 %M = Sub_evaporation(h,HH2,M0); M = Sub_FourParaRefPro(M0,h,HH1,HH2,slopeElvated1,slopeElvated2); %M = Sub_StanRefPro(h); % M = hunhebodao(M0,h,HH1,HH2,HH3,slopeElvated1,slopeElvated2); for i = 1:1:length(h)-1 g(i) = (M(i+1)-M(i))/(h(i+1)-h(i)); end g(i+1) = g(i); ae = 6371*1000; %% 计算 sitax = zeros(1,2); iii = 2; sitax(1,:) = [0,0.0000001]; sitamax = sitamax*180/pi; for sita0 =-0.6:0.05:0.6 id = find(ht == h); %%返回ht=h时的i值 gi = g(id); gi1 = g(id-1); dhi = dh; sita1 = sita0*pi/180; x1 = 0; h1 = ht; u = []; u(1,1) = x1; u(1,2) = h1; for ii = 1:1:hmax if gi == 0 %% 梯度为0时 dhi = dh; sita2 = sita1; h2 = h1+dhi; x2 = x1+cot(sita2)*dh; id = find(round(h2) == h); gi = g(id); gi1 = g(id-1); else %% 梯度不为0时 if sita1>0 && sita1<= (3*pi/180) %% 仰角大于0小于3°时，射线向上传播 r = sita1*sita1+2*gi*dhi; if r>=0 dhi = min(h(id+1)-h1,dh); sita2 = sqrt(sita1*sita1+2*gi*dhi); h2 = h1+dhi; x2 = x1+(sita2-sita1)/gi; else sita2 = 0; dhi = (sita2*sita2-sita1*sita1)/2/gi; h2 = h1+dhi; x2 = x1+(sita2-sita1)/gi; end id = min(find( h >= h2)); gi = g(id); gi1 = g(id-1); elseif sita1<0 && sita1>= -3*pi/180 %% 仰角小于0大于-3°时，射线向下传播 r = sita1*sita1-2*gi1*dhi; if r>=0 if h1 == 0 dhi = 0; sita2 = -sita1; %% 擦地角 h2 = h1; x2 = x1; id = min(find(h >= h2)); gi = g(id); sitax(iii,:) = [x1/1000,sita2]; %%保存擦地角 iii = iii+1; else dhi = min(h1-h(id-1),dh); sita2 = -sqrt(sita1*sita1-2*gi1*dhi); h2 = h1-dhi; x2 = x1+(sita2-sita1)/gi1; id = min(find(h >= h2)); if id == 1 gi = g(id); gi1 = g(id); else gi = g(id); gi1 = g(id-1); end end else sita2 = 0; dhi = (sita2*sita2-sita1*sita1)/2/gi1; h2 = h1+dhi; x2 = x1+(sita2-sita1)/gi1; id = min(find(h >= h2)); gi = g(id); end elseif sita1 == 0 %% 仰角为0 if gi>0 %% 射线向上传播 r = sita1*sita1+2*gi*dhi; if r>=0 dhi = min(h(id+1)-h1,dh); sita2 = sqrt(sita1*sita1+2*gi*dhi); h2 = h1+dhi; x2 = x1+(sita2-sita1)/gi; else sita2 = 0; dhi = (sita2*sita2-sita1*sita1)/2/gi; h2 = h1+dhi; x2 = x1+(sita2-sita1)/gi; end id = min(find(h >= h2)); if id == 1 gi = g(id); gi1 = g(id); else gi = g(id); gi1 = g(id-1); end elseif gi<0 %% 射线向下传播 r = sita1*sita1-2*gi1*dhi; if r>=0 if h1 == 0 dhi = 0; sita2 = -sita1; h2 = h1; x2 = x1; id = min(find(h >= h2)); gi = g(id+1); sitax(iii,:) = [x1/1000;sita2]; %%保存擦地角 iii = iii+1; else dhi = min(h1-h(id-1),dh); sita2 = -sqrt(sita1*sita1-2*gi1*dhi); h2 = h1-dhi; x2 = x1+(sita2-sita1)/gi1; id = min(find(h >= h2)); if id == 1 gi = g(id); gi1 = g(id); else gi = g(id); gi1 = g(id-1); end end else sita2 = 0; dhi = (sita2*sita2-sita1*sita1)/2/gi1; h2 = h1+dhi; x2 = x1+(sita2-sita1)/gi1; id = min(find(h >= h2)); gi = g(id); end end elseif sita1>(3*pi/180) %% 仰角大于3°时 dhi = min(h(id+1)-h1,dh); h2 = h1+dhi; sq = sqrt(M(id+1)^2-(M(1)*cos(sita1))^2); x2 = x1+(ae/(ae+h2))*(M(1)*cos(sita1)*dhi)/sq; sita2 = atan(dhi/(x2-x1)); id = min(find(h >= h2)); gi = g(id); gi1 = g(id-1); elseif sita1<(-3*pi/180) %% 仰角小于3°时 if h1 == 0 dhi = 0; sita2 = -sita1; h2 = h1; x2 = x1; id = min(find(h >= h2)); gi = g(id); else dhi = min(h1-h(id-1),dh); h2 = h1-dhi; sq = sqrt(M(id+1)^2-(M(1)*cos(sita1))^2); x2 = x1+(ae/(ae+h2))*(M(1)*cos(sita1)*dhi)/sq; sita2 = -atan(dhi/(x2-x1)); id = min(find(h >= h2)); if id == 1 gi = g(id); gi1 = g(id); else gi = g(id); gi1 = g(id-1); end end end end u(ii,1) = x2; u(ii,2) = h2; x1 = x2; h1 = h2; sita1 = sita2; if x1>xmaxkm*1000 break; end if h1 >= hmax-1 break; end end figure(1); plot(u(:,1)/1000,u(:,2),'r'); hold on; xlabel('距离(km)'); ylabel('高度(m)'); end xsitax = sortrows(sitax,1); if xsitax(end,1)<xmaxkm xsitax(end+1,1) = xmaxkm; xsitax(end,2) = 0.000001; end