stcylinder_rcs.m

这是计算stchlinder的散射程序

% stcylinder_rcs.m
% =========================================================================
% MATLAB code for the calculation of
%	Frequency Response of the Polarization Scattering Matrix 
% of short and thin cylinder at broadside, by converting 2-D (infinite bodies) 
% solutions of the wave equation in cylindrical coordinates to the 3-D 
% (finite bodies) cases. Simplified formulas are used. 
% Syntax:
%		[thh,tvv,thv]=stcylinder_rcs(a,l,tau,fmin,fmax,nf)
% Input:
%		a -- radius of the cylinder (cm)
%		l -- length of the cylinder (cm)
% 		tau-- angle tilted above the local horizontal (deg.)
%		fmin -- minimum frequency (MHz)
%		fmax -- maximum frequency (MHz)
%		nf -- No. of frequency samples
% Output:
%		thh -- complex RCS for H-H polarization 
%		tvv -- complex RCS for V-V polarization
%		thv -- complex RCS for H-V polarization (thv=tvh)
% Conditions when holding accuracy: 
%		both the radius and the length of the cylinder should be 
%   much less than the radar wavelength
% ==========================================================================

function [thh,tvv,thv]=stcylinder_rcs(a,l,tau,fmin,fmax,nf)

l=l/100;    % (m)
a=a/100.;

tau=tau*pi/180.;  % (rad.)

j=sqrt(-1);
c=30000.0;  % (cm/s)

if nf~=1
	fstep=(fmax-fmin)/(nf-1);
else 
   fstep=0.0;
end

for i=1:nf
   freq(i)=fmin+(i-1)*fstep;
   wavelength=c/freq(i)/100.;  % (m)
   k=2*pi/wavelength;
   
   te(i)=k^2*l^3/(3*(log(4*l/a)-1));
end

thh=te.*sin(tau)^2;
tvv=te.*cos(tau)^2;
thv=te.*(-sin(tau)*cos(tau));

gresult=stcylinder_graph(thh,tvv,thv,freq);

% Graphics
% --------------------------------------------------
function gresult=stcylinder_graph(thh,tvv,thv,freq)   
% --------------------------------------------------
nf=max(size(freq));

if nf~=1
   figure
   
   thha=20*log10(abs(thh)+eps);
   tvva=20*log10(abs(tvv)+eps);
   thva=20*log10(abs(thv)+eps);
   amax=max([max(thha) max(tvva) max(thva)]);
   amin=min([min(thha) min(tvva) min(thva)]);
   if amin<-65
      amin=-65;
   end
   
   thhp=atan2(imag(thh),real(thh))*180/pi;
   tvvp=atan2(imag(tvv),real(tvv))*180/pi;
   thvp=atan2(imag(thv),real(thv))*180/pi;
   
   subplot(321)
   plot(freq,thha)
   axis([freq(1) freq(nf) amin-5 amax+5])
   grid on
   title('Magnitude of RCS (HH)')
  % xlabel('Frequency (MHz)')
   ylabel('RCS (dBsm)')
   
   subplot(322)
   plot(freq,thhp)
   axis([freq(1) freq(nf) -180 180])
   grid on
   title('Phase of RCS (HH)')
  % xlabel('Frequency (MHz)')
   ylabel('Phase (deg.)')
   
   subplot(323)
   plot(freq,tvva)
   axis([freq(1) freq(nf) amin-5 amax+5])
   grid on
   title('Magnitude of RCS (VV)')
  % xlabel('Frequency (MHz)')
   ylabel('RCS (dBsm)')
   
   subplot(324)
   plot(freq,tvvp)
   axis([freq(1) freq(nf) -180 180])
   grid on
   title('Phase of RCS (VV)')
  % xlabel('Frequency (MHz)')
   ylabel('Phase (deg.)')
   
   subplot(325)
   plot(freq,thva)
   axis([freq(1) freq(nf) amin-5 amax+5])
   grid on
   title('Magnitude of RCS (HV)')
   xlabel('Frequency (MHz)')
   ylabel('RCS (dBsm)')
   
   subplot(326)
   plot(freq,thvp)
   axis([freq(1) freq(nf) -180 180])
   grid on
   title('Phase of RCS (HV)')
   xlabel('Frequency (MHz)')
   ylabel('Phase (deg.)')
	gresult=1;   
else
   thha=20*log10(abs(thh)+eps)
   thhp=atan2(imag(thh),real(thh))*180/pi
   tvva=20*log10(abs(tvv)+eps)
   tvvp=atan2(imag(tvv),real(thh))*180/pi
   thva=20*log10(abs(thv)+eps)
   thvp=atan2(imag(thv),real(thv))*180/pi
   gresult=0;
end