function [p,R] = anderson(k,a,c1,c0,rho1,rho0,theta,n,r); 
% [p,R] = anderson(k,a,c1,c0,rho1,rho0,theta,n,r)
%
% 2004.07
% Ralph E. Baddour
% rbaddour@uhnres.utoronto.ca
% Dept. of Medical Biophysics, University of Toronto, Canada
%
% Acoustic scattering solution for a fluid sphere (i.e., does not support
% shear) in a medium that does not support shear (i.e., a fluid as well).
% (based on the solution of Anderson, V. C., JASA 22:426-431, 1950)
%
% Calculates the normalised scattered steady state pressure (normalised to the
% incident pressure, i.e., incident pressure amplitude, Pin = 1) that arrives
% at an observation point r [mm] away from the sphere.
%
% Also calculates reflectivity factor (as defined in Anderson, 1950):
% Rtheta = |p| / |pgeom|, where |pgeom| = Pin[(pi*a^2)/(4*pi*r^2)]^0.5 = Pin*a/(2*r)
% (Pin is considered to be scattered uniformly through a solid angle 4pi sterad)
%
% INPUT
% Incident wave wavenumber                      (k)     (may be a vector) 
% Sphere radius                                 (a)     [mm] 
% Longitudinal speed of sound in sphere fluid   (c1)    [mm/s]
% Longitudinal speed of sound in medium         (c0)    [mm/s]
% Density of the sphere fluid                   (rho1)  [gm/cm^3]
% Density of surrounding medium                 (rho0)  [gm/cm^3]
% Angle, in radians, of observation point       (theta) (0 is backscatter, unlike Faran!)
% Order of the solution (number of terms)       (n) 
% Distance of observation point                 (r)     [mm] 
%                                                       
% OUTPUT
% normalised scattered pressure (complex)       (p)
% reflectivity factor                           (R)

warning off;
i = sqrt(-1);
p = zeros(size(k));
R = zeros(size(k));

g = rho1/rho0;
h = c1/c0;
k1 = k/h;
% omega = k*c0;

for m = 0:n 
    Cm_top1 = (alpham(m,k1*a)./alpham(m,k*a)) .* (nm(m,k*a)./jm(m,k1*a));
    Cm_top2 = (betam(m,k*a)./alpham(m,k*a)) .* (g*h);
    Cm_top = Cm_top1 - Cm_top2;
    Cm_bottom = ( (alpham(m,k1*a)./alpham(m,k*a)) .* (jm(m,k*a)./jm(m,k1*a)) ) - (g*h);
    Cm = Cm_top ./ Cm_bottom;
    
    Pm = legendre(m,cos(theta));
    Pm = Pm(1,1);
    
    p = p + ( (((-i)^m)*((2*m)+1)./(1+(i*Cm))) .* (Pm.*(jm(m,k*r) + (i*nm(m,k*r)))) );
    R = R + ( Pm*((-1)^m)*((2*m)+1)./(1+(i*Cm)) );
end
R = (2./(k*a)) .* abs(R);
warning on;

%----------------------------
%     REQUIRED FUNCTIONS
%----------------------------

function y = alpham(m,x) 
% Function to return alpham, as defined in Anderson 1950.
%
% INPUT
% Function parameter        (x)     (may be a vector) 
% Order of the term         (m)
y = (m * jm(m-1,x)) - ((m+1) * jm(m+1,x));

function y = betam(m,x) 
% Function to return betam, as defined in Anderson 1950.
%
% INPUT
% Function parameter        (x)     (may be a vector) 
% Order of the term         (m)
y = (m * nm(m-1,x)) - ((m+1) * nm(m+1,x));

function y = jm(m,x) 
% Function to return the spherical Bessel function of the first kind. 
%
% INPUT
% Function parameter        (x)     (may be a vector) 
% Order of the solution     (m)
y = sqrt(pi ./ (2*x)) .* besselj(m+0.5, x); 
  
function y = nm(m,x) 
% Function to return the spherical Bessel function of the second kind. 
% (also called the Neumann function, denoted as nm in Anderson 1950)
%
% INPUT
% Function parameter        (x)     (may be a vector) 
% Order of the solution     (m)
y = sqrt(pi ./ (2*x)) .* bessely(m+0.5, x); 
