mie_bistaic_far_scattering_dg.m

本程序为matlab版的Mie级数展开程序

function Mie_bistaic_far_scattering()
%%本程序计算导体球的远场后向收发分置散射场


%%%%%%%%%%%%%%%%%%%%%%%%%%%%系统参数设置%%%%%%%%%%%%%%%%%%%%%%%%%
m=800+800*i;%%导体球的反射率
M=800;%%采样点数

% a=0.0745;%%球半径
a=0.0505;%%球半径


r=4.0/a;%%等效距离
c=3.0e+8;%光速

sita_temp=10;%%收发分置角
sita=pi-sita_temp*pi/180;%%接收天线位置

Sample_Frequency=M-1;%频域采样点数；
T=25e-012;%%时域采样间隔
Max_Frequency=(M-1)*( 1/(T*M));%最高频率
x_max=a*2*pi*Max_Frequency/c;
nmax=round(x_max+4*(x_max)^(1/3)+2)
%%%注：以下覆盖40G的范围
frequency=zeros(Sample_Frequency,1);

E_scatter=[];
jishuqi=0;

%%%%%%%%%%%%%%%%%%%%%%%%扫频计算各频点的散射场%%%%%%%%%%%%%%%%%
for ii=1:Sample_Frequency
    jishuqi=jishuqi+1
     
    frequency=ii*Max_Frequency/Sample_Frequency;
     x=a*2*pi.*frequency./c;%%把频率转化为波数
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     [an,bn]=mie_ab(m,x,nmax);%%计算给定频率点的系数
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     u=cos(sita);
     [p,t]=Mie_pt(u,nmax);%%计算给定频率点的Pi_n和Tao_n
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %%计算S2
     S2=0;
     %%ss=0;
     for nn=1:nmax
         S2=S2+(2*nn+1)*(an(nn)*t(nn)+bn(nn)*p(nn))/(nn*nn+nn);
         %%ss=ss+(2*nn+1)*(-1)^nn*(an(nn)-bn(nn));
     end
     
    E_scatter=[E_scatter;cos(pi)*S2*exp(j*x*r)/(-j*x*r)];
    %% E_scatter=[E_scatter;(abs(ss))^2/x/x];%%
end
E_Far_Scatter1=[0;E_scatter];
E_Far_Scatter=abs(E_Far_Scatter1).*cos(angle(E_Far_Scatter1))-j*abs(E_Far_Scatter1).*sin(angle(E_Far_Scatter1));%%相位反转
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算双高斯散射场%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
time=(0:M-1)*T;

t0=30*T;
arfa1=80.08*1e+18;
arfa2=180.35*1e+18;

A1=(arfa1)^0.5/((arfa1)^0.5-(arfa2)^0.5);
A2=(arfa2)^0.5/((arfa1)^0.5-(arfa2)^0.5);

%%产生双极高斯脉冲
Dgausi_signal=A1*exp(-arfa1*(time-t0).^2)-A2*exp(-arfa2*(time-t0).^2);
Dgausi_signal_analy=hilbert(Dgausi_signal');
input=Dgausi_signal_analy;
DoubGauss=fft(input);

%%%%%计算散射场
response_DG=E_Far_Scatter.*DoubGauss;

save sphere_theory_Far_scattering_5_05cm_10du E_Far_Scatter
f=0:40/800:40-40/800;
figure(1)
plot(f,abs(E_Far_Scatter1))
title('RCS')
xlabel('GHz')
%%axis([0 800 0 0.014])
figure(2)
plot(angle(E_Far_Scatter))
title('相位')

figure(3)
plot(time,real(ifft(E_Far_Scatter)))
title('脉冲响应')
xlabel('s')
figure(4)
plot(abs(response_DG))
title('双高斯激励下散射场的频谱')
figure(5)
plot(real(ifft(response_DG)))
title('双高斯激励下散射场时域波形')


% E_far_scatter2=abs(E_Far_Scatter).*cos(angle(E_Far_Scatter))-j*abs(E_Far_Scatter).*sin(angle(E_Far_Scatter));
% figure(6)
% plot(real(ifft(E_far_scatter2)))


%%计算an和bn
function [an,bn] = mie_ab(m,x,nmax)

% Computes a matrix of Mie Coefficients, an, bn, 
% of orders n=1 to nmax, for given complex refractive-index
% ratio m=m'+im" and size parameter x=k0*a where k0= wave number in ambient 
% medium for spheres of radius a;
% Eq. (4.88) of Bohren and Huffman (1983), BEWI:TDD122
% using the recurrence relation (4.89) for Dn on p. 127 and 
% starting conditions as described in Appendix A.


z=m.*x;  
nmx=round(max(nmax,abs(z))+16);
n=(1:nmax); nu = (n+0.5); 

sx=sqrt(0.5*pi*x);
px=sx.*besselj(nu,x);
p1x=[sin(x), px(1:nmax-1)];
chx=-sx.*bessely(nu,x);
ch1x=[cos(x), chx(1:nmax-1)];
gsx=px-i*chx; gs1x=p1x-i*ch1x;
dnx(nmx)=0+0i;
for j=nmx:-1:2      % Computation of Dn(z) according to (4.89) of B+H (1983)
    dnx(j-1)=j./z-1/(dnx(j)+j./z);
end;
dn=dnx(n);          % Dn(z), n=1 to nmax
da=dn./m+n./x; 
db=m.*dn+n./x;

an=(da.*px-p1x)./(da.*gsx-gs1x);
bn=(db.*px-p1x)./(db.*gsx-gs1x);

%%计算Pi_n和Tao_n
function [p,t]=Mie_pt(u,nmax)
% pi_n and tau_n, -1 <= u=cos(sita) <= 1, n1 integer from 1 to nmax 
% angular functions used in Mie Theory
% Bohren and Huffman (1983), p. 94 - 95

p(1)=1; 
t(1)=u;
p(2)=3*u; 
t(2)=3*cos(2*acos(u));
for n1=3:nmax,
    p1=(2*n1-1)./(n1-1).*p(n1-1).*u;
    p2=n1./(n1-1).*p(n1-2);
    p(n1)=p1-p2;
    t1=n1*u.*p(n1);
    t2=(n1+1).*p(n1-1);
    t(n1)=t1-t2;
end;
result=[p;t];