function bscat_ave = copepod_DWBA_bscat_analy(freq, Lmf, Lml, Wmf, Wml, th_p, L_p, cw, g_f, g_lip, h_f, h_lip)
% compute the averaged analytical differential backscattering cross section as a function of
% angle of orientation and length (keeping the aspect ratio the same for different lengths)
% The shapes of flesh and lipid are concentric prolate spheroids of the same length
% 
% INPUTS:
%        freq = frequency in kHz
%         Lmf = mean length of the copepods in mm
%         Lml = mean length of the lipid inside the body in mm (assumed the same as its body flesh length)
%         Wmf = mean width of the copepods in mm
%         Wml = mean width of the lipid inside the body in mm 
%        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_f = density contrast of the copepod flesh
%       g_lip = density contrast of the copepod lipid
%         h_f = sound speed contrast of the copepod flesh
%       h_lip = sound speed contrast of the copepod lipid

%% 
% kw = 2*pi*freq/cw;
% ks = -2*kw;  % vec(ki) - vec(ks) = -2 vec(ks)
ar_f = Lmf/Wmf;  % aspect ratio of flesh outer boundary
ar_lip = Lml/Wml; % aspect ratio for both inner flesh and lipid

%% distorted wave numbers in flesh and lipid
% k_f = ks/h_f;
% k_lip = ks/h_lip;

a_f = Wmf/2;
a_lip = Wml/2;
bscat_f_outer = bscat_ave_ellipsoid(freq, cw, a_f, 1, ar_f, th_p, L_p, g_f, h_f);
bscat_f_inner = bscat_ave_ellipsoid(freq, cw, a_lip, 1,ar_lip, th_p, L_p, g_f, h_f);
bscat_lip = bscat_ave_ellipsoid(freq, cw, a_lip, 1, ar_lip, th_p, L_p, g_lip, h_lip);

bscat_ave = (bscat_f_outer - bscat_f_inner + bscat_lip);





