;+
;NAME:
; curl_xyz
;PURPOSE:
; Does the curl for a 3d vector, using the DERIV routine,
; in Cartesian coordinates
;CALLING SEQUENCE:
; curl_xyz, ax, ay, az, x, y, z, curl_x, curl_y, curl_z
;INPUT:
; ax, ay, az = are 3d arrays defined on (x,y,z)
; x, y, z, = the coordinates on which the arrays is defined.
;OUTPUT:
; curl_x, curl_y, curl_z = the curl of a
;HISTORY:
; 18-aug-2000, jmm, jimm@ssl.berkeley.edu
; switched to use deriv_dx, _dy, and _dz, 
; 1-apr-2003, jmm, added _extra
; 29-dec-2004, jmm, added spherical keyword, if this is set, 
;              call curl_rtp.pro
;-
PRO curl_xyz, ax, ay, az, x, y, z, curl_x, curl_y, curl_z, $
              spherical = spherical, _extra = _extra

   IF(keyword_set(spherical)) THEN BEGIN ;spherical q's
       curl_rtp, ax, ay, az, x, y, z, curl_x, curl_y, curl_z, _extra=_extra
       RETURN
   ENDIF
   sizax = size(ax)
   IF(sizax[0] NE 3) THEN BEGIN
      message, /info, 'Array must be 3-d'
      RETURN
   ENDIF
   sizay = size(ay)
   IF(sizay[0] NE 3) THEN BEGIN
      message, /info, 'Array must be 3-d'
      RETURN
   ENDIF
   sizaz = size(az)
   IF(sizaz[0] NE 3) THEN BEGIN
      message, /info, 'Array must be 3-d'
      RETURN
   ENDIF
   IF(total(sizax[0:3]) NE total(sizay[0:3]) OR $
      total(sizax[0:3]) NE total(sizaz[0:3])) THEN BEGIN
      message, /info, 'Arrays have mismatched sizes'
      RETURN
   ENDIF
   
   nx = N_ELEMENTS(x)
   IF(Nx NE sizax[1]) THEN BEGIN
      message, /info, 'Bad X array'
      RETURN
   ENDIF
   ny = N_ELEMENTS(y)
   IF(Ny NE sizax[2]) THEN BEGIN
      message, /info, 'Bad Y array'
      RETURN
   ENDIF
   nz = N_ELEMENTS(z)
   IF(Nz NE sizax[3]) THEN BEGIN
      message, /info, 'Bad Z array'
      RETURN
   ENDIF

;X component
;   curl_x = dazdy - daydz
   curl_x = deriv_dy(y, az, _extra = _extra)-deriv_dz(z, ay, _extra = _extra)
;Y component
;   curl_y = daxdz - dazdx
   curl_y = deriv_dz(z, ax, _extra = _extra)-deriv_dx(x, az, _extra = _extra)
;Z component
;   curl_z = daydx - daxdy
   curl_z = deriv_dx(x, ay, _extra = _extra)-deriv_dy(y, ax, _extra = _extra)
   RETURN
END