;+
; NAME:
;	vu_timeseries
; PURPOSE:
;	Extract a time series at Earth from a 3D heliospheric matrix obtained with
;	one of the tomography programs
; CATEGORY:
;	WWW
; CALLING SEQUENCE:
	FUNCTION vu_timeseries, XC3D, R3D, dR3D, F3D, $
		ut		= ut		, $
		range	= range		, $
		delt	= delt		, $
		T		= T			, $
		periodic= periodic	, $
		sequence= sequence	, $
		radius	= radius	, $
		nearL1	= nearL1	, $
		power	= power		, $
		body	= body		, $
		_extra	= _extra
; INPUTS:
;	The 3D heliospheric matrix is described by the following input:
;	XC3D			array[2] or array[2,ntime]; type: float
;						Carrington range of reconstruction
;	R3D				scalar or array[ntime]; type: float
;						heliocentric distance at inner boundary of matrix
;	dR3D			scalar or array[ntime]; type: float
;						heliocentric distance resolution of the heliospheric matrix
;	F3D				array[n,m,l,ntype,ntime]; float
;						3D heliocentric matrix;
;						1st dim: Carrington variable, covering range XC3D
;						2nd dim: Heliographic latitude, covering [-90,90] degrees
;						3rd dim: Heliocentric distance, with inner boundary at R3D and resolution dR3D
;						4th dim: different data types (usually 2: velocity/density, Br/Bt
;						5th dim: # members in a (time-dependent) sequence
; OPTIONAL INPUT PARAMETERS:
;	sequence=sequence
;					array[ntim]; type: time structure of float; default: none
;						if specified, these are the times corresponding to the
;						ntime data arrays specified in argument F3D.
;	T=T				array; type: time structure or float
;						times at which time series is evaluated. A float array
;						is interpreted as an array of Carrington times. If T is
;						not set then the times are set up using the 'range' and
;						'ut' keywords.
;	range=range		array[2]; type: float; default: XC3D[*,0]
;						(NOT used if explicit times are specified using the T
;						keyword) the Carrington range (start and end Carrington
;						variable of time series to be extracted)
;	ut=ut			scalar; type: integer
;					array[1]; type: time structure
;						(NOT used if explicit times are specified using the T
;						keyword). Defines the time (UT) used to center the time
;						series. A scalar is interpreted as the difference
;						between local time and UT, and is used to get a UT from
;						the system time. The time series corresponds to
;						UT +/- half a the range specified in keyword 'range'.
;	delt=delt		scalar; type: float
;						time steps in time series (days)
;						if not specified then the time step is determined from
;						the resolution in Carrington variable of F3D
;	radius=radius
;
;	body=body		scalar; type: string
;						usually set to return value of function href=jpl_body=.
;						if SET then the JPL ephemeris is used to calculate the
;						heliospheric locations where the time series is extracted.
;						if NOT set then function href=big_eph= is used to extract
;						a time series at Earth (i.e. NOT setting 'body' should give the
;						save results as setting body=jpl_body(/earth)). 
;						Note the keyword /nearL1 only makes sense if body=jpl_body(/earth)
;
;	/nearL1			indicates whether time series is for spacecraft at L1;
;					if set the heliocentric distance calculated from
;					big_eph is multiplied with 0.99)
;
;	power=power		array[ntype] or array[ntype,ntime]
;						power of heliocentric distance to be removed from the time series
;						(0 for velocity, 2 for normalized density;
;						 2 for radial component of the magnetic field; 1 for tangential component)
; OUTPUTS:
;	F				array[*]; float
;						the function values for the time series
; OPTIONAL OUTPUT PARAMETERS:
;	T				array[*]; time structure
;						times at which function values where calculated
;						1st and last time correspond to the selected time range
; INCLUDE:
	@compile_opt.pro					; On error, return to caller
