;+
; NAME:
;	InsituTimeSeries
; PURPOSE:
;	Extract in-situ data in specified interval, and average if necessary
; CALLING SEQUENCE:
	FUNCTION InsituTimeseries,$
		data_request		, $
		t_request			, $
		t_given 			, $
		data_given			, $
		delt	= delt		, $
		weight	= weight	, $
		traceback=traceback , $
		source	= source_	, $
		nocvsky	= nocvsky	, $
		silent	= silent	, $
		cvlast	= cvlast
; INPUTS:
;	data_request	array, type: string
;						determines which data are extracted.
;						Identifiers are:
;						'lng','lat','dis','den','vel','br','bt','bn','bb','pdyn','vlng','vlat'
;						'lng'	= heliographic longitude (deg)
;						'lat'	= heliographic latitude (deg)
;						'dis'	= heliocentric distance (AU)
;						'den'	= density (cm^-3)
;						'vel'	= velocity (km/s)
;						'br'	= radial magnetic field (nT)
;						'bt'	= tangential magnetic field (nT)
;						'bn'	= normal magnetic field (nT)
;						'bb'	= field strength (nT)
;						'pdyn'	= dynamic pressure (nPa)
;						'vlng'  = azimuthal flow angle
;						'vlat'  = elevation flow angle
;	t_request		array[2],array[N] type: standard time structure
;						If array[2]: data are extracted between t_request[0] and t_request[1]
;							at times extracted from s/c insitu files
;						If array[N]: data are extracted at times t_request
;						Should be in chronological order
; OPTIONAL INPUT PARAMETERS:
;	source=source	scalar; type: integer or string
;						identifies insitu instrument either by name or integer
;						as provided by function href=Instrument=.
;	delt=delt		scalar; type: time difference structure, or float; default: 0
;						if present, and non-zero, time window used for averaging.
;						in-situ data inside [t-delt/2,t+delt/2] are averaged.
;						If a scalar (float or integer) is specified, then it is
;						assumed to be the time window in days.
;	/weight			if set, standard deviations are used to weight the averages (if available
;						for selected instrument and data)
;						if not set, an unweighted average is calculated
;	traceback=traceback
;					scalar; type: any; default:0
;						trace back distance (AU)
;					if set and nonzero then the in situ data are traced back to a distance of
;					'traceback' AU, i.e. the time and heliographic longitude reflect the time of
;					origin at the traceback distance.
;					This option explicitly assumes that the distances in the in situ files are
;					in AU, the velocities in km/s, and the heliographic longitudes in degrees.
;					If this option is used date_request MUST contain 'lng','dis','vel'
; OUTPUTS:
;	Result			scalar; type: integer
;					0: if something went wrong
;					1: time series succesfully set up
;					(first check this status indicator before using the output arrays!!!)
;	t_given			array[*]; type: time structure
;						times for output time series
;						If delt is set to a positive value and t_request contains more than two
;						elements then t_given=t_request and the time series contains data
;						averaged over delt days.
;						If delt is zero, or t_request contains only 2 elements then
;						t_given contains all times for which spacecraft data were found
;						(averaged over delt days if delt > zero).
;						(if /traceback is used then the time range is shifted relative to
;						the input range)
;	data_given		array[*,n_elements(data_request)], type: float
;						data_given[*,i] are the data for the quantity indicated by data_request[i]
;						missing data are marked by the value !values.f_nan
; OPTIONAL OUTPUT PARAMETERS:
;	source=source	scalar; type: integer
;						if not set on input, the return value is set to whatever default
;						is set by the function href=Instrument= (currently /somni).
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, BadValue, CvSky, TimeFixYear, Instrument, IsTime
;	TimeOp, TimeGet, TimeLimits, TimeSet, TimeUnit, flt_read
;	big_body, big_eph, say, boost
; RESTRICTIONS:
;	Needs system variables: !sun.spiral and !sun.au
;	The input time array should be in chronological order
; PROCEDURE:
; >	Input files for density and velocity are hourly-averages. The filename is:
;		fil = filepath(root=getenv('SSW_SMEI_DAT'),sub='insitu','sw'+spacecraft+'_'+strcompress(year,/rem)+'.hravg')
;		e.g. file = 'imp8_2000.hravg'
; >	The files are assumed to be chronologically ordered.
; >	The file structure is described internally as column numbers for all quantities.
; > Heliographic longitudes are returned in monotonic decreasing order by adding a multiple of 360
;	degrees when the timeseries crosses from one Carrington rotation to the next.
; MODIFICATION HISTORY:
;	AUG-1998, Paul Hick (UCSD/CASS)
;	JUL-1999, Paul Hick (UCSD/CASS); documentation added
;	SEP-2001, Paul Hick (UCSD/CASS)
;		Added GSE & GSM to RTN conversion; currently this is only used
;		for the real-time ACE data.
;	JUL-2002, Paul Hick (UCSD/CASS)
;		Fixed a bug with the delt keyword. When delt was input as a
;		difference time structure it was converted to hours, instead of days.
;	DEC-2003, Paul Hick (UCSD/CASS)
;		Changed /omni to /somni (accesses the files with B in RTN)
;		Added /eomni (accesses the files with B in GSM)
;		The insitu data files can now also be located in a subdirectory
;		of $SSW_SMEI_DATA/insitu, and the filenames can be more general
;		than the default template described in PROCEDURE.
;	SEP-2006, Paul Hick (UCSD/CASS)
;		Added Dana Crider Mars Orbiter data
;	OCT-2006, Paul Hick (UCSD/CASS)
;		All ephemeris calculations now done by big_eph.
;		Added 'vlng','vlat' to pull the off-radial flow angles of the solar
;		wind speed out of the OMNI databases. Note that ONMI2 (EOMNI) and
;		modified OMNI (SOMNI) use different coordinates systems and different
;		definitions for the flow angles.
;	FEB-2007, Paul Hick (UCSD/CASS)
;	    Replaced "from" keyword by "source" keyword. "from" is still accepted
;		but is now considered obsolete.
;	AUG-2007, Paul Hick (UCSD/CASS)
;		Input keyword can now also be a string (instead of just an integer)
;		In addition made the order of the various spacecraft independent
;		of the order in which function Instruments sets them up.
;		Removed keyword 'from'.
;	SEP-2007, Paul Hick (UCSD/CASS)
;		Added instrument vhm to plot Ulysses magnetic field data.
;		Added instrument windb to plot Wind magnetic field data
;		Rename instrument wind to windsw (plots n,v data)
;	FEB-2012, Paul Hick (UCSD/CASS)
;		Bugfix. For magnetic field data requiring a conversion from GSM
;		or GSE, the conversion would be silently skipped if the data request
;		did not include all three magnetic field components, br, bt, and bn.
;		There now is a check which results in an empty data array (and
;		return value 0) if only one or two magnetic field component are
;		specified.
;	JAN-2013, Paul Hick (UCSD/CASS)
;		Added instrument 'rcelias' (real-time SOHO/Celias data)
;		Bugfix. Test for presence of doy instead of month,day of month
;		did not always work.
;	JUN-2013, Paul Hick (UCSD/CASS)
;		Bumped 'atleast' value for rcelias from 5 to 6 after
;		minor change in header lines of insitu data file.
;	MAY-2015, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		For magnetic fields the coordinate transformation from GSM/GSE
;		to RTN is now done prior to the averaging over time delt,
;		instead of after.
;-

