一、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;