;+
;NAME:
; curl_rtp
;PURPOSE:
; Does the curl for a 3d vector, using the DERIV routine,
; in Cartesian coordinates
;CALLING SEQUENCE:
; curl_rtp, 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
; 21-jan-2004, spherical q's x -> r, y -> theta, z -> phi
;-
PRO curl_rtp, ax, ay, az, x, y, z, curl_x, curl_y, curl_z, $
              use_saved_grids = use_saved_grids, _extra = _extra

   Common saved_rtp, rx, sintheta
   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

   If(n_elements(rx) Eq 0 Or (Not keyword_set(use_saved_grids))) Then Begin
       rx = rebin(x, nx, ny, nz)
       sintheta = rebin(sin(y), ny, nz, nx)
       sintheta = transpose(sintheta, [2, 0, 1])
   Endif
   rsintheta = rx*sintheta
;X component
;   curl_x = d(sintheta*az)dy/(rsintheta) - daydz/rsintheta
   curl_x = deriv_dy(y, sintheta*az, _extra = _extra)/rsintheta - $
     deriv_dz(z, ay, _extra = _extra)/rsintheta
;Y component
;   curl_y = daxdz/rsintheta - d(r*az)dx/r
   curl_y = deriv_dz(z, ax, _extra = _extra)/rsintheta - $
     deriv_dx(x, rx*az, _extra = _extra)/rx
;Z component
;   curl_z = d(r*ay)dx/r - daxdy/r
   curl_z = deriv_dx(x, rx*ay, _extra = _extra)/rx - $
     deriv_dy(y, ax, _extra = _extra)/rx
   
   RETURN
END