function [fe,fr] = faran(k,a,c1,c3,rho1,rho3,sigma,theta,n); 
% MATLAB function to solve the far-field form function for 
% the echo from an elastic/inelastic sphere in a fluid. 
% Martin Anderson 
% Created: 3-25-96      Last modified: 10-29-99 
% 
% Accepts the fluid isonification wavenumber    (k) (may be a vector) 
% the sphere radius                             (a)     [mm] 
% the longitudinal velocities in the sphere     (c1)    [mm/s] 
% and the fluid                                 (c3)    [mm/s] 
% the densities of the sphere material          (rho1)  [gm/cm^3] 
% and the fluid                                 (rho3)  [gm/cm^3] 
% the Poisson ratio of the sphere material      (sigma) 
% the polar angle in radians                    (theta) 
% and the required order of the solution.       (n) 
% Calculates the far-field form function for 
% the echo from an elastic sphere in a fluid as defined 
% in Faran, J. J., JASA 23:405-418, 1951, as well 
% as the same function for an inelastic sphere, accomplished 
% as described in Hickling, R. H., JASA 34:1582-1592, 1962, 
% by setting Fn = 0. 
% 
% Syntax: 
% [fe,fr] = f_faran(k,a,c1,c3,rho1,rho3,sigma,theta,n); 
% Dep: enun.m rnun.m 
x3 = k * a; 
x1 = x3 * c3 / c1; 
i = sqrt(-1); 

fe = zeros(size(k)); 
fr = zeros(size(k)); 

for ord = 0:n 
        enu = enun(ord,x1,x3,rho1,rho3,sigma); 
        rnu = rnun(ord,x1,x3,rho1,rho3,sigma); 
        Pn = legendre(ord,cos(theta)); 
        Pn = Pn(1,1); 

% The following formula solves for the complex pressure in the far-field: 
% (After Eq. 10 in Hickling, J., JASA 34:1582-1592, 1962) 
        einc = 2 * (-1)^n * (2 * ord + 1) .* sin(enu) .* exp(-i * enu) * Pn ./ (a*k); 
        rinc = 2 * (-1)^n * (2 * ord + 1) .* sin(rnu) .* exp(-i * rnu) * Pn ./ (a*k); 

%       disp([num2str(round(ord*100/n)) '% done...']) 
        fe = fe + einc; 
        fr = fr + rinc; 
end 

%----------------------------------------------------------------- 

function nu = enun(n,x1,x3,rho1,rho3,sigma) 
% MATLAB function to find the nu-sub-n angle required for 
% calculating the far_field pressure of 
% the echo from an elastic sphere in a fluid. 
% Martin Anderson 
% Created: 3-25-96      Last modified: 6-26-97 
% Accepts the required order of the solution (n), 
% the longitudinal wavennumbers in the sphere (x1), 
% the longitudinal wavennumbers in the fluid (x3), 
% the density of the sphere (rho1), 
% the density of the fluid (rho3), 
% and the Poisson's ratio of the sphere material (sigma). 
% Calculates the nu-sub-n angle described in 
% equation 30 of Faran, J. J., "Sound scattering by 
% solid cylincers and spheres," JASA 23:405-418, 1951. 
% 
% Syntax: 
% nu = enun(n,x1,x3,rho1,rho3,sigma) 

% Set Fn = 0 to solve for inelastic sphere: 
% Fn = 0; 

% Calculate Fn to solve for elastic sphere: 
Fn = -(rho3/rho1) * gamman(n,x1,sigma); 

tnu = -(jn(n,x3) .* Fn - x3 .* djn(n,x3)); 
tnu = tnu ./ (yn(n,x3) .*Fn - x3 .* dyn(n,x3)); 
nu = atan(tnu); 
  

%----------------------------------------------------------------- 

function nu = rnun(n,x1,x3,rho1,rho3,sigma) 
% MATLAB function to find the nu-sub-n angle required for 
% calculating the far_field pressure of 
% the echo from a rigid sphere in a fluid. 
% Martin Anderson 
% Created: 3-25-96      Last modified: 6-27-97 
% Accepts the required order of the solution (n), 
% the longitudinal wavennumbers in the sphere (x1), 
% the longitudinal wavennumbers in the fluid (x3), 
% the density of the sphere (rho1), 
% the density of the fluid (rho3), 
% and the Poisson's ratio of the sphere material (sigma). 
% Calculates the nu-sub-n angle described in 
% equation 30 of Faran, J. J., "Sound scattering by 
% solid cylincers and spheres," JASA 23:405-418, 1951. 
% 
% Syntax: 
% nu = rnun(n,x1,x3,rho1,rho3,sigma) 

% Set Fn = 0 to solve for inelastic sphere: 
Fn = 0; 

% Calculate Fn to solve for elastic sphere: 
% Fn = -(rho3/rho1) * gamman(n,x1,sigma); 

