Procedure: Grad Purpose: Calculates the gradient of a 2d or 3d grid in one of two ways. In 2d: Method1(default): gradientX = (grid[x+1,y] - grid[x,y] + grid[x+1,y+1] - grid[x,y+1]) / (2*dx) gradientY = (grid[x,y+1] - grid[x,y] + grid[x+1,y+1] - grid[x+1,y]) / (2*dy) Method2(leftright): gradientX = (grid[x+1,y] - grid[x,y] + grid[x,y] - grid[x-1,y]) / (2*dx) gradientY = (grid[x,y+1] - grid[x,y] + grid[x,y] - grid[x,y-1]) / (2*dy) This method is actually equivalent to: gradientX = (grid[x+1,y] - grid[x-1,y]) / (2*dx) gradientY = (grid[x,y+1] - grid[x,y-1]) / (2*dy) In 3d: Method1(default): gradientX = (grid[x+1,y,z] - grid[x,y,z] + grid[x+1,y+1,z] - grid[x,y+1,z] + grid[x+1,y,z+1] - grid[x,y,z+1] + grid[x+1,y+1,z+1] - grid[x,y+1,z+1]) / (4*dx) gradientY = (grid[x,y+1,z] - grid[x,y,z] + grid[x+1,y+1,z] - grid[x+1,y,z] + grid[x,y+1,z+1] - grid[x,y,z+1] + grid[x+1,y+1,z+1] - grid[x+1,y,z+1]) / (4*dy) gradientZ = (grid[x,y,z+1] - grid[x,y,z] + grid[x+1,y,z+1] - grid[x+1,y,z] + grid[x,y+1,z+1] - grid[x,y+1,z] + grid[x+1,y+1,z+1] - grid[x+1,y+1,z]) / (4*dz) Method2(leftright): gradientX = grid[x+1,y,z] - grid[x-1,y,z] / (2*dx) gradientY = grid[x,y+1,z] - grid[x,y-1,z] / (2*dy) gradientZ = grid[x,y,z+1] - grid[x,y,z-1] / (2*dz) Method1 will produce an output that is one element smaller in each dimension and whose element centers are offset by half the nominal spacing of the grid. Method2 will have the same centers and same number of elements as the original grid(if the original grid had regular spacing). Example: Inputs: grid: an NxM grid of points, if it contains NaNs the output may be unpredictable.(or an NxMxP) x(optional): An N length array specifying the positions of the grid points on the x-axis xc should be monotonic and should contain no NaNs. If unset this routine will assume dx = 1.0 y(optional): An M length array specifying the positions of the grid points on the y-axis yc should be monotonic and should contain no NaNs. If unset this routine will assume dy = 1.0 z:(optional) a P length array specifying the positions of the grid points on the z-axis. zc should be monotonic and should contain no NaNs. If unset this routine will assume dz = 1.0 Keywords: grad: The gradient is output through this keyword as an NxMx2 array of points. grad[*,*,0] is the x gradient & grad[*,*,1] is the y gradient xout: The positions of the gradient outputs on the x axis are output through this keyword as an N length array yout: The positions of the gradient outputs on the y axis are output through this keyword as an M length array xy: The positions for each output point are passed out as pairs through this keyword. The output array will have dimensions N*Mx2,(N times M by 2) dxy: The gradient for each point is passed out as pairs through this keyword. The output array will have dimensions N*Mx2,(N times M by 2) leftright: Set this keyword if you want to use the second method of gradient calculation. Notes: 1. This procedure is not particularly tolerant of NaNs in the input, so you should remove them before passing them into this routine. 2. The output may have slightly different centers/ dimensions as the input. This is will definitely be the case if the input array had irregular dimensions. 3. xy,dxy are useful output format keywords for the plotxyvec routine While grad,xout, & yout may be easier for other tasks.