FUNCTION sgp4_tle, tt_obs, body_, source=source, tle=tle, tt_tle=tt_tle, use_tle=use_tle
;+
; NAME:
;	sgp4_tle
; PURPOSE:
;	Extract orbital elements covering bracketing tt_obs from
;	tle data base file
; CATEGORY:
;	gen/idl/ephem
; CALLING SEQUENCE:
;	status = sgp4_tle(tt_obs,body_,source=source,tle=tle,tt_tle=tt_tle)
; INPUTS:
;	tt_obs			array[n]; type: time structure
;						UT time
; OPTIONAL INPUTS:
;	body_			scalar; type: string; default: 'sat27640'
;						spacecraft designation
;						The default sat27640 is Coriolis.
;	source=source	scalar; type: string; default: who_am_i(/dir)/sgp
;						directory where TLE file is located (by default the
;						subdir sgp of the dir where this file is located.
; OUTPUTS:
;	status			scalar; type: integer
;						0: the ephemeris find is not found or there
;							is a read error. In this case tle and tt_tle
;							won't exist
;						1: tles found
;	tle=tle			array[2,ntle]; type: string
;						tle (each tle consists of two records)
;						the tles will bracket the input time period tt_obs
;	tt_tle			array[ntle]; type: time structure
;						times extracted from tle records
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	txt_read, InitVar, who_am_i, sgp_alias, BadValue, hide_env
;	TimeSet, TimeOp, TimeUnit, TimeLimits, TimeFixYear, TimeOp
;	SuperArray, ArrayLocation
; PROCEDURE:
;	The file with orbital elements is 'body.txt' i.e. sat27640.txt for 
;	Coriolis. The file is looked for in the subdirectory sgp of the
;	directory where this source code is located. Used keyword 'source'
;	to point to another subdirectory.
;	Orbital elements can be retrieved from www.space-track.org
; MODIFICATION HISTORY:
;	DEC-2005, Paul Hick (UCSD/CASS)
;		Extracted from sgp4_eph
;-

InitVar, source, filepath(root=who_am_i(/dir),'sgp') 
InitVar, body_, 'sat27640', set=body
InitVar, use_tle, /key

; Translate aliases

body = sgp_alias(body)

tmp = filepath(root=source,body+'.txt')
IF NOT txt_read(tmp,ch_tle,/silent) THEN BEGIN
	message, /info, 'no ephemeris file, '+hide_env(tmp)
	RETURN, 0
ENDIF

ch_tle 	= ch_tle[5:n_elements(ch_tle)-2]

ntle	= n_elements(ch_tle)/2
ch_tle 	= reform(ch_tle,2,ntle,/overwrite)
tt_tle	= TimeSet(yr=TimeFixYear(long( reform(strmid(ch_tle[0,*],18,2)))),doy=double(reform(strmid(ch_tle[0,*],20,12))))

InitVar, tt_obs, tt_tle
nobs	= n_elements(tt_obs)
szobs	= size(tt_obs, /dim)
tt_obs	= reform(tt_obs,nobs,/overwrite)

umin = TimeUnit(/minute)

; Extract the TLE vectors bracketing the requested time period
; (strictly speaking not necessary, but reduces the size of uu_obs below
; usually by a lot if the ephemeris file gets big)

i = where( TimeOp(/subtract,tt_tle,TimeLimits(tt_obs,/min),umin) 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_tle,TimeLimits(tt_obs,/max),umin) GT 0 )
IF j[0] EQ -1 THEN j = ntle-1 ELSE j = j[0] 			; Follows time range

ch_tle = ch_tle[*,i:j]
tt_tle = tt_tle[  i:j]
ntle   = n_elements(tt_tle)

CASE 1 OF

use_tle: BEGIN

	nobs   = ntle
	tt_obs = tt_tle
	dtt    = TimeOp(/subtract,tt_obs,tt_tle,umin)

END

ntle EQ 1: BEGIN					; Happens if all tt_obs are outside the
									; time range covered by the ephemeris
	ch_tle = SuperArray(ch_tle,nobs,/trail)
	dtt = TimeOp(/subtract,tt_obs,tt_tle,umin)

END

ELSE: BEGIN

	; Create array [nobs,ntle] for obs and tle time to analyze time differences

	uu_obs	= replicate(!TheTime,nobs,ntle)
	uu_tle	= replicate(!TheTime,nobs,ntle)

	uu_obs.(0) = tt_obs.(0)#replicate(1,ntle)
	uu_obs.(1) = tt_obs.(1)#replicate(1,ntle)

	uu_tle.(0) = replicate(1,nobs)#tt_tle.(0)
	uu_tle.(1) = replicate(1,nobs)#tt_tle.(1)
										; Time diff obs - tle in minutes
	dtt 	= TimeOp(/subtract,uu_obs,uu_tle,umin)	; Time diff obs - tle in minutes

	tmp 	= min(abs(dtt), dim=2, i)
	i		= ArrayLocation(i,size=size(dtt))
	i		= reform(i[1,*])

	ch_tle 	= ch_tle[*,i]				; tle nearest to obs time
	dtt 	= dtt[lindgen(nobs),i]		; Time diff of obs with nearest tle

END

ENDCASE

; rr(2) and rr(4) are read with format F6.5. rr(8) with F7.7
; The corresponding numbers in the TLE record do not contain an explicit
; decimal point. E.g. rr(2) and rr(4) could be 54147. rr(8) could be 0013639.
; IDL 6.1 reads this as 54157.0 instead of 0.54157 (it works correctly in
; Fortran). Rather than messing with the format the numbers are explicitly
; corrected below.

format = '18X,D14.8,1X,F10.8,2(1X,F6.5,I2),/,7X,2(1X,F8.4),1X,F7.7,2(1X,F8.4),1X,F11.8'

tmp = {tle_struct, EPOCH: 0d0, XNDT2O: 0d0, XNDD6O: 0d0, IEXP: 0d0, BSTAR: 0d0, 	$
	IBEXP: 0d0, XINCL: 0d0, XNODEO: 0d0, EO: 0d0, OMEGAO: 0d0, XMO: 0d0, XNO: 0d0}

tle = replicate(tmp,nobs)
FOR i=0,nobs-1 DO BEGIN
	reads, ch_tle[*,i], format='('+format+')', tmp
	tle[i] = tmp
ENDFOR
;reads, ch_tle, format='('+format+')', tle

tle.XNDD6O	= tle.XNDD6O*1d-5		; Fix read error
tle.BSTAR	= tle.BSTAR *1d-5		; Fix read error
tle.EO		= tle.EO    *1d-7		; Fix read error

RETURN, 1  &  END