function    f_bscat_ave = bscat_ave_ellipsoid(freq, cw, a, e_ba, e_ca, th_p, L_p, g, h)
% function bscat_ave_ellipsoid returns the backscattering by a homogeneous ellipsoid using  
% Distorted Wave Born Approximation (DWBA) as afunction of frequency, and the scattering is 
% averaged over angle of orientation (th_pdf) and length (Ln)
% INPUTS:
%    freq = frequency (kHz)
%      cw = sound speed in water
%       a = semi-minor axis in x-axis (mm)
%    e_ba = aspect ratio of the two semi-minor axes  = b/a
%    e_ca = aspect ratio of semi-major to semi-minor axis (x-axis) = c/a
%    th_p = angle parameter: 
%           if th_p is a scaler or an array, it implies a single orientation, 
%               or uniform average over a given range of angle of oreintation 
%           if th_p is struct:  th_p.thm and th_p.th_std is the mean and
%               standard deviation of a Gaussian distribution N(thm,
%               th_std)and the average is taken ove thm +- 3*th_std
%    L_p = Length parameter: 
%           if L_p is a scaler or an array, it implies a single normalized length (unity), 
%               or uniform average over a given range of normalized length  
%           if L_p is struct:  L_p.L_std is the standard deviation of a 
%               Gaussian distribution N(Lm, L_p.L_std) and the average is 
%               taken ove Lm +- 3*L_std (note Lm = 2*a*e_ca)
%       g = density contrast = rho_scat/rho_water
%       h = sound speed contrast = c_scat/cw
% OUTPUT:
%   f_bscat_ave = averged backscattering amplitude based on DWBA
%   
% 
%% Dezhang Chu, NOAA Fisheries, NWFSC
%% Created on March 5th, 2017
%% Revised on 5/23/2017 add dth and dL to the coresponding PDFs

nf = length(freq);
freq = reshape(freq, 1, nf);   % a row vector
kw = 2*pi*freq*1000/cw;     % kHz --> Hz
a = a*1e-3;                 % mm --> m
b = a*e_ba;
c = a*e_ca;
Cb = (1 - g*h^2)/(g*h^2) - (g - 1)/g;

if isstruct(th_p)
    nth = 30;
    th = linspace(th_p.thm - 3*th_p.th_std, th_p.thm + 3*th_p.th_std, nth);
    th_pdf = pdf('norm', th, th_p.thm, th_p.th_std);
else
    nth = 100;
    th_pdf = ones(nth,1)/nth;
end
th = th(:);   % column vector
th_pdf = th_pdf(:);         % column vector
dth = th(2)-th(1);

if isstruct(L_p)
    nL = max(10,5*max(round(kw*a*e_ca)));
    L = linspace(1 - 3*L_p.L_std, 1 + 3*L_p.L_std, nL);
    L_pdf = pdf('norm', L, 1, L_p.L_std);
    L_pdf = L_pdf(:);          % a column vector
else
    L_pdf = ones(nL,1)/nL;     % a column vector
end
L = reshape(L, 1, length(L));  % a row vector
dL = L(2)-L(1);

Th_pdf = repmat(th_pdf, 1, nL);  % 2D pdf for angle of orientation
coef0 = Cb*a*b*c*L.^3;
arg0 = sqrt(sind(th).^2 + e_ca^2*cosd(th).^2)*(2*kw/h);  % 
Coef = repmat(coef0,nth,1);     % 2D coefficient matrix (th, L)

for i = 1:nf
     arg_f = arg0(:,i)*a*L;
     j1 = sqrt(pi./(2*arg_f)).*besselj(1.5,arg_f);
     j1_arg = j1./arg_f;
     %% backscattering amplitude
     f_bscat_i  = kw(i)^2*Coef.*j1_arg;    
     if nf > 1
         % average over angle of orientation
         f_bscat_ave_th = sum(f_bscat_i.*Th_pdf,1)*dth;
         % average over length
         f_bscat_ave(i) = f_bscat_ave_th*L_pdf*dL;
     end
end
if nf == 1
    f_bscat_ave = f_bscat_i;
end

