;+
; NAME:
;	CarringtonT0
; PURPOSE:
;	Calculate list of start times for subsequent 'Carrington rotations'
; CATEGORY:
;	Celestial mechanics
; CALLING SEQUENCE:
	FUNCTION CarringtonT0, T, ScName, ddoy = dDoy, ncar = nCar,		$
		eastlimb = EastLimb, westlimb = WestLimb, n_c = N_C, f_c = F_C
; INPUTS:
;	T		array; type: time structure
;				times
;	ScName	scalar; type: string; default: 'ScEarth'
;				function name; identifies spacecraft
;	dDoy	scalar; type: float
;				required accuracy (in days) for the Carrington start times JDs
;	nCar	scalar; type: integer
;				# start times to be calculated (default=1)
; OUTPUTS:
;	Ts		array; type: standard time structure
;					Carrington start times
; OPTIONAL OUTPUTS:
;	n_c=N_C	scalar; type: integer
;				Rotation for which Ts[0] is the start time
;				(this is the rotation containing T)
;	f_c=F_C	fraction of rotation between input T and Ts[0]
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, TimeOp, TimeUnit, TimeSet, AngleRange, ScEarth
;	ArrayLocation
; PROCEDURE:
; >	ScName is a user-written function which calculates the heliographic
;		longitude of a spacecraft for a given time T
;	The call to ScName has the form:
;		LNG = ScName(T, degrees=degrees)
;	where T, a standard time structure is input and LNG (float) is output.
;	The function should not modify T in any way.
; >	The spacecraft is supposed to move in the ecliptic, circling the Sun in
;	the same direction as Earth (direct motion)
; >	The start time of a new Carrington rotation is defined as the time for
;	which the heliographic longitude of the spacecraft is zero.
; >	the input T is the time where the search for start times begins.
;	the first nCar Carrington start times after T, including
;	the one containing T, are calculated)
; MODIFICATION HISTORY:
;	JAN-1992, Paul Hick (UCSD)
;	NOV-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu); all times are now
;		time structures (used to be yr, doy or Julian days).
;-

InitVar, ScName , ''
InitVar, nCar	, 1
InitVar, dDoy	, 0.0001d0

dDoy = double(dDoy[0])

mCar = nCar+2				; Calculate extra start times as safety margin

; Reform to one-dim array. The structure is restored before returning.
; Note that T is a structure, so it is never a scalar.

szT = size(T)
nT  = szT[szT[0]+2]
T   = reform(T,nT)

; The standard Earth-based Carrington rotations can be calulated faster.

IF ScName EQ '' THEN BEGIN
	PSun = !sun.synodicp

	Tref = TimeSet(jd=2445871, sec=78624)		; JD 2445871.91d0)

	N_C = [round(1750.0d0+TimeOp(/subtract,T,Tref,TimeUnit(/day))/PSun)]; Guess at rotation # containing T
	N_C = N_C-1							; The actual rotation could be one earlier,
										; 	so step back one rotation
	; Set up 2D array N_C with rotation numbers. The first dimension (nDat)
	; corresponds to the input dates. The 2nd dimension (mCar) contains successive
	; rotation numbers, for nDat=3 and mCar=4 the array would look like:
	;	1800	1859	1888
	;	1801	1860	1889
	;	1802	1861	1890
	;	1803	1862	1891

	N_C = N_C#replicate(1,mCar)+replicate(1,nT)#indgen(1,mCar)
	Ts = TimeOp(/add, Tref, TimeSet(day=(N_C-1750)*PSun))	; Approximate start times

	n  = indgen(nT*mCar)
	nl = n
	dl = lonarr(nT*mCar)

	REPEAT BEGIN

		n = n[nl]
		Ts[n] = TimeOp(/add, Ts[n], TimeArray(dl[n],/linear))

		dl = ScEarth(Ts[n], eastlimb=EastLimb, westlimb=WestLimb)
		dl = AngleRange(dl, /pi)/(2*!dpi)*!sun.synodicp
		dl = TimeArray(TimeSet(day=dl),/linear)

		nl = where(dl NE 0)

	ENDREP UNTIL nl[0] EQ -1

