% dihedral_rcs_f.m% ==============================================================% MATLAB code for the calculation of%  frequency Response of the Polarization Scattering Matrix % of dihedral corners near boresight, by using PO approximation,% with consideration of single-bounce and double-bounce % contributions.% Syntax:% 		[dhh,dvv,dhv]=dihedral_rcs_f(a,b,l,beta,tau,fmin,fmax,phi,nf)% Input:%		a,b -- lengths of the sides of the dihedral corner (cm)%		l -- Height of the dihedral corner (cm)% 		tau-- angle tilted above the local horizontal (deg.)%		beta -- the angle between the two faces (deg.)%		fmin -- minimum frequency (MHz)%		fmax -- maximum frequency (MHz)%	 	phi -- azimuthal angle (-50=<phi<=50 deg.)%		nf -- No. of frequency samples% Output:%		dhh -- complex RCS for H-H polarization %		dvv -- complex RCS for V-V polarization%		dhv -- complex RCS for H-V polarization (dhv=dvh)% Conditions when holding accuracy: %		a,b,l>= several wavelengths% ==============================================================function [dhh,dvv,dhv]=dihedral_rcs_f(a,b,l,beta,tau,fmin,fmax,phi,nf)c=300.0;  % (m/s)a=a/100.; %(m)b=b/100.; %(m)l=l/100.; %(m)beta0=beta/2.; %(deg.)phi=phi*pi/180.; %(rad.)fstep=(fmax-fmin)/(nf-1);j=sqrt(-1);for i=1:nf   freq=fmin+fstep*(i-1);	wavelength=c/freq;  % (m)	k=2*pi/wavelength;   % calculation of sa,sb   ka=k*a;   kb=k*b;   beta=beta0*pi/180.;   beta3=beta*3.;      x=ka*cos(beta+phi)/pi;   sa=-j*ka*l/wavelength*sin(beta+phi)*exp(-j*ka*cos(beta+phi))*sinc(x);      x=kb*cos(beta-phi)/pi;   sb=-j*kb*l/wavelength*sin(beta-phi)*exp(-j*kb*cos(beta-phi))*sinc(x);   % calculation of sab,sba   gama=atan2(b*sin(2.*beta),a-b*cos(2*beta));   alpha=pi-beta3;      if phi<=-alpha      aa=0.0;   elseif phi>=gama-alpha      aa=b*sin(beta-phi)/sin(beta3-phi);   else      aa=a;   end      if phi>=alpha      bb=0;   elseif phi<=gama-beta      bb=a*sin(beta+phi)/sin(beta3+phi);   else      bb=b;   end    %e-plane   x=k*bb*cos(2*beta)*cos(beta+phi)/pi;   sin3=sin(beta3+phi);   sab=-j*k*bb*l/wavelength*exp(-j*k*bb*cos(beta*2)*cos(beta+phi))*sinc(x);   sabe=sab*sin3;      x=k*aa*cos(2*beta)*cos(beta-phi)/pi;   sin3=sin(beta3-phi);   sba=-j*k*aa*l/wavelength*exp(-j*k*aa*cos(beta*2)*cos(beta-phi))*sinc(x);   sbae=sba*sin3;    %h-plane   sin3=-sin(beta-phi);   sabh=sab*sin3;      sin3=-sin(beta+phi);   sbah=sba*sin3;   %Total e&h-plane fields   tvv(i)=(sa+sb+sabh+sbah)*wavelength/pi;   thh(i)=(sa+sb+sabe+sbae)*wavelength/pi;enddhh=thh*cos(2*tau*pi/180.);dvv=tvv*cos(2*tau*pi/180.);dhv=-thh*sin(2*tau*pi/180.);%Graphicsgresult=dihedral_graph(dhh,dvv,dhv,fmin,fmax);% Graphics for dihedral_rcs% -------------------------------------------------------------function gresult=dihedral_graph(dhh,dvv,dhv,fmin,fmax)np=max(size(dhh));pstep=(fmax-fmin)/(np-1);figurei=1:np;p=((i-1)*pstep+fmin);thha=20*log10(abs(dhh)+eps);thhp=atan2(imag(dhh),real(dhh))*180/pi;thhamin=max([min(thha),-75]);thhamax=max([max(thha),-75]);subplot(321)plot(p,thha,'b-')axis([fmin fmax thhamin-5 thhamax+5])grid ontitle('RCS Magnitude of the Dihedral (HH)')%xlabel('Frequency (MHz)')ylabel('RCS (dBsm)')subplot(322)plot(p,thhp);axis([fmin fmax -180 180]);grid ontitle('RCS Phase of the Dihedral (HH)')%xlabel('Frequency (MHz)')ylabel('Phase (deg.)')tvva=20*log10(abs(dvv)+eps);tvvp=atan2(imag(dvv),real(dvv))*180/pi;tvvamin=max([min(tvva),-75]);tvvamax=max([max(tvva),-75]);subplot(323)plot(p,tvva,'b-')axis([fmin fmax tvvamin-5 tvvamax+5])grid ontitle('RCS Magnitude of the Dihedral (VV)')%xlabel('Frequency (MHz)')ylabel('RCS (dBsm)')subplot(324)plot(p,tvvp);axis([fmin fmax -180 180]);grid ontitle('RCS Phase of the Dihedral (VV)')%xlabel('Frequency (MHz)')ylabel('Phase (deg.)')thva=20*log10(abs(dhv)+eps);thvp=atan2(imag(dhv),real(dhv))*180/pi;thvamin=max([min(thva),-75]);thvamax=max([max(thva),-75]);subplot(325)plot(p,thva,'b-')axis([fmin fmax thvamin-5 thvamax+5])grid ontitle('RCS Magnitude of the Dihedral (HV)')xlabel('Frequency (MHz)')ylabel('RCS (dBsm)')subplot(326)plot(p,thvp);axis([fmin fmax -180 180]);grid ontitle('RCS Phase of the Dihedral (HV)')xlabel('Frequency (MHz)')ylabel('Phase (deg.)')gresult=1;