function del_sq, x, stride = stride
;+
; NAME:
;	DEL_SQ
;
;
; PURPOSE:
;	Return an approximation to Del^2 of the supplied array.
;
;
; CATEGORY:
;	Utilities (Maths)
;
;
; CALLING SEQUENCE:
;	d2 = del_sq(x)
;
;
; INPUTS:
;	x	(float)	The array of 1, 2 or 3 dimensions to be
;			differentiated. 
;
;
; KEYWORD PARAMETERS:
;	stride	(float)	The spacing of the points in the array (must
;			be the same in all dimensions (if not given
;			unit spacing is assumed).
;
;
; OUTPUTS:
;	d2	Float	The values of del^2 of the data. A border of
;			zeroes will be left around the array.
;
;
; MODIFICATION HISTORY:
;	Original: 2/6/05; SJT
;-

on_error, 2

sz = size(x)

if sz[0] lt 1 or sz[0] gt 3 then $
  message, "The array passed must have 1, 2 or 3 dimensions"

case sz[0] of
    1: begin
        dx = (shift(x, -1)+shift(x, 1))/2.
        dx -= x
        dx[0] = 0.
        dx[sz[1]-1] = 0.
    end
    2: begin
        dx = (shift(x, -1, 0) + shift(x, 1, 0)+ $
              shift(x, 0, -1) + shift(x, 0, 1))/2.
        dx -= 2.*x
        dx[0, *] = 0.
        dx[*, 0] = 0.
        dx[sz[1]-1, *] = 0.
        dx[*, sz[2]-1] = 0.
    end
    3: begin
        dx = (shift(x, -1, 0, 0) + shift(x, 1, 0, 0) + $
              shift(x, 0, -1, 0) + shift(x, 0, 1, 0) + $
              shift(x, 0, 0, -1) + shift(x, 0, 0, 1))/2.
        dx -= 3.*x
        dx[0, *, *] = 0.
        dx[*, 0, *] = 0.
        dx[*, *, 0] = 0.
        dx[sz[1]-1, *, *] = 0.
        dx[*, sz[2]-1, *] = 0.
        dx[*, *, sz[3]-1] = 0.
    end
endcase

if keyword_set(stride) then dx /= stride

return, dx

end