function [Tkxky,wn_x,wn_y,JPsiZeta] = psi2trans2d_v2(dx,dy,psi)
% calculates the spectral transfer of the velocity field given by
% the streamfunction
%-------
% INPUT:
%-------
% dx 	grid spacing, x direction
% dy    grid spacing, y direction
% psi 	streamfunction, not necessarily square, psi(y,x)
%-------
% OUTPUT:
%-------
% T 		
% wn_T 	is the wavenumbers corresponding to T
% JPsiZeta	J(psi,zeta), where zeta=nable^2 psi
%-------------
% History:
% Built by Rob Scott on Feb. 27th, (TPA, 5:30 AM)
% from psi2flux_v2.m
%------------------------------------------------------------------------------

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


% 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)));

  ddxq = real(ifft2(j*kx.*fft2(q)));
  ddyq = real(ifft2(j*ky.*fft2(q)));

  JPsiZeta = u.*ddxq+v.*ddyq;

  Tkxky = real(conj(fft2(psi)).*fft2(JPsiZeta))/(nx^2 * ny^2);

return