function [ds,ts] = rectbody(x,y,z,sb,ts,sigbody, ...                lx,ux,ly,uy,lz,uz)% ----------------------------------------------% Filling delta sigma into a rectangular body% ----------------------------------------------%% ts      : total sigma in each cell (being updated)% ds      : delta sigma in each cell%% x       : vector of the x coordinates of the mesh% y       : vector of the y coordinates of the mesh% z       : vector of the z coordinates of the mesh% sb      : Nxyz,1 vector of background conductivity% sigbody : total conductivity of the body% lx      : lower x boundary of the body% ux      : upper x boundary of the body% ly      : lower y boundary of the body% uy      : upper y boundary of the body% lz      : lower z boundary of the body% uz      : upper z boundary of the body[xx,yy,zz] = ndgrid(x,y,z);xyz = [xx(:) yy(:) zz(:)];ind = find(xyz(:,1) >= lx & xyz(:,1) <= ux & ...           xyz(:,2) >= ly & xyz(:,2) <= uy & ...           xyz(:,3) >= lz & xyz(:,3) <= uz);ts(ind) = sigbody*ones(size(ind));ds = ts - sb;