function [flux,wn_flux,T,wn_T] = vel2flux_v2(dx,dy,u_in,v_in)
% calculates the spectral ke flux of the given velocity field
%-------
% INPUT:
%-------
% dx 	is the grid spacing in x direction
% dy 	is the grid spacing in y direction
% u_in	zonal velocity
% v_in	meridional velocity
%-------
% OUTPUT:
%-------
% flux 		is the spectral KE flux
% wn_flux 	is the wavenumbers corresponding to flux
% T 		is the spectral transfer of KE into wavenumber, wn_T
%		(i.e. it's the divergence of flux)
% wn_T	 	is the wavenumbers corresponding to T
%-------
% History:
% Written by Rob Scott
% Built from vel2flux
% Modified Jan 20, 2008:
% 	allows for non-square domain
%-------
% Notes:
% Where do wavenumbers come from?
% The DFT projects the field onto
%    exp(-j 2 pi (k-1) (n-1)/N), 1<= k <= N, 1<= j <= N
% To find wn's in rad/m, multiply top and bottom by dx [m],
%    exp(-j 2 pi (k-1) (n-1) dx/ N dx), note (n-1)dx = x, and Length = N dx.
%  = exp( -j wn x), so 
%    wn = 2 pi (k-1)/ N dx, with 1<= k <= N
% But note the wn is ambiguous because of the shift rule:
%    k = N-p  --> k = -p
%    exp(-j 2 pi ((N-p)-1) (n-1)/N) =  exp(-j 2 pi (-p-1) (n-1)/N)
% It's convenient for interpretation to represent the wn as the smallest number
% possible. This is done by replacing k=(N-p)>N/2 with -p when N is even
% or k=(N-p)>(N-1)/2 with -p when N is odd
%------------------------------------------------------------------------------


nx = size(u_in,2);
ny = size(u_in,1);

if size(v_in,1) ~= ny | size(v_in,2)~= nx
   error('ERROR: size u, v must be same')
end

% size of domain:
Lx = dx * nx;
Ly = dy * ny;


% here's the heart of the directional transfer calculation:

% check if nx is odd or even
if round(nx/2)*2 == nx
   Nx = nx/2;
   wn_x = 2*pi*[0:Nx,-(Nx-1):-1]' / Lx;
else
   Nx = (nx-1)/2;
   wn_x = 2*pi*[0:Nx,-Nx:-1]' / Lx;
end

% check if ny is odd or even
if round(ny/2)*2 == ny
   Ny = ny/2;
   wn_y = 2*pi*[0:Ny,-(Ny-1):-1]' / Ly;
else
   Ny = (ny-1)/2;
   wn_y = 2*pi*[0:Ny,-Ny:-1]' / Ly;
end


[kx,ky] = meshgrid(wn_x,wn_y);
wvsq = kx.*kx+ky.*ky;

%  u = real(ifft2(-j*ky.*fft2(psi)));
%  v = real(ifft2(j*kx.*fft2(psi)));
%  q = real(ifft2((-wvsq).*fft2(psi)));

  ddx_u = real(ifft2(j*kx.*fft2(u_in)));
  ddy_u = real(ifft2(j*ky.*fft2(u_in)));

  ddx_v = real(ifft2(j*kx.*fft2(v_in)));
  ddy_v = real(ifft2(j*ky.*fft2(v_in)));

  adv_u = u_in.*ddx_u + v_in.*ddy_u;
  adv_v = u_in.*ddx_v + v_in.*ddy_v;

  Tkxky = real( - conj(fft2(u_in)).*fft2(adv_u) ...
                - conj(fft2(v_in)).*fft2(adv_v)  )/(nx^2 * ny^2);

% convert this to an isotropic flux

wn_max = sqrt(max(wn_y)^2 + max(wn_x)^2)/sqrt(2);
% somewhat arbitrary increment in wavenumber, used for derivative below:
dk = min((2*pi/Lx),(2*pi/Ly));

%step_max = ceil(wn_max/dk);

step_max = min(Nx,Ny)+1;

% get the "cumulative" KE spectral flux 
% (analogous to cumulative density function, CDF)
for step = 1:step_max
   k = dk*step;
   index=find(wvsq >= k^2);
   flux(step) = sum(sum(Tkxky(index)));
   wn_step(step) = k;
   clear index
end

% make a column vector:
if size(flux,1) < size(flux,2)
   flux = flux';
end


% Isotropic transfer:
% centred difference applied to [cumulative] flux:
% (analogous to geting PDF by differentiating CDF)

T = -(flux(2:end) - flux(1:end-1)) /dk;

% wavenumber of psd is the centre point of the wavenumber interval used
% in the centred difference:
wn_T = (wn_step(2:end)+wn_step(1:end-1))/2;
wn_flux = wn_step;

return