function [FFv,a_s]=FormFunctionAtten(k,a,rho_s)
%
%  function [FFv,a_s] = FormFunctionAtten(k,a,[rho_s])
%
%  Input: k = acoustic wavenumber (=2*pi/lambda, lambda=c/f) in 1/m
%         a = sediment radius (=d/2) in m
%         rho_s = sediment density [optional, default=2650]
%
%  Output: FFV = form function volume 
%          a_s = sediment attenuation constant
%  
% calculates the Form Function (FFv) and the sediment attenuation constant (a_s)
% for sediment of radius a with sound of acoustic wavenumber k. It 
% assumes a sediment density 2650. From Thorne, P.D.& D.M. Hanes, CSR, 2003
% equations 9, 10a and 10b.
%
%  The sediment attenuation is smaller with this method espaecially
%  at the very small values (Vincent's method, shows a secondary peak
%  at radius 0.0002 m)
%
%   George Voulgaris, USC, Feb 15th, 2004
%
if nargin<3;
    rho_s=2650;
end
k=k(:);
a=a(:);
Na=length(a);
Nk=length(k);
%
kf=1.1;
ka=0.18;
%
for i=1:Na;
    x=k.*a(i);
    chi1=1.1*(4/3)*ka*x.^4;
    chi2=1 + 1.3*x.^2 + (4/3)*ka*x.^4;
    chi=chi1./chi2;
    a_s(:,i)=3*chi/(4*a(i)*rho_s);   % sediment attenuation constant
    %
    c1=1-0.25*exp(-((x-1.4)/0.5).^2);
    c2=1+0.37*exp(-((x-2.8)/2.2).^2);
    Co=1.1*c1.*c2;
    FFv(:,i)=Co.*(kf*x.^2)./(1+kf*x.^2); % form function
end