% Backscattering coefficient as a function of ka for
% a cylinder.
%
% Code developoed by Jonathan Mamou with the help of 
% Dr Michael Oelze at the
%
% Bioacoustics Research Lab at the University of Illinois
% at Urbana-Champaign
% 
% in March 2003.
%
% This code also displays the usual Form Factor functions.
% The code is based on Faran's theory.
%
% Please feel free to contact the author at jonathan@brl.uiuc.edu

clear all
ka = (0:499)*.01;

glass=0;
poly=1;
fat=2;
tung=3;
Nume=10;


material=glass;
shear=1

opt=0;



if material==glass
  c1=5570;
  pwire=2.38;
  Poisson=0.21;
elseif material==poly
  c1=1350;
  pwire=1.06;
  Poisson=0.35;
elseif material==fat 
  c1=1460;
  pwire=0.94;
  Poisson=0.4993;
elseif material==tung 
  c1=5200;
  pwire=19.3;
  Poisson=0.284;
end

if shear==0
  Poisson=0.5-10^(-13);
end



c2 = sqrt((.5*c1^2-Poisson*c1^2)/(1-Poisson));

ph2o = 1.0;

c = 1540;
im = sqrt(-1);
a = 25e-6;
x1 = ka.*c/c1;
x2 = ka.*c/c2;
eta = zeros(500,11);
eta_el = zeros(500,11);


for i = 0:10
   J = bessel(i,ka);
   N=  bessely(i,ka);
   
   Jx1= bessel(i,x1);
   Jx2= bessel(i,x2);
   
 
   
   xJp = ka.*bessel(i-1,ka)-i*J;

   xNp = ka.*bessely(i-1,ka)-i*N;
	
   x1Jp= x1.*bessel(i-1,x1)-i*Jx1;
   x2Jp= x2.*bessel(i-1,x2)-i*Jx2;
   

 
   del = -J./N;
   alpha = -xJp./J;
   beta = -xNp./N;
   

   alphax1 = -x1Jp./Jx1;
 
   
   
 
   alphax2 = -x2Jp./Jx2;
  
   

   Num1= alphax1./(alphax1+1);
   Num2=(i^2)./(alphax2+i*i-x2.^2/2);
   
  Num= Num1-Num2;
 
  Den1=(alphax1+i^2-x2.^2/2)./(alphax1+1);
  Den2=(i^2*(alphax2+1))./(alphax2+i^2-x2.^2/2);
  
  
   Denom= Den1-Den2;
   
   

   etc = (-x2.^2/2).*Num./Denom;
   
  
   
   phi=-etc.*ph2o/pwire;
   
   eta_el(:,i+1)= atan(del'.*(alpha'+phi')./(phi'+beta'));
   
   eta(:,i+1) = atan((del'.*alpha')./beta');


end
sum = zeros(500,1);
%sum = sum+sin(eta(:,1)).^2;

sum_el = zeros(500,1);
%sum_el = sum_el+sin(eta_el(:,1)).^2;

i=0;
sum = sum+sin(eta(:,i+1)).*exp(-im.*eta(:,i+1)).*(-1)^i;
sum_el = sum_el+sin(eta_el(:,i+1)).*exp(-im.*eta_el(:,i+1)).*(-1)^i;  

for i = 1:Nume
  sum = sum+2*sin(eta(:,i+1)).*exp(-im.*eta(:,i+1)).*(-1)^i;
  sum_el = sum_el+2*sin(eta_el(:,i+1)).*exp(-im.*eta_el(:,i+1)).*(-1)^i;  
end

sum = 2*abs(sum).^2./(ka)';
sum_el=2*abs(sum_el).^2./(ka)';


figure
plot(ka,sum);
hold on
plot(ka,sum_el,'--');
xlabel('ka');

sum = sum./ka.^3';
sum_el=sum_el./ka.^3';
sum=sum/(sum(2));
sum_el=sum_el/sum_el(2);


figure
plot(ka,sum);
hold on
plot(ka,sum_el,'--');
xlabel('ka');

% Gaussian Match:



Min=10^100;
Min2=Min;
Min3=Min;
Min4=Min;


ka_min_index=min(find(ka>=0.5))
%ka_max_index=min(find(ka>1.5))
ka_max_index=min(find(ka>=1.2))

if opt==1
for dd=1:200
  d=dd/100;
  
  G=exp(-(d*ka).^2);
%  E=1./(1+(d*2*ka).^2);
E=(bessel(1,d*2*ka)./(d*ka)).^2;

F=bessel(0,d*2*ka).^2;

H=(1./(1+(d*ka).^2)).^2;

  sumg=sum_el(ka_min_index:ka_max_index);
  Gg=G(ka_min_index:ka_max_index);
  Eg=E(ka_min_index:ka_max_index);
  Fg=F(ka_min_index:ka_max_index);
  Hg=H(ka_min_index:ka_max_index);
 
  Err=mean((sumg'-Gg).^2);
  Err2=mean((sumg'-Eg).^2);
  Err3=mean((sumg'-Fg).^2);
  Err4=mean((sumg'-Hg).^2);
  
  if Err<Min;
    d_star=d;
    Min=Err;
  end
  
  if Err2<Min2;
    d_star2=d;
    Min2=Err2;
  end
  
  if Err3<Min3;
    d_star3=d;
    Min3=Err3;
  end
  
  if Err4<Min4;
    d_star4=d;
    Min4=Err4;
  end
  
  
  
end

d_star
d_star2
d_star3
d_star4

Min
Min2
Min3
Min4

else
  d_star=1;
  d_star2=1;
  d_star3=1;
  d_star4=1;
end

xlabel('ka');
plot(ka,exp(-(d_star*ka).^2),'r')
%plot(ka,1./(1+(d_star2*2*ka).^2),'k')
plot(ka,(bessel(1,d_star2*2*ka)./(d_star2*ka)).^2,'k')
plot(ka,bessel(0,2*d_star3*ka).^2,'g')

%plot(ka,1./(1+(d_star4*ka).^2).^2,'c')