% Backscattering coefficient as a function of ka for
% a sphere.
%
% 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
nnn=500;
ka = (0:nnn-1)*.01;

Nume=20;
eta = zeros(nnn,Nume+1);
eta_el = zeros(nnn,Nume+1);

glass=0;
poly=1;
fat=2;
tung=3;
res=4;


material=glass;
shear=1;

opt=0;

if material==glass 
  c1=5570;
  pwire=2.38;
  Poisson=0.21;

elseif material==poly
  c1=2000;
  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;
elseif material==res
  c1=1500;
  pwire=0.90;
  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;



opt=0;

for i = 0:Nume
  
  
   J = sbessel(i,ka);
   N = sbessely(i,ka);
   
   Jx1= sbessel(i,x1);
   Jx2= sbessel(i,x2);
   
 
   xJp=-0.5*J+0.5*ka.*(sbessel(i-1,ka)-sbessel(i+1,ka));
   xNp=-0.5*N+0.5*ka.*(sbessely(i-1,ka)-sbessely(i+1,ka));

   x1Jp=-0.5*Jx1+0.5*x1.*(sbessel(i-1,x1)-sbessel(i+1,x1));
 
     x2Jp=-0.5*Jx2+0.5*x2.*(sbessel(i-1,x2)-sbessel(i+1,x2));
  
 
   del = -J./N;
   alpha = -xJp./J;
   beta = -xNp./N;
   
   alphax1 = -x1Jp./Jx1;
   alphax2 = -x2Jp./Jx2;
 
   
   
   Num= alphax1./(alphax1+1)-(i*i+i)./(alphax2+i*i+i-1-x2.*x2/2);
   Denom= (2*alphax1+i*i+i-x2.*x2/2)./(alphax1+1)-((i*i+i)*(alphax2+1))./(alphax2+i*i+i-1-x2.*x2/2);
 
   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(nnn,1);
sum_el = zeros(nnn,1);

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


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

sume=((sum)./ka.^2').^2;
sumeN=sume/sume(2);

sume_el=((sum_el)./ka.^2').^2;
sume_elN=sume_el/sume_el(2);



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



figure
xlabel('ka')
plot(ka,sumeN);
hold on
plot(ka,sume_elN,'--');
 

% Gaussian Match:

if opt==1;

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))


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

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

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

  sumg=sum_elN(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
end

%d_star=sqrt(0.827)
%d_star2=1
%d_star3=1
%d_star4

d_star=0.8271
d_star2=1
d_star3=1
d_star4=1.2114

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