InitVar, delt	  , 0
InitVar, traceback, 0
InitVar, weight   , /key
InitVar, silent	  , 0
InitVar, nocvsky  , /key
InitVar, cvlast   , /key

InitVar, source_, Instrument(), set=source

; Convert string input to integer (if the string is not
; set up in Instrument then it is ignored) and the
; default defined in Instrument is returned.

IF IsType(source,/string) THEN source = Instrument(source)

; Spacecraft selection.
; 'nearearth' identifies s/c in earth orbit or near Earth (e.g. L1)

spacecraft = Instrument(source, /name)
nearearth  = Instrument(source, /near)

;====== ADDING SPACECRAFT REQUIRES MODIFICATION HERE =======

; List of spacecraft for which in situ data are made available here.
; This list should be a subset of the list of instruments
; set up in function Instrument (i.e. in instrument_list)

names = [	'imp8'   ,'celias' ,'rcelias','windsw','windb','somni'	, $
			'helios1','helios2',							  $
			'swoops' ,'acesw'  ,'aceb'  ,'eomni','mgs','vhm', $
			'stereoa','stereob',							  $
			'swace']

; Convert 'source' to an index into 'names' (this is the order
; in which the spacecraft are stored in arrays below).

tmp = (where(spacecraft EQ names))[0]
IF tmp EQ -1 THEN BEGIN
	say, silent=silent, threshold=1, 'no in situ data available for '+source
	RETURN, 0
ENDIF
source = tmp

