pem.m

计算电波传播的一个MATLAB程序,用抛物线法求解的..还是初始阶段,共享一下.

function PEM(hs,frequency)
epsilon0=1e-9/(36*pi);
mu0=4*pi*1e-7;
const=3e8;
k0=2*pi*frequency*1e6/const;

x=0.1:180:200e3;
height=0.1:.5:500;
for k=1:length(x)  %
    tantheta(:,k)=abs(height' - hs)./(x(k));
    dist(:,k) = sqrt((height' - hs).^2 + x(k)^2);
end
theta=atan(tantheta);

thetabw=15/180*pi;
direct=(4*pi)/(thetabw)^2;
freep = sqrt(60*direct) ./ dist;

p=k0*sin(theta(:,1));
omega=sqrt(2*log(2))/(k0*sin(thetabw/2));
fp=exp(-p.^2*omega^2/4);
% psibar=fp.*exp(-i*p*hs) - conj(fp).*exp(i*p*hs);
psibar=2*real(fp.*exp(-i*p*hs) ); % .*exp(-i*k0*dist(:,1))./dist(:,1) );
psi0=ifft(psibar);
q=psi0;
for m=2:length(x)
    q=exp(i*(k0).*.1*height'/6400000).*ifft(exp(i*.1*(sqrt((k0)^2-(k0.*(sin(theta(:,m)))).^2)-k0)).*fft(q));
    psi0=[psi0 q];
end

% figure;
% plot(height,theta(:,1)'*180/pi,'r');
% ylabel('theta');
% figure;
% plot(height,20*log10(freep(:,1)'),'b');
% ylabel('freep');
% figure;
% plot(height,p,'r');
% ylabel('p');
% figure;
% plot(height,fp);
% ylabel('fp');
% figure;
plot(height,20*log10(abs(psibar)));
ylabel('psibar');
figure;
plot(height,20*log10(abs(psi0)));
ylabel('psi0');

figure;
h1=40;
plot(x,-20*log10(abs(psi0(h1,:))),'r');
xlabel('range in km');
ylabel('power in dB');
titletext = strcat('场强顺距离的变化，高度为(m)： ',num2str(h1*0.5));
title(texlabel(titletext));
hold on;
h2=50;
plot(x,-20*log10(abs(psi0(h2,:))),'g');
h3=60;
plot(x,-20*log10(abs(psi0(h3,:))),'b');
legend(num2str(0.5*h1),num2str(0.5*h2),num2str(0.5*h3));
hold off;