tnu = -(jn(n,x3) .* Fn - x3 .* djn(n,x3)); 
tnu = tnu ./ (yn(n,x3) .*Fn - x3 .* dyn(n,x3)); 
nu = atan(tnu); 

%----------------------------------------------------------------- 

function res = gamman(n,x1,sigma) 
% MATLAB function to find the gamma-sub-n factor required 
% for calculating the far_field pressure of 
% the echo from an elastic sphere in a fluid. 
% Martin Anderson 
% Created: 3-25-96      Last modified: 3-27-96 
% Accepts the required order of the solution (n), 
% the longitudinal wavennumbers in the sphere (x1), 
% and the Poisson's ratio of the sphere material (sigma). 
% Calculates the gamma-sub-n factor described in 
% equation 30 of Faran, J. J., "Sound scattering by 
% solid cylincers and spheres," JASA 23:405-418, 1951. 
% 
% Syntax: 
% res = gamman(n,x1,sigma); 

x2  = sqrt((sigma-1)/(sigma-0.5)) * x1; 

tax1 = tan(alphn(n,x1)); 
tax2 = tan(alphn(n,x2)); 

nterm = n^2 + n; 
x2term = (nterm -1 -0.5*x2.^2 +tax2); 

num1 = tax1./(tax1+1); 
num2 = nterm./x2term; 
num = -x2.^2 .* (num1 - num2); 

denom1 = (nterm -0.5*x2.^2 + 2*tax1) ./ (tax1+1); 
denom2 = (nterm * (tax2+1)) ./ x2term; 
denom = 2 * (denom1 - denom2); 
res = num ./ denom; 

%----------------------------------------------------------------- 

function y = djn(n,x) 
% MATLAB function to solve the first derivative of 
% the spherical bessel function of the first kind. 
% Martin Anderson 
% Created: 3-25-96      Last modified: 3-26-96 
% Accepts the argument(s)                       (x) (may be a vector) 
% and the required order of the solution.       (n) 
% Returns the first derivative at the argument(s) 
% using the central difference formula with Richardson 
% extrapolation. 
% 
% Syntax: 
%  y = djn(n,x); 

h = .001; 
h2 = 2 * h; 
len = length(x); 
y = zeros(size(x)); 
for c = 1:len 
        ord = x(c)-h2:h:x(c)+h2; 
        f = sqrt(pi ./ (2 * ord)) .* besselj(n+0.5, ord); 
        F = (f(4) - f(2)) / h2; 
        F2 = (f(5) - f(1))/ (2 * h2); 
        y(c) = F + (F - F2)/3; 
end 
  

%----------------------------------------------------------------- 

function y = dyn(n,x) 
% MATLAB function to solve the first derivative of 
% the spherical bessel function of the second kind. 
% Martin Anderson 
% Created: 3-25-96      Last modified: 3-26-96 
% Accepts the argument(s)                       (x) (may be a vector) 
% and the required order of the solution.       (n) 
% Returns the first derivative at the argument(s) 
% using the central difference formula with Richardson 
% extrapolation. 
% 
% Syntax: 
%  y = dyn(n,x); 

h = .001; 
h2 = 2 * h; 
len = length(x); 
y = zeros(size(x)); 
for c = 1:len 
        ord = x(c)-h2:h:x(c)+h2; 
        f = sqrt(pi ./ (2 * ord)) .* bessely(n+0.5, ord); 
        F = (f(4) - f(2)) / h2; 
        F2 = (f(5) - f(1))/ (2 * h2); 
        y(c) = F + (F - F2)/3; 
end 
  

%----------------------------------------------------------------- 

function y = alphn(n,x) 
% MATLAB function to solve the intermediate angle alpha-sub-n 
% Martin Anderson 
% Created: 3-25-96      Last modified: 3-26-96 
% Accepts the wavenumbers * a                   (x) (may be a vector) 
% and the required order of the solution.       (n) 
% 
% Syntax: 
%  y = alphn(n,x); 

y = atan(( - x .*  djn(n,x)) ./ jn(n,x)); 

%----------------------------------------------------------------- 

function y = jn(n,x) 
% MATLAB function to solve the spherical bessel function 
% of the first kind. 
% Martin Anderson 
% Created: 3-25-96      Last modified: 3-26-96 
% Accepts the argument(s)                       (x) (may be a vector) 
% and the required order of the solution.       (n) 
% 
% Syntax: 
%  y = jn(n,x); 

y = sqrt(pi ./ (2 * x)) .* besselj(n+0.5, x); 
  

%----------------------------------------------------------------- 

function y = yn(n,x) 
% MATLAB function to solve the spherical bessel function 
% of the second kind. 
% Martin Anderson 
% Created: 3-25-96      Last modified: 3-26-96 
% Accepts the argument(s)                       (x) (may be a vector) 
% and the required order of the solution.       (n) 
% 
% Syntax: 
%  y = yn(n,x); 

y = sqrt(pi ./ (2 * x)) .* bessely(n+0.5, x); 
  

%----------------------------------------------------------------- 
  
  
