;+
; NAME:
;	vu_insitu
; PURPOSE:
;	Produces a plot, comparing a time series from a tomographic reconstruction
;	and insitu observations from instrument 'source'.
; CATEGORY:
;	vupack
; CALLING SEQUENCE:
	PRO vu_insitu, hdr, ff, ut0=ut0, ihdr=ihdr,		$

		ts_all 		= ts_all	,	$	; Defines times in time series
		ts_range	= ts_range	,	$
	 	ts_step		= ts_step	,	$

		ti_range	= ti_range	,	$	; Controls s/c insitu data
		nofuture	= nofuture	,	$
		delt		= delt		,	$
		source		= source	,	$

		periodic	= periodic	,	$
		hardcopy	= hardcopy	,	$	; to be phased out
		printer		= printer	,	$
		noplot		= noplot	,	$
		tm_label	= tm_label	,	$
		power		= power		,	$

	 	type		= type		,	$

		nosequence	= nosequence,	$
		use_forecast_time = use_forecast_time,	$
		cast_offset = cast_offset,	$
		silent		= silent	,	$

		body		= body		,	$
		module		= module	,	$

		ascii_file	= ascii_file,	$

 		_extra      = _extra	,	$

		f1=f1, f2=f2, t1=t1, t2=t2
; INPUTS:
;	hdr				array[k]; type: structure
;						information about the tomographic reconstruction
;						(single file or time-dependent sequence)
;						if 'hdr' does not exist, then user is prompted for
;						a data file using the 'pickfile' dialog.
;	ff				array[n,m,l,2,k]; type: float
;						3D volume data associated with 'hdr'
;						if not present then the files
;						specified in the 'hdr' argument are read.
; OPTIONAL INPUT PARAMETERS:
;	ut0=ut0			scalar or array[1]; type: time structure, float or integer
;						Specifies the time on which the time series is centered.
;						If not specified then the time from 'hdr[ihdr]' is used.
;						See href=Carrington= for more details.
;	ihdr=ihdr		scalar; type: integer; default: 0
;						index of header to be used for display
;
;	/nosequence		Only used if n_elements(hdr) > 1
;					if n_elements(hdr) > 1 then an attempt is made to determine whether
;					the input data can be treated as a time-ordered sequence, e.g. output
;					from the time-dependent tomography (see href=vu_is_sequence=).
;					Override this check by setting /nosequence.
;
;	/forecast		if set, a dashed line at the plot center is added, and the word
;						'FORECAST' is plotted to the right of the dashed line.
;	/use_forecast_time
;					(only for internal use by vu_rawtimeseries)
;
;	power=power 	scalar, array[2] or array[2,k]; type: float
;						normalization applied to the data if 'ff' by vu_timeseries.
;						If NOT set then then the values are extracted from 'hdr'.
;						Override these values by setting this keyword. If setting these
;						to an array be careful to be consistent with the array structure
;						of 'hdr' and 'ff'.
;
;	body=body		scalar; type: string; default: none
;						selects a body for which we have an ephemeris available
;						the timeseries for the selected body is plotted.
;						E.g. body=jpl_body(/mars,/string) plots a time series for Mars.
;
;	Keywords to define the time series (i.e. the times at which data are evaluated).
;	Either ts_all is used, or the pair ts_range and ts_step. If neither of these is set then
;	the times in 'hdr' are used).
;
;	ts_all=ts_all	array; type: float or time structures
;						times at which data are to be evaluated, as Carrington times
;						or UT times.
;	ts_step=ts_step	scalar; type: float or difference time structure: default: 0.25
;						time step in time series in days.
;	ts_range=ts_range
;					array[2]; type: float; default: [-15,15]
;						range of times centered on ut0 in days. The time series covers
;						[ut0-ts_range[0],ut0+ts_range[1]] in steps of ts_step.
;
;	Keywords to limit the range of in situ data plotted.
;
;	ti_range=ti_range
;					array[1] or array[2]; type: time_structure or Carrington variable
;						if array[1] only in situ data upto ti_range[0] are plotted
;						if array[2] only in situ data inside time range ti_range are plotted.
;	/nofuture		if set, ti_range is set to ut0, i.e. only in situ data up to time ut0
;						(usually the forecast time) are plotted.
;
;	Keywords to select in situ data (see href=InsituTimeSeries=):
;
;	source=source	scalar; type: integer
;						identifies insitu instrument (see function href=Instrument=
;	delt=delt		scalar; type: time difference structure, or float; default: 0.0
;						if present, and non-zero, time window used for averaging (days)
;						the in situ date from instrument 'source'
;
;	For keywords controlling plot layout see href=vu_type_insitu=.
;
;	ascii_file=ascii_file
;					scalar; type: string
;						name of file used to write timeseries data in ascii
;						If ascii_file='terminal" output is written to terminal.
; OUTPUTS:
;	(graphics)
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
	@vu_fnc_include.pro 		;
; CALLS:
;	InitVar, Carrington, vu_getdata, vu_get, vu_type_insitu, vu_get_page, vu_timeseries
;	vu_prefix, InsituTimeSeries, vu_InsituCurve, TimeUnit, TimeSet, TimeOp, TimeGet
;	gridgen, BadValue, Instrument, reset_colors, set_page, vu_is_sequence, destroyvar
;	SetFileSpec, IsType, IsTime, Carrington, jpl_body, hide_env
; PROCEDURE:
; >	The input file(s) come from a tomographic reconstruction, and contain 3D
;	heliospheric matrices for density and velocity (nv3d* file) or magnetic
;	field (wson* file). The file(s) are read using vu_getdata.
; >	If a time-dependent sequence is processed the times from the headers are
;	used for the time series at Earth. For a single matrix ut0, ts_range and
;	ts_step determine the time series.
; MODIFICATION HISTORY:
;	MAR-2000, Paul Hick (UCSD/CASS)
;	SEP-2002, Paul Hick (UCSD/CASS)
;		Added ts_all and /user_forecast_time keywords.
;	OCT-2002, Paul Hick (UCSD/CASS)
;		Dropped /earth keyword and associated section of code that
;		was used to compare time series from nv3d* and ea* files.
;	FEB-2007, Paul Hick (UCSD/.CASS; pphick@ucsd.edu)
;		Added option to print to terminal
;	AUG-2007, George Megally(UCSD/CASS; gmegally@ucsd.edu
;		Expanded program to include STEREO A/B satellites.
;	AUG-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Renamed keyword 'from' to 'source' to avoid conflicts
;		with 'from' keyword 'from' Time* functions
;-

InitVar, silent  , 0

InitVar, nofuture, /key
InitVar, printer , /key
InitVar, noplot	 , /key
InitVar, use_forecast_time, /key

IF IsType(hardcopy,/defined) THEN printer = hardcopy ; to be phased out

utt = vu_getdata(hdr, ff, ut0=ut0, ihdr=ihdr, silent=silent, _extra=_extra, count=count)

IF count EQ 0 THEN RETURN

utt  = Carrington(utt,/get_time)
uday = TimeUnit(/day)

;message, /info, TimeGet(utt,/ymd,upto=TimeUnit(/minute))+strcompress(ihdr)

CASE use_forecast_time OF
0: BEGIN
	InitVar, nosequence, /key
	nosequence = nosequence OR n_elements(hdr) EQ 1
	IF NOT nosequence THEN nosequence = 1-vu_is_sequence(hdr,silent=silent)
END
1: nosequence = 0+2*IsType(ts_all, /defined)
ENDCASE

IF use_forecast_time THEN BEGIN

	; vu_rawtimeseries uses this.
	; If ts_all is not defined then a fore/aftcast timeseries for
	; time-dependent tomography is required (nosequence = 0)
	; If ts_all is set then we are dealing with a fore/aftcast timeseries for
	; corotating tomography or a fixed-time time series for time-dependent
	; tomography (nosequence=2).

 	InitVar, ts_all, vu_get(hdr, /uttime)

ENDIF ELSE IF IsType(ts_range, /undef) AND IsType(ts_step, /undef) THEN BEGIN

	; It neither ts_range nor ts_step are set, use the times from the headers
	; (unles ts_all is input)

	IF nosequence EQ 0 THEN InitVar, ts_all, vu_get(hdr, /uttime)

ENDIF


nts_all = n_elements(ts_all)

IF nts_all EQ 0 THEN BEGIN

	InitVar, ts_range, [-15.,15.]			; Range +/- 15 days
	InitVar, ts_step , 0.25

	IF IsTime(ts_range) THEN ts_range = TimeSet(day=ts_range)
	IF IsTime(ts_step ) THEN ts_step  = TimeSet(day=ts_step )

	nts_all = round( (ts_range[1]-ts_range[0])/float(ts_step)+1 )

	; The next roundoff is not really necessary, but makes the timeseries in
	; movies look less jumpy, because identical times are used.

	;ts_all = TimeGet(utt, /roundt, day=ts_step)
	ts_all = utt
	ts_all = TimeOp(/add, ts_all, TimeSet(day=gridgen(nts_all, range=ts_range)) )
											; utt +/- 15 days with time every 6 hours
ENDIF

nts_all = n_elements(ts_all)


InitVar, power, vu_get(hdr, /power_norm)		; Needed by vu_timeseries (default is zero)

CASE IsType(body, /defined) OF

0: BEGIN

	InitVar, source, Instrument()				; Set default if not defined yet

	IF NOT Instrument(source,/nearearth) THEN BEGIN
		message, /info, '"'+source+'" is not a near-Earth instrument'
		RETURN
	ENDIF

	sc_label = Instrument(source, /label)
	nearL1   = Instrument(source, /L1)			; Needed by vu_timeseries

	IF IsType(tm_label, /defined) THEN tm_label_ = tm_label

	type = vu_type_insitu('earth_insitu', _extra=_extra, hdr=hdr, type=type)

	no_insitu = 0

END

1: BEGIN

	_mercury = jpl_body(/mercury,/string)
	_venus   = jpl_body(/venus  ,/string)
	_earth   = jpl_body(/earth  ,/string)
	_mars    = jpl_body(/mars   ,/string)
	_ulysses = big_body('Ulysses')
	_stereoa = big_body('Stereo A')
	_stereob = big_body('Stereo B')

	_body    = strcompress(strlowcase(body), /rem)

	CASE 1 OF
	strpos(_body,strlowcase(_mercury)) EQ 0: BEGIN
		type = vu_type_insitu('mercury_insitu' , _extra=_extra, hdr=hdr, type=type)
		body = _mercury
		no_insitu = 1
	END	
	strpos(_body,strlowcase(_venus  )) EQ 0: BEGIN
		type = vu_type_insitu('venus_insitu', _extra=_extra, hdr=hdr, type=type)
		body = _venus
		no_insitu = 1
	END	
	strpos(_body,strlowcase(_earth  )) EQ 0: BEGIN
		InitVar, source, Instrument()
		sc_label = Instrument(source, /label)
		nearL1   = Instrument(source, /L1)			; Needed by vu_timeseries
		type = vu_type_insitu('earth_insitu', _extra=_extra, hdr=hdr, type=type)
		body = _earth
		no_insitu = 0
	END
	strpos(_body,strlowcase(_mars   )) EQ 0: BEGIN
		InitVar, source, Instrument(/mgs)
		sc_label = Instrument(source, /label)
		type = vu_type_insitu('mars_insitu' , _extra=_extra, hdr=hdr, type=type)
		body = _mars
		no_insitu = Instrument(source) NE Instrument(/mgs)
	END
	strpos(_body,strlowcase(_ulysses)) EQ 0: BEGIN
		InitVar, source, Instrument(/swoops)
		sc_label = Instrument(source, /label)
		type = vu_type_insitu('ulysses_insitu' , _extra=_extra, hdr=hdr, type=type)
		body = _ulysses
		no_insitu = Instrument(source) NE Instrument(/swoops)
	END
;8/29 insert begin 
	strpos(_body,strlowcase(strcompress(_stereoa,/rem))) EQ 0: BEGIN
		InitVar, source, Instrument(/stereoa)
		sc_label = Instrument(source, /label)
		type = vu_type_insitu('stereoa_insitu' , _extra=_extra, hdr=hdr, type=type)
		body = _stereoa
		;Currently no insitu data is available, change later.
		no_insitu = Instrument(source) NE Instrument(/stereoa)
	END
	strpos(_body,strlowcase(strcompress(_stereob,/rem))) EQ 0: BEGIN
		InitVar, source, Instrument(/stereob)
		sc_label = Instrument(source, /label)
		type = vu_type_insitu('stereob_insitu' , _extra=_extra, hdr=hdr, type=type)
		body = _stereob
		no_insitu = Instrument(source) NE Instrument(/stereob)
	END
;end 8/29

	ELSE: BEGIN
		type = vu_type_insitu('earth_insitu', _extra=_extra, hdr=hdr, type=type)
		no_insitu = 1
	END
	ENDCASE

	CASE IsType(tm_label, /defined) OF
	0: tm_label_ = '3D @ '+body
	1: tm_label_ = tm_label +' @ '+body
	ENDCASE

END

ENDCASE

CASE nosequence OF

; Time-dependent sequence of tomography data processed as a whole.
; If /use_forecast_time and /fixed_time is NOT set the timeseries
; is calculated here also.

0: f1 = vu_timeseries( vu_get(hdr, /array_range), vu_get(hdr, /start_distance),	$
		vu_get(hdr, /distance_step), ff, T=ts_all,					$
		periodic=periodic, nearL1=nearL1, power=power, body=body, 	$
		sequence=vu_get(hdr, /uttime), _extra=_extra)


; Only one matrix (ihdr) of a collection of tomography data is processed. The collection
; could be a group of corotating tomography data, or a time-dependent sequence.

1: f1 = vu_timeseries( vu_get(hdr[ihdr], /array_range), vu_get(hdr[ihdr], /start_distance),	$
		vu_get(hdr[ihdr], /distance_step), ff[*,*,*,*,ihdr], T=ts_all,	$
		periodic=periodic, nearL1=nearL1, power=power, body=body, _extra=_extra)


; Used by vu_rawtimeseries (if /use_forecast_time and 'ts_all' are set).
; hdr is a collection of raw tomography data. For a time-dependent data set 'ts_all' is
; an array containing only one value. For a corotating tomography 'ts_all' is the
; forecast time with a fore/aftcast offset added/subtracted.

2: BEGIN
	f1 = make_array(type=IsType(ff), dim=[nts_all,(size(ff,/dim))[3]])
	FOR i=0,nts_all-1 DO		$
		f1[i,*] = vu_timeseries( vu_get(hdr[i], /array_range), vu_get(hdr[i], /start_distance),	$
			vu_get(hdr[i], /distance_step), ff[*,*,*,*,i], T=ts_all[i],	$
			periodic=periodic, nearL1=nearL1, power=power, body=body, _extra=_extra)
END

ENDCASE

t1 = Carrington(ts_all,/get_time)	; InsituTimeSeries expects time structures

Badf1 = BadValue(f1)
IF nts_all EQ 1 THEN f1 = reform(f1, [nts_all,size(f1,/dim)])

; Kinda kludgy.

tmp = f1			; Contains n,V, or Br,Bt

dim = size(f1,/dim)
dim[1] = vu_fnc(/count)

; Create space for derived quantities by making a 2nd dim of vu_fnc(/count)

f1 = make_array(type=IsType(f1),dim=dim, value=Badf1)

; Put n,V or Br,Bt in the right spot

CASE type[0].prefix OF
vu_prefix (/silent			): f1[*,[vu_fnc_speed,vu_fnc_density],*] = tmp[*,[vu_t3d_speed,vu_t3d_density],*]
vu_prefix (/silent,/magnetic): f1[*,[vu_fnc_brad ,vu_fnc_btang	],*] = tmp[*,[vu_t3d_brad ,vu_t3d_btang  ],*]
ENDCASE

destroyvar, tmp

; Fill in derived quantities for tomographic reconstructions

f1[*,vu_fnc_brtmag  ,*] = vu_fnc(vu_fnc_brtmag  , f1); Br^2+Bt^2 = B^2
f1[*,vu_fnc_bmag    ,*] = vu_fnc(vu_fnc_bmag    , f1); Br^2+Bt^2+Bn^2 = B^2
f1[*,vu_fnc_pdynamic,*] = vu_fnc(vu_fnc_pdynamic, f1); Dynamic pressure

; Pick up the in situ data

cdata = strarr( vu_fnc (/count) )

cdata[ vu_fnc_speed   ] = vu_fnc(vu_fnc_speed	 , /shortid)
cdata[ vu_fnc_density ] = vu_fnc(vu_fnc_density  , /shortid)
cdata[ vu_fnc_brad    ] = vu_fnc(vu_fnc_brad 	 , /shortid)
cdata[ vu_fnc_btang   ] = vu_fnc(vu_fnc_btang	 , /shortid)
cdata[ vu_fnc_bnorm   ] = vu_fnc(vu_fnc_bnorm	 , /shortid)
cdata[ vu_fnc_brtmag  ] = vu_fnc(vu_fnc_brtmag	 , /shortid)
cdata[ vu_fnc_bmag	  ] = vu_fnc(vu_fnc_bmag	 , /shortid)
cdata[ vu_fnc_pdynamic] = vu_fnc(vu_fnc_pdynamic , /shortid)


IF NOT no_insitu THEN	$
	no_insitu = 1-InsituTimeSeries(cdata, t1, t2, f2, delt=delt, source=source, silent=2)

CASE no_insitu OF
0: f2[*,vu_fnc_brtmag ] = vu_fnc (vu_fnc_brtmag, f2); Br^2+Bt^2 (Bn is NOT zero)
1: BEGIN
	t2 = t1
	f2 = replicate(Badf1, nts_all, n_elements(cdata))
END
ENDCASE

IF use_forecast_time THEN BEGIN				; Use forecast times to plot time axis.

	; For a fixed time ts_all (and hence t1 and t2) contains only a single time.
	; Instead we use the forecast times to plot the time axis.
	; For the insitu data we replicate the first data point.

	IF n_elements(uniq(TimeOp(/subtract,ts_all,ts_all[0],uday))) LE 1 THEN BEGIN ; For fixed-time option only
		t1 = Carrington( vu_get(hdr, /forecast_time) )
		t2 = t1
		f2 = replicate(1,nts_all)#f2[0,*]
	ENDIF

ENDIF

; Check for limits on time.
; If /nofuture is specified, then in situ data with times later than
; the forecast time ut are flagged as bad.

IF nofuture THEN ti_range = utt

; If ti_range is specified than all in situ data outside the time range
; specified by ti_range are flagged as bad.

IF IsType(ti_range, /defined) THEN BEGIN
	ti_range = Carrington(ti_range, /get_time)
	InitVar, ti_range, [t1[0], ti_range[0]], count=1
	tt_diff  = TimeOp(/subtract, ti_range[1], ti_range[0], uday, signdiff=tt_sign)
	tt_limit = tt_sign*TimeOp(/subtract, t2, ti_range[0], uday)
	tt_diff  = tt_sign*tt_diff
	tt_limit = where(tt_limit LT 0 OR tt_limit GT tt_diff)
	IF tt_limit[0] NE -1 THEN f2[tt_limit,*] = BadValue(f2)
ENDIF



IF IsType(ascii_file, /string) THEN BEGIN
	CASE ascii_file EQ 'terminal' OF
	0: BEGIN
		SetFileSpec, ascii_file, /parse, status=tmp
		IF tmp THEN BEGIN
			message, /info, 'time series, '+hide_env(ascii_file)
			openw, /get_lun, iu, ascii_file
			FOR i=0,n_elements(t1)-1 DO $
				printf, iu, TimeGet(t1[i], /ymd), TimeGet(t1[i], /doy),	$
					reform(f1[i,*]), format='(A,F12.5,'+strcompress(dim[1],/rem)+'F12.5)'
			free_lun, iu
		ENDIF
	END
	1:	FOR i=0,n_elements(t1)-1 DO $
			print, TimeGet(t1[i], /ymd), TimeGet(t1[i], /doy),	$
				reform(f1[i,*]), format='(A,F12.5,'+strcompress(dim[1],/rem)+'F12.5)'
	ENDCASE
ENDIF



old_device = !D.name

CASE printer OF
0: BEGIN
	IF noplot THEN RETURN
	set_page, /zbuffer, xysize=type[0].xysize, /silent
	reset_colors
END
1: set_page, /printer, aspect=0.89, _extra=_extra
ENDCASE

CASE n_elements(type) EQ 1 OF
0: data_type = round(total(type.data,2))
1: data_type = type.data
ENDCASE

vu_InsituCurve, t1, f1, t2, f2, ut0=utt, sc_label=sc_label,	$
	tm_label=tm_label_, data_type=data_type,				$
	fmin=type.Fmin, fmax=type.Fmax, _extra=_extra,			$
	plot_type=plot_type, silent=silent

display = type[0].display
label   = type[0].label

module = display

IF IsType(cast_offset, /defined) THEN	$
	module = module+'_'+(['a','f'])[cast_offset ge 0]+string(format='(I3.3)',abs(round(24*cast_offset)))

module = module+label

IF NOT plot_type[0] AND plot_type[1] THEN module = module+'_c'

vu_get_page, hdr[ihdr], title='', module, ut=utt, device=old_device, silent=silent, _extra=_extra

RETURN  &  END