; CALLS:
;	InitVar, IsTime, IsType, Carrington, SuperArray
;	InterpolateHeliosphere, gridgen, BadValue, TimeOp, TimeUnit, big_eph
; PROCEDURE:
; >	The time associated with a given Carrington variable is the time at which the
;	corresponding heliographic longitude crossed the center of the solar disk.
; >	The times (or equivalent Carrington variables) to be extracted can be set in several ways:
;	1. Explicitly set the times array T
;	2. Set 'range' to a two-element array with start and end time or Carrington variable.
;		Not setting 'range' is the same as setting 'range' to XC3D
;	3. Setting 'range' to a scalar fraction of a rotation is the same as setting
;		range = mean(XC3D)+0.5*range*[-1,1], i.e. range=0.5 picks one rotation centered
;		on the mean of XC3D.
; > Setting 'ut' shifts 'range' and centers it on 'ut'
; > A time array with the specified step size is set up, and converted to an array of
;	Carrington variable. The function values are obtained by interpolating on the 3D matrix.
;
; > The heliocentric distance by default is the Sun-Earth distance. This can be modified
;	using the 'radius' and 'nearL1' keywords. 'radius' can be used to explicitly
;	set the distance (usually this will be a scalar). Setting /nearL1 multiplies
;	the heliocentric Sun-Earth distance with 0.99.
;
; MODIFICATION HISTORY:
;	OCT-1999, Paul Hick (UCSD/CASS)
;	MAY-2001, Paul Hick (UCSD/CASS
;		added /nearL1 keyword)
;	APR-2004, Paul Hick (UCSD/CASS)
;		activated keyword 'body' to enable use of JPL ephemeris.
;	OCT-2006, Paul Hick (UCSD/CASS)
;		Replaced jp_eph by big_eph and made keyword body a string.
;	    This allows access to all bodies in the list returned by
;	    the function big_body()
;	FEB-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Replaced longitudes keyword by cv2carrington keyword in
;		InterpolateHeliosphere calls.
;-

InitVar, nearL1, /key

sz   = size(F3D)
uday = TimeUnit(/day)

; First set up the times to be used for the time series (TT = UT times)

CASE IsType(T, /defined) OF

0: BEGIN

	CASE n_elements(range) OF
	0   : dxc = XC3D[*,0]
	1   : dxc = mean(XC3D[*,0])+0.5*range[0]*[-1,1]
	ELSE: dxc = range[0:1]
	ENDCASE

	IF IsType(ut, /defined) THEN	$
		dxc = Carrington(ut, /get_variable)+(dxc-mean(dxc))

	CASE IsType(delt, /defined) OF
	0: nt = sz[1]
	1: nt = round((dxc[1]-dxc[0])*!sun.synodicp/delt)
	ENDCASE

	TT = Carrington(gridgen(nt, range=dxc))
	T  = TT 						; Set return argument
END

1: TT = Carrington(T,/get_time)

ENDCASE


CASE sz[0] OF
3: BEGIN
	ntype = 1					; One matrix: one data type, one time
	ntime = 1
END
4: BEGIN
	IF sz[4] EQ 2 THEN BEGIN
		ntype = sz[4]
		ntime = 1
	ENDIF ELSE BEGIN
		ntype = 1
		ntime = sz[4]
	ENDELSE
END
5: BEGIN
	ntype = sz[4]				; Multiple data types; multiple times
	ntime = sz[5]
END
ENDCASE

nT = n_elements(TT)

F3D = reform(F3D, [sz[1:3], ntype, ntime], /overwrite)


xc = XC3D
IF n_elements(xc) EQ 2 THEN xc = SuperArray(xc, ntime, /trail)

rr = R3D
IF n_elements(rr) EQ 1 THEN rr = replicate(rr[0], ntime)

drr = dR3D
IF n_elements(drr) EQ 1 THEN drr = replicate(drr[0], ntime)

nseq = n_elements(sequence)		; Should be the same as ntime if specified
at_earth = IsType(body, /undefined)

CASE at_earth OF

0: BEGIN

	bodies = big_body(count=count)
	FOR i=0,count-1 DO IF strpos(body,bodies[i]) EQ 0 THEN break

	CASE i EQ count OF
	0: _body = bodies[i]
	1: BEGIN
		_body = flt_string(body,numfmt=i)
		CASE i OF
		1: _body = [_body,0,0]
		2: _body = [_body  ,1]
		3:
		ELSE: message, 'no valid big_eph body, and no location specified: '+body
		ENDCASE
		message, /info, 'location: Lng,Lat,Rad='+strjoin(_body)
		_body[0:1] = AngleUnits(from_degrees=_body[0:1], /to_radians)
	END
	ENDCASE

END

1: _body = jpl_body(/earth,/string)

ENDCASE


