;+
; NAME:
;	Carrington
; PURPOSE:
;	Get times when Earth was at Carrington variable XC, and v.v
; CALLING SEQUENCE:
	FUNCTION Carrington, xc_or_t		, $
		fraction		= fraction		, $
		longitude		= longitude 	, $
		rotation		= rotation		, $
		get_fraction	= get_fraction	, $
		get_longitude	= get_longitude , $
		get_rotation	= get_rotation	, $
		get_time		= get_time		, $
		get_variable	= get_variable	, $
		near_longitude	= near_longitude, $
		degrees 		= degrees
; INPUTS:
;	xc_or_t			array; type: any numerical type
;						if input is a Carrington variable, the
;						return value is the corresponding UT time
;					array; type: standard time structure
;						if input is a UT time, the return value
;						is a Carrington variable.
; OPTIONAL INPUTS:
;	/get_fraction	return fraction of rotation
;	/get_longitude	return heliographic longitude
;	/get_rotation	return integer rotation number
;	/get_time		always return UT time (even if input already
;						is a UT time)
;	/get_variable	always return Carrington variable (even if
;						input already is a Carrington variable)
;
;	near_longitude=near_longitude
;					heliographic longitude in [0,360]
;
;					If specified this heliographic longitude is
;					translated into a Carrington variable within half
;					a rotation from input Carrington variable/time
;					xc_or_t. Note that the rules for the return value
;					are the same as without the near_longitude, i.e.
;					if xc_or_t is a Carrington variable then the
;					return value is a time structure, unless one of the
;					get_* keywords is set.
;
; OUTPUTS:
;	Result			If no /get_* keywords are set:
;						array; type: standard time structure
;						array; type: double
;					/get_fraction : fraction of rotation
;					/get_longitude: heliographic longitude
;					/get_rotation : integer Carrington rotation nr.
;					/get_time	  : UT time
;					/get_variable : Carrington variable
; OPTIONAL OUTPUTS:
;	fraction=fraction
;					array; type: double
;					positive fraction of a rotation
;	longitude=longitude
;					array; type: double
;					heliographic longitude
;	rotation=rotation
;					array; type: integer
;					integer Carrington rotation number
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	TimeOrigin, TimeOp, TimeArray, TimeSet, TimeUnit
;	AngleRange, jpl_eph, CvSky, SubArray
; PROCEDURE:
;	An estimate for the Carrington start time is set up first.
;	This is iteratively refined.
;
;	If keyword near_longitude is defined then
;	the difference in longitude between near_longitude and
;	xc_or_t is calculated. If the difference is less than
;	-180 deg then 360 deg is added; if greater than 180 deg
;	then 360 deg is subtracted.
; 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, and avoid calls to
;		CarringtonT0
;	APR-2004, Paul Hick (UCSD/CASS)
;		Added code to prevent iteration loop to get stuck in
;		infinite loop
;	JUN=2006, Paul Hick (UCSD/CASS)
;		Bug fix (du was subscripted with n instead of nn)
;	JUL-2007, Paul Hick (UCSD/CASS)
;		Merged carringtonvar, carringtont, carringtonlng,
;		carringtonnr and arg_time into this procedure.
;		Merged CarringtonNear by adding keyword near_longitude.
;	AUG-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Removed calls to scearth.
;-

OnePi = !dpi
TwoPi = 2.0d0*!dpi

rpm = ToRadians(degrees=degrees)


IF IsType(near_longitude,/defined) THEN BEGIN

	xc = Carrington(xc_or_t, /get_variable, longitude=lng)

	;-------
	; The decimal fraction of XC determines the heliographic longitude.
	;
	; Get the difference in heliographic longitude.
	; The point of observation is located east or west relative to the observer.
	; If it lies in the eastern hemisphere -180 <  DLNG <=   0;
	; If it lies in the western hemisphere    0 <= DLNG <  180.

	lng = near_longitude*rpm-lng		; Radians !!

	; If abs(lng) > 180 than add or subtract 360 degrees.
	; If lng >  180 subtract 360 making -180 <= lng < 0
	; If lng < -180 add      360 making    0 <= lng < 180

	lng += (abs(lng) GT OnePi)*(1-2*(lng GT 0))*TwoPi

	; Carrington at near_longitude

	xc -= lng/TwoPi

	RETURN, Carrington(xc, $
		fraction		= fraction		, $
		longitude		= longitude 	, $
		rotation		= rotation		, $
		get_fraction	= get_fraction	, $
		get_longitude	= get_longitude , $
		get_rotation	= get_rotation	, $
		get_time		= get_time		, $
		get_variable	= get_variable	, $
		degrees 		= degrees)

ENDIF


CASE 1 OF

IsTime(xc_or_t): BEGIN

	longitude = jpl_eph(xc_or_t,/location,/to_sphere,/precess,/to_ecliptic)
	longitude = CvSky(xc_or_t, from_ecliptic=longitude,/to_heliographic,/silent)
	longitude = SubArray(longitude)

	longitude = AngleRange(longitude)	; Heliographic longitude in [0,360]
	fraction  = (TwoPi-longitude)/TwoPi	; Exact positive fraction of one

	dd = TimeSet(jd=2445871, sec=78624)	; JD 2445871.91d0)
	dd = TimeOp(/subtract, xc_or_t, dd, TimeUnit(/day) )/!sun.synodicp
	rotation = floor(dd)
	dd -= rotation						; Approx positive fraction of one
	rotation += 1750

	dd = dd-fraction					; Error in fraction

	n = where(dd LT -0.5d0)
	IF n[0] NE -1 THEN rotation[n] = rotation[n]-1

	n = where(dd GT 0.5d0)
	IF n[0] NE -1 THEN rotation[n] = rotation[n]+1

	t_or_xc = rotation+fraction

	CASE 1 OF

	IsType(get_time 	,/defined): t_or_xc = xc_or_t

	IsType(get_rotation ,/defined): t_or_xc = rotation

	IsType(get_fraction ,/defined): t_or_xc = fraction

	IsType(get_longitude,/defined): t_or_xc = longitude/rpm

	ELSE						  : t_or_xc = rotation+fraction

	ENDCASE

END

ELSE: BEGIN

	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_or_t-1750)*!sun.synodicp))

	rotation  = floor(xc_or_t)
	fraction  = xc_or_t-rotation			; Positive fraction of one

	longitude = TwoPi*(1.0d0-fraction)		; Heliographic longitude in [0,360]

	CASE 1 OF

	IsType(get_variable ,/defined): t_or_xc = xc_or_t

	IsType(get_rotation ,/defined): t_or_xc = rotation

	IsType(get_longitude,/defined): t_or_xc = longitude/rpm

	IsType(get_fraction ,/defined): t_or_xc = fraction

	ELSE: BEGIN

		; 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 = jpl_eph(T[n],/location,/to_sphere,/precess,/to_ecliptic)
			dt = CvSky(T[n], from_ecliptic=dt,/to_heliographic,/silent)
			dt = SubArray(dt)

			dt = AngleRange(dt-longitude[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

		t_or_xc = T

	END

	ENDCASE

END

ENDCASE

longitude /= rpm

RETURN, t_or_xc  &  END