How does intgrad2 function?

Dear all,

I am no expert in programming and can not understand what is the basis of this function intgrad2. Someone can explain, for example, which translates into the following extension of the algorithm?

Thanks

% do the trailing edge in x, backward difference

indx = nx;

indy = (1:ny)';

ind = indy + (indx-1)*ny;

rind = repmat(L+(1:ny)',1,2);

cind = [ind-ny,ind];

dfdx = repmat([-1 1]./dx(end),ny,1);

Af(L+(1:ny),:) = [rind,cind,dfdx];

rhs(L+(1:ny)) = fx(:,end);

L = L+ny;

----->Algorith

function fhat = intgrad2(fx,fy,dx,dy,f11)

% intgrad: generates a 2-d surface, integrating gradient information.

% usage: fhat = intgrad2(fx,fy)

% usage: fhat = intgrad2(fx,fy,dx,dy)

% usage: fhat = intgrad2(fx,fy,dx,dy,f11) %

% arguments: (input)

% fx,fy - (ny by nx) arrays, as gradient would have produced. fx and

% fy must both be the same size. Note that x is assumed to

% be the column dimension of f, in the meshgrid convention.

% % nx and ny must both be at least 2.

% % fx and fy will be assumed to contain consistent gradient

% information. If they are inconsistent, then the generated

% gradient will be solved for in a least squares sense.

% % Central differences will be used where possible.

% % dx - (OPTIONAL) scalar or vector - denotes the spacing in x

% if dx is a scalar, then spacing in x (the column index

% of fx and fy) will be assumed to be constant = dx.

% if dx is a vector, it denotes the actual coordinates

% of the points in x (i.e., the column dimension of fx

% and fy.) length(dx) == nx

% % DEFAULT: dx = 1

% % dy - (OPTIONAL) scalar or vector - denotes the spacing in y

% if dy is a scalar, then the spacing in x (the row index

% of fx and fy) will be assumed to be constant = dy.

% if dy is a vector, it denotes the actual coordinates

% of the points in y (i.e., the row dimension of fx

% and fy.) length(dy) == ny

% % DEFAULT: dy = 1 % % f11 - (OPTIONAL) scalar - defines the (1,1) eleemnt of fhat

% after integration. This is just the constant of integration.

% % DEFAULT: f11 = 0

% % arguments: (output)

% fhat - (nx by ny) array containing the integrated gradient

% % Example usage 1: (Note x is uniform in spacing, y is not.)

% xp = 0:.1:1;

% yp = [0 .1 .2 .4 .8 1];

% [x,y]=meshgrid(xp,yp);

% f = exp(x+y) + sin((x-2*y)*3);

% [fx,fy]=gradient(f,.1,yp);

% tic,fhat = intgrad2(fx,fy,.1,yp,1);toc

% % Time required was 0.06 seconds

% % Example usage 2: Large grid, 101x101

% xp = 0:.01:1;

% yp = 0:.01:1;

% [x,y]=meshgrid(xp,yp);

% f = exp(x+y) + sin((x-2*y)*3);

% [fx,fy]=gradient(f,.01);

% tic,fhat = intgrad2(fx,fy,.01,.01,1);toc

% % Time required was 4 seconds

% Author; John D'Errico

% Current release: 2

% Date of release: 1/27/06

% size

if (length(size(fx))>2) (length(size(fy))>2)

error 'fx and fy must be 2d arrays'

end

[ny,nx] = size(fx);

if (nx~=size(fy,2)) (ny~=size(fy,1))

error 'fx and fy must be the same sizes.'

end

if (nx<2) (ny<2)

error 'fx and fy must be at least 2x2 arrays'

end

% supply defaults if needed

if (nargin<3) isempty(dx)

% default x spacing is 1

dx = 1;

end

if (nargin<4) isempty(dy)

% default y spacing is 1

dy = 1;

end

if (nargin<5) isempty(f11)

% default integration constant is 0

f11 = 0;

end

% if scalar spacings, expand them to be vectors

dx=dx(:);

if length(dx) == 1

dx = repmat(dx,nx-1,1);

elseif length(dx)==nx

% dx was a vector, use diff to get the spacing

dx = diff(dx);

else

error 'dx is not a scalar or of length == nx'

end

dy=dy(:);

if length(dy) == 1

dy = repmat(dy,ny-1,1);

elseif length(dy)==ny

% dy was a vector, use diff to get the spacing

dy = diff(dy);

else

error 'dy is not a scalar or of length == ny'

end

if (length(f11) > 1) ~isnumeric(f11) isnan(f11) ~isfinite(f11)

error 'f11 must be a finite scalar numeric variable.'

end

% build gradient design matrix, sparsely. Use a central difference

% in the body of the array, and forward/backward differences along

% the edges.

% A will be the final design matrix. it will be sparse.

% The unrolling of F will be with row index running most rapidly.

rhs = zeros(2*nx*ny,1);

% but build the array elements in Af

Af = zeros(2*nx*ny,6);

L = 0;

% do the leading edge in x, forward difference

indx = 1;

indy = (1:ny)';

ind = indy + (indx-1)*ny;

rind = repmat(L+(1:ny)',1,2);

cind = [ind,ind+ny];

dfdx = repmat([-1 1]./dx(1),ny,1);

Af(L+(1:ny),:) = [rind,cind,dfdx];

rhs(L+(1:ny)) = fx(:,1);

L = L+ny;

% interior partials in x, central difference

if nx>2

[indx,indy] = meshgrid(2:(nx-1),1:ny);

indx = indx(:);

indy = indy(:);

ind = indy + (indx-1)*ny;
m = ny*(nx-2);

rind = repmat(L+(1:m)',1,2);

cind = [ind-ny,ind+ny];

dfdx = 1./(dx(indx-1)+dx(indx));

dfdx = dfdx*[-1 1];

Af(L+(1:m),:) = [rind,cind,dfdx];

rhs(L+(1:m)) = fx(ind);

L = L+m;

end

% do the trailing edge in x, backward difference

indx = nx;

indy = (1:ny)';

ind = indy + (indx-1)*ny;

rind = repmat(L+(1:ny)',1,2);

cind = [ind-ny,ind];

dfdx = repmat([-1 1]./dx(end),ny,1);

Af(L+(1:ny),:) = [rind,cind,dfdx];

rhs(L+(1:ny)) = fx(:,end);

L = L+ny;

% do the leading edge in y, forward difference

indx = (1:nx)';

indy = 1;

ind = indy + (indx-1)*ny;

rind = repmat(L+(1:nx)',1,2);

cind = [ind,ind+1];

dfdy = repmat([-1 1]./dy(1),nx,1);

Af(L+(1:nx),:) = [rind,cind,dfdy];

rhs(L+(1:nx)) = fy(1,:).';

L = L+nx;

% interior partials in y, use a central difference

if ny>2

[indx,indy] = meshgrid(1:nx,2:(ny-1));

indx = indx(:);

indy = indy(:);

ind = indy + (indx-1)*ny;

m = nx*(ny-2);

rind = repmat(L+(1:m)',1,2);

cind = [ind-1,ind+1];

dfdy = 1./(dy(indy-1)+dy(indy));
dfdy = dfdy*[-1 1];

Af(L+(1:m),:) = [rind,cind,dfdy];

rhs(L+(1:m)) = fy(ind);

L = L+m;

end

% do the trailing edge in y, backward diffeence

indx = (1:nx)';

indy = ny;

ind = indy + (indx-1)*ny;

rind = repmat(L+(1:nx)',1,2);

cind = [ind-1,ind];

dfdy = repmat([-1 1]./dy(end),nx,1);

Af(L+(1:nx),:) = [rind,cind,dfdy];

rhs(L+(1:nx)) = fy(end,:).';

% finally, we can build the rest of A itself, in its sparse form.

A = sparse(Af(:,1:2),Af(:,3:4),Af(:,5:6),2*nx*ny,nx*ny);

% Finish up with f11, the constant of integration.

% eliminate the first unknown, as f11 is given.

rhs = rhs - A(:,1)*f11;

% Solve the final system of equations. They will be of

% full rank, due to the explicit integration constant.

% Just use sparse \

fhat = A(:,2:end)\rhs;

fhat = reshape([f11;fhat],ny,nx);