一维介质粗糙面双站散射系数计算

%%%%........矩量法计算一维介质粗糙面 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