% 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;endfor 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));endthh=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