function rad6fgh6, rad6, radec_star, dradec_star

;+
; NAME:
;  radrec6
; PURPOSE: (one line)
;  Convert from [r, ra, dec, r', ra', dec'] to [f,g,h,f',g',h']
; DESCRIPTION:
;  Convert from [r, ra, dec, r', ra', dec'] to [f,g,h,f',g',h']
; CATEGORY:
;  Astronomy
; CALLING SEQUENCE:
;  fghstate = rad6fgh6(rad6, radec_star, dradec_star)
; INPUTS:
;  rad6 - vector of 
;   r - range in LENGTH UNITS (eg km)
;   ra - right ascention in radians
;   dec - declination in radians
;   r' - range derivative in LENGTH/TIME (km/s)
;   ra' - right ascention derivative in radians/TIME (NOT cosdec dra/dt)
;   dec' - declination derivative in radians/TIME
;  radec_star - vector of;
;   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']
;
;  f - distance on F axis in units of LENGTH
;  g - distance on G axis in units of LENGTH
;  h - distance on H axis in units of LENGTH
;  f' - f derivative in LENGTH/TIME
;  g' - g derivative in LENGTH/TIME
;  h' - h derivative in LENGTH/TIME
;
;   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 r is in km, r' in km/s, and ra' and dec' in radian/s, 
;  then f, g, and h are in km and f', g' and h' are in km/s
;
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  None
; PROCEDURE:
;  This method is slightly preferred over rec6fgh6 when the
;  fundemental calculated variables are ra, dec
; MODIFICATION HISTORY:
;  Written 2006 Jan 02 Leslie Young
;  2006 May 11 LAY corrected documentation
;-
    
    if n_params() lt 3 then dradec_star = [0.d,0.d]

    r = double(rad6[0])
    ra = double(rad6[1])
    dec = double(rad6[2])
    rdot = double(rad6[3])
    radot = double(rad6[4])
    decdot = double(rad6[5])
    A = double(radec_star[0])
    D = double(radec_star[1])
    Adot = double(dradec_star[0])
    Ddot = double(dradec_star[1])
    cosdec = cos(dec)
    sindec = sin(dec)
    cosD = cos(D)
    sinD = sin(D)
    cos_dra = cos(ra-A)
    sin_dra = sin(ra-A)
    cos_ddec = cos(dec-D)
    sin_ddec = sin(dec-D)
    eps = 1-cos_dra

    f = r * sin_dra * cosdec
    g = r * (sin_ddec + eps * sinD * cosdec)
    h = r * (cos_ddec - eps * cosD * cosdec)
    
    ; calculate dfdr = sin_dra * cosdec rather than f/r etc to
    ; avoid having to check for zero denominators
    
    dfdr    = sin_dra * cosdec
    dfdra   = r * cos_dra * cosdec
    dfddec  = - r * sin_dra * sindec
    dfdA    = - dfdra
    dfdD    = 0.d

    dgdr    = sin_ddec + eps * sinD * cosdec
    dgdra   = r * (sin_dra * sinD * cosdec)
    dgddec  = r * (cos_ddec - eps * sinD * sindec)
    dgdA    = - dgdra
    dgdD    = r * (-cos_ddec + eps * cosD * cosdec)

    dhdr    = cos_ddec - eps * cosD * cosdec
    dhdra   = r * (- sin_dra * cosD * cosdec)
    dhddec  = r * (-sin_ddec + eps * cosD * sindec)
    dhdA    = - dhdra
    dhdD    = r * (sin_ddec + eps * sinD * cosdec)

    d = [ [dfdr,   dfdra,  dfddec, dfdA, dfdD], $
          [dgdr,   dgdra,  dgddec, dgdA, dgdD], $
          [dhdr,   dhdra,  dhddec, dhdA, dhdD] ]
    drad = [rdot, radot, decdot, Adot, Ddot]
    dfgh = reform(d ## drad)

    fgh6 = reform([f,g,h,dfgh])

    return, fgh6
	
end