;+
; NAME:
;  oc_lla2xyz
; PURPOSE: (one line)
;  Return the geocentric J2000 XYZ of a site on Earth at a given time
; DESCRIPTION:
;  Return the geocentric J2000 state vector of a site on Earth at a given time
; CATEGORY:
;  Occultations
; CALLING SEQUENCE:
;  rec6 = oc_lla2xyz(lonlatalt, et)
; INPUTS:
;  lonlatalt = lon,lat,alt vector
;   lonlatalt[0] = East longitude in radian
;   lonlatalt[1] = latitude in radian
;   lonlatalt[2] = altitude in km
;  et = TBD seconds after J2000 ('ephemeris time' in NAIF lingo)
;   a scalar or a vector
; OPTIONAL INPUT PARAMETERS:
;  frame (default 'J2000')
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  rec6 = [X,Y,Z, dX/dt, dY/dt, dZ/dt]
;   This is a J2000 state vector as descibed in the NAIF documents
;   XYZ are in km, derivatives in km/s
;   Coordinate system is J2000. 
;   Origin is the Earth center, head is the site at the specified lon,lat,alt
;  If et is a vector of n elements, then
;   rec6 will be a matrix 6 x n
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  None
; PROCEDURE:
;
; MODIFICATION HISTORY:
;  2005 Dec 29 Leslie A Young SwRI
;-

function oc_lla2xyz, lonlatalt, et, debug=debug

    naifinit

    lon = lonlatalt[0]  ; E lon in radian
    lat = lonlatalt[1]  ; lat in radian
    alt = lonlatalt[2]  ; alt in km

    ;;;;;;;; 
    ;; Code fragment from cspice_pxform.html

    cspice_bodvar, 399L, 'RADII', abc
    equatr =  abc[0]
    polar  =  abc[2]
    
    ;;
    ;; Calculate the flattening factor for Earth.
    ;;
    f =  ( equatr - polar ) / equatr
 
    ;;
    ;; Calculate the Cartesian coordinates on Earth for the
    ;; location at 'lon', 'lat', 'alt'.
    ;;
    cspice_georec, lon, lat, alt, equatr, f, epos
    if keyword_set(debug) then begin
        print, 'oc_lla2xyz: epos'
        print, epos
        print, vabs(epos)
    end

    if not keyword_set(frame) then frame = 'J2000'

    if n_elements(et) eq 1 then begin
    
    ;;
    ;; Retrieve the transformation matrix from "IAU_EARTH"
    ;; to "J2000" at epoch 'et'.
    ;;
        cspice_sxform, "IAU_EARTH", frame, et, rotate6
    
    ;;
    ;; Transform the Cartesian position vector from "IAU_EARTH"
    ;; to "J2000."
    ;;
        rec6 = transpose(rotate6) # [epos, dblarr(3) ]

    endif else begin

        n = n_elements(et)
        rec6 = dblarr(6,n)
        for i=0,n-1 do begin
            cspice_sxform, "IAU_EARTH", frame, et[i], rotate6
            rec6[*,i] = transpose(rotate6) # [epos, dblarr(3) ]
        end
               
    end

    return, rec6

end