longley-rice模型的matlab实现

以下是一个简单的matlab代码实现longley-rice模型：

function [pl] = longley_rice(f, d, h1, h2, pol, T, P, eps_r)

% f：频率（MHz）
% d：距离（km）
% h1：发射天线高度（m）
% h2：接收天线高度（m）
% pol：偏振（0-垂直，1-水平）
% T：温度（K）
% P：大气压力（hPa）
% eps_r：相对介电常数

% 常数
c = 299792458; % 光速
k = 1.38064852e-23; % 玻尔兹曼常数

% 计算传播路径损耗
lambda = c / (f * 1e6);
gamma = 1.333 - 0.267*log10(f);
gamma_v = 7.56 - 7.52*exp(-(h1/1000)^0.726);
gamma_h = 9.07 - 3.44*log10(h1);
delta_h = h2 - h1;
theta_e = atan(delta_h / (d * 1000));
theta_e = max(theta_e, -0.785398163);
theta_e = min(theta_e, 0.785398163);
d_e = d * 1000 / cos(theta_e);
k_e = 0.00838 * (f/1000)^0.588 + 0.01715*exp(0.0896*(T/273-1)) * (f/1000)^0.7;
k_v = k_e * (1.607 - 1.6088 / (1 + (f/1000)^2.61)) * (P/1013)^(-0.753);
k_h = k_e * (1.333 - 1.3373 / (1 + (f/1000)^2.61)) * (P/1013)^(-0.742);
alpha = (1.607 - 1.6088 / (1 + (f/1000)^2.61)) * (P/1013)^(-0.753) * (d_e/1000)^(-1.607) * (f/1000)^0.912;
beta = (1.333 - 1.3373 / (1 + (f/1000)^2.61)) * (P/1013)^(-0.742) * (d_e/1000)^(-1.333) * (f/1000)^0.912;
gamma_p = 1.444 * log10(f) - 0.02 * (h1-h2) * log10(f) - 0.267;
gamma_r = gamma_p - 0.1 * (gamma_p - 6.9) * exp(-0.12 * (h2/1000-5));

if (pol == 0)
    a_v = -20*log10(1 + exp(-1.7858+1.448*log10(f)-0.157*f^0.407)-exp(-3.8052+1.3296*log10(f)-0.1414*f^0.407));
    a_h = -20*log10(1 + exp(-1.7858+1.448*log10(f)-0.157*f^0.407)-exp(-3.8052+1.3296*log10(f)-0.1414*f^0.407)) - gamma_h;
    b_v = -22.2*log10(d_e) - 3.2*log10(d_e)^2 + 0.1*log10(f)*log10(d_e) + gamma_v;
    b_h = -20.8*log10(d_e) - 4.2*log10(d_e)^2 + 0.1*log10(f)*log10(d_e) + gamma_h;
    c_hv = 0.03 * (h2-h1) * exp(-sqrt((d_e/1000)/(2*h1/1000)));
    c_vh = c_hv;
elseif (pol == 1)
    a_h = -20*log10(1 + exp(-1.7858+1.448*log10(f)-0.157*f^0.407)-exp(-3.8052+1.3296*log10(f)-0.1414*f^0.407));
    a_v = -20*log10(1 + exp(-1.7858+1.448*log10(f)-0.157*f^0.407)-exp(-3.8052+1.3296*log10(f)-0.1414*f^0.407)) - gamma_v;
    b_h = -22.2*log10(d_e) - 3.2*log10(d_e)^2 + 0.1*log10(f)*log10(d_e) + gamma_h;
    b_v = -20.8*log10(d_e) - 4.2*log10(d_e)^2 + 0.1*log10(f)*log10(d_e) + gamma_v;
    c_hv = 0.03 * (h2-h1) * exp(-sqrt((d_e/1000)/(2*h1/1000)));
    c_vh = c_hv;
end

a = a_v + (a_h - a_v) / (1 + (f/(38.3*sqrt(eps_r)))^2);
b = b_v + (b_h - b_v) / (1 + (f/20.6)^2);
c = c_hv + (c_vh - c_hv) / (1 + (f/39.8)^1.8);

pl = a + b*log10(d_e) + c*log10(d_e)*log10(d_e) - alpha - beta*log10(f) + gamma_p - gamma_r;
pl = max(pl, 0);

end

该函数的输入参数包括：频率（MHz）、距离（km）、发射天线高度（m）、接收天线高度（m）、偏振（0-垂直，1-水平）、温度（K）、大气压力（hPa）、相对介电常数。输出为传播路径损耗（dB）。