IF nseq EQ 0 THEN BEGIN

	; A time series at times TT is extracted from each of the ntype*ntime
	; data arrays in F3D

	CASE IsType(_body,/string) OF

	0: R = CvSky(TT, _body, _extra=_extra, /to_heliographic)

	1: BEGIN
		R = big_eph(TT		, $
			body=_body		, $
			/to_sphere		, $
			/to_heliographic, $
			/precess		, $
			/onebody		, $
			/silent			)
		IF at_earth THEN R[0,*] = Carrington(TT)
	END

	ENDCASE

	CASE IsType(radius, /defined) OF
	0: IF nearL1 THEN R[2,*] = 0.99*R[2,*]		; L1 is at 0.99 x Sun-Earth distance
	1: R[2,*] = radius
	ENDCASE

	F = make_array( type=IsType(F3D), dim=[nT, ntype, ntime] )

	FOR n=0,ntime-1 DO				$
		FOR i=0,ntype-1 DO			$
			F[*,i,n] = InterpolateHeliosphere(R, F3D[*,*,*,i,n], rr[n], drr[n], $
				xcrange=xc[*,n], periodic=periodic, $
				cv2carrington=1-at_earth, is_carrington=at_earth)

	; Here we convert from 'normalized' quantities (read from file) to 'real' quantities:
	;	nr^2 --> n;		r^2Br --> Br;	rBt --> Bt
	; The 'power' values are usually extracted from a file header.

	IF IsType(power, /defined) THEN BEGIN

		CASE size(power, /n_dim) OF
		0:	FOR n=0,ntime-1 DO			$
				FOR i=0,ntype-1 DO		$
					IF power      NE 0 THEN F[*,i,n] = F[*,i,n]/R[2,*]^power
		1:	FOR n=0,ntime-1 DO			$
				FOR i=0,ntype-1 DO		$
					IF power[i  ] NE 0 THEN F[*,i,n] = F[*,i,n]/R[2,*]^power[i]
		2:	FOR n=0,ntime-1 DO			$
				FOR i=0,ntype-1 DO		$
					IF power[i,n] NE 0 THEN F[*,i,n] = F[*,i,n]/R[2,*]^power[i,n]
		ENDCASE

	ENDIF

ENDIF ELSE BEGIN

	IF ntime NE nseq THEN message, 'oops'

	; When processing a time sequence of ntime data arrays, first extract 'ntype' time 
	; series at the sequence times. Then interpolate on those time series to get the time
	; series at the requested times TT.

	tseq = Carrington(sequence, /get_time)

	CASE IsType(_body,/string) OF

	0: R = CvSky(tseq, _body, _extra=_extra, /to_heliographic)

	1: BEGIN
		R = big_eph(tseq	, $
			body=_body		, $
			/to_sphere		, $
			/to_heliographic, $
			/precess		, $
			/onebody		, $
			/silent			)
		IF at_earth THEN R[0,*] = Carrington(sequence, /get_variable)
	END

	ENDCASE

	CASE IsType(radius, /defined) OF
	0: IF nearL1 THEN R[2,*] = 0.99*R[2,*]			; L1 is at 0.99 x Sun-Earth distance
	1: R[2,*] = radius
	ENDCASE

	G = make_array( type=IsType(F3D), dim=[ntype, ntime] )

	FOR n=0,ntime-1 DO			$
 		FOR i=0,ntype-1 DO		$
			G[i,n] = InterpolateHeliosphere(R[*,n], F3D[*,*,*,i,n], rr[n], drr[n],	$
				xcrange=xc[*,n], periodic=periodic, $
				cv2carrington=1-at_earth, is_carrington=at_earth)

	; Here we convert from 'normalized' quantities (read from file) to 'real' quantities:
	;	nr^2 --> n;		r^2Br --> Br;	rBt --> Bt
	; The 'power' values are usually extracted from a file header.

	IF IsType(power, /defined) THEN BEGIN

		CASE size(power, /n_dim) OF
		0:	FOR i=0,ntype-1 DO		$
				IF power    NE 0 THEN G[i,*] = G[i,*]/R[2,*]^power
		1:	FOR i=0,ntype-1 DO		$
				IF power[i] NE 0 THEN G[i,*] = G[i,*]/R[2,*]^power[i]

		; Here we only use power[0..ntype-1] i.e. the same power is used for all sequence times.

		2:	FOR i=0,ntype-1 DO		$
				IF power[i] NE 0 THEN G[i,*] = G[i,*]/R[2,*]^power[i]
		ENDCASE

	ENDIF

	; Now interpolate at times TT to get the requested time series.

	F = replicate(BadValue(F3D), nT, ntype)

	good = where(TimeOp(/subtract,TT,tseq[0],uday) GE 0 AND TimeOp(/subtract,TT,tseq[ntime-1],uday) LE 0)

	IF good[0] NE -1 THEN BEGIN
		TT = TT[good] 
		FOR i=0,ntype-1 DO F[good,i] = TimeInterpol(reform(G[i,*]), tseq, TT)
	ENDIF


ENDELSE

F3D = reform(F3D, sz[1:sz[0]], /overwrite)

RETURN, reform(F)  &  END