;+
; NAME:
;	CarringtonT
; PURPOSE:
;	Get times when Earth was at Carrington variable XC
; CALLING SEQUENCE:
	FUNCTION CarringtonT, xc
; INPUTS:
;	XC			array; type: float
; OUTPUTS:
;	T			array; type: standard time structure
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	TimeOrigin, TimeOp, TimeArray, TimeSet, AngleRange
; PROCEDURE:
;	An estimate for the Carrington start time is set up first.
;	This is refined by a call to CarringtonT0.
; MODIFICATION HISTORY:
;	SEP-1999, Paul Hick (UCSD/CASS)
;	JAN-2004, Paul Hick (UCSD/CASS)
;		Substantial rewrite to improve precision from several minutes
;		to about a milli-second
;	APR-2004, Paul Hick (UCSD/CASS)
;		Added code to prevent iteration loop to get stuck in infinite loop
;	JUN=2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Bug fix (du was subscripted with n instead of nn)
;-

units = !TimeUnits.in_day
TimeOrigin, /time, units_in_day=86400000L

T = TimeSet(jd=2445871, sec=78624)		; JD 2445871.91d0)
T = TimeOp(/add, T, TimeSet(day=(xc-1750)*!sun.synodicp))

TwoPi = 2.0d0*!dpi
HLng = TwoPi*(1.0d0-(xc-floor(xc)))

; n is a list of times that have not converged yet.
; Start out with the full list of all times.
; dt are the time differences that are pushed to zero.

n  = lindgen(n_elements(T))
dt = 0*n

; After each iteration nn is a subset of indices into n that have
; not converged yet. Initialize to full list.

nn = n
count = 12

REPEAT BEGIN

	n  = n [nn]
	du = dt[nn]

	T[n] = TimeOp(/add, T[n], TimeArray(du,/linear))

	dt = AngleRange(ScEarth(T[n])-HLng[n],/pi)/TwoPi*!sun.synodicp
	dt = TimeArray( TimeSet(day=dt), /linear )

	; If dt = 0 we found the correct time.
	; If du = -dt the iteration loop will start bouncing between
	; two times. If abs(dt) = 1 and du = -dt we also accept the time
	; as the correct one 

	nn = dt EQ 0 OR (abs(dt) EQ 1 AND dt EQ -du)
	nn = where(1B-nn)

ENDREP UNTIL nn[0] EQ -1 OR count-- LT 0

; A typical count should be 3 or 4. Abort if we actually need 12
; (i.e. three times the expected count).
			
IF count LT 0 THEN message, 'iteration count exceeded'

T = TimeOp(/units, T, old_units=86400000L, new_units=units)

TimeOrigin, /time, units_in_day=units

RETURN, T  &  END