;+
;       NAME:
;               oc_ephicy_arclen.Pro
;
;       PURPOSE:
;		Returns a structure with the following information:
;                 {et, ra, dec, radot, decdot, ephinfo} where ephinfo
;                 contains the observer range and observer range rate.
;
;	DESCRIPTION:
;              Calculates the ra, dec, radot, decdot, observer range
;              and range rate using the NAIF ICY routines.
;
;       CALLING SEQUENCE:
;               oc_ephicy_arclen, target, observer, utcStrStart, utcStrEnd,
;               arclen, eph
;
;       INPUTS:
;              target - the target for the ephemeris, i.e. "Pluto"
;              observer - the observer location, i.e. "Earth"
;              utcStrStart - the UTC string of the start time,
;                            i.e. "January 11, 2006 00:00:00"
;              utcStrEnd - the UTC string of the end time
;              arclen - the increment in sky-plane distance traveled
;                       by the target, radians
;       OUTPUTS:
;              eph - a structure with the following elements
;              eph.ra - the Pluto RA in radians
;              eph.dec - the Pluto Dec in radians
;              eph.radot - the RA rate in radians/s
;                          (d(ra)*cos(dec)/dt)
;              eph.decdot -- the Dec rate in radians/s
;              eph.ephinfo - the observer range (in km)
;                                          and range rate (km/s).
;
;
;       REVISON HISTORY:
;       17-January-2005 CBO SwRI
;       26-January-2005 LAY SwRI added newton-raphson step
;       19-November-2005 LAY SwRI Changed from procedure to function
;       20-November-2005 LAY SwRI added call to rd2xieta to
;                            avoid xi=2pi when ra crosses zero 
;       02-July-2006 LAY SwRI add kernel keyword to single_eph  to avoid extra calls
;                             to naif_laodedkernels
;-

function oc_ephicy_arclen, target, observer, utcStrStart, utcStrEnd, arclen, $
                          neph, debug=debug

       naifinit

; Generate the time list
       CSPICE_STR2ET, utcStrStart, et0
       CSPICE_STR2ET, utcStrEnd, et1

       eph1 = single_eph( target, et0, neph1, observer=observer)
       kernels = eph1.naif.kernels

       eph = eph1
       neph = neph1

       arcmin = !dpi/(180.d*60.)   ; radian/arcmin

       time = et0
       
       WHILE time LT et1 DO BEGIN
           
           ra0=eph1.ra          ; LAY 2005/01/26 for iteration step
           dec0=eph1.dec        ; LAY 2005/01/26 for iteration step
           cosdec0 = cos(dec0)  ; LAY 2005/01/26 for iteration step
           
           time0 = eph1.et

           dt = arclen/SQRT(eph1.radot^2 + eph1.decdot^2)

           i = 0
           test = -1.

           WHILE test LT 0. DO BEGIN
               i = i + 1
               t = time + dt*i/10.
               IF t GT et1 THEN RETURN, eph
               eph1 = single_eph( target, t, neph1, observer=observer,ker=kernels )
               rd2xieta, ra0, dec0, eph1.ra, eph1.dec, xi, eta
               test = (xi^2 + eta^2 - arclen^2)
               if keyword_set(debug) then begin
                   print, t-time0, $
                     xi/arcmin, eta/arcmin, $
                     test, $
                     form='("OC_EPHICY_ARCLEN: (1) ", F8.0, 2E10.2, E10.2)'
               end
           ENDWHILE

           time = time + dt*i/10
           if keyword_set(debug) then begin
               print, time-time0, $
                 form='("OC_EPHICY_ARCLEN: (2) ", F8.0)'
           end

           iter = 0             ; LAY 2005/01/26 for iteration step
           eps = 1e-6           ; LAY 2005/01/26 for iteration step
           repeat begin ; BEGIN NEWTON-RAPHSON STEP LAY 2005/01/26
               IF time GT et1 THEN RETURN, eph
               eph1 = single_eph( target, time, neph1, observer=observer, ker=kernel)
               rd2xieta, ra0, dec0, eph1.ra, eph1.dec, xi, eta
               if keyword_set(debug) then begin
                   print, time-time0, $
                   xi/arcmin, eta/arcmin, eph1.radot/arcmin, $
                     eph1.decdot/arcmin, $
                     form='("OC_EPHICY_ARCLEN: (3) ", F8.0,2F6.2,2E10.2)'
               end
               f = (xi^2 + eta^2 - arclen^2)
               df = 2.*(xi*eph1.radot+eta*eph1.decdot) 
               time = time - f/df


               iter=iter+1
               if time lt time0 then begin
                   cspice_et2utc, time0,"C",3,utcstr
                   if keyword_set(debug) then print, utcstr
                   stop
               end
           endrep until abs(f)/arclen^2 lt eps or iter ge 20 
                                ; END NEWTON-RAPHSON STEP LAY 2005/01/26
           time = eph1.et       ; LAY 2005/01/26 cleanup after iteration step
           IF keyword_set(debug) then print, "OC_EPHICY_ARCLEN: (4) ",eph1.et
           eph=[eph,eph1]
           neph = neph + neph1
           
       ENDWHILE

       return, eph

 END

;+