Select Git revision
lsqgradient2d.m
jml1g18 authored
lsqgradient2d.m 1.52 KiB
function [dFdx, dFdy] = lsqgradient2d(F, dx, dy)
%% [dFdx, dFdy] = lsqgradient(F, dx, dy)
% 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 2D array of size Nx x Ny
% dFdx_i an Nx x Ny array
% 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;
elseif nargin < 3
dy = 1;
end
%% allocate memory for dFdx
Nx = size(F, 1);
Ny = size(F, 2);
dFdx = nan(size(F), class(F));
dFdy = nan(size(F), class(F));
%% forwards and backwards difference at edges
dFdx(1, :) = (F(2,:) - F(1,:) ) / dx;
dFdy(:, 1) = (F(:,2) - F(:,1) ) / dy;
dFdx(end,:) = (F(end,:) - F(end-1,:)) / dx;
dFdy(:,end) = (F(:,end) - F(:,end-1)) / dy;
%% 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
%% second order least squares near middle
% use an expansion of f(x)' which is accurate to O(dx^4)
for i = 3:(Nx-2)
dFdx(i, :) = (2*F(i+2,:) + F(i+1,:) - F(i-1,:) - 2*F(i-2,:)) / (10*dx);
end
for j = 3:(Ny-2)
dFdy(:, j) = (2*F(:,j+2) + F(:,j+1) - F(:,j-1) - 2*F(:,j-2)) / (10*dy);
end
end