function rec6fgh6, rec6, radec_star, dradec_star

;+
; NAME:
;  radrec6
; PURPOSE: (one line)
;  Convert from [x,y,z,x',y',z'] to [f,g,h,f',g',h']
; DESCRIPTION:
;  Convert from [x,y,z,x',y',z'] to [f,g,h,f',g',h']
; CATEGORY:
;  Astronomy
; CALLING SEQUENCE:
;  fgh6 = rec6fgh6(rec6, radec_star, dradec_star)
; 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
;  radec_star - vector of [A,D]
;   A - right ascention in radians
;   D - declination in radians
; OPTIONAL INPUT PARAMETERS:
;  dradec_star - vector of [dA, dD]
;   dA = d(A)/dt (NOT cosdec dA/dt)
;   dD = dD/dt
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  fgh6 - vector of [f,g,h,f',g',h']
;   fgh are distances along unit vectors IJK, as defined in 
;   Green Spherical Astronomy, sec 13.2.
;  K - in the direction of star (A,D)
;  I - perpendicular to K and celestial N pole ("East on sky")
;  J - perpendicular to I and K ("North on sky")
;
;  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 f, g, and h are km, and f', g', and h' are in km/s
;
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  None
; PROCEDURE:
; MODIFICATION HISTORY:
;  Written 2006 Jan 02, by Leslie Young, SwRI
;  2006 Mar 11 LAY. Corrected decumentation
;-

    if n_params() lt 3 then dradec_star = [0.d,0.d]

    A = double(radec_star[0])
    D = double(radec_star[1])
    dA = double(dradec_star[0])
    dD = double(dradec_star[1])

    m = oc_x2f_rotmat(A, D)

    ; IDL index is (x, y) = (col, row)
    ; m00 = -sinA
    ; m10 =  cosA
    ; m20 =   0
    ; m01 = -cosA sinD
    ; m11 = -sinA cosD
    ; m21 = cosD
    ; m02 = cosA cosD
    ; m12 = sinA cosD
    ; m22 = sinD

    dm = [ [ -m[1,0] * dA              ,   $
              m[0,0] * dA              ,   $
              0.                       ],  $
           [ -m[1,1] * dA - m[0,2] * dD,   $
              m[0,1] * dA - m[1,2] * dD,   $
                           -m[2,2] * dD ], $
           [ -m[1,2] * dA + m[0,1] * dD,   $
              m[0,2] * dA + m[1,1] * dD,   $
                            m[2,1] * dD ] ]
    
    xyz = double(rec6[0:2])
    dxyz = double(rec6[3:5])
    fgh = transpose(m) # xyz ; like m ## xyz, but for fgh, xyz row vectors
    dfgh = transpose(m) # dxyz  + transpose(dm) # xyz
    
    fgh6 = [fgh, dfgh]

    return, fgh6

end







;    dm00 = -m[1,0] * dA
;    dm10 =  m[0,0] * dA
;    dm11 =  0.

;    dm01 = -m[1,1] * dA - m[0,2] * dD
;    dm11 =  m[0,1] * dA - m[1,2] * dD
;    dm12 =               -m[2,2] * dD
    
;    dm02 = -m[1,2] * dA + m[0,1] * dD
;    dm12 =  m[0,2] * dA + m[1,1] * dD
;    dm22 =                m[2,1] * dD
    
;    dm22_da = 0.
;    dm22_dd = m[2,1]
;    dm00_da = - m[1,0]
;    dm00_dd = 0.d
;    dm10_da = m[0,0]
;    dm10_dd =  0.d
;    dm11_da = 0.d
;    dm11_dd = 0.d
;    dm01_da = -m[1,1]
;    dm01_dd = -m[0,2]
;    dm11_da = m[0,1]
;    dm11_dd = -m[1,2]
;    dm12_da = 0. d
;    dm12_dd = -m[2,2]
;    dm02_da = -m[1,2]
;    dm02_dd = m[0,1]
;    dm12_da = m[0,2]
;    dm12_dd = m[1,1]
;    dm22_da = 0.
;    dm22_dd = m[2,1]