% =================================================================
%	MATLAB Function for Sphere-tiped-Cone RCS Calculation
% =================================================================
%
% Syntax:
%    [sver,shor]=sphcone_rcs(iangle,freq,alfa,a,b,h)
% 
% Input parameters:  
%    iangle: Incident angle (to the axis of symmetry of the sphere-
%		     tipped cone, [-pi,pi] radians;
%	 freq:  frequency, GHz;
%	 alfa:  the half cone angle, radians;
%	 a:  Base radius, cm;
%	 b:  Sphere radius, cm;
%	 h:  Cone length, cm.
%	
% Output items: 
%    sver   complex RCS for vertical polarization (voltage)
%    shor   Complex RCS for horizon polarization (voltage)
%
% Authored by Xiaojian Xu, March 10, 2003

function [sver,shor]=sphcone_rcs(iangle,freq,alfa,a,b,h)
% Check input parameters
	if a==0.0  % base radius a=0 illegal
	  disp('Sphere-tipped Cone structure unreasonable!')
	  return
    end
	if iangle<0 
       iangle=-iangle;
    end
    if iangle>pi
        iangle=2*pi-iangle;
    end

% Common parameters preparation
    wavelength=30.0/freq; %cm
    k=2*pi/wavelength; % wave number
    n=1.5+alfa/pi;
	x=2.*k*a*sin(iangle);
	if iangle~=0.0&iangle<pi 
	   t0=sin(pi/n)/n*sqrt(a/(sin(iangle)*k));
    end
	t1=1./(cos(pi/n)-1.);
	t=pi/n;
	ca=asin(1.22/(k*a));

% NOTE: Singularities may arise in RCS equations of the scattering
% centers S1 & S2 , at axial aspects(iangle=0,pi) and the broadside
% aspect (iangle=pi/2-alfa).

%	****   The case around iangle = 0	 ****
	if iangle<=alfa     %0<=iangle<=alfa
	   v0=2.*sqrt(pi)*a*sin(t)/n;
	   v1=1./(cos(t)-cos(3.*t));
	   v2=2.*tan(iangle)*sin(3.*t)/n;
	   v3=(cos(t)-cos(3.*t))^2;
	   v4=1./(cos(t)-1.);
       if abs(x)<=1.e-5 
	      bj0=1.;
	      bj1=0.;
	      bj2=0.;
	   else
          bj0=besselj(0,x); %call bessel(x,bj0,0)
	      bj1=besselj(1,x); %call bessel(x,bj1,1)
	      bj2=besselj(2,x); %call bessel(x,bj2,2)
       end 
       realh=v0*(v1*bj0-v4*bj2);
	   realv=v0*(v1*bj0+v4*bj2);
	   simage=-1.*v0*v2/v3*bj1;
       sver=realv+sqrt(-1)*simage;
	   shor=realh+sqrt(-1)*simage;
       
%	****   The case around iangle = pi	  ****
%	To avoid occurance of singularity,calculate RCS of
%	the Circle-Plate (radius=a).
    elseif iangle>pi-ca  %iangle>pi-ca
	   if abs(x)<=1.e-5
		  sever=k*k*pi*a^4;
		  sehor=sever;
		  sver=complex(0.,-1.*sqrt(sever));
		  shor=sver;
	   else
		  bj1=besselj(1,x);  %call bessel(x,bj1,1)
		  v=2.*sqrt(pi)*k*a*a*bj1/x;
		  sver=sqrt(-1)*v;
		  shor=sver;
       end
% ****  Other incidental angle cases  ****          
    else    

% scattering center-2
	   [ctp1,ctp2]=cones1(iangle,k,a,h,alfa,n,t0,t1);
%  At the broadside aspect, iangle=pi/2-alfa, we modify RCS equation
%  with the Physical-Optic Method.
	  if abs(iangle-(pi/2-alfa))<=1.75e-2
		  v0=4./9.*k*cos(alfa)/tan(alfa)^2;
		  v1=sqrt(a*a*a)-sqrt((b*cos(alfa))^3);
		  v1=sqrt(v0*v1*v1);
		  phase=pi/4.-2.*k*(a*sin(iangle)+h/2*cos(iangle));
		  s1ver=complex(v1*cos(phase),v1*sin(phase));
		  s1hor=s1ver;
		  if abs(s1ver)>abs(ctp1)
		     s1ver=ctp1;
		     s1hor=ctp2;
          end
      else
          s1ver=ctp1;
		  s1hor=ctp2;
      end
          
% scattering center-3          
	   [s2ver,s2hor]=cones2(iangle,k,a,h,alfa,n,t0,t1);
       sver=s1ver+s2ver;
       shor=s1hor+s2hor;
   end   
      
% sphere tip scattering  
   s0=sph_rcs(iangle,k,b,a,h,alfa);

% total scattering
   sver=sver+s0;
   shor=shor+s0;
% scale to (m) from (cm)
   sver=sver/100.0;
   shor=shor/100.0;

%------------------------------------------------------------------
% Calculate RCS of the single scattering centers
%------------------------------------------------------------------
% RCS of the spherical part of sphere-tipped cone  
function s=sph_rcs(iangle,k,b,a,h,alfa)

    if b==0.0 % radius=0
	   s=complex(0.0,0,0);
       return
    end
	angle=pi/2.-alfa;
	if iangle>=angle  % exceeding broadside in azimuth  
	   s=complex(0.0,0.0);
	else
	   if iangle>0.0
	      rcs=sqrt(pi)*b;
	   else  %iangle=0
	      t1=2.*k*b*(1.-sin(alfa));
	      t2=k*b*cos(alfa)*cos(alfa);
	      temp=1.-sin(t1)/t2;
	      rcs=sqrt(abs(pi*temp))*b;
       end
        phase=-2*k*cos(iangle)*(h/2+b);
	   s=complex(rcs*cos(phase),rcs*sin(phase));
    end

% RCS of the second scattering center: s1ver,s1hor
function [sver,shor]=cones1(iangle,k,a,h,alfa,n,t0,t1)
	t2=cos(pi/n)-cos((3.*pi-2.*iangle)/n);
	t2=1./t2;
	rcsver=t0*(t1-t2);
	rcshor=t0*(t1+t2);
     phase=pi/4.-2*k*(a*sin(iangle)-h/2*cos(iangle));
	
    sver=complex(rcsver*cos(phase),rcsver*sin(phase));
	shor=complex(rcshor*cos(phase),rcshor*sin(phase));
 
% RCS of the third scattering center
function [sver,shor]=cones2(iangle,k,a,h,alfa,n,t0,t1)
    if iangle>alfa&iangle<pi/2.
	   sver=0;
	   shor=0;
       return
    elseif iangle<=alfa&iangle>0.0
	   t2=cos(pi/n)-cos((3.*pi+2.*iangle)/n);
	else
	   t2=cos(pi/n)-cos((pi-2.*iangle)/n);
    end
    t2=1.0/t2;
	rcsver=t0*(t1-t2);
	rcshor=t0*(t1+t2);	
 	phase=-pi/4+2*k*(a*sin(iangle)-h/2*cos(iangle));
	sver=complex(rcsver*cos(phase),rcsver*sin(phase));
	shor=complex(rcshor*cos(phase),rcshor*sin(phase));