; All insitu data are located in $SSW_SMEI_DAT/insitu
; ($SSW_SMEI_DAT is set to $DAT in the IDL startup script)

insitu_dir = filepath(root=getenv('SSW_SMEI_DAT'),'insitu')

; Data files can be located in 'insitu_dir' itself, or in a subdirectory of 'insitu_dir'.
; The array 'dirs' and 'name' control several different options:

; OPTION 1:  dir='' and name=''
;	data are organized in yearly files that are all located in 'insitu_dir'
;	The filenames match the template <spacecraft>_<yr>.hravg
;	with <spacecraft> the name of the spacecraft or instrument
;	(as extracted from the 'names' array above)
;	and <yr> the year of the data in the file.
;
; OPTION 2: dir='' and name=<filename>
;	all data are located in a single file <filename> located in 'insitu_dir'
;
; The following two options were introduced to handle 'real-time' data
; that typically are kept up to date by periodic cronjobs.

; OPTION 3: dir=<subdir> and name=<prefixYEARpostfix>
;	files are organized in yearly files that are all located in 'insitu_dir/<subdir>'
;	The year is coded into the filename at the location of the string 'YEAR'.
;
; OPTION 4: dir=<subdir> and name=<filename>
;	all data are located in a single file <filename> located in 'insitu_dir/<subdir>'
;	Note that <filename> should NOT contain the string 'YEAR'
;	(this is option 3; use option 2 in this case).

dir 	= [ '','','realcelias'           ,'wind_swe'                 ,'','omni'          ,'','','swoops'         ,'','','omni'          ,'mars'            ,'vhm'        ,'stereo'         ,'stereo'         ,'']; Subdirectory
name	= [ '','','celias_realtime.hravg','YEAR_WIND_hourly_averages','','omni_mYEAR.dat','','','hourav_YEAR.asc','','','omni2_YEAR.dat','swp_mgs_YEAR.dat','rtnYEAR.dat','stereoaYEAR.dat','stereobYEAR.dat','']; File name template

; In situ data files: columns where times, velocities and densities are stored

nyr		= [  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]	; Column # for year
ndoy	= [  1, 1, 2, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1, 1, 1, 1, 1]	; Column # for doy of year
nmon	= [ -1,-1,-1,-1, 1,-1,-1,-1,-1, 1, 1,-1,-1,-1,-1,-1,-1]	; Column # for month
nday	= [ -1,-1,-1,-1, 2,-1,-1,-1,-1, 2, 2,-1,-1,-1,-1,-1,-1]	; Column # for day of month

nhr		= [  2, 2, 3, 2, 3, 2, 2, 2, 2, 3, 3, 2,-1,-1, 2, 2, 2]	; Column # for hour
nmin	= [ -1, 3, 4,-1, 4,-1,-1,-1, 3,-1,-1,-1,-1,-1, 3, 3, 3]	; Column # for minutes

lng		= [ -1,-1,-1,-1,-1, 4,-1,-1, 7,-1,-1,-1,-1,-1,-1,-1,-1]	; Column # for heliolongitude
lat		= [ -1,-1,-1,-1,-1, 3,-1,-1, 6,-1,-1,-1,-1,-1,-1,-1,-1]	; Column # for heliolatitude
dis		= [ -1,-1,-1,-1,-1,-1,-1,-1, 5,-1,-1,-1,-1,-1,-1,-1,-1]	; Column # for heliodistance

; 'Bad data' flags for quantities used in input files

