function [var, Dsx, Dsy] = err_s2p (phi, Sx, Sy)
%function [var,Dsx,Dsy]= err_s2p(phi,Sx,Sy,Nx,Ny)
% [variance, Dsx, Dsy]=err_s2p(phi,Sx,Sy,Nx,Ny) returns errors in solution
%  Dsx= phi(y,x+1)-phi(y,x)-Sx(y,x)
%  Dsy= phi(y+1,x)-phi(y,x)-Sy(y,x)
%   variance=sum(Dsx^2+Dsy^2)

[Ny, Nx] = size (Sx);
Dsx=zeros(Ny,Nx); Dsy=Dsx;
Dsx(1:Ny,1:Nx-1)=phi(1:Ny,2:Nx)-phi(1:Ny,1:Nx-1)-Sx(1:Ny,1:Nx-1);
Dsy(1:Ny-1,1:Nx)=phi(2:Ny,1:Nx)-phi(1:Ny-1,1:Nx)-Sy(1:Ny-1,1:Nx);
var=trace(Dsx*Dsx') +trace(Dsy*Dsy');