function [psd,wn_out] = vel2psd(dx,dy,u_in,v_in)
% calculates the power spectral density PSD of the velocity given
% the two orthogonal components, u_in, v_in
%-------
% 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:
%-------
% psd 		is the PSD of velocity
% wn_out 	is the wavenumbers corresponding to PSD
%------------------------------------------------------------------------------
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;

% wavenumbers: (units are e.g. [rad/m] if Lx is in [m])
%Brian (Glenn Flierl) use: k256=(1/10)*[0:128,-127:-1];
wn_x = 2*pi*[0:nx/2,-nx/2+1:-1]' / Lx;
wn_y = 2*pi*[0:ny/2,-ny/2+1:-1]' / Ly;

% 2d-space of wavenumbers:
[kx,ky] = meshgrid(wn_x,wn_y);
wvsq = kx.*kx+ky.*ky;

% Discrete Fourier transforms of velocities:
uhat = fft2(u_in);
vhat = fft2(v_in);

% 2-dimensional power spectral density:
% psd(ky,kx)
psd_2d = 1/2 * (abs(uhat).^2 + abs(vhat).^2 )/(nx^2 * ny^2);

% Integrate 2d PSD around lines of constant k
%---------------------------------------

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

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

% centred difference applied to cumulative spectrum:
% (analogous to geting PDF by differentiating CDF)

psd = (psc(2:end) - psc(1:end-1)) /dk;

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

if size(psd,1) < size(psd,2)
   psd = psd';
end
return