function [eflux_lin,eflux_nonlin,invfr_simple,mflux_east,mflux_north] = ...
          get_flux_pgoff_spc_v2(rho_0,u0,v0,f0,N0,aa,ksq,k_2d,l_2d,kiso,...
          flag_spc, pgoff_2d_q1,pgoff_2d_q2)
%--------------------------------------------------------------------------
% History:
%---------
% "Regularized" i.e. the get_flux_pgoff_v4 adds 1mm/s to U0 in
% 1/Fr def.

% Nov 21, confirmed that quadrants 1 and 3, and 2 and 4 were identical for
% all 3 integrands: eflux_integrand, mflux_integrand. Removed the loop over
% these, and confirmed that: E_lin_map, Taux_map, Tauy_map, InvFr_map
%-----------
% Output:
% eflux_lin     generation rate of Lee waves found using simple linear
%               theory. Units = [W].
%------------------------------------------------
% u0 = U0*sind(Theta0);   % zonal velocity
% v0 = U0*cosd(Theta0);   % meridional velocity
%--------------------------------------------------------------------------
% phi = atan2(v0,u0)*180/pi;    % angle in degrees anticlockwise must rotate x-y axis to align (u0,v0) along x-axis

noctave = length(kiso)-1;

% initialize storage 
eflux_lin = NaN(1,noctave);

% FIRST QUADRANT

s_2d = u0*k_2d + v0*l_2d;

if flag_spc == 1
        ksm2 = ksq;           % total horizontal wavenumber squared
    elseif flag_spc == 2
        ksm2 = (N0^2 - s_2d.^2)./(s_2d.^2 -f0^2);      % m^2/K^2 = tan^2(wave angle to horizontal)
    elseif flag_spc == 3
        ksm2 = ksq.*(N0^2 - s_2d.^2) ./ (s_2d.^2 - f0^2);   % m^2
    else
        error('define flag_spc 1 , 2, or 3')
end




hms         = 0;
mflux_east  = 0;
mflux_north = 0;

for octave = 1:noctave
    
    idx  = find(abs(s_2d)>=f0 & abs(s_2d)<=N0 & ...
                ksm2>=kiso(octave)^2 & ksm2<kiso(octave+1)^2);
    
    % just at relevant points
    wsq_vec   = s_2d(idx).^2;
    a_vec     = sqrt(wsq_vec-f0^2);
    b_vec     = sqrt(N0^2-wsq_vec);
    k_vec     = sqrt(k_2d(idx).^2+l_2d(idx).^2);
    c_vec     = 1./(4 .* k_vec *pi^2);
    pgoff     = pgoff_2d_q1(idx);
    d_vec     = rho_0 * a_vec .*b_vec .* pgoff .* c_vec;
    eflux_vec = abs(s_2d(idx)) .* d_vec;
    mflux_vec_x = k_2d(idx).*sign(s_2d(idx)) .*d_vec;
    mflux_vec_y = l_2d(idx).*sign(s_2d(idx)) .*d_vec;
    aa_vec      = aa(idx);
    
    eflux_lin(octave) = 2*(eflux_vec'   * aa_vec);
    mflux_east        = mflux_east + 2*(mflux_vec_x' * aa_vec);
    mflux_north       = mflux_north + 2*(mflux_vec_y' * aa_vec);
    hms               = hms + 2*(pgoff' * aa_vec);
end
% hrms_u =   sqrt( sum(sum(aa(ip,:).*squeeze(sum(p_2d(ip,:,1:2),3)),2),1) ...
%           + sum(sum(aa(in,:).*squeeze(sum(p_2d(in,:,3:4),3)),2),1) )/2/pi;

% SECOND QUADRANT
% k --> -k
% just at relevant points, all same formulas as above, except:
% pgoff quad different, and 
% adding results to previous
% idx is different because of the change in s_2d
s_2d = -u0*k_2d + v0*l_2d;   % % k --> -k

for octave = 1:noctave
    
    idx  = find(abs(s_2d)>=f0 & abs(s_2d)<=N0 & ...
                ksm2>=kiso(octave)^2 & ksm2<kiso(octave+1)^2);
    %idx  = find(abs(s_2d)>=f0 & abs(s_2d)<=N0);
    
    % just at relevant points
    wsq_vec   = s_2d(idx).^2;
    a_vec     = sqrt(wsq_vec-f0^2);
    b_vec     = sqrt(N0^2-wsq_vec);
    k_vec     = sqrt(k_2d(idx).^2+l_2d(idx).^2);
    c_vec     = 1./(4 .* k_vec *pi^2);
    pgoff     = pgoff_2d_q2(idx);        % second quadrant , where k --> -k
    d_vec     = rho_0 * a_vec .*b_vec .* pgoff .* c_vec;
    eflux_vec = abs(s_2d(idx)) .* d_vec;
    mflux_vec_x = -k_2d(idx).*sign(s_2d(idx)) .*d_vec;   % k --> -k
    mflux_vec_y =  l_2d(idx).*sign(s_2d(idx)) .*d_vec;
    aa_vec      = aa(idx);
    
    eflux_lin(octave)  = eflux_lin(octave)   + 2*(eflux_vec'   * aa_vec); % SUM OVER ALL QUADRANTS  , [W/m^2]
    mflux_east         = mflux_east  + 2*(mflux_vec_x' * aa_vec);
    mflux_north        = mflux_north + 2*(mflux_vec_y' * aa_vec);
    
    hms                = hms + 2*(pgoff' * aa_vec);
end
% hms is the integral over spectrum, 
hrms_u       = sqrt(hms)/(2*pi);
% Inverse Froude:
U0           = sqrt(u0^2+v0^2);    % current speed
invfr_simple = hrms_u *N0/(U0+0.001);

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

%------------------------------------------------
% Compare this to the definition of "Long Number" Lo = NH/U , Aguilar &
% Sutherland (2006, Eqn 4).

%------------------------------------------------
if invfr_simple > 0.7
    eflux_nonlin = eflux_lin * (0.7./invfr_simple).^2;
else 
    eflux_nonlin = eflux_lin;
end


return