function [dFdx, dFdy, dFdz] = lsqgradient(F, dx, dy, dz)
	%% [dFdx, dFdy, dFdz] = lsqgradient(F, dx, dy, dz)
	% computes numerical gradient based on second order,
	% least squares approximation as detailed in PIV: A Practical Guide
	%
	% dF/dx_i ~=  2 f_(i+2) + f_(i+1) - f(i-1) - 2f_(i-2)
	%              ------------------------------------
	%                          10 dx_i
	%
	% F must be a 3D array of size Nx x Ny x Nz
	% dFdx_i an Nx x Ny x Nz arrays
	% at the edges, dFdx_i is calculated with either first order central 
	% differences or at the very edges, forwards/backwards differences 
	% 

	%% parse input
	if nargin < 2
		dx = 1;
		dy = 1;
		dz = 1;
	elseif nargin < 3
		dy = 1;
		dz = 1;
	elseif nargin < 4
		dz = 1;
	end
	
	%% allocate memory for dFdx
	Nx		= size(F, 1);
	Ny    = size(F, 2);
	Nz		= size(F, 3);
	
	%% generate kernel
	K						= [2 1 0 -1 -2] / 10;	% least squares kernel, third order accurate
	kernel_x				= reshape(K(1:5) / dx, [5 1 1]);
	kernel_y				= reshape(K(1:5) / dy, [1 5 1]);
	kernel_z				= reshape(K(1:5) / dz, [1 1 5]);
	
	%% convolve over given dimension
	% but values near edges will be nonsense
	dFdx					= convn(F, kernel_x, 'same');
	dFdy					= convn(F, kernel_y, 'same');
	dFdz					= convn(F, kernel_z, 'same');
	
	%% forwards and backwards difference at edges
	dFdx(1, :, :)		= (F(2,:,:)		- F(1,:,:)		) / dx;
	dFdy(:, 1, :)		= (F(:,2,:)		- F(:,1,:)		) / dy;
	dFdz(:, :, 1)		= (F(:,:,2)		- F(:,:,1)		) / dz;
	dFdx(end,:,:)		= (F(end,:,:)  - F(end-1,:,:)	) / dx; 
	dFdy(:,end,:)		= (F(:,end,:)	- F(:,end-1,:)	) / dy;
	dFdz(:,:,end)		= (F(:,:,end)	- F(:,:,end-1) ) / dz;	
	
	%% first order central difference near edges
	if Nx >= 3
		for i = [2 Nx-1]
			dFdx(i, :, :) = (F(i+1,:,:) - F(i-1,:,:)) / (2*dx);
		end
	end
	if Ny >= 3
		for j = [2 Ny-1]
			dFdy(:, j, :) = (F(:,j+1,:) - F(:,j-1,:)) / (2*dy);
		end
	end
	if Nz >= 3
		for k = [2 Nz-1]
			dFdz(:, :, k) = (F(:,:,k+1) - F(:,:,k-1)) / (2*dz);
		end
	end
end
