电磁波模型
一、1、插值函数
clear all x=0:2*pi; y=sin(x); xx=0:0.5:2*pi; %interp1对sin函数进行分段线性插值,调用interp1的时候,默认的是分段线性插值 y1=interp1(x,y,xx); figure plot(x,y,'o',xx,y1,'*') title('分段线性插值') legend('初始值','分段线性插值') %临近插值 y2=interp1(x,y,xx,'nearest'); figure plot(x,y,'o',xx,y2,'*'); title('临近插值') legend('初始值','临近插值') %球面线性插值 y3=interp1(x,y,xx,'spline'); figure plot(x,y,'o',xx,y3,'*') title('球面插值') legend('初始值','球面插值') %三次多项式插值法 y4=interp1(x,y,xx,'cubic'); figure plot(x,y,'o',xx,y4,'*'); title('三次多项式插值') legend('初始值','三次多项式插值')
二、
1、射线模型
%program:Ray tracing routines 蒸发波导射线模型 clc;close;clear %% %************************************初始化模型************************** %-----------基本参数----------- tx_height=27;%天线高度(m) rmax_user=50*1e3;%传输范围(m) zmax_user=60;%最大高度(m)??? zdelta=0.5;%步进值(m) height=(0:floor(zmax_user*2/zdelta))*zdelta;%离散高度 %-----------蒸发波导-------- M1=[0 315;40 295;100 302];%修正折射率剖面????线性的修正折射率剖面。波导仿真的修改条件 m1=interp1(M1(:,1),M1(:,2),height,'linear','extrap')*1e-6+1;%插值法,从修正折射率剖面到修正折射指数(1.00025~1.0004) g1=(m1(2:end)-m1(1:end-1))/zdelta;%计算梯度 g值 %----------输入期望角度-------- angle=[4.6/1000;3.7/1000;3.6/1000;-3.6/1000;-3.7/1000;-4.6/1000;];%逆时针为正,角度为弧度制(mrad) N_angle=6;%角度个数 %% %**************************开始射线模型计算************************** x=zeros(N_angle,10000); z=zeros(N_angle,10000); steps=zeros(N_angle,1); for ii=1:N_angle StepNum=1; x(ii,StepNum)=0.0;%初始值 z(ii,StepNum)=tx_height;%天线高度 theta_i=angle(ii);%角度 h_i=abs(zdelta); while z(ii,StepNum)<=zmax_user && x(ii,StepNum)<=rmax_user %停止条件 if theta_i >0 %up going deltam=interp1(m1,(z(ii,StepNum)+zdelta)/zdelta+1)-interp1(m1,z(ii,StepNum)/zdelta+1);%折射指数差值 r=theta_i^2+2*deltam; if r>=0 h_i=min(zdelta,(floor(z(ii,StepNum)/zdelta)+1)*zdelta-z(ii,StepNum));%????? h_i=max(h_i,0.1*zdelta);%???? deltam=interp1(m1,(z(ii,StepNum)+h_i)/zdelta+1)-interp1(m1,z(ii,StepNum)/zdelta+1);%折射指数差值 theta_ip1=sqrt(theta_i^2+2*deltam);%出射角 StepNum=StepNum+1; z(ii,StepNum)=z(ii,StepNum-1)+h_i;%高度 x(ii,StepNum)=x(ii,StepNum-1)+(theta_ip1-theta_i)./(deltam/h_i);%距离 theta_i=theta_ip1; else theta_ip1=0; %到达上表面折射点 h_i=(theta_ip1^2-theta_i^2)./(2*g1(ceil(z(ii,StepNum)/zdelta)));%正 StepNum=StepNum+1 z(ii,StepNum)=z(ii,StepNum-1)+h_i; x(ii,StepNum)=x(ii,StepNum-1)+(theta_ip1-theta_i)/g1(ceil(z(ii,StepNum-1)/zdelta)); theta_i=theta_ip1; end elseif theta_i<0 % down going deltam=interp1(m1,(z(ii,StepNum)-zdelta)./zdelta+1)-interp1(m1,z(ii,StepNum)./zdelta+1); r=theta_i.^2+2*deltam; if z(ii,StepNum)<=0 %到达下表面折射点 theta_i=-theta_i; StepNum=StepNum+1; z(ii,StepNum)=z(ii,StepNum-1); x(ii,StepNum)=x(ii,StepNum-1); elseif r>=0 h_i=min(zdelta,z(ii,StepNum)-(ceil(z(ii,StepNum)/zdelta)-1)*zdelta); h_i=max(h_i,0.1*zdelta); deltam=interp1(m1,(z(ii,StepNum)-h_i)/zdelta+1)-interp1(m1,z(ii,StepNum)/zdelta+1); theta_ip1=-sqrt(theta_i^2+2*deltam); StepNum=StepNum+1; z(ii,StepNum)=z(ii,StepNum-1)-h_i; x(ii,StepNum)=x(ii,StepNum-1)+(theta_ip1-theta_i)./(deltam/(-h_i)); theta_i=theta_ip1; else theta_ip1=0; h_i=(theta_ip1^2-theta_i^2)./(2*g1(ceil(z(ii,StepNum)/zdelta)));%正 StepNum=StepNum+1; z(ii,StepNum)=z(ii,StepNum-1)-h_i; x(ii,StepNum)=x(ii,StepNum-1)+(theta_ip1-theta_i)/g1(ceil(z(ii,StepNum-1)/zdelta)); theta_i=theta_ip1; end else if g1(ceil(z(ii,StepNum)/zdelta))>0 h_i=zdelta; deltam=interp1(m1,(z(ii,StepNum)+h_i)/zdelta+1)-interp1(m1,z(ii,StepNum)/zdelta+1); theta_ip1=sqrt(theta_i^2+2*deltam); StepNum=StepNum+1; z(ii,StepNum)=z(ii,StepNum-1)+h_i; x(ii,StepNum)=x(ii,StepNum-1)+(theta_ip1-theta_i)./(deltam/h_i); theta_i=theta_ip1; else h_i=zdelta; deltam=interp1(m1,(z(ii,StepNum)-h_i)/zdelta+1)-interp1(m1,z(ii,StepNum)/zdelta+1); theta_ip1=-sqrt(theta_i^2+2*deltam); StepNum=StepNum+1; z(ii,StepNum)=z(ii,StepNum-1)-h_i; x(ii,StepNum)=x(ii,StepNum-1)+(theta_ip1-theta_i)./(deltam/(-h_i)); end end end steps(ii)=StepNum;%保存步数 end %% %*********************************轨迹展示**************************** plot(x(1,1:steps(1)-1)./1000,z(1,1:steps(1)-1),'-ro','LineWidth',1.5);%angle= 4.6 hold on plot(x(2,1:steps(2)-1)./1000,z(2,1:steps(2)-1),'-cd','LineWidth',1.5);% angle= 3.7 plot(x(3,1:steps(3)-1)./1000,z(3,1:steps(3)-1),'-bs','LineWidth',1.5);%angle=3.6 plot(x(4,1:steps(4)-1)./1000,z(4,1:steps(4)-1),'-gs','LineWidth',1.5);%angle=3.6 plot(x(5,1:steps(5)-1)./1000,z(5,1:steps(5)-1),'-kd','LineWidth',1.5);%angle=3.7 plot(x(6,1:steps(6)-1)./1000,z(6,1:steps(6)-1),'-mo','LineWidth',1.5);%angle=4.6 mrad axis([0 rmax_user/1000 0 zmax_user ]) set(gca,'fontsize',15,'fontname','Time New Roman') set(gca,'xtick',0:5:rmax_user/1000) L=legend(' 4.6 mrad',' 3.7 mrad',' 3.6 mrad',' -3.6 mrad',' -3.7 mrad','-4.6 mrad'); set(L,'fontsize',14,'fontname','Times New Roman','Fontweight','bold') xlabel('Range(km)','FontName','Times New Roman','fontsize',16) ylabel('Height"(m)','FontName','Times New Roman','fontsize',16)
3、
%% ray tracing routine %发射天线高度6m,发射处波导高度14m,50km处波导高度7m,100公里处波导高度14m。 close all clear tic %% 输入计算区域参数 tx_height = 6; % 天线高度m rmax_user1=40*1e3;%传输范围(m) rmax_user2=80*1e3; rmax_user3=125*1e3; % 传播最远距离km zmax_user = 100; % 传播最大高度m zdelta = 0.25;% 将射线的递进高度步长与剖面的步长一致 height = (0:floor(zmax_user*2/zdelta))*zdelta;%高度和折射剖面对应 %% 蒸发波导剖面,0km~50km 伪折射率方法,已知波导高度,求剖面 evap_const = 1.5e-4; % meter, evaporation duct constant, typical value 空气动力学粗糙度因子 M0 = 330;%不影响,求的是变化率 Mdelta = 10; % 波导高度 M1 = 0.125*(height - Mdelta*log(height./evap_const+1)); M1 = M0 + M1;% 大气修正折射率 m1 = M1*1e-6+1;% 大气修正折射指数 g1 = (m1(2:end)-m1(1:end-1))/zdelta; %梯度 %% 第二个折射率剖面50~100km剖面 evap_const = 1.5e-4; % evaporation duct constant, typical value 空气动力学粗糙度因子 M0 = 330; Mdelta = 14; % 波导高度 M2 = 0.125*(height - Mdelta*log(height./evap_const+1)); M2 = M0 + M2;% 大气修正折射率 m2 = M2*1e-6+1;% 大气修正折射指数 g2 = (m2(2:end)-m2(1:end-1))/zdelta; %% 第三个折射率剖面 evap_const = 1.5e-4; % evaporation duct constant, typical value 空气动力学粗糙度因子 M0 = 330; Mdelta = 20; % 波导高度 M3 = 0.125*(height - Mdelta*log(height./evap_const+1)); M3 = M0 + M3;% 大气修正折射率 m3 = M3*1e-6+1;% 大气修正折射指数 g3 = (m3(2:end)-m3(1:end-1))/zdelta; %% 角度参数,输入角度是与水平面之间的俯仰角,与grazing angle(入射余角)一样都是与水平面的夹角。内部计算的时候用的是与z轴顺时针的夹角 angle_start = 0.1; %anticlockwise 逆时针正 angle_end = -0.1; tetadel = 0.002; % Ntheta = floor(abs(angle_start-angle_end)/tetadel);% 步进数,角度个数 angle=[angle_start:-tetadel:angle_end]*pi/180;%度数 N_angle=length(angle);%角度个数 %% %****** 射线方程开始 *********** x = zeros(N_angle,10000);%距离 z = zeros(N_angle,10000);%高度 for ii = 1:N_angle StepNum = 1; x(ii,StepNum) = 0.0; % 第0步 z(ii,StepNum) = tx_height; theta_i = angle(ii);% 初始出射角,公式里面关于角度的计算都是用的是弧度 h_i = abs(zdelta); % 初始步长,取参考值 % 这里假设出射点不在折射率转折点附近 %% %******************************** while z(ii,StepNum) <= zmax_user && x(ii,StepNum)<= rmax_user1% 终止条件,超过计算空间范围 if theta_i > 0 && theta_i<= (3*pi/180) % = 0.0524 向上传播(怎么来的)3度 deltam = interp1(m1,(z(ii,StepNum)+zdelta)/zdelta+1)-interp1(m1,z(ii,StepNum)/zdelta+1); r = theta_i^2 + 2*deltam;%出射角 if r>= 0 % 向上传播 h_i = min(zdelta,(floor(z(ii,StepNum)/zdelta)+1)*zdelta-z(ii,StepNum)); h_i = max(h_i,0.1*zdelta); deltam = interp1(m1,(z(ii,StepNum)+h_i)/zdelta+1)-interp1(m1,z(ii,StepNum)/zdelta+1); theta_ip1 = sqrt(theta_i^2 +2*deltam); StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)+ h_i ;%加高度 x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1-theta_i)./(deltam/h_i); theta_i = theta_ip1; else % 向上传播的时候发生反转 theta_ip1 = 0;% 射线在未到达下一剖面层时已经发生反转,现在需要确定是在多大的高度发生反转 h_i = (theta_ip1^2-theta_i^2)./(2*g1(ceil(z(ii,StepNum)/zdelta))); %positive StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)+ h_i; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1 - theta_i)/ g1(ceil(z(ii,StepNum-1)/zdelta)) ; theta_i = theta_ip1; end elseif theta_i<0 && theta_i >= -(3*pi/180) %向下传播 deltam = interp1(m1,(z(ii,StepNum)-zdelta)./zdelta+1)-interp1(m1,z(ii,StepNum)./zdelta+1); r = theta_i.^2 + 2*deltam; if z(ii,StepNum) <= 0 % zdelta % 到达下界面反转 theta_i = -theta_i; StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1); x(ii,StepNum) = x(ii,StepNum-1); elseif r>=0 h_i = min(zdelta,z(ii,StepNum)-(ceil(z(ii,StepNum)/zdelta)-1)*zdelta); h_i = max(h_i,0.1*zdelta); deltam = interp1(m1,(z(ii,StepNum)-h_i)/zdelta+1)-interp1(m1,z(ii,StepNum)/zdelta+1); theta_ip1 = - sqrt(theta_i^2 +2*deltam); StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)- h_i ; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1-theta_i)./(deltam/(-h_i)); theta_i = theta_ip1; elseif r < 0 theta_ip1 = 0; % 射线在未到达下一剖面层时已经发生反转,现在需要确定是在多大的高度发生反转 h_i = (theta_ip1^2-theta_i^2)./(2*g1(ceil(z(ii,StepNum)/zdelta))); %positive StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)- h_i; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1 - theta_i)/ g1(ceil(z(ii,StepNum-1)/zdelta)) ; theta_i = theta_ip1; end elseif theta_i > (3*pi/180) % 0.0524 h_i = zdelta; m_i = interp1(m1,(z(ii,StepNum))/zdelta+1); m_ip1 = interp1(m1,(z(ii,StepNum)+h_i)/zdelta+1); re = 6370*1e3; % 地球半径 StepNum = StepNum + 1; z(ii,StepNum) = z(ii,StepNum-1) + h_i; x(ii,StepNum) = x(ii,StepNum-1) + re/(re+z(ii,StepNum))*m_i*cos(theta_i)*h_i/sqrt(m_ip1^2-(m_i*cos(theta_i))^2); theta_ip1 = atan(h_i/(x(ii,StepNum)-x(ii,StepNum-1))); theta_i = theta_ip1; elseif theta_i < -(3*pi/180) % 0.0524 if z(ii,StepNum)>= zdelta h_i = zdelta; m_i = interp1(m1,(z(ii,StepNum))/zdelta+1); m_ip1 = interp1(m1,(z(ii,StepNum)-h_i)/zdelta+1); re = 6370*1e3; % 地球半径 StepNum = StepNum + 1; z(ii,StepNum) = z(ii,StepNum-1) - h_i; x(ii,StepNum) = x(ii,StepNum-1) + re/(re+z(ii,StepNum))*m_i*cos(theta_i)*h_i/sqrt(m_ip1^2-(m_i*cos(theta_i))^2); theta_ip1 = -atan(h_i/(x(ii,StepNum)-x(ii,StepNum-1))); theta_i = theta_ip1; else theta_i = -theta_i; StepNum = StepNum+1; x(ii,StepNum) = x(ii,StepNum-1); z(ii,StepNum) = z(ii,StepNum-1); end elseif theta_i == 0 if g1(ceil(z(ii,StepNum)/zdelta))>0 % 继续向上走一步 h_i = zdelta; deltam = interp1(m1,(z(ii,StepNum)+h_i)/zdelta+1)-interp1(m1,z(ii,StepNum)/zdelta+1); theta_ip1 = sqrt(theta_i^2 +2*deltam); StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)+ h_i ; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1-theta_i)./(deltam/h_i); theta_i = theta_ip1; else % 继续向下走一步 h_i = zdelta; deltam = interp1(m1,(z(ii,StepNum)-h_i)/zdelta+1)-interp1(m1,z(ii,StepNum)/zdelta+1); theta_ip1 = -sqrt(theta_i^2 +2*deltam); StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)- h_i ; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1-theta_i)./(deltam/(-h_i)); theta_i = theta_ip1; end end end %% % *********************** while z(ii,StepNum)<=zmax_user && x(ii,StepNum)>rmax_user1 && x(ii,StepNum)<=rmax_user2 if theta_i > 0 && theta_i<= (3*pi/180) deltam = interp1(m2,(z(ii,StepNum)+zdelta)/zdelta+1)-interp1(m2,z(ii,StepNum)/zdelta+1); r = theta_i^2 + 2*deltam; if r>= 0 h_i = min(zdelta,(floor(z(ii,StepNum)/zdelta)+1)*zdelta-z(ii,StepNum)); h_i = max(h_i,0.1*zdelta); deltam = interp1(m2,(z(ii,StepNum)+h_i)/zdelta+1)-interp1(m2,z(ii,StepNum)/zdelta+1); theta_ip1 = sqrt(theta_i^2 +2*deltam); StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)+ h_i ; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1-theta_i)./(deltam/h_i); theta_i = theta_ip1; else % 向上传播的时候发生反转 theta_ip1 = 0;% 射线在未到达下一剖面层时已经发生反转,现在需要确定是在多大的高度发生反转 h_i = (theta_ip1^2-theta_i^2)./(2*g2(ceil(z(ii,StepNum)/zdelta+0.01))); %positive StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)+ h_i; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1 - theta_i)/ g2(ceil(z(ii,StepNum-1)/zdelta+0.01)) ; theta_i = theta_ip1; end elseif theta_i<0 && theta_i >= -(3*pi/180) %向下传播 deltam = interp1(m2,(z(ii,StepNum)-zdelta)./zdelta+1)-interp1(m2,z(ii,StepNum)./zdelta+1); r = theta_i.^2 + 2*deltam; if z(ii,StepNum) <= 0 %zdelta % 到达下界面反转 theta_i = -theta_i; StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1); x(ii,StepNum) = x(ii,StepNum-1); elseif r>=0 h_i = min(zdelta,z(ii,StepNum)-(ceil(z(ii,StepNum)/zdelta)-1)*zdelta); h_i = max(h_i,0.1*zdelta); deltam = interp1(m2,(z(ii,StepNum)-h_i)/zdelta+1)-interp1(m2,z(ii,StepNum)/zdelta+1); theta_ip1 = - sqrt(theta_i^2 +2*deltam); StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)- h_i ; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1-theta_i)./(deltam/(-h_i)); theta_i = theta_ip1; elseif r < 0 theta_ip1 = 0;% 射线在未到达下一剖面层时已经发生反转,现在需要确定是在多大的高度发生反转 h_i = (theta_ip1^2-theta_i^2)./(2*g2(ceil(z(ii,StepNum)/zdelta))); %positive StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)- h_i; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1 - theta_i)/ g2(ceil(z(ii,StepNum-1)/zdelta)) ; theta_i = theta_ip1; end elseif theta_i > (3*pi/180) % 0.0524 h_i = zdelta; m_i = interp1(m2,(z(ii,StepNum))/zdelta+1); m_ip1 = interp1(m2,(z(ii,StepNum)+h_i)/zdelta+1); re = 6370*1e3; % 地球半径 StepNum = StepNum + 1; z(ii,StepNum) = z(ii,StepNum-1) + h_i; x(ii,StepNum) = x(ii,StepNum-1) + re/(re+z(ii,StepNum))*m_i*cos(theta_i)*h_i/sqrt(m_ip1^2-(m_i*cos(theta_i))^2); theta_ip1 = atan(h_i/(x(ii,StepNum)-x(ii,StepNum-1))); theta_i = theta_ip1; elseif theta_i < -(3*pi/180) % 0.0524 if z(ii,StepNum)>= zdelta h_i = zdelta; m_i = interp1(m2,(z(ii,StepNum))/zdelta+1); m_ip1 = interp1(m2,(z(ii,StepNum)-h_i)/zdelta+1); re = 6370*1e3; % 地球半径 StepNum = StepNum + 1; z(ii,StepNum) = z(ii,StepNum-1) - h_i; x(ii,StepNum) = x(ii,StepNum-1) + re/(re+z(ii,StepNum))*m_i*cos(theta_i)*h_i/sqrt(m_ip1^2-(m_i*cos(theta_i))^2); theta_ip1 = -atan(h_i/(x(ii,StepNum)-x(ii,StepNum-1))); theta_i = theta_ip1; else theta_i = -theta_i; StepNum = StepNum+1; x(ii,StepNum) = x(ii,StepNum-1); z(ii,StepNum) = z(ii,StepNum-1); end elseif theta_i == 0 if g2(ceil(z(ii,StepNum)/zdelta))>0 % 继续向上走一步 h_i = zdelta; deltam = interp1(m2,(z(ii,StepNum)+h_i)/zdelta+1)-interp1(m2,z(ii,StepNum)/zdelta+1); theta_ip1 = sqrt(theta_i^2 +2*deltam); StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)+ h_i ; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1-theta_i)./(deltam/h_i); theta_i = theta_ip1; else % 继续向下走一步 h_i = zdelta; deltam = interp1(m2,(z(ii,StepNum)-h_i)/zdelta+1)-interp1(m2,z(ii,StepNum)/zdelta+1); theta_ip1 = -sqrt(theta_i^2 +2*deltam); StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)- h_i ; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1-theta_i)./(deltam/(-h_i)); theta_i = theta_ip1; end end end %********************************************** %% % *********************** while z(ii,StepNum)<=zmax_user && x(ii,StepNum)>rmax_user2 && x(ii,StepNum)<=rmax_user3 if (theta_i > 0 && theta_i<= (3*pi/180)) deltam = interp1(m3,(z(ii,StepNum)+zdelta)/zdelta+1)-interp1(m3,z(ii,StepNum)/zdelta+1); r = theta_i^2 + 2*deltam; if r>= 0 h_i = min(zdelta,(floor(z(ii,StepNum)/zdelta)+1)*zdelta-z(ii,StepNum)); h_i = max(h_i,0.1*zdelta); deltam = interp1(m3,(z(ii,StepNum)+h_i)/zdelta+1)-interp1(m3,z(ii,StepNum)/zdelta+1); theta_ip1 = sqrt(theta_i^2 +2*deltam); StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)+ h_i ; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1-theta_i)./(deltam/h_i); theta_i = theta_ip1; else % 向上传播的时候发生反转 theta_ip1 = 0;% 射线在未到达下一剖面层时已经发生反转,现在需要确定是在多大的高度发生反转 h_i = (theta_ip1^2-theta_i^2)./(2*g2(ceil(z(ii,StepNum)/zdelta+0.01))); %positive StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)+ h_i; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1 - theta_i)/ g2(ceil(z(ii,StepNum-1)/zdelta+0.01)) ; theta_i = theta_ip1; end elseif (theta_i<0 && theta_i >= -(3*pi/180)) %向下传播 deltam = interp1(m3,(z(ii,StepNum)-zdelta)./zdelta+1)-interp1(m3,z(ii,StepNum)./zdelta+1); r = theta_i.^2 + 2*deltam; if z(ii,StepNum) <= 0 %zdelta % 到达下界面反转 theta_i = -theta_i; StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1); x(ii,StepNum) = x(ii,StepNum-1); elseif r>=0 h_i = min(zdelta,z(ii,StepNum)-(ceil(z(ii,StepNum)/zdelta)-1)*zdelta); h_i = max(h_i,0.1*zdelta); deltam = interp1(m3,(z(ii,StepNum)-h_i)/zdelta+1)-interp1(m3,z(ii,StepNum)/zdelta+1); theta_ip1 = - sqrt(theta_i^2 +2*deltam); StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)- h_i ; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1-theta_i)./(deltam/(-h_i)); theta_i = theta_ip1; elseif r < 0 theta_ip1 = 0;% 射线在未到达下一剖面层时已经发生反转,现在需要确定是在多大的高度发生反转 h_i = (theta_ip1^2-theta_i^2)./(2*g2(ceil(z(ii,StepNum)/zdelta))); %positive StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)- h_i; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1 - theta_i)/ g2(ceil(z(ii,StepNum-1)/zdelta)) ; theta_i = theta_ip1; end elseif theta_i > (3*pi/180) % 0.0524 h_i = zdelta; m_i = interp1(m3,(z(ii,StepNum))/zdelta+1); m_ip1 = interp1(m3,(z(ii,StepNum)+h_i)/zdelta+1); re = 6370*1e3; % 地球半径 StepNum = StepNum + 1; z(ii,StepNum) = z(ii,StepNum-1) + h_i; x(ii,StepNum) = x(ii,StepNum-1) + re/(re+z(ii,StepNum))*m_i*cos(theta_i)*h_i/sqrt(m_ip1^2-(m_i*cos(theta_i))^2); theta_ip1 = atan(h_i/(x(ii,StepNum)-x(ii,StepNum-1))); theta_i = theta_ip1; elseif theta_i < -(3*pi/180) % 0.0524 if z(ii,StepNum)>= zdelta h_i = zdelta; m_i = interp1(m3,(z(ii,StepNum))/zdelta+1); m_ip1 = interp1(m3,(z(ii,StepNum)-h_i)/zdelta+1); re = 6370*1e3; % 地球半径 StepNum = StepNum + 1; z(ii,StepNum) = z(ii,StepNum-1) - h_i; x(ii,StepNum) = x(ii,StepNum-1) + re/(re+z(ii,StepNum))*m_i*cos(theta_i)*h_i/sqrt(m_ip1^2-(m_i*cos(theta_i))^2); theta_ip1 = -atan(h_i/(x(ii,StepNum)-x(ii,StepNum-1))); theta_i = theta_ip1; else theta_i = -theta_i; StepNum = StepNum+1; x(ii,StepNum) = x(ii,StepNum-1); z(ii,StepNum) = z(ii,StepNum-1); end elseif theta_i == 0 if g2(ceil(z(ii,StepNum)/zdelta))>0 % 继续向上走一步 h_i = zdelta; deltam = interp1(m3,(z(ii,StepNum)+h_i)/zdelta+1)-interp1(m3,z(ii,StepNum)/zdelta+1); theta_ip1 = sqrt(theta_i^2 +2*deltam); StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)+ h_i ; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1-theta_i)./(deltam/h_i); theta_i = theta_ip1; else % 继续向下走一步 h_i = zdelta; deltam = interp1(m3,(z(ii,StepNum)-h_i)/zdelta+1)-interp1(m3,z(ii,StepNum)/zdelta+1); theta_ip1 = -sqrt(theta_i^2 +2*deltam); StepNum = StepNum+1; z(ii,StepNum) = z(ii,StepNum-1)- h_i ; x(ii,StepNum) = x(ii,StepNum-1) + (theta_ip1-theta_i)./(deltam/(-h_i)); theta_i = theta_ip1; end else break end end %********************************************** %% xx = zeros(1,StepNum-1); zz = zeros(1,StepNum-1); for jj=1:StepNum-1 xx(jj) = x(ii,jj)./1e3; zz(jj) = z(ii,jj); end plot(xx,zz,'-r') % axis([0 337 0 100]) hold on end set(gca,'fontsize',18, 'FontName','Times') set(gca,'xtick',0:20:200) set(gca,'ytick',0:10:100) axis([0 125 0 50]) xlabel('Range (km)') ylabel('Height (m)') toc % dis = 0:200; % edh = [ones(101,1)*16; ones(100,1)*30]; % plot(dis,edh,'--b','linewidth',3) % print(gcf,'-r600','-djpeg','Ray-tracing_EDH=30,16_AT=20.jpg')
三、PE模型
1、主函数
%PE_model函数 % SSPE_function(freq, ... % thetabw, thetae, polrz, tx_height, range, zmax_user, edge_range, edge_height, ... % duct_type, duct_M, duct_height, duct_range, terrain_type, interp_type, backward,... % ground_type, epsilon, sigma, delx, delz) clear all %% *********************************************************** %初始赋值 %*************************************************** %d大气波导数据 duct_M_array1=[300,385.400000000000,0,0,0; 350,300,350,0,0; 340,356,340,358,0; 300,330,310,350,0; 300,0,0,0,0; 450,250,300,0,0];%大气波导折射率值 duct_height_array1=[0,300,0,0,0; 0,200,300,0,0; 0,135,150,300,0; 0,50,150,300,0; 20,0,0,0,0; 0,100,250,0,0];%大气波导高度值 duct_type_array1=3;%大气波导类型选择 duct_M = duct_M_array1(duct_type_array1,:); %折射率值 duct_height = duct_height_array1(duct_type_array1,:);%高度值 duct_type = duct_type_array1; duct_range = 0; freq =3000;%频率 thetabw=2;%三分贝波束宽度 thetae=0;% polrz=1;%极化方式 tx_height=10;%天线高度 range=100;%水平距离 zmax_user=300;%竖直高度 edge_range=[5.13636000000000,7.31818000000000,9.50000000000000,11.1364000000000,... 12.9545000000000,15.2273000000000,16.9545000000000,19.4091000000000,21.4091000000000,... 23.6818000000000,25.2273000000000,26.8636000000000,28.9545000000000,31.1364000000000,... 33.7727000000000,36.0455000000000,37.9545000000000,39.5909000000000,41.1364000000000,... 42.7727000000000,44.5000000000000,46.5909000000000,48.4091000000000,49.3182000000000,... 49.6818000000000];%地形x轴(km) edge_height=[2.88461500000000,10.5769200000000,18.2692300000000,33.6538500000000,... 41.3461500000000,50.9615400000000,60.5769200000000,74.0384600000000,97.1153800000000,... 106.730800000000,75.9615400000000,50.9615400000000,29.8076900000000,22.1153800000000,... 14.4230800000000,14.4230800000000,20.1923100000000,45.1923100000000,89.4230800000000,... 108.653800000000,135.576900000000,129.807700000000,139.423100000000,81.7307700000000,... 14.4230800000000];%地形y轴(m) terrain_type=2;%地形类型 interp_type=3;%插值方式 backward=1;%后项模式 ground_type=1;%表面模式 epsilon=15;%介电常数 sigma=0.2293;%电磁导电率 delz=0.2900;%高度步进 delx=200;%水平距离步进(m) %% %函数主体 tic % warning_gui;%调用提醒函数 % pause(.1) % f1 = findobj('tag','figurew'); stopflag = 0; path_loss = [];%预留空间 prop_fact = []; free_space_loss = []; %% ************************************************************************ % CONSTANTS and DEFAULT INPUT PARAMETERS(常数) % ************************************************************************* inputs.c = 3*1e8; % meter/second, velocity of light in free space光速 inputs.eps0 = (1e-9)/(36*pi); % F/m, permittivity 介电常数 inputs.mu0 = 4*pi*1e-7; % H/m, permeability of free space自由场介电常数 inputs.nu0 = 120*pi; % ohm, intrinsic impedance of the free space自由场固有阻抗 inputs.ae = 6378*1e3; % meter, earth radius 地球半径 inputs.sgrad = 0.118; % 1/meter, standart atmosphere gradient of modified refractivity标准大气修正折射率梯度 inputs.evap_const = 1.5e-4; % meter, evaporation duct constant, typical value波导常数类型值 inputs.N = 1024*1; % FFT size, # of points btw [0, zmax_user] fft点数 inputs.maxangle = 15; % degree, max allowable angle 最大允许角度 inputs.wind_ratio = 1/2; % the ratio of height extension (wrt N) to apply window function窗函数高度延伸比 % N used in the program = N*(wind_ratio+1) minf = 2001;%频率 if (freq <= 401) inputs.wind_ratio = 8; % N used in the program = N*wind_ratio elseif (freq <= 1001) inputs.wind_ratio = 4; % N used in the program = N*wind_ratio elseif (freq < minf) inputs.wind_ratio = 2; % N used in the program = N*wind_ratio end inputs.win_type = 1; % win_type=1 ==> hanning window窗函数选择 % win_type=2 ==> hamming window inputs.eps = 1e-10; inputs.threshold = 0.025; % Threshold value of difference matrix in backward propagation反向传播中差分矩阵的阈值 inputs.max_iter = 1000; % Max. number of iterations 迭代次数 %% ************************************************************************ % INPUT PARAMETERS 输入参数 % ************************************************************************* inputs.freq = freq*1e6; % Hz, frequency (conversion from MHz to Hz)频率 inputs.range = range*1e3; % m, output range (conversion from km to m) 水平距离 inputs.zmax_user = zmax_user; % meter, maximum height (max desired calculation height)高度 inputs.polrz = polrz; % polrz=1 ==> odd, horizontal polarization极化方式 % polrz=2 ==> even, vertical polarization inputs.thetabw = thetabw*pi/180;%三分贝波束宽度(单位rad) inputs.thetae = thetae*pi/180;%仰角度数(单位rad) inputs.tx_height = tx_height;%天线高度 %% ************************************************************************ % PARAMETERS CALCULATED FROM INPUT PARAMETERS 由输入值得到的常数 % delx, delz, zmax, z, p, pmax, etc. %************************************************************************** lamda = inputs.c/inputs.freq; % meter, wavelength 波长 beta = 2*pi/lamda; % 1/meter, wavenumber 理论物理中定义为:k=2π/λ,波数 zmax = inputs.zmax_user*(1+inputs.wind_ratio); %????什么高度 %if (ground_type == 2) && ((inputs.freq*1e-6) < 400) % mixed if ((inputs.freq*1e-6) < minf) % mixed zmax = inputs.zmax_user*inputs.wind_ratio; end %%delz = lamda / (2*sin(inputs.maxangle*pi/180)); %delz = inputs.zmax_user/(inputs.N-1); %if ground_type == 2 % mixed % delz = 19.53125; %end z = 0:delz:zmax;%高度方向递进 inputs.N = length(z); pmax = (pi/delz);%角度个数? p = linspace(0, pmax, inputs.N);%??? %delx = 2*beta*delz^2; %delx = min(delx, 1000); %delx = max(delx, 30); if inputs.tx_height == 0 %天线高度 inputs.tx_height = inputs.tx_height + delz; end; tx_range_index = 1; inputs.delx = delx;%水平步进 inputs.delz = delz;%高度步进 inputs.z = z;%??? inputs.beta = beta;%波数 inputs.delp = pi/zmax;%0.007?? %% ************************************************************************ % OUTPUT domain points %************************************************************************** range_vec = [0:delx:inputs.range]; % output.x,画图x轴 [val, savez] = min(abs(z-inputs.zmax_user));%值与位置,0.14,1035 z_user = z(1:savez); % output.z ,画图y轴 zae = z / inputs.ae;%高度除以半径 Nx = length(range_vec);%水平方向个数 N = inputs.N;%高度方向个数 %% ************************************************************************ % DEFINITION OF TERRAIN PROFILE 地形参数定义 %************************************************************************** if ~isempty(edge_range) temp(:,1) = edge_range.'; temp(:,2) = edge_height.'; ind = find(edge_range <= range+1e-8);%单位不统一???km ind2 = find(edge_height <= zmax+1e-8);%(m) ind = cat(2, ind, ind2);%连接两个矩阵 ind = unique(sort(ind));%排序 temp = temp(ind,:);%???? temp = sortrows(temp);%按列排序 edge_range = temp(:,1).'; edge_height = temp(:,2).';%大小重新排列 clear temp; num_of_edges = length(edge_range); edge_range = edge_range.*1e3;%变换单位 edge_range_index = floor(edge_range/delx) + 1; edge_height_index = floor(edge_height/delz) + 1; %地形下标数 ind = find(diff(edge_range_index) == 0);%将相邻差值为0的数找出来 edge_range_index(ind) = []; edge_height_index(ind) = [];%去掉相同的值 if (terrain_type == 2) && (num_of_edges > 1) && (interp_type > 1) % terrain case地形模式、插值模式 if interp_type == 2 % linear线性插值 aa = fit(edge_range_index.', edge_height_index.', 'linear'); elseif interp_type == 3 % cubic spline aa = fit(edge_range_index.', edge_height_index.', 'cubicspline'); %拟合曲线 end edge_range_index = edge_range_index(1):edge_range_index(end);%按顺序排列 edge_height_index = aa(edge_range_index); edge_height_index = floor(edge_height_index) + 1;%扩充到224个数据,怎么用??? end end %% ************************************************************************ % WINDOW FUNCTION CALCULATION 窗函数计算 %************************************************************************** %if (ground_type == 2) && ((inputs.freq*1e-6) < 400) % mixed if ((inputs.freq*1e-6) < minf) % mixed wind = fun_window2(inputs); % WINDOW FUNCTION DEFINITION else wind = fun_window(inputs); % WINDOW FUNCTION DEFINITION end wind = wind.'; %% ************************************************************************ % DEFINITION OF REFRACTIVITY PROFILE 折射剖面 %************************************************************************** temp(:,1) = duct_range'; temp(:,2) = duct_type'; num = size(duct_height,2); temp(:,3:3+num-1) = duct_height; temp(:,3+num:3+2*num-1) = duct_M; ind = find(duct_range <= range+1e-8); temp = temp(ind,:); temp = sortrows(temp); duct_range = temp(:,1)'; duct_type = temp(:,2)'; duct_height = temp(:,3:3+num-1); duct_M = temp(:,3+num:3+2*num-1); clear temp; np = length(duct_type); duct_range = duct_range.*1e3; %duct_range_index = [1 round(duct_range/delx) Nx]; duct_range_index = floor(duct_range/delx)+1; nprofile = zeros(np, N); mprofile = zeros(np, N); for i = 1:np zero_vec = find(duct_height(i,:) == 0); if length(zero_vec) == 1 duct_heighti = duct_height(i,:); duct_Mi = duct_M(i,:); else duct_heighti = duct_height(i, 1:zero_vec(2)-1); duct_Mi = duct_M(i, 1:zero_vec(2)-1); end; [nprofile(i,:), mprofile(i,:)] = fun_refrac(duct_Mi, duct_type(i), duct_heighti, inputs); end if np == 1 %if ground_type == 1 % PEC % mn = nprofile.'; %else mn = mprofile.'; %end else % Interpolate (Linear) mn = zeros(N, Nx); for i = 1:np-1 indi = duct_range_index(i); indi1 = duct_range_index(i+1); for j = 1:N %if ground_type == 1 % PEC % mn(j,indi:indi1) = linspace(nprofile(i,j), nprofile(i+1,j), indi1-indi+1); %else mn(j,indi:indi1) = linspace(mprofile(i,j), mprofile(i+1,j), indi1-indi+1); %end end; end; wind = repmat(wind,1,Nx); end %clear nprofile; clear mprofile; %% ************************************************************************ % WINDOWED ENVIRONMENT TERM 窗口环境 %************************************************************************** %if ground_type == 1 % PEC %wide-angle prop = exp(1j*beta*delx*(sqrt(1-p.^2./beta^2)-1)); envpr = exp(1j * beta * delx * (mn-1)); % ENVIRONMENT TERM %prop = exp(-1j * p.^2 .* delx / (2*beta)); %envpr = exp(1j*0.5 * beta * delx * mn); % ENVIRONMENT TERM prop = prop.' .* wind(:,1); envpr = envpr .* wind; %else %envpr = exp(1j*0.5 * beta * delx * mn); % ENVIRONMENT TERM % envpr = exp(1j*beta * delx * mn); % ENVIRONMENT TERM % envpr = envpr .* wind; %end; %clear mn; clear wind; %% ************************************************************************ % INITIAL FIELD DEFINITION 初始场定义 %************************************************************************** u0z = fun_initial_field(inputs); % INITIAL FIELD DEFINITION u0z = u0z.'; %% ************************************************************************ % CALCULATE PARAMETERS if GROUND is MIXED TYPE 计算地面混合类型时的参数 %************************************************************************** % Reflection Coefficient ref_coef = -1; alg = 0; % if alg==0, PEC solution % if alg==1, central difference DMFT % if alg==2, backward difference DMFT if ground_type == 2 % mixed %u0z = u0z./max(abs(u0z)); % DMFT Calculations %if (inputs.polrz==2) %vertical polrz alg = 2; %backward difference DMFT if (inputs.freq*1e-6) < 400 alg = 1; %central difference DMFT end %end er = epsilon + 1i*sigma*60*lamda; if inputs.polrz == 2 %vertical polarization alfa = 1i*beta*sqrt(er-1)./er; else %horizontal polarization alfa = 1i*beta*sqrt(er-1); end if alg == 1 %central difference DMFT if inputs.polrz == 2 %vertical polarization r = sqrt(1+(alfa*delz).^2)-alfa*delz; % eq.20 else %horizontal polarization r = -sqrt(1+(alfa*delz).^2)-alfa*delz; % eq.21 end else %backward difference DMFT r = 1 / (1+(alfa*delz)); end NN = inputs.N-1; rn = r.^(0:NN).'; ad = (log(r)/delz).^2; if alg == 1 %central difference DMFT dmft.rk = 2*(1-r^2) / (1+r^2) / (1-r^(2*N)); dmft.c1x = exp(1j*delx*sqrt(beta^2+ad)); ad = ((log(r)-1j*pi)/delz).^2; dmft.c2x = exp(1j*delx*sqrt(beta^2+ad)); else dmft.c3x = exp(1j*delx*sqrt(beta^2+ad)); end %********** if alg == 1 %central difference DMFT dmft.ck1 = 0.5*(u0z(1) + u0z(end)*rn(end)); dmft.ck2 = 0.5*(u0z(1)*rn(end) + u0z(end)); dmft.ck1 = dmft.ck1 + sum(u0z(2:end-1) .* rn(2:end-1)); uzf = fliplr(u0z); dmft.ck2 = dmft.ck2 + sum(uzf(2:N-1) .* rn(2:N-1) .* (-1).^(1:N-2)'); dmft.ck1 = dmft.ck1*dmft.rk; dmft.ck2 = dmft.ck2*dmft.rk; else %backward difference DMFT dmft.ck3 = sum(u0z(1:end-1) .* rn(1:end-1)); end dmft.alfa = alfa; dmft.r = r; dmft.rn = rn; clear rn; %********** % END OF DMFT Calculations % Reflection Coefficient for TERRAIN if inputs.polrz == 2 %vertical polarization nu2 = inputs.nu0 * sqrt(er-1) ./ er; else %horizontal polarization nu2 = inputs.nu0 * sqrt(er-1); end nu1 = inputs.nu0; ref_coef = (nu2-nu1) ./ (nu2+nu1); end %% ************************************************************************ % ITERATION STARTS (FIELD CALCULATION)迭代计算(场计算) %************************************************************************** u_matrix = zeros(inputs.N, Nx); % total matrix u_matrix(:,tx_range_index) = u0z; if terrain_type == 1, % NO TERRAIN CASE (FLAT) uz = u0z; % initial field for i = tx_range_index:Nx-1, if np == 1 % range-indep. refrac. envpri = envpr; else envpri = envpr(:,i); end if ground_type == 1 % PEC uz = fun_field_pec(uz, inputs, prop, envpri); else % mixed if (er == 0) || (alg == 0) uz = fun_field_pec(uz, inputs, prop, envpri); else [uz, dmft] = fun_field_mixed(uz, inputs, prop, envpri, alg, dmft); end end u_matrix(:, i+1) = abs(uz); set (findobj(f1,'tag','text2'),'String',' '); pause(.00001) if isempty(findobj('tag','figurew')) stopflag = 1; break; end; end if stopflag == 1 return; end elseif terrain_type == 2, % terrain case 1WAY if backward == 1 % 1-way SSPE uz = u0z; % initial field edgeh = sparse(1, Nx); edgeh(edge_range_index) = edge_height_index; for i = tx_range_index:Nx-1, if np == 1 % range-indep. refrac. envpri = envpr; else envpri = envpr(:,i); end if ground_type == 1 % PEC uz = fun_field_pec(uz, inputs, prop, envpri); else if (er == 0) || (alg == 0) uz = fun_field_pec(uz, inputs, prop, envpri); else [uz, dmft] = fun_field_mixed(uz, inputs, prop, envpri, alg, dmft); end end uz(1:full(edgeh(i+1)),1) = 0; u_matrix(:, i+1) = uz; % set (findobj(f1,'tag','text2'),'String',' '); % pause(.00001) % % if isempty(findobj('tag','figurew')) % stopflag = 1; % break; % end; end if stopflag == 1 return; end else % backward == 2 2-way SSPE RECURSIVE uz = u0z; % initial field %u_matrix(:,tx_range_index) = uz; iter_user = 1; dif = 10^10; %ref_coef = -1.0; threshold = inputs.threshold; max_iter = inputs.max_iter; output.bounce = sparse(1, Nx); %% isedge = sparse(1, Nx); isedge(edge_range_index) = 1; edger = sparse(1, Nx); edger(edge_range_index) = edge_range_index; edgeh = sparse(1, Nx); edgeh(edge_range_index) = edge_height_index; ind = find(edgeh == 1); edgeh(ind) = 0; isedge(ind) = 0; edger(ind) = 0; if tx_range_index == 1 move.now = tx_range_index; way.now = 'F'; %FORWARD matrix.now(:,1) = uz; else move.now = [tx_range_index tx_range_index]; way.now(1) = 'F'; way.now(2) = 'B'; % BACKWARD matrix.now(:,1) = uz; matrix.now(:,2) = uz; end matrix.next = []; move.next = []; iter = 1; while dif > threshold, if isempty(findobj('tag','figurew')) stopflag = 1; return end; set (findobj(f1,'tag','text2'),'String',cat(2,'Threshold = ',num2str(dif), ' > ', num2str(threshold))); pause(.0001) if (iter ~= 1) move.now = move.next; move.next = []; way.now = way.next; way.next = 'T'; matrix.now = matrix.next; matrix.next = []; ind = find(way.now ~= 'T'); % Not TERMINATION move.now = move.now(ind); way.now = way.now(ind); matrix.now = matrix.now(:,ind); end len = length(move.now); ic = 1; for i = 1:len, uz = matrix.now(:,i); %uznew = fun_field_pec(uz, inputs, prop, envpr); if np == 1 % range-indep. refrac. envpri = envpr; else envpri = envpr(:,move.now(i)); end if ground_type == 1 % PEC uznew = fun_field_pec(uz, inputs, prop, envpri); else if (er == 0) || (alg == 0) uznew = fun_field_pec(uz, inputs, prop, envpri); else [uznew, dmft] = fun_field_mixed(uz, inputs, prop, envpri, alg, dmft); end end %% FORWARD PROPAGATION if (way.now(i) == 'F') && (move.now(i) ~= Nx) ind = find(move.next == move.now(i)+1, 1); add_flag = 1; if ~isempty(ind) if (way.next(ind) == 'F') matrix.next(:,ind) = matrix.next(:,ind) + uznew; if isedge(move.now(i)+1) % meet edge matrix.next(1:edgeh(move.now(i)+1), ind) = 0; end add_flag = 0; end end if add_flag move.next(ic) = move.now(i)+1; matrix.next(:,ic) = uznew; if isedge(move.now(i)+1) % meet edge matrix.next(1:edgeh(move.now(i)+1), ic) = 0; end way.next(ic) = 'F'; if ((move.now(i)+1) == Nx) way.next(ic) = 'T'; end ic = ic+1; end if isedge(move.now(i)+1) % meet edge if isedge(move.now(i)) if (edgeh(move.now(i)) < edgeh(move.now(i)+1)) uznew = uznew*ref_coef *exp(1j*beta*2*(edger(move.now(i)+1)-1)*delx); uznew(edgeh(move.now(i)+1)+1:end) = 0; %uznew = fun_field_pec(uznew, inputs, prop, envpr); if np == 1 % range-indep. refrac. envpri = envpr; else envpri = envpr(:,move.now(i)); end if ground_type == 1 % PEC uznew = fun_field_pec(uznew, inputs, prop, envpri); else if (er == 0) || (alg == 0) uznew = fun_field_pec(uznew, inputs, prop, envpri); else [uznew, dmft] = fun_field_mixed(uznew, inputs, prop, envpri, alg, dmft); end end ind = find(move.next == move.now(i), 1); add_flag = 1; if ~isempty(ind) if (way.next(ind) == 'B') matrix.next(:,ind) = matrix.next(:,ind) + uznew; matrix.next(1:edgeh(move.now(i)), ind) = 0; add_flag = 0; end end if add_flag move.next(ic) = move.now(i); matrix.next(:,ic) = uznew; matrix.next(1:edgeh(move.now(i)), ic) = 0; way.next(ic) = 'B'; ic = ic+1; end end else uznew = uznew*ref_coef *exp(1j*beta*2*(edger(move.now(i)+1)-1)*delx); uznew(edgeh(move.now(i)+1)+1:end) = 0; %uznew = fun_field_pec(uznew, inputs, prop, envpr); if np == 1 % range-indep. refrac. envpri = envpr; else envpri = envpr(:,move.now(i)); end if ground_type == 1 % PEC uznew = fun_field_pec(uznew, inputs, prop, envpri); else if (er == 0) || (alg == 0) uznew = fun_field_pec(uznew, inputs, prop, envpri); else [uznew, dmft] = fun_field_mixed(uznew, inputs, prop, envpri, alg, dmft); end end ind = find(move.next == move.now(i), 1); add_flag = 1; if ~isempty(ind) if (way.next(ind) == 'B') matrix.next(:,ind) = matrix.next(:,ind) + uznew; add_flag = 0; end end if add_flag move.next(ic) = move.now(i); matrix.next(:,ic) = uznew; way.next(ic) = 'B'; ic = ic+1; end end output.bounce(move.now(i)+1) = output.bounce(move.now(i)+1) + 1; end %% BACKWARD PROPAGATION elseif (way.now(i) == 'B') && (move.now(i) ~= 1) ind = find(move.next == move.now(i)-1, 1); add_flag = 1; if ~isempty(ind) if (way.next(ind) == 'B') matrix.next(:,ind) = matrix.next(:,ind) + uznew; if isedge(move.now(i)-1) % meet edge matrix.next(1:edgeh(move.now(i)-1), ind) = 0; end add_flag = 0; end end if add_flag move.next(ic) = move.now(i)-1; matrix.next(:,ic) = uznew; if isedge(move.now(i)-1) % meet edge matrix.next(1:edgeh(move.now(i)-1), ic) = 0; end way.next(ic) = 'B'; if (move.now(i)-1) == 1 way.next(ic) = 'T'; end ic = ic+1; end if isedge(move.now(i)-1) % meet edge if isedge(move.now(i)) if (edgeh(move.now(i)) < edgeh(move.now(i)-1)) uznew = uznew*ref_coef *exp(1j*beta*2*(edger(move.now(i)+1)-1)*delx); uznew(edgeh(move.now(i)-1)+1:end) = 0; %uznew = fun_field_pec(uznew, inputs, prop, envpr); if np == 1 % range-indep. refrac. envpri = envpr; else envpri = envpr(:,move.now(i)); end if ground_type == 1 % PEC uznew = fun_field_pec(uznew, inputs, prop, envpri); else if (er == 0) || (alg == 0) uznew = fun_field_pec(uznew, inputs, prop, envpri); else [uznew, dmft] = fun_field_mixed(uznew, inputs, prop, envpri, alg, dmft); end end ind = find(move.next == move.now(i), 1); add_flag = 1; if ~isempty(ind) if (way.next(ind) == 'F') matrix.next(:,ind) = matrix.next(:,ind) + uznew; matrix.next(1:edgeh(move.now(i)), ind) = 0; add_flag = 0; end end if add_flag move.next(ic) = move.now(i); matrix.next(:,ic) = uznew; matrix.next(1:edgeh(move.now(i)), ic) = 0; way.next(ic) = 'F'; ic = ic+1; end end else uznew = uznew*ref_coef *exp(1j*beta*2*(edger(move.now(i)+1)-1)*delx); uznew(edgeh(move.now(i)-1)+1:end) = 0; %uznew = fun_field_pec(uznew, inputs, prop, envpr); if np == 1 % range-indep. refrac. envpri = envpr; else envpri = envpr(:,move.now(i)); end if ground_type == 1 % PEC uznew = fun_field_pec(uznew, inputs, prop, envpri); else if (er == 0) || (alg == 0) uznew = fun_field_pec(uznew, inputs, prop, envpri); else [uznew, dmft] = fun_field_mixed(uznew, inputs, prop, envpri, alg, dmft); end end ind = find(move.next == move.now(i), 1); add_flag = 1; if ~isempty(ind) if (way.next(ind) == 'F') matrix.next(:,ind) = matrix.next(:,ind) + uznew; add_flag = 0; end end if add_flag move.next(ic) = move.now(i); matrix.next(:,ic) = uznew; way.next(ic) = 'F'; ic = ic+1; end end output.bounce(move.now(i)-1) = output.bounce(move.now(i)-1) + 1; end else disp('Error'); end end utemp = zeros(inputs.N, Nx); lenn = length(move.next); dif = 0; for i = 1:lenn, utemp(:,move.next(i)) = utemp(:,move.next(i)) + matrix.next(:,i); u_matrix(:,move.next(i)) = u_matrix(:,move.next(i)) + matrix.next(:,i); dif = dif + norm(matrix.next(:,i)); end %dif = norm(utemp); output.convergence(iter) = dif; iter = iter+1; %if mod(iter, 1)==0, % disp(way.next); % disp(num2str(move.next)); % disp(['Difference (< ' sprintf('%.2e', threshold) ') = ' sprintf('%.5f', dif)]); % disp(['Iteration = ' sprintf('%d', iter)]); % disp(' '); %disp(output.bounce(find(output.bounce~=0))); % disp('*****************************************************'); %end if (dif < threshold) && (iter < Nx) while dif < threshold threshold = threshold*(0.9); end end if (iter > iter_user*max_iter), iter_user = iter_user+1; %disp('Max number of iterations is exceeded..'); %disp('Do you want to continue for iteration?'); %answer = input('If yes, press enter; otherwise press ''0'' ==> '); %if answer == 0, % disp('The iteration is stopped by the user without satisfying the predefined threshold!'); % break; %end; end; end end; end; range_matrix = repmat(range_vec, inputs.N, 1); %u_matrix = u_matrix.*exp(-1j*beta*range_matrix); u_matrix = abs(u_matrix); u_matrix = u_matrix(1:savez,:); log10lamda = log10(lamda); log10umatrix = log10(u_matrix(1:savez,:)); %log10umatrix2 = log10(u_matrix(1:savez,:) ./ sqrt(range_matrix(1:savez,:))); %% ************************************************************************ % CALCULATE PROPAGATON FACTOR %************************************************************************** prop_fact = 20*log10umatrix + 10*log10(range_matrix(1:savez,:)) + 10*log10lamda; %dB %% ************************************************************************ % CALCULATE PATH LOSS %************************************************************************* path_loss = -20*log10umatrix + 20*log10(4*pi)+ ... 20*log10(range_matrix(1:savez,:)) - 20*log10lamda; %dB %% ************************************************************************ % CALCULATE FREE-SPACE LOSS %************************************************************************** free_space_loss = 32.45+20*log10(inputs.freq*1e-6)+20*log10(range_vec*1e-3); %dB toc %% %画图表示 plot_array = prop_fact; plot_flag = 2; figure axes; new_x = range_vec./1e3;%水平距离 new_y = z_user; %高度 xmin = new_x(1);%画图范围 xmax = new_x(length(new_x)); ymin = z_user(1); ymax = z_user(length(z_user)); hhimage = imagesc(new_x, new_y, plot_array); %画图 % shading interp; %阴影 % view(0,90);%视角选择 % set(gcf,'WindowButtonMotionFcn','guimousevalue'); % imageflag = 1; set(gca,'Ydir','normal');%y轴设置 hx=xlabel('Range (km)','Linewidth',2,'fontsize',9);%坐标标注 hy=ylabel('Altitude (m)','Linewidth',2,'fontsize',9); % set(hx,'pos',get(hx,'pos') + [0 10.0 0]); %set(hy,'pos',get(hy,'pos') + [0.3 0.0 0]); axis tight; colormap jet cb = colorbar('eastoutside');%东边外部 %%
2、窗函数
function [wind] = fun_window(inputs) N = inputs.N; wind_ratio = inputs.wind_ratio/2; wind = ones(1,N); if inputs.win_type == 1, ha = hanning(2*round(N*wind_ratio)+1).'; ha = ha(round(N*wind_ratio)+1:end); elseif inputs.win_type == 2, ha = hamming(2*round(N*wind_ratio)+1).'; ha = ha(round(N*wind_ratio)+1:end); end; wind(N-length(ha)+1:N) = ha;
3、窗函数
% WINDOWING FUNCTION窗函数 function [wind] = fun_window2(inputs) N = round(inputs.N/2); wind = ones(1,inputs.N); if inputs.win_type == 1, ha = hanning(2*N+1).'; ha = ha(N+1:end); elseif inputs.win_type == 2, ha = hamming(2*N+1).'; ha = ha(N+1:end); end; wind(inputs.N-length(ha)+1:inputs.N) = ha;
4、折射剖面函数
%% ************************************************************************ % REFRACTIVITY PROFILE FUNCTION function [n, m] = fun_refrac(duct_M, duct_type, duct_height, inputs) M(1) = duct_M(1); % standard atmosphere if duct_type == 1, term = inputs.delz*inputs.sgrad; M = M(1)*ones(1,inputs.N); M = M + (0:term:(inputs.N-1)*term); % evaporation duct elseif duct_type == 5, term = 0.125*(inputs.z(2:inputs.N)-duct_height(1)*log10(inputs.z(2:inputs.N)/inputs.evap_const)); M(2:inputs.N) = M(1) + term; % linear profile % (user defined, surface duct, surface_based duct, elevated duct) else num = length(duct_height); rg1 = 2; for i = 1:num-1, if i == (num-1), rg2 = inputs.N; else rg2 = round(duct_height(i+1)/inputs.delz); end; duct_grad = (duct_M(i+1)-duct_M(i))/(duct_height(i+1)-duct_height(i)); % 1/meter, gradient of the profile term = inputs.delz*duct_grad; for ii = rg1:rg2, M(ii) = M(ii-1)+term; end; rg1 = rg2+1; end; end; n = M*(1e-6) - inputs.z / inputs.ae + 1; m = n.*n - 1 + 2.*inputs.z / inputs.ae; %N = M-(1e6).*inputs.z / inputs.ae; %n = M*(1e-6) + 1; %m = n.*n - 1; %n = M*(1e-6) + 1; %m = n.*n - 1;
5、天线初始场计算
%% *********************************************************************** % INITIAL FIELD CALCULATION FROM ANTENNA PATTERN function [u0z] = fun_initial_field(inputs) u0z = 0; for i = 1:length(inputs.tx_height) ww = sqrt(2*log(2))./(inputs.beta*sin(inputs.thetabw(i)/2)); ufsp = 1/(sqrt(pi)*ww)*exp(-1i*inputs.beta*sin(inputs.thetae(i))*inputs.z).*exp(-((inputs.z-inputs.tx_height(i))/ww).^2); ufsn = 1/(sqrt(pi)*ww)*exp(-1i*inputs.beta*sin(inputs.thetae(i))*(-inputs.z)).*exp(-((-inputs.z-inputs.tx_height(i))/ww).^2); if inputs.polrz == 1, % H polrz u0z = u0z + (ufsp-ufsn); else % V polrz u0z = u0z + (ufsp+ufsn); end; end
6、
function [uz] = fun_field_pec(uz, inputs, prop, envpr) N = length(uz); if inputs.polrz == 1, % H polrz up = dst(uz(2:end-1)); % If Partial Differential Toolbox doesn't exist in your Matlab, dst % and dct functions may not work. In this case, disable above comand % and enable the following commands: %up = cat(2, -fliplr(uz(2:end)), uz); %up = fftshift(fft(ifftshift(up))); else % V polrz up = dct(uz); % If Partial Differential Toolbox doesn't exist in your Matlab, dst % and dct functions may not work. In this case, disable above comand % and enable the following commands: %up = cat(2, fliplr(uz(2:end)), uz); %up = fftshift(fft(ifftshift(up))); end; if inputs.polrz == 1, up = up.*prop(2:end-1); uz = [0; idst(up); 0]; % If Partial Differential Toolbox doesn't exist in your Matlab, dst % and dct functions may not work. In this case, disable above comand % and enable the following commands: %up = up.*cat(2, fliplr(prop(2:end)), prop); %uz = fftshift(ifft(ifftshift(up))); %uz = uz(N:end); %uz(1)=0; else up = up.*prop; uz = idct(up); % If Partial Differential Toolbox doesn't exist in your Matlab, dst % and dct functions may not work. In this case, disable above comand % and enable the following commands: %up = up.*cat(2, fliplr(prop(2:end)), prop); %uz = fftshift(ifft(ifftshift(up))); %uz = uz(N:end); end; uz = uz.*envpr;
7、
function [uz, dmft] = fun_field_mixed(uz, inputs, prop, envpr, alg, dmft) %tic N = length(uz); if alg == 1 %central difference DMFT w = (uz(3:end)-uz(1:end-2)) / (2*inputs.delz) + dmft.alfa*uz(2:end-1); w = [0; w; 0]; %********** if inputs.polrz == 1, % H polrz wp = dst(w(2:end-1)); else % V polrz wp = dct(w); end; if inputs.polrz == 1, wp = wp.*prop(2:end-1); w = [0; idst(wp); 0]; else wp = wp.*prop; w = idct(wp); end; %********** dmft.ck1 = dmft.ck1 * dmft.c1x; dmft.ck2 = dmft.ck2 * dmft.c2x; ym = zeros(1,N-1); for i = 2:N-1 ym(i) = 2*inputs.delz*w(i) + dmft.r * ym(i-1); end %********** uz(N) = 0; for i = 2:N, nmi = N-i+1; uz(nmi) = dmft.r*(ym(nmi)-uz(nmi+1)); end; %********** %ar = dmft.ck1 - dmft.rk*(0.5*uz(1) + 0.5*uz(end)*dmft.rn(end) + sum(uz(2:N-1).*dmft.rn(2:N-1))); ar = dmft.ck1 - dmft.rk*(0.5*uz(1) + 0.5*uz(end)*dmft.rn(end) + sum(uz.*dmft.rn)); uzf = fliplr(uz); br = dmft.ck2 - dmft.rk*(0.5*uz(1)*dmft.rn(end) + 0.5*uz(end) + sum(uzf(2:N-1) .* dmft.rn(2:N-1) .* (-1).^(1:N-2)')); %********** rnf = fliplr(dmft.rn); uz = uz + ar*dmft.rn + br*rnf .* (-1).^(N:-1:1)'; else %backward difference DMFT w = uz(2:N-1) - dmft.r * uz(1:N-2); w = [0; w; 0]; %********** if inputs.polrz == 1, % H polrz wp = dst(w(2:end-1)); else % V polrz wp = dct(w); end; if inputs.polrz == 1, wp = wp.*prop(2:end-1); w = [0; idst(wp); 0]; else wp = wp.*prop; w = idct(wp); end; %********** dmft.ck3 = dmft.ck3 * dmft.c3x; uz(1) = 0; for i = 2:N, uz(i) = w(i) + dmft.r * uz(i-1); end; ar = dmft.ck3 - sum(uz(1:N-1).*dmft.rn(1:N-1)); uz = uz + ar*dmft.rn; end uz = uz.*envpr;