;+
; NAME:
;   oc_tracklla
;
; PURPOSE:
;   Calculate the longitude and latitude for occultation tracks
;
; DESCRIPTION:
;   Calculate the longitude and latitude for occultation tracks
;
; CALLING SEQUENCE:
;   oc_tracklonlat(et,targ,radec_star, r, lon, lat, found, fg0=fg0, $
;      obs=obs,abcorr=abcorr,frame=frame, deg=deg)
;
; INPUTS:
;   et - the ET at which you want to know the coordinates (may be an
;        array).  Units are seconds of TDB (barycentric dynamical
;        time) after J2000.
;   targ = target code
;   radec_star - vector of;
;    A - right ascention in radians
;    D - declination in radians
;   r - km north (south if r < 0) of track centerline
; OPTIONAL INPUT PARAMETERS:
;
; KEYWORDS
;   fg0 - final fg = nominal - fg0
;   deg - if set, return lon, lat in decimal degrees
;   obs - the keyword argument that indicates the observer location, 
;         default to Earth center (399).
;   abcorr  - aberation correction.  Default 'LT'
;   frame   - frame of state, ra, dec.  Default 'J2000'
;
; OUTPUTS:
;   lon - longitude
;   lat - latitude
;   found - 1 if lat, lon on Earth
; 
;  REVISON HISTORY:
;   2006 Jan 11 Leslie Young.  Based on single_eph
;-

pro oc_tracklonlat, et,targ,radec_star, r, lon, lat, found, fg0=fg0, $
               obs=obs,abcorr=abcorr,frame=frame, deg=deg
   
   ; set keywords -------------------------------------------------------
   if not keyword_set(obs) then obs = 399 ; Earth
   if not keyword_set(abcorr) then abcorr = 'LT' ; light-time only
   if not keyword_set(frame) then frame = 'J2000' ; 
   if not keyword_set(fg0) then fg0 = [0.d, 0.d]
   
   n = n_elements(et)
   lon = fltarr(n)
   lat = fltarr(n)
   found = bytarr(n)
   for i=0, n-1 do begin
       llactr = oc_fg2lla(et[i],targ, radec_star, [0.d,0.d]+fg0, foundi)
       found[i] = foundi
       if foundi then begin
           fgh6 = oc_naiffgh(et[i],targ,radec_star, lonlatalt=llactr)
           fgv = normv( fgh6[3:4] )
           dfg = [-fgv[1], fgv[0] ] ; delta f,g
           if dfg[1] lt 0.d then dfg = -dfg
           fg = r * dfg
           lla = oc_fg2lla(et[i],targ, radec_star, fg+fg0, foundi)
           if foundi then begin
               lon[i] = lla[0]
               lat[i] = lla[1]
           endif
       endif
   endfor

   if keyword_set(deg) then begin
       lon = lon * 180.d / !dpi
       lat = lat * 180.d / !dpi
   end
   
end