ENDIF ELSE BEGIN

	; The following calculation is much slower than the previous one, but it
	; should work for any spacecraft in a direct motion, heliocentric orbit
	; (i.e. the Helios spacecraft) if the appropriate ScName function is available.

	; PSun is sideric solar rotation period at latitude 16 deg (in days)
	; (Does not have to be very accurate. Make sure to underestimate it.)

	PSun = !sun.siderealp-0.4

	Del	 = 0.1*PSun						; Time step size

	Lng0 = fltarr(nT)
	N_C  = intarr(nT)

	Ts   = replicate(!TheTime, nT, mCar); nT x mCar array

	Ts[*,0] = TimeOp(/subtract, T, TimeSet(day=PSun))	; Safety margin

	; Collect mCar start times for subsequent Carrington rotations

	FOR nn=0,mCar-1 DO BEGIN

		; In the next statement the assumption that the spacecraft is in a direct
		; orbit is used. The new Ts must be set earlier than the start time we need.

		if nn ne 0 then Ts[*,nn] = TimeOp(/add, Ts[*,nn-1], TimeSet(day=PSun))

		T1   = Ts[*,nn]
		Lng1 = call_function(ScName, T1)

		; The spacecraft heliographic longitude is a decreasing function of time
		; within one Carrington rotation (the Sun catches up with the spacecraft).

		; Step forward in time from Ts with step size Del until the heliographic
		; longitude jumps from (almost) 0 at time T0 to (almost) 360 at time Ts.

		n  = indgen(nT)
		nl = n
		REPEAT BEGIN
		    n = n[nl]
		    Lng0[n] = Lng1[n]
		    T1  [n] = TimeOp(/add, T1[n], TimeSet(day=Del))
		    Lng1[n] = call_function(ScName, T1[n])
		    nl = where(Lng1[n] LT Lng0[n])
		ENDREP UNTIL nl[0] EQ -1

		; A new Carrington rotation starts less then Del days before Ts. Refine the
		; determination of the start time (until it is better than dDoy) by repeated
		; bisection on the interval [T0,T1] = [Ts-Del,Ts]

		T0 = TimeOp(/subtract, T1, TimeSet(day=Del))

		n = indgen(nT)
		nl = n
		REPEAT BEGIN
			n = n[nl]
			Ts[n,nn] = TimeOp(/meant, T0[n], T1[n])
			Lng = call_function(ScName, Ts[n,nn])

			nl = where(Lng LE Lng0[n])
			IF nl[0] NE -1 THEN begin
				T0  [n[nl]] = Ts [n[nl],nn]
				Lng0[n[nl]] = Lng[nl]
			ENDIF

			nl = where(Lng GE Lng1[n])
			IF nl[0] NE -1 THEN BEGIN
				T1  [n[nl]] = Ts [n[nl],nn]
				Lng1[n[nl]] = Lng[nl]
			ENDIF

			nl = where(TimeOp(/subtract, T1[n], T0[n], TimeUnit(/day)) GT dDoy)
		ENDREP UNTIL nl[0] EQ -1

		; Take the center of the final interval [T0,T1] as the final estimate for
		; the start time of a new Carrington rotation.

		Ts[*,nn] = TimeOp(/meant, T0, T1)

	ENDFOR

ENDELSE

; JD0 are the Julian days corresponding to the input T (array structure S)
; JDs are start times for Carrington rotations N_C (structure [S,mCar])

JDs = TimeGet(Ts, /njd)
JD0 = TimeGet(T , /njd)

; Find the rotations with start times > JD0. Collect the rotation numbers in N_C[S]

i  = -1
nl = replicate(i,nT)

REPEAT BEGIN
	i = i+1
	n = where( JDs[*,i] LE JD0)
	IF n[0] NE -1 THEN nl[n] = i
ENDREP UNTIL n[0] EQ -1 OR i EQ mCar-1

IF (where(nl EQ -1))[0] NE -1 THEN message, 'Something wrong in paradise

; nl is an array with structure S, each element s has a value nl[s] with
; 0<=nl[s]<=mCar-1 such that JDs[s,nl[s]] <= JD0[s] and JDs[s,nl[s]+1] > JD0[s]

N_C = N_C[*,0]+nl					; Rotation # for 1st of the nCar rotations
									; x is 1-dim index into Ts array
nl = ArrayLocation( transpose([[indgen(nT)],[nl]]), size=size(Ts), /onedim)
				 					; [S,2] array bracketing rotations containing T
IF ScName EQ '' THEN BEGIN			; Lng[x] has structure [S,2]
	LngIn = call_function('ScEarth', T, eastlimb=EastLimb, westlimb=WestLimb)
	Lng   = call_function('ScEarth', Ts[[[nl],[nl+nT]]],eastlimb=EastLimb,westlimb=WestLimb)
ENDIF ELSE BEGIN
	LngIn = call_function(ScName, T)
	Lng   = call_function(ScName, Ts[[[nl],[nl+nT]]])
ENDELSE

Lng = AngleRange(Lng, /pi)			; Longitudes at start of rotation  should be close to zero
									; Fraction of rotation from T to start of 1st rotation
F_C = AngleRange(Lng[*,0]-LngIn)/(Lng[*,0]-Lng[*,1]+2*!dpi)

nl  = nl#replicate(1,nCar)+replicate(1,nT)#(indgen(nCar)*nT)
Ts  = Ts[nl]

szT = szT[1:szT[0]]
T   = reform(T  , szT, /overwrite)
N_C = reform(N_C, szT, /overwrite)
F_C = reform(F_C, szT, /overwrite)

IF nCar EQ 1 THEN TS = reform(Ts,szT, /overwrite) ELSE TS = reform(Ts,[szT,nCar], /overwrite)

IF nT EQ 1 THEN BEGIN			; Only one time: remove leading dimension of one
	N_C = N_C[0]
	F_C = F_C[0]
	Ts  = reform(Ts, /overwrite)
ENDIF

RETURN, Ts  &  END