;+
; NAME:
;	TimeLInterpol
; PURPOSE:
;	Linear interpolation on a function of time
; CATEGORY:
;	gen/idl/toolbox/time
; CALLING SEQUENCE:
	FUNCTION TimeLInterpol, ff, tt, t, unit, $
		scalar	= scalar	, $
		inear	= inear		, $
		dtnear	= dtnear	, $
		ilow	= ilow		, $
		dtlow	= dtlow		, $
		fraction= fraction	, $
		ihigh	= ihigh		, $
		dthigh	= dthigh
; INPUTS:
;	ff			array[n]; type: any
;					function values
;					if ff does not exist then the array index
;						lindgen(n) is used
;	tt			array[n]; type: time structure
;					UT times
;	t			array[m]; type: time structure
;					times were interpolated values are needed
; OPTIONAL INPUTS:
;	/scalar 	if n=1 then return values will be scalars instead
;				of array[1]
; OUTPUTS:
;	inear	= inear
;	dtnear	= dtnear
;	ilow	= ilow
;	dtlow	= dtlow
;	ihigh	= ihigh
;	dthigh	= dthigh
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, TimeOp, TimeUnit, TimeLimits, ArrayLocation
; PROCEDURE:
; MODIFICATION HISTORY:
;	DEC-2005, Paul Hick (UCSD/CASS)
;-

InitVar, scalar, /key
InitVar, unit, TimeUnit(/day)

nt = n_elements(t)
sz = size(t,/dim)
t  = reform(t,nt,/overwrite)

nl = n_elements(tt)

; Extract the time from the list bracketing the requested time period
; (strictly speaking not necessary, but reduces the size of uu_obs below
; usually by a lot if the list is big)

i = where( TimeOp(/subtract,tt,TimeLimits(t,/min),unit) LT 0 )
IF i[0] EQ -1 THEN i = 0 ELSE i = i[n_elements(i)-1]	; Preceeds time range

j = where( TimeOp(/subtract,tt,TimeLimits(t,/max),unit) GT 0 )
IF j[0] EQ -1 THEN j = nl-1 ELSE j = j[0]		; Follows time range

IF i EQ nl-1 THEN i = (i-1) > 0  
IF j EQ 0    THEN j = (j+1) < (nl-1)

; i:j brackets the whole time range t

offset = i

tl = tt[i:j]
nl = n_elements(tl)

CASE nl EQ 1 OF

0: BEGIN

	; Create array [nobs,ntlm] for obs and tlm time to analyze time differences

	vtt = replicate(!TheTime,nt,nl)
	vtl = replicate(!TheTime,nt,nl)

	vtt.(0) = t.(0)#replicate(1,nl)
	vtt.(1) = t.(1)#replicate(1,nl)

	vtl.(0) = replicate(1,nt)#tl.(0)
	vtl.(1) = replicate(1,nt)#tl.(1)

	dt = TimeOp(/subtract,vtt,vtl,unit); Time diff tt-tl

	a = min(abs(dt), dim=2, inear)
	inear = ArrayLocation(inear,size=size(dt))
	inear = reform(inear[1,*])		; Indices in tl for the nearest index

	dtnear = dt[lindgen(nt),inear]	; Time diff with nearest tl

END

1: BEGIN							; Happens if all tt are outside range of tt_list

	dt = TimeOp(/subtract,t,tl,unit)

	inear = replicate(0,nt)

	dtnear = dt

END

ENDCASE

a = 2*(dtnear GE 0)-1

ilow  = inear
ihigh = ((ilow+a) > 0) < (nl-1)

a = where(ilow GT ihigh)
IF a[0] NE -1 THEN BEGIN
	k = ilow[a]
	ilow [a] = ihigh[a]
	ihigh[a] = k
ENDIF

a = where(ilow EQ ihigh)
IF a[0] NE -1 THEN BEGIN
	k = where(ilow [a] EQ 0)
	IF k[0] NE -1 THEN ihigh[a[k]] = (ihigh[a[k]]+1) < (nl-1)
	k = where(ihigh[a] EQ nl-1)
	IF k[0] NE -1 THEN ilow [a[k]] = (ilow [a[k]]-1) > 0
ENDIF

a = lindgen(nt)

dtlow  = dt[a,ilow ]			; Time diff with nearest tl
dthigh = dt[a,ihigh]

fnc = dtlow
fraction = dtlow

a = where(ilow NE ihigh)

IF a[0] NE -1 THEN BEGIN
	CASE IsType(ff,/defined) OF
	0: fnc[a] = offset+(dtlow[a]*ihigh[a]-dthigh[a]*ilow[a] )/(dtlow[a]-dthigh[a])
	1: fnc[a] = (dtlow[a]*ff[offset+ihigh[a]]-dthigh[a]*ff[offset+ilow[a]])/(dtlow[a]-dthigh[a])
	ENDCASE
	fraction[a] = dtlow[a]*(ihigh[a]-ilow[a])/(dtlow[a]-dthigh[a])
ENDIF

a = where(ilow EQ ihigh)

IF a[0] NE -1 THEN BEGIN
	CASE IsType(ff,/defined) OF
	0: fnc[a] =    offset+ilow[a]
	1: fnc[a] = ff[offset+ilow[a]]
	ENDCASE
	fraction[a] = 0
ENDIF

inear = offset+inear
ilow  = offset+ilow
ihigh = offset+ihigh

inear  = reform( inear,sz,/overwrite)
dtnear = reform(dtnear,sz,/overwrite)

ilow   = reform( ilow ,sz,/overwrite)
dtlow  = reform(dtlow ,sz,/overwrite)

ihigh  = reform( ihigh,sz,/overwrite)
dthigh = reform(dthigh,sz,/overwrite)

fnc    = reform(fnc   ,sz,/overwrite)

IF scalar AND n_elements(t) EQ 1 THEN BEGIN
	inear  =  inear[0]
	dtnear = dtnear[0]

	ilow   =  ilow [0]
	dtlow  = dtlow [0]

	ihigh  =  ihigh[0]
	dthigh = dthigh[0]

	fnc    = fnc   [0]
ENDIF

RETURN, fnc  &  END