lngflag = [0.0,0.0,0.0,0.0,0.0,9999.9,0.0,0.0,999.99,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
latflag = [0.0,0.0,0.0,0.0,0.0,9999.9,0.0,0.0,999.99,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
disflag = [0.0,0.0,0.0,0.0,0.0,   0.0,0.0,0.0,9.9999,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]

; Make sure avoid setting both gse2rtn and gsm2rtn for the same source

gse2rtn = [  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]	; Trigger for B-conversion from GSE --> RTN
gsm2rtn = [  0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0]	; Trigger for B-conversion from GSM --> RTN

; This is passed to flt_read, and is used to skip past header records
; that contain less than the indicated number).

atleast = [  1, 1, 6, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]	; Minimum nr of values on each record


;====== ADDING OTHER QUANTITIES REQUIRES MODIFICATIONS HERE =======

posin = [ [ 19, 6, 7, 6,-1,12,11,11, 8, 7,-1,23,-1,-1, 5, 5, 5], $	; Column # for in situ densities
		  [  4, 5, 6, 3,-1, 9,10,10,12, 8,-1,24,-1,-1, 6, 6, 8], $	; Column # for in situ velocities
		  [ -1,-1,-1,-1, 5, 5, 6, 6,-1,-1, 7,12,-1, 3,-1,-1,-1], $	; Column # for in situ RTN Br, or GSE/GSM Bx
		  [ -1,-1,-1,-1, 8, 6, 7, 7,-1,-1, 8,15,-1, 4,-1,-1,-1], $	; Column # for in situ RTN Bt, or GSE/GSM By
		  [ -1,-1,-1,-1, 9, 7, 8, 8,-1,-1, 9,16,-1, 5,-1,-1,-1], $	; Column # for in situ RTN Bn, or GSE/GSM Bz
		  [ -1,-1,-1,-1,10, 8, 9, 9,-1,-1,10, 9,-1, 6,-1,-1,-1], $	; Column # for in situ magnetic field magnitude
		  [ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,38,-1,-1, 7, 7,-1], $	; Column # for Kp
		  [ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,40,-1,-1,-1,-1,-1], $	; Column # for DST
		  [ -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 2,-1,-1,-1,-1], $	; Column # for dynamic pressure
		  [ -1,-1,-1,-1,-1,11,-1,-1,-1,-1,-1,25, 2,-1,-1,-1,-1], $	; Column # for vel azimuth
		  [ -1, 8,-1,-1,-1,10,-1,-1,-1,-1,-1,26, 2,-1,-1,-1,-1]] 	; Column # for vel elevation


flag = [[ 0.0,-1.0,-1.0,9999.0,     0.0, 999.9,  999.90,  999.90, 0.0, -9999.9,   0.0, 999.9, 0.0, 0.000, 0.0, 0.0,-9999.9], $	; Flag for in situ densities
		[ 0.0,-1.0,-1.0,9999.0,     0.0,9999.0, 9999.90, 9999.90, 0.0, -9999.9,   0.0,9999.0, 0.0, 0.000, 0.0, 0.0,-9999.9], $	; Flag for in situ velocities
		[ 0.0, 0.0, 0.0,   0.0,-1.0E+31, 999.9,99999.99,99999.99, 0.0,     0.0,-999.9, 999.9, 0.0,99.999, 0.0, 0.0,    0.0], $	; Flag for in situ radial magnetic field
		[ 0.0, 0.0, 0.0,   0.0,-1.0E+31, 999.9,99999.99,99999.99, 0.0,     0.0,-999.9, 999.9, 0.0,99.999, 0.0, 0.0,    0.0], $	; Flag for in situ tangential magnetic field
		[ 0.0, 0.0, 0.0,   0.0,-1.0E+31, 999.9,99999.99,99999.99, 0.0,     0.0,-999.9, 999.9, 0.0,99.999, 0.0, 0.0,    0.0], $	; Flag for in situ normal field
		[ 0.0, 0.0, 0.0,   0.0,-1.0E+31, 999.9,99999.99,99999.99, 0.0,     0.0,-999.9, 999.9, 0.0,99.999, 0.0, 0.0,    0.0], $	; Flag for in situ magnetic field strength
		[ 0.0, 0.0, 0.0,   0.0,     0.0,   0.0,    0.00,    0.00, 0.0,     0.0,   0.0,    99, 0.0, 0.000, 0.0, 0.0,    0.0], $	; Flag for Kp
		[ 0.0, 0.0, 0.0,   0.0,     0.0,   0.0,    0.00,    0.00, 0.0,     0.0,   0.0, 99999, 0.0, 0.000, 0.0, 0.0,    0.0], $	; Flag for DST
		[ 0.0, 0.0, 0.0,   0.0,     0.0,   0.0,    0.00,    0.00, 0.0,     0.0,   0.0,   0.0, 0.0, 0.000, 0.0, 0.0,    0.0], $	; Flag for dynamic pressure
		[ 0.0, 0.0, 0.0,   0.0,     0.0, 999.9,    0.00,    0.00, 0.0,     0.0,   0.0, 999.9, 0.0, 0.000, 0.0, 0.0,    0.0], $	; Flag for velocity azimuth
		[ 0.0,99.0, 0.0,   0.0,     0.0, 999.9,    0.00,    0.00, 0.0,     0.0,   0.0, 999.9, 0.0, 0.000, 0.0, 0.0,    0.0]	]	; Flag for velocity elevation

;======= END MODIFICATIONS

n_data_request = n_elements(data_request)
extended_data_request = data_request

; If magnetic field data are requested for insitu data requiring conversion
; from GSE/GSM to RTN (e.g. aceb), then all three components need to be present
; so make sure that extended_data_request contains all three by tagging missing components
; at the end.

rtnout	= [	(where(extended_data_request EQ 'br' ))[0],	$
			(where(extended_data_request EQ 'bt' ))[0],	$
			(where(extended_data_request EQ 'bn' ))[0]	]

nB = round(total(rtnout NE -1))			; Nr of requested B components
gse2rtn = nB GT 0 AND gse2rtn[source]	; nB = 0: no magnetic field requested
gsm2rtn = nB GT 0 AND gsm2rtn[source]

atleast = atleast[source]

IF (gse2rtn OR gsm2rtn) AND nB NE 3 THEN BEGIN

	say, silent=silent, threshold=1, 'request only includes '+strcompress(nB,/rem)+' B component(s):  '+strjoin(extended_data_request[where(rtnout NE -1)],',')

	; Conversion for magnetic field required, but not all components
	; are present. Add the missing components to extended_data_request

	IF rtnout[0] EQ -1 THEN extended_data_request = [extended_data_request,'br']
	IF rtnout[1] EQ -1 THEN extended_data_request = [extended_data_request,'bt']
	IF rtnout[2] EQ -1 THEN extended_data_request = [extended_data_request,'bn']

	rtnout	= [	(where(extended_data_request EQ 'br' ))[0],	$
				(where(extended_data_request EQ 'bt' ))[0],	$
				(where(extended_data_request EQ 'bn' ))[0]	]

	nB = round(total(rtnout NE -1))			; Nr of requested B components (now is 3)

ENDIF

lngout	=	(where(extended_data_request EQ 'lng'))[0]
disout	=	(where(extended_data_request EQ 'dis'))[0]
velout	=	(where(extended_data_request EQ 'vel'))[0]

IF traceback NE 0 THEN BEGIN
	IF lngout EQ -1 THEN BEGIN
		extended_data_request = [extended_data_request,'lng']
		lngout	=	(where(extended_data_request EQ 'lng'))[0]
	ENDIF
	IF disout EQ -1 THEN BEGIN
		extended_data_request = [extended_data_request,'dis']
		disout	=	(where(extended_data_request EQ 'dis'))[0]
	ENDIF
	IF velout EQ -1 THEN BEGIN
		extended_data_request = [extended_data_request,'vel']
		velout	=	(where(extended_data_request EQ 'vel'))[0]
	ENDIF
ENDIF

n_extended_data_request = n_elements(extended_data_request)

latout	=	(where(extended_data_request EQ 'lat' ))[0]

posout	= [	(where(extended_data_request EQ 'den' ))[0],	$
			(where(extended_data_request EQ 'vel' ))[0],	$
			(where(extended_data_request EQ 'br'  ))[0],	$
			(where(extended_data_request EQ 'bt'  ))[0],	$
			(where(extended_data_request EQ 'bn'  ))[0],	$
			(where(extended_data_request EQ 'bb'  ))[0],	$
			(where(extended_data_request EQ 'kp'  ))[0],	$
			(where(extended_data_request EQ 'dst' ))[0],	$
			(where(extended_data_request EQ 'pdyn'))[0],	$
			(where(extended_data_request EQ 'vlng'))[0],	$
			(where(extended_data_request EQ 'vlat'))[0]		]

dir		= dir [source]
name	= name[source]

nyr		= nyr [source]
ndoy	= ndoy[source]
nmon	= nmon[source]
nday	= nday[source]
nhr		= nhr [source]
nmin	= nmin[source]

lng		= lng [source]
lat		= lat [source]
dis		= dis [source]

lngflag	= lngflag[source]
latflag	= latflag[source]
disflag	= disflag[source]

posin	= reform(posin[source,*])
flag	= reform(flag [source,*])

rtnin   = posin[[2,3,4]]
rtnflag = flag [[2,3,4]]

trange	= TimeLimits(t_request,/bounds)

djd = (IsTime(delt) ? TimeGet(delt,/full,/day) : delt)/2.0d0	; Time window is t+[-djd,+djd] days

tmp = strjoin(TimeGet(trange, /ymd, upto=TimeUnit(/day)),'-')
IF djd NE 0 THEN tmp += ', '+(['unweighted','weighted'])[weight]+' '+ $
	strcompress(2*djd,/rem)+'-day mean'+(['','using standard deviations'])[weight]
say, silent=silent, threshold=1, tmp 

yy		= TimeGet(trange, /year)
trange	= TimeGet(trange, /njd )

datafound = 0
destroyvar, errors
threshold = 0

FOR year=yy[0],yy[1] DO BEGIN					; Loop over all years

	CASE dir EQ '' OF
	0: BEGIN
		fil = strpos(name,'YEAR')
		CASE fil NE -1 OF
		0: fil = filepath(root=insitu_dir,sub=dir,name)
		1: fil = filepath(root=insitu_dir,sub=dir,strmid(name,0,fil)+strcompress(year,/rem)+strmid(name,fil+4))
		ENDCASE
	END
	1: BEGIN
		CASE name EQ '' OF
		0: fil = filepath(root=insitu_dir,name)
		1: fil = filepath(root=insitu_dir,spacecraft+'_'+strcompress(year,/rem)+'.hravg')
		ENDCASE
	END
	ENDCASE

	status = flt_read(fil, tmp, /exp, error=error, silent=silent+1, atleast=atleast)
	IF status THEN BEGIN
		status = max(posin) LT (size(tmp,/dim))[0]
		IF NOT status THEN BEGIN
			boost, errors, ['only'+strcompress((size(tmp,/dim))[0])+' columns found; need at least'+strcompress(max(posin)+1),hide_env(fil)]
			threshold = 2
		ENDIF
	ENDIF ELSE	$
		boost, errors, error

	IF status THEN BEGIN

		;======= DEAL WITH PECULARITIES IN IN SITU DATA FILES HERE
		; Ace has hours and minutes merged as 4-char string: 1200 is 12 h, 00 m

		CASE spacecraft OF
		'acesw'  : tmp[nhr,*] = round(tmp[nhr,*]/100.0)+(tmp[nhr,*] mod 100)/60.0
		'aceb'   : tmp[nhr,*] = round(tmp[nhr,*]/100.0)+(tmp[nhr,*] mod 100)/60.0
		'rcelias': tmp[nhr,*] = tmp[nhr,*]+0.5
		ELSE:
		ENDCASE
		;=======

		tmp[nyr,*] = TimeFixYear(tmp[nyr,*])

		CASE ndoy NE -1 OF
		0: t_tmp = TimeSet(yr=tmp[nyr,*],mon=tmp[nmon,*],day=tmp[nday,*],hour=(nhr NE -1)*tmp[nhr > 0,*],minute=(nmin NE -1)*tmp[nmin > 0,*])
		1: t_tmp = TimeSet(yr=tmp[nyr,*],doy=tmp[ndoy,*]				,hour=(nhr NE -1)*tmp[nhr > 0,*],minute=(nmin NE -1)*tmp[nmin > 0,*])
		ENDCASE
		
		t_tmp = reform(t_tmp,/overwrite)

		;IF nmin NE -1 THEN t_tmp = TimeOp(/add,t_tmp,TimeSet(/diff,minute=tmp[nmin,*]))

		; The IMP data sometimes contain more than one entry for the same time

		jdtmp = TimeGet(t_tmp, /njd)

		list = uniq(jdtmp, sort(jdtmp))
		IF n_elements(list) NE n_elements(jdtmp) THEN BEGIN
			say, silent=silent, threshold=2, $
				'# duplicate times:'+strcompress(n_elements(jdtmp)-n_elements(list))+'; first occurrence used only'
			t_tmp = t_tmp[list]
			jdtmp = jdtmp[list]
			tmp	  = tmp[*,list]
		ENDIF

		;list = where(trange[0]-1.1*djd le jdtmp and jdtmp le trange[1]+1.1*djd, nlist)
		list = where(trange[0]-djd LE jdtmp AND jdtmp LE trange[1]+djd, nlist)

		IF nlist NE 0 THEN BEGIN				; Extract all data between start & end time
			CASE datafound OF
			0: BEGIN
				t_given = t_tmp[list]
				jd   = jdtmp[list]				; Accumulate Julian dates
				data = tmp[*,list]				; Accumulate in situ data
				datafound = 1
			END
			1: BEGIN
				t_given = [t_given, t_tmp[list]]
				jd   = [jd  ,jdtmp[list]]
				data = [[data],[tmp[*,list]]]
			END
			ENDCASE
		ENDIF

	ENDIF

ENDFOR

IF datafound THEN BEGIN

	sctimes = n_elements(t_request) EQ 2 OR djd EQ 0

	; Array t_given contains the full time structure corresponding to
	; jd (times from insitu files)

	CASE sctimes OF
	0: BEGIN							; Use requested input times
		t_given = t_request
		jdt = TimeGet(t_given, /njd)
	END
	1: jdt = jd							; Use s/c times read from file
	ENDCASE

	; Now t_given corresponds to the requested times jdt
	; jd contains the times for the insitu data

	list = where(trange[0] LE jdt AND jdt LE trange[1], nlist)
	IF nlist EQ 0 THEN datafound = 0

ENDIF

IF NOT datafound THEN BEGIN
	boost, errors, spacecraft+'_*.hravg, no data found in requested time range'
	say, silent=silent, threshold=threshold, errors
	RETURN, 0
ENDIF

; If magnetic fields are present in GSM or GSE, then these need to coverted
; to RTN prior to smoothing

nB = round(total(rtnin NE -1))			; Nr of B components in input (should be 0 or 3)

IF NOT cvlast AND (gse2rtn OR gsm2rtn) AND nB EQ 3 THEN BEGIN	; Magnetic field requested, and present in input

	b_found = data[rtnin,*]
	good = where(b_found[0,*] NE rtnflag[0] AND b_found[1,*] NE rtnflag[1] AND b_found[2,*] NE rtnflag[2])

	IF good[0] EQ -1 THEN BEGIN
		say, silent=silent, threshold=2, 'no valid magnetic data'
	ENDIF ELSE IF nocvsky THEN BEGIN
		say, silent=silent, threshold=1, 'skipping magnetic field conversion from '+(['GSM','GSE'])[gse2rtn]+' to RTN'
	ENDIF ELSE BEGIN
		b_found[*,good] = CvSky(TimeSet(njd=jd[good]), b_found[*,good], from_GSE=gse2rtn, from_GSM=gsm2rtn, /to_RTN, /rectangular, /silent)
		say, silent=silent, threshold=1, 'magnetic field converted from '+(['GSM','GSE'])[gse2rtn]+' to RTN before averaging'
		data[rtnin,*] = b_found
	ENDELSE

ENDIF

njd = n_elements(jdt)
data_given = replicate(BadValue(0.0),njd,n_extended_data_request)

IF sctimes THEN BEGIN
	IF lng NE -1 THEN BEGIN			; Heliolongitudes present
		list = where(data[lng,*] NE lngflag, nlist)
		IF nlist NE 0 AND lngout NE -1 THEN data_given[list,lngout] = data[lng,list]
	ENDIF

	IF lat NE -1 THEN BEGIN			; Heliolatitudes present
		list = where(data[lat,*] NE latflag, nlist)
		IF nlist NE 0 AND latout NE -1 THEN data_given[list,latout] = data[lat,list]
	ENDIF

	IF dis NE -1 THEN BEGIN			; Heliocentric distances present
		list = where(data[dis,*] NE disflag, nlist)
		IF nlist NE 0 AND disout NE -1 THEN data_given[list,disout] = data[dis,list]
	ENDIF
ENDIF

IF NOT sctimes OR lng EQ -1 OR lat EQ -1 OR  dis EQ -1 THEN BEGIN

	IF nearearth THEN spacecraft = jpl_body(/earth,/string)

	body = big_body()
	j = (where(strlowcase(big_body(spacecraft)) EQ strlowcase(body)))[0]

	CASE j NE -1 OF
	0: say, silent=silent, threshold=1, 'no  heliographic coordinates for '+spacecraft
	1: BEGIN
	    body = body[j]
		tmp = big_eph(t_given	, $
				body = body	, $
				/to_ecliptic, $
				/to_sphere	, $
				/precess	, $
				/degrees	)
		tmp = reform(tmp,/overwrite)
		say, silent=silent, threshold=1, 'calculating heliographic coordinates for '+body
		tmp = CvSky(t_given, from_ecliptic=tmp, /to_heliographic, /degrees,/silent)
		IF lngout NE -1 THEN data_given[*,lngout] = tmp[0,*]
		IF latout NE -1 THEN data_given[*,latout] = tmp[1,*]
		IF disout NE -1 THEN data_given[*,disout] = tmp[2,*]
		tmp = 0
	END
	ENDCASE

ENDIF

IF lngout NE -1 THEN BEGIN
	list = data_given[*,lngout]					; Make sure longitudes decrease monotonically
	list = list-shift(list,1)					; .. by adding multiple of 360 deg where needed
	list[0] = 1
	list = where(list GT 0)
	nlist = n_elements(list)
	FOR i=nlist-1,1,-1 DO data_given[list[i-1]:list[i]-1,lngout] = data_given[list[i-1]:list[i]-1,lngout]+(nlist-i)*360
ENDIF

IF djd EQ 0 THEN BEGIN							; No smoothing needed

	FOR j=0,n_elements(posin)-1 DO BEGIN		; Loop over all selected quantities
		IF posin[j] NE -1 THEN BEGIN			; Data present
			list = where(data[posin[j],*] NE flag[j], nlist)
			IF nlist NE 0 AND posout[j] NE -1 THEN data_given[list,posout[j]] = data[posin[j],list]
		ENDIF
	ENDFOR

ENDIF ELSE BEGIN								; Smooth selected quantities

	FOR i=0,njd-1 DO BEGIN						; Loop over all requested times
		tlist = abs(jd-jdt[i]) LE djd

		FOR j=0,n_elements(posin)-1 DO BEGIN	; Loop over den and vel

			IF posin[j] NE -1 THEN BEGIN
				list = where(tlist AND data[posin[j],*] NE flag[j], nlist)
				IF nlist NE 0 THEN BEGIN
					flist = data[posin[j],list]

					CASE weight OF
					0: wlist = replicate(1.,nlist); Unweighted average
					1: BEGIN					; Weighted average
						wlist = data[posin[j]+1,list]; Standard deviations

						dev = where(wlist EQ 0, ndev)

						CASE nlist OF
						ndev: wlist = replicate(1.,nlist) 	; Unweighted average
						ELSE: wlist[dev] = total(wlist)/(nlist-ndev) ; Fill in s.d. = 0.
						ENDCASE

						wlist = 1./(wlist*wlist); Weights 1./(s.d)^2
					END
					ENDCASE

					IF posout[j] NE -1 THEN data_given[i,posout[j]] = total(wlist*flist)/total(wlist)

				ENDIF
			ENDIF

		ENDFOR
	ENDFOR

	list = where(trange[0] LE jdt AND jdt LE trange[1])
	t_given = t_given[list]
	data_given = data_given[list,*]

ENDELSE

IF cvlast AND (gse2rtn OR gsm2rtn) THEN BEGIN

	b_given = transpose(data_given[*,rtnout])
	IF total( finite(b_given) ) EQ 0 THEN BEGIN
		say, silent=silent, threshold=2, 'no valid magnetic data'
	ENDIF ELSE IF nocvsky THEN BEGIN
		say, silent=silent, threshold=1, 'skipping magnetic field conversion from '+(['GSM','GSE'])[gse2rtn]+' to RTN'
	ENDIF ELSE BEGIN
		b_given = CvSky(t_given, b_given, from_GSE=gse2rtn, from_GSM=gsm2rtn, /to_RTN, /rectangular, /silent)
		say, silent=silent, threshold=1, 'magnetic field converted from '+(['GSM','GSE'])[gse2rtn]+' to RTN after averaging'
		data_given[*,rtnout] = transpose(b_given)
	ENDELSE

ENDIF

; Traceback time and longitude to distance specified in variable traceback (in AU)

IF traceback NE 0 THEN BEGIN				; Traceback

	list = where( finite(data_given[*,disout]) AND finite(data_given[*,velout]) )
	IF list[0] EQ -1 THEN BEGIN
		say, silent=silent, threshold=1, 'no data available for traceback'
		RETURN, 0
	ENDIF
	t_given = t_given[list]
	data_given = data_given[list,*]						; Only keep data where traceback is done
	dT  = (data_given[*,disout]-traceback)/data_given[*,velout]
	data_given[*,lngout] = data_given[*,lngout]+!sun.spiral*dT*360
	t_given = TimeOp(/add, t_given, TimeSet(/diff, -dT*(!sun.au/0.00086400d0), TimeUnit(/day)))

ENDIF

IF traceback NE 0 THEN BEGIN

	trange = TimeLimits(t_given, /bounds)		; Exact time range
	say, silent=silent, threshold=2, $
		['original time range:'+strjoin(TimeGet(trange, /ymd, upto=TimeUnit(/day)),'-'), $
		 'original heliolongitude:'+strcompress(min(data_given[*,lngout]))+' to'+strcompress(max(data_given[*,lngout]))]

ENDIF

IF n_extended_data_request GT n_data_request THEN BEGIN
	data_given = data_given[*,0:n_data_request-1]
ENDIF

RETURN, 1  &  END