sphere_rcs.m

这是计算各种形状的散射程序

%   sphere_rcs.m
%================================================================
%  MATLAB code for the calculation of
%   Frequency Response 
% of the backscattering from conducting spheres. Exact solution
% of wave equation.  
% Syntax:
%		[sigma,phi]=sphere_rcs(radius,fs,fe,nf)
% Input: 
%		fs--start frequency (MHz)
%		fe--Stop (End) frequency (MHz) 
%		nf--No. of frequency samples
%		radius -- Radius of the sphere (cm)
% Output:
%		sigma--Amplitude of the RCS (dBsm)  
%		phi--Phase of the RCS (deg.)
%================================================================
function [sigma,phi]=sphere_rcs(radius,fs,fe,nf)

  df=(fe-fs)/real(nf-1);  %MHz
  
  for i=1:nf
     f=fs+(i-1)*df; %frequency (MHz)
     wavlen=30000.0/f; % wavelen -- wavelenth (cm)      
     ka=2*pi*radius/wavlen;
     
     [sigma0,phi0]=sphrcs(ka);
     
     sigma(i)=10*log10(sigma0*pi*radius*radius/10000.0+eps); %(dBsm)
     phi(i)=phi0*180/pi; %(deg.)  		   
  end   
  
% Graphics
  
  hs=sphere_graph(fs,fe,sigma,phi,radius);
    
% ----------------------------------------------------------------
% Normalized RCS amplitude and phase for a conducting sphere
%
%  Syntax: [sigma,phi]=sphrcs(ka)
%
%  Input: 
%        ka -- a=the radius of the sphere, k=2*pi/wavelength
%  Output:
%        sigma -- RCS Amplitude normalized by pi*a^2/4.0
%        phi -- RCS phase (deg.)
% ----------------------------------------------------------------

function [sigma,phi]=sphrcs(ka)
   sg=0.0+j*0.0;
   sg0=100.+j*100.;
   n=1;
   eps0=1e-12;  %accuracy control constant
   
   while abs(sg0)>eps0 & n<1000 
 	    [aj0,y0]=shankl(ka,n-1);
	    [aj1,y1]=shankl(ka,n);
	    han0=aj0+j*y0;
	    han1=aj1+j*y1;
	    an=aj1/han1;
	    bn=(ka*aj0-n*aj1)/(ka*han0-n*han1);
	    sg0=(-1)^n*(n+.5)*(bn-an);
	    sg=sg+sg0;
       n=n+1;
   end
   sigma=(2./ka*abs(sg))^2;
   phi=atan2(real(sg),imag(sg));
	
% ------------------------------------------------------------
	function [aj,y]=shankl(x,n)
% ------------------------------------------------------------
   i=1:1:1000;
   b=0.0*i;
   if x>=250.0
      m=round(1.03*x+65.0);
   elseif x>=45.0
      m=round(1.1*x+40.0);
   elseif x>=10.0    
      m=round(1.4*x+25.0);
   else 
      m=round(2.5*x+13.0);
   end  
   am=m;
   a=sin(x);
   if abs(a)>1e-6
      m1=0;
   else
      a=cos(x);
      m1=-1;
   end   
      y1=0.0;
      y2=1.e-12;
   while m1<m    
      if m<1000 
         b(m+1)=y2;
      end
      y3=(2.*am+1.0)*y2/x-y1;
      m=m-1;
      am=m;
      y1=y2;
      y2=y3;
   end      
   if m1==0 
      b(1)=y2;
   end   
   b=b*a/(y2*x);
   y=-cos(x)/x;
   y2=y;
   if n~=0 
     m=0;
     am=0.0;
     y1=sin(x)/x;
     while m<n
       y=y2*(2.*am+1.)/x-y1;
       m=m+1;
       am=m;
       y1=y2;
       y2=y;
     end
   end   
   aj=b(n+1);
   
% Graphics  
%-----------------------------------------------------
function hs=sphere_graph(fs,fe,sigma,phi,radius)
%-----------------------------------------------------
m0=max(size(sigma));

ar=10.^(sigma./20.);
ai=phi*pi/180.;
dat=ar.*cos(ai)+j*ar.*sin(ai);
  
% calculate the impulse response

hamm=hamming(m0); %hamming window for IFFT
dat=dat.*hamm';
fdat=fftshift(ifft(dat,m0*2));
hs=real(fdat)/.2682;  %impulse response,.2682 is a modification factor for IFFT & Hamming
amp2=20*log10(abs(fdat)/.2682); %range profile

i=1:1:m0*2;
r(i)=(i-m0-1)*150.0/(fe-fs)/2;

fstep=(fe-fs)/(m0-1);
i=1:m0;
f(i)=fs+(i-1)*fstep;

figure
subplot(221),plot(f,sigma)
axis([fs fe min(sigma)-1 max(sigma)+1])
grid on
xlabel('Frequency (MHz)')
ylabel('RCS (dBsm)')
title('RCS Magnitude of the Sphere')

subplot(222),plot(f,phi)
axis([fs fe -180 180])
grid on
xlabel('Frequency (MHz)')
ylabel('Phase (deg.)')
title('RCS Phase of the Sphere')

tt=r*100/15;
rscale=fix(6*radius/100/r(m0*2)*m0)

subplot(223),plot(tt,hs)
axis([tt(m0-rscale),tt(m0+rscale),min(hs)*1.1 max(hs)*1.1]),grid on
xlabel('Time delay (ns)')
ylabel('Amplitude')
title('Impulse Response')

subplot(224),plot(r,amp2)
axis([r(m0-rscale),r(m0+rscale),max(amp2)-55 max(amp2)+1]),grid on
xlabel('Down range (m)')
ylabel('RCS (dBsm)')
title('Range profile')

% END of the code