function interp2d, A, x0, y0, x1, y1, nxny, missing=missing,   $
    grid=grid, quintic=quintic, regular=regular, cubic=cubic

; Check that A is a 2d array

sz = size(a)

if sz(0) ne 2 then begin
  message,'A must be a 2-d array',/cont
  return,-1
endif

; Check that the dimensions of x0 and y0 match A

if (sz(1) ne n_elements(X0)) or     $
   (sz(2) ne n_elements(Y0)) then begin
  message,'Dimensions of A, X0, Y0 are not consistent',/cont
  return,-1
endif

if not keyword_set(regular) then begin
  if n_elements(nxny) eq 0 then nxny = [51,51]

; Call triangulate and trigrid to get a regularly spaced grid

  nx = n_elements(X0) & ny = n_elements(Y0)
  gs = [(max(X0)-min(X0))/(nxny(0)-1), (max(Y0)-min(Y0))/(nxny(1)-1)]

  xx = rebin(reform(x0),nx,ny)
  yy = rebin(transpose(reform(y0)),nx,ny)
  triangulate, xx, yy, tr
  if n_elements(quintic) eq 0 then quintic = 0  ; Make sure quintic is define
  zz = trigrid(xx, yy, A, tr, gs, quintic=quintic)

  zz = zz(0:nxny(0)-1,0:nxny(1)-1)  ; Make sure the dimensions are matched
endif else zz = A           ; /regular was set -- x0 and y0 are linear
sz = size(zz)

xslope = (max(X0)-min(X0)) / (sz(1)-1)
yslope = (max(Y0)-min(Y0)) / (sz(2)-1)
xx = min(X0) + sz(1)*xslope
yy = min(Y0) + sz(2)*yslope

; Map the coordinates

x2 = (x1 - min(x0)) / xslope
y2 = (y1 - min(y0)) / yslope

; Now interpolate

if n_elements(grid)    eq 0 then grid = 0
if n_elements(cubic)   eq 0 then cubic= 0
if n_elements(missing) eq 0 then $
        return,interpolate(zz,x2,y2,grid=grid,cubic=cubic) else $
        return,interpolate(zz,x2,y2,grid=grid,missing=missing,cubic=cubic)
end