%%%%........矩量法计算一维介质粗糙面 HH 极化和 VV 极化的函数........%%%%
function [sca,bsc]=MoM_media(lamda,N,l,h,L,g,in,eps,p)
% 输入
% wave: 波长
% N: 表面点数
% l; 相关长度
% h: 均方根高度
% L: 表面长度
% g: 锥形因子
% in: 入射角
% epsilon: 相对介电常数
% p: rho A矩阵前的系数
% 输出
% sca: 散射角
% bsc: 双站散射系数
gamma = 1.78107; % 欧拉常数
inc = in*pi/180; % 入射角
k = 2*pi/lamda; % 自由空间波束
k1 = k*sqrt(eps); % 介质空间波束
dx = L/N; % 每段长度
% 计算双站散射系数分母部分
denom = 8*pi*k*g*sqrt(pi/2)*cos(inc)*... % 连接符
(1-(1+2*(tan(inc)).^2)/(2*(k*g*cos(inc)).^2)); % 双站散射系数分母
[x,f,f1,f2] = surface(N,l,h,L); % 调用粗糙面子程序
b = incid(k,x,f,inc,g); % 调用入射锥形波子程序，生成矩阵b
% for m = 1:N;
% for n = 1:N
% dl = sqrt(1+f1(n)^2)*dx;
% if m==n;
% A(m,n) = (1i/4)*dx*(1+(2*1i/pi)*(log(gamma*k*dl/4)-1));
% B(m,n) = 0.5;
% A1(m,n) = (1i/4)*dx*(1+(2*1i/pi)*(log(gamma*k1*dl/4)-1));
% B1(m,n) = -0.5;
% else
% r = sqrt((x(m)-x(n))^2+(f(m)-f(n))^2);
% A(m,n) = dx*1i*besselh(0,1,k*r)/4;
% B(m,n) = -dx*(1i*k/4)*(-f1(n)*(x(m)-x(n))+(f(m)-f(n)))/r*besselh(1,1,k*r);
% A1(m,n) = dx*1i*besselh(0,1,k1*r)/4;
% B1(m,n) = -dx*(1i*k1/4)*(-f1(n)*(x(m)-x(n))+(f(m)-f(n)))/r*besselh(1,1,k1*r);
% end
% end
% end
% 32-48行与51-79行一样，meshgrid可以快速运算
% 计算矩阵A
[xm,xn] = meshgrid(x);
[fm,fn] = meshgrid(f);
gap = k*sqrt((xm-xn).^2+(fm-fn).^2); % 导体问题第一类零阶和一阶汉克尔函数的参数
gap1 = k1*sqrt((xm-xn).^2+(fm-fn).^2); % 介质问题第一类零阶和一阶汉克尔函数的参数
for m = 1:N;
gap(m,m) = 1;
gap1(m,m) = 1;
end
A = dx*1i*besselh(0,1,gap)/4; % 利用汉克尔函数生成m~=n时的矩阵A
A1 = dx*1i*besselh(0,1,gap1)/4;
K = (1i*k/4)*besselh(1,1,gap); % 第一类一阶汉克尔函数
B = -dx*K; % 生成m~=n时的矩阵A
K1 = (1i*k1/4)*besselh(1,1,gap1);
B1 = -dx*K1;
for m = 1:N;
dl = sqrt(1+f1(m)*f1(m))*dx;
A(m,m) = (1i/4)*dx*(1+(2*1i/pi)*(log(gamma*k*dl/4)-1)); % m==n时的A矩阵
A1(m,m) = (1i/4)*dx*(1+(2*1i/pi)*(log(gamma*k1*dl/4)-1));
I = dx/(4*pi)*f2(m)/(1+(f1(m))^2);
B(m,m) = 0.5-I; % m==n时的B矩阵
B1(m,m) = -0.5-I;
end
M = [A B;p*A1 B1]; % 总阻抗矩阵
z = zeros(N,1); % 下方无入射激励，激励矩阵为零
Z = [b';z]; % 总激励矩阵
H = M\Z; % 为M*H=Z形式，求解H矩阵
H = H';
U = H(1,1:N);
psin = H(1,N+1:2*N);
for n = 1:181;
sca(n) = -90+(n-1); % 散射角度
sc = sca(n)*pi/180;
integ = exp(-1i*k*(sin(sc)*x+f*cos(sc)));
fac = 1i*k*(f1*sin(sc)-cos(sc));
integ = integ.*(-U+psin.*fac)*dx;
psis = sum(integ); % 散射波
bsc(n) = abs(psis)^2/denom; % 双站散射系数BSC
end
% 入射锥形波函数
function b = incid(k,x,z,inc,g)
% 输入：
% k: 空间波束
% x: 横坐标
% z; 高度函数
% inc: 入射角
% g: 锥形因子
% 输出：
% b: 入射锥形波
ki = k*(x*sin(inc)-z*cos(inc)); % 入射波矢量ki
r = ((x+z*tan(inc))/g).^2;
w = (2*r-1)/((k*g*cos(inc)).^2); % 附加相位项w
b = exp(1i*ki.*(1+w)-r);
% 高斯随机粗糙面函数
function [x,f,f1,f2] = surface(N,l,h,L)
% 输入：
% N： 表面点数
% l： 相关长度
% h： 均方根高度
% L： 表面长度
% 输出
% x: 横坐标
% f: 高度
% f1: 一阶导数
% f2: 二阶导数
format long;
x = linspace(-L/2, L/2, N);
Z = h.*randn(1,N);
F = exp(-x.^2/(l^2/2));
f = sqrt(2/sqrt(pi))*sqrt(L/N/l)*ifft(fft(Z).*fft(F));
dx = L/N;
for n = 1:N;
if n==1;
f1(1) = (f(2)-f(N))/(2*dx);
elseif n==N;
f1(N) = (f(1)-f(N-1))/(2*dx);
else
f1(n) = (f(n+1)-f(n-1))/(2*dx);
end
end
for n = 1:N;
if n==1;
f2(1) = (f1(2)-f1(N))/(2*dx);
elseif n==N;
f2(N) = (f1(1)-f1(N-1))/(2*dx);
else
f2(n) = (f1(n+1)-f1(n-1))/(2*dx);
end
end
评论5