function rec6rad6, rec6

;+
; NAME:
;  radrec6
; PURPOSE: (one line)
;  Convert from [x,y,z,x',y',z'] to [r, ra, dec, r', ra', dec']
; DESCRIPTION:
;  Convert from [x,y,z,x',y',z'] to [r, ra, dec, r', ra', dec']
; CATEGORY:
;  Astronomy
; CALLING SEQUENCE:
;  rad6 = rec6rad6(rad6)
; INPUTS:
;  rec6 - vector of [x,y,z,x',y',z']
;  x - distance on X axis in units of LENGTH
;  y - distance on Y axis in units of LENGTH
;  z - distance on Z axis in units of LENGTH
;  x' - x derivative in LENGTH/TIME
;  y' - y derivative in LENGTH/TIME
;  z' - z derivative in LENGTH/TIME
; OPTIONAL INPUT PARAMETERS:
;  none
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  rad6 - vector of 
;   r - range in LENGTH UNITS
;   ra - right ascention in radians
;   dec - declination in radians
;   r' - range derivative in LENGTH/TIME
;   ra' - right ascention derivative in radians/TIME
;   dec' - declination derivative in radians/TIME
;
;  If x, y, and z are in km, and x', y', z' in km/s
;  (e.g., the same format as states returned by spice routines)
;  then r will be in km, r' in km/s, and ra' and dec' in radian/s, 
;
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  None
; PROCEDURE:
; MODIFICATION HISTORY:
;  Written 2005 Apr 26, by Leslie Young, SwRI
;  2006 Mar 11 LAY, corrected documentation
;-

    x = double(rec6[0])
    y = double(rec6[1])
    z = double(rec6[2])

    r2 = x^2d + y^2d + z^2
    if r2 le 0 then begin
        r = 0.d
        ra = 0.d
        dec = 0.d
    endif else begin
      r = sqrt(r2)
      dec = asin(z/r)
      if x eq 0.d and y eq 0.d then begin
          ra = 0.d
      endif else begin
        ra = atan(y,x)	; Find angles.
        if ra lt 0 then ra = ra+2*!dpi ; add 2 pi to angles < 0.
      endelse
    endelse
    
    r2 = r^2.d
    
    polar = ( x^2 + y^2 eq 0.d)
    origin = (r eq 0.d)
    if origin then drdx = 0. else drdx = x/r
    if origin then drdy = 0. else drdy = y/r
    if origin then drdz = 0. else drdz = z/r
    if polar then dradx = 0. else dradx = -(y/(x^2 + y^2))
    if polar then drady = 0. else drady =  (x/(x^2 + y^2))
    dradz = 0.
    if polar then ddecdx = 0. else ddecdx = -((x*z)/(r2*Sqrt(x^2+y^2)))
    if polar then ddecdy = 0. else ddecdy = -((y*z)/(r2*Sqrt(x^2+y^2)))
    if origin then ddecdz = 0. else ddecdz = sqrt(x^2+y^2)/r2
    
    d = [ [drdx,     drdy,     drdz], $
          [dradx,    drady,    dradz], $
          [ddecdx,   ddecdy,   ddecdz] ]

	drec = double(rec6[3:5])
	drad = reform(d ## drec)
	
 	rad6 = reform([r,ra,dec,drad])
	return, rad6
	
end