function rad6rec6, rad6

;+
; NAME:
;  radrec6
; PURPOSE: (one line)
;  Convert from [r, ra, dec, r', ra', dec'] to [x,y,z,x',y',z']
; DESCRIPTION:
;  Convert from [r, ra, dec, r', ra', dec'] to [x,y,z,x',y',z']
; CATEGORY:
;  Astronomy
; CALLING SEQUENCE:
;  state = rad6rec6(rad6)
; INPUTS:
;  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
; OPTIONAL INPUT PARAMETERS:
;  none
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  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
;
;  If r is in km, r' in km/s, and ra' and dec' in radian/s, 
;  then x, y, and z will be in km, and x', y', z' in km/s
;  (e.g., the same format as states returned by spice routines)
;
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  None
; PROCEDURE:
; MODIFICATION HISTORY:
;  Written 2005 Apr 26, by Leslie Young, SwRI
;  2006 Mar 11 LAY correct documentation
;-

    r = double(rad6[0])
    ra = double(rad6[1])
    dec = double(rad6[2])
    cosra = cos(ra)
    sinra = sin(ra)
    cosdec = cos(dec)
    sindec = sin(dec)
    x = r*cosdec*cosra
    y = r*cosdec*sinra
    z = r*sindec
    
    ; calculate dxdr = cosdec cosra rather than x/r etc to
    ; avoid having to check for zero denominators
    dxdr   =  cosdec * cosra
    dxdra  = -r * cosdec * sinra
    dxddec = -r * cosra * sindec
    dydr   =  cosdec * sinra
    dydra  =  r * cosdec * cosra
    dyddec = -r * sinra * sindec
    dzdr   =  sindec
    dzdra  =  0
    dzddec =  r * cosdec

    d = [ [dxdr,   dxdra,  dxddec], $
          [dydr,   dydra,  dyddec], $
          [dzdr,   dzdra,  dzddec] ]
    drad = double(rad6[3:5])
    drec = reform(d ## drad)

    rec6 = reform([x,y,z,drec])
    return, rec6
	
end