function out = PE_ssft_duct
clear all;clc;close all;
%% 参数设置
freq = 5e9;
ht = 10;
bw = 3;
elv = 0;
pol = 0;
Nz = 8192;
Nx = 1500;
dx = 100;
dz = 3e8/freq;
c0 = 3e8; %光速m/s
waveL = c0/freq;
k0 = 2*pi/waveL;
x = (1:Nx)*dx;
z = (1:Nz)*dz;
%% 分配内存
u = zeros(1,2*Nz); %进行迭代的场;
uxz = zeros(Nx,Nz); %距离-高度剖面;
E = zeros(Nx,Nz);
PL = zeros(Nx,Nz);
PF = zeros(Nx,Nz);
%% 吸收层
abs_thick = round(Nz/5); sig0 = 0.0;%abs_thick 吸收层厚度
hamming = @(vx, sig0) 0.54 - 0.46*cos( pi*(vx * (acos((0.54-sig0)/0.46)/pi-1) +1) );
wind1D_up = @(NX, wind) [ones(NX-length(wind), 1); wind(:)];
abs_wind = wind1D_up( Nz, hamming(linspace(0, 1, abs_thick), sig0) );
abs_wind = [abs_wind; flipud(abs_wind)];
%% 地表处理
epr_ground = 80; %地面相对介电常数;
sigma_ground = 4; %地面电导率;
%% 高斯方向图
ww = @(k1) sqrt(2*log(2))/(k1*sind(bw/2));
ufs = @(vx,k1) 1/(sqrt(pi)*ww(k1))*exp(-1i*k1*sind(elv)*vx).*exp(-((vx-ht)/ww(k1)).^2); %PETOOL
ufss = @(vx,k1) ufs(vx,k1)+(-1)^pol*ufs(-vx,k1); %[幅值归一化]
%% PE算法核心程序
fft_kernel = @(f1, f2) ifft( f1 .* fft(f2) );
fft_k = @(Nfft,dt) fftshift((0:Nfft-1)*2*pi/Nfft/dt - pi/dt);
kz = fft_k(2*Nz, dz);
%n=Evap(30,z);% 蒸发波导
n=surf_nbase(40,-0.5,z); %表面波导
%n=elev(0.1,-0.5,80,30,z); %悬空波导
C1 = exp(-1i*k0*dx*(n-2)); %折射项
C2 = exp(-1i*dx*conj(sqrt(k0^2-kz.^2))); %绕射项
u(1:Nz) = ufss(z,k0); uxz(1,:)= u(1:Nz); %初始场
% 函数预定义
f_PL = @(vu,vx,vf) -20*log10(abs(vu)) + 20*log10(4*pi) + 10*log10(vx) - 30*log10(c0/vf); %传播损耗PL函数
f_PF = @(vx,vu,vf) 20*log10(sqrt(vx)*abs(vu)) + 10*log10(c0/vf); %传播因子PF函数
f_FreeLoss = @(vf,vx) 32.45 + 20*log10(vf/1e6) + 20*log10(vx/1e3); %FSPL函数 Hz、m
% 开始计算
tic,
for ix = 1:Nx %PE迭代计算
thetaR = atan(ht./(ix*dx)); %掠入射角度
R = RCoeffi(thetaR, epr_ground, sigma_ground, freq, pol); %反射系数
R0 = 1;
u(Nz+1:1:2*Nz) = R.*R0.*u(Nz:-1:1);
u = fft_kernel(C2,u).*abs_wind';
u(1:Nz) = u(1:Nz).*C1;
uxz(ix,:) = u(1:Nz); %注意是否除以sqrt(r)!!!
E(ix,:) = u(1:Nz)/sqrt(x(ix));
PL(ix,:) = f_PL(u(1:Nz),x(ix),freq);
PF(ix,:) = f_PF(x(ix),uxz(ix,:),freq);
end
out.ctime = toc;
FreeLoss = f_FreeLoss(freq,x)'; % 自由空间损耗/dB
%% 绘图
figure
todB = @(vx)20*log10(abs(vx)); %unit:dB
[X,Z]=meshgrid(x,z);mesh(X*1e-3,Z,todB(abs(E')));%view()
colorbar; colormap(jet); view(2);%colorbar('SouthOutside');
axis([x(1)*1e-3 x(end)*1e-3 0 120]); caxis([-120,-40]);
xlabel('range (km)');ylabel('height (m)'); title('X-Z PL(dBV/m)');
end
%%
function n = Evap(ductH,z)%蒸发波导
%输入参数:蒸发波导高度-ductH,高度向量-z;
z0 = 1.5e-4;
modifiedM = @(vx)(330 + 0.125*(vx-ductH*log((vx+z0)/z0)));
M = modifiedM(z); n = M*1e-6+1; %n = 1;
end
%%
function n = surf_nbase(d, c1, z)%无基础层表面波导
%波导高度d,陷获层斜率c1
M0 = 315; %地表修正折射率指数;
M = M0 + (z<=d).*z*c1+(z>d).*(0.125*(z-d)+c1*d); n = M*1e-6+1;
end
%%
function n = elev(c1, c2, z1, z2, z)%含基础层表面/悬空波导
M0 = 330;%地表修正折射率指数;
M1 = M0 + c1*z1 + c2*z2 + 0.118*(z - z1 - z2);
M2 = M0 + c1*z1 + c2*(z - z1);
M3 = M0 + c1*z;
M = (z>=z1+z2).*M1 + (z<z1+z2&z>z1).*M2 + (z<=z1).*M3; n = M*1e-6+1;
end
%%
function R = RCoeffi(theta, epsr, sigma, freq, pol)
%极化方式,1-水平,0-垂直:
%theta:掠入射角
%%
eps0 = 1/36/pi*1e-9; %真空介电常数
omega = 2*pi*freq;%角频率
epsrc = epsr - 1i*sigma/omega/eps0;%建筑物,复相对介电常数
if pol == 1 %极化方式,1-水平
Zg = sqrt(epsrc-cos(theta).*cos(theta));
else %0-垂直极化
Zg = sqrt(epsrc-cos(theta).*cos(theta))/epsrc;
end
R = (sin(theta)-Zg)./(sin(theta)+Zg);
end
在只改变参数名字和注释下修改代码,使最终呈现的功能和源代码一致。
最新发布