;+
; NAME:
;	vu_insitucurve
; PURPOSE:
;	Plots time series for tomography and in situ data;
;	and/or correlation plots of tomography vs. in situ data
; CALLING SEQUENCE:
	PRO vu_InsituCurve, ttomo, ytomo, tsitu, ysitu,	$
		data_type	= data_type		, $

		timeseries	= timeseries	, $
		correlation	= correlation	, $

		ut0			= ut0			, $
		tx_range	= tx_range		, $
		forecast	= forecast		, $

		fmin		= fmin			, $
		fmax		= fmax			, $

		upto		= upto			, $
		sc_label	= sc_label		, $
		tm_label	= tm_label		, $
		digits		= digits		, $

		_extra		= _extra		, $

		plot_type	= plot_type 	, $
		vert_ut 	= vert_ut		, $
		colored		= colored	
; INPUTS:
;	ttomo		array[N]; type: time structure
;					times of available data points for tomographic data
;	ytomo		array[N1,M,K]; type: float
;					tomographic time series.with 'bad data' indicated by
;					!values.f_nan. The second dimension (m) identifies the
;					physical quantity
;
;					m=0: density
;					m=1: velocity
;					m=2: radial magnetic field component
;					m=3: tangential magnetic field component
;					m=4: normal magnetic field component
;					m=5: bt^2+bn^2
;					m=6: total field strength
;					m=7: dynamic pressure
;
;					The 3rd dimension identifies multiple sets of tomographic
;					data (this is currently only used to compare timeseries from
;					the nv3d files with the ea files).
;
;	tsitu		array[N2]; type: time structure
;					times of available data point for in situ data
;					(always the same as ttomo???)
;	ysitu		array[N2,M]; type: float
;					spacecraft in situ time series with 'bad data' indicated
;					by !values.f_nan containing same quantities as ytomo.
; OPTIONAL INPUT PARAMETERS:
;	ut0=ut0		scalar; type: integer or time structure
;					If absent, then array 'ttomo' defines the timeaxis
;					If present, 'ttomo' determines the range of the time axis,
;					but the time axis will be centered on 'ut0'.
;					If 'ut0' is specified as an integer then it is interpreted
;					as the difference between local time and UT in hours
;	tx_range=tx_range
;				array[2]; float or time structure
;					defines the range of the time axis. A float array is
;					interpreted as a time range in days centered on 'ut0'
;	/forecast	only used if keyword 'ut0' is set.
;					Adds a vertical dashed line at the center of the time
;					series (marking time 'ut0'), and adds the word 'FORECAST'
;					to the right of the dashed line.
;
;	If 'ut0' is absent then the time series reflects the time range of the
;	'ttomo' array.
;
;	fmin	= fmin	scalar, array[2]; type: any
;					scalar is same as [scalar,scalar]
;	fmax	= fmax	scalar, array[2]; type: any
;					scalar is same as [scalar,scalar]
;					sets minimum and maximum on the two timeseries plotted
;
;	Speed/density plots cannot be mixed with magnetic field plots.
;
;	/speed			select velocity plots
;	/density		select density plots
;	/bradial		select radial component of magnetic field
;	/btang			select tangential component of magnetic field
;	/bnorm			select normal component of magnetic field
;	/brtmag
;	/bmag			select magnetic field strength
;	/pdynamic		select dynamic (ram) pressure
;
;	data_type=data_type
;					array[n] (n=8 currently); type: integer
;					array with 0 (data type not plotted) or 1 (data type plotted)
;					If not specified the array is constructed from the data type
;					keywords listed above.
;					If specified it usually is the tag entry with the same name
;					from the structure defined by href=vu_type_insitu=.
;
;	/timeseries		plot time series for selected quantities
;	/correlation	plot correlations only
;					(if both are omitted than time series and correlation are
;					plotted).
;
;	sc_label=sc_label
;					scalar; type: string
;					label used to indicate origin of in situ data
;					(usually this is obtained from href=Instrument= with
;					the keyword /label set).
;	tm_label=tm_label
;					scalar or array; type: string; default: ['corotating', 'time dependent']
;					labels used for the tomography time series
;
;	_extra=_extra	IDL graphics keyword passed to plot statement
;					(charsize, thick, etc.)
; OPTIONAL OUTPUT PARAMETERS:
;	tx_range=tx_range
;					array[2]; time structure
;					actual time range used to plot time axis
;	plot_type=plot_type
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
	@vu_fnc_include.pro 	;
; CALLS:
;	InitVar, lsqNormalFit, TimeXAxis, PlotCurve, Carrington, Instrument
;	gridgen, IsTime, TimeUnit, TimeGet, TimeSet, TimeOp, TimeLimits
;	boost, anicenum
; PROCEDURE:
;	The 'tomo' arrays are usually extracted from a tomography output file
;	The 'insitu' arrays is set up by a call to InsituTimeseries
; MODIFICATION HISTORY:
;	NOV-1999, Paul Hick (UCSD/CASS)
;	JUL-2002, Paul Hick (UCSD/CASS)
;		Added keyword /forecast (before this was functionality was
;		included with keyword 'ut0'
;	AUG-2002, Paul Hick (UCSD/CASS)
;		Added keyword tx_range. Also modified the calculation of tx_range
;		if it is not specified.
;	SEP-2002, Paul Hick (UCSD/CASS)
;		Renamed PlotCompareInsitu to PlotInsituCurve and revamped.
;	APR-2004, Paul Hick (UCSD/CASS)
;		Bug fix (call to TimeAdd instead of TimeOp)
;		Added keyword digits to pass info to function 'anicenum'.
;	JUL-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Restructured the plotting of the timeseries to provide better
;		plots if no tomography data are available. It used to make a
;		black window; now it plots 'no data'.
;		Removed keyword noplot.
;		Removed keywords result_str and result_fit returning infro from
;		correlation between tomographic and in situ timeseries (now done
;		by vu_correlate.
;		Renamed PlotInsituCurve to vu_InsituCurve.
;-

InitVar, timeseries	, /key
InitVar, correlation, /key

InitVar, tm_label, 'corotating'
InitVar, tm_label, [tm_label,'time dependent'], count=1
InitVar, digits, 4

; The order in this list must be coordinated with vu_type_insitu.

;IF IsType(data_type, /undefined) THEN BEGIN

;	data_type = intarr(vu_fnc(/count))

;	data_type[ vu_fnc_speed   ] = IsType( speed	  , /defined )
;	data_type[ vu_fnc_density ] = IsType( density , /defined )
;	data_type[ vu_fnc_bradial ] = IsType( bradial , /defined )
;	data_type[ vu_fnc_btang   ] = IsType( btang   , /defined )
;	data_type[ vu_fnc_bnorm   ] = IsType( bnorm   , /defined )
;	data_type[ vu_fnc_brtmag  ] = IsType( brtmag  , /defined )
;	data_type[ vu_fnc_bmag	  ] = IsType( bmag	  , /defined )
;	data_type[ vu_fnc_pdynamic] = IsType( pdynamic, /defined )

;ENDIF

ndata = n_elements(data_type)
IF total(data_type) EQ 0 THEN data_type = replicate(1,ndata)

plot_type = [timeseries, correlation]
IF total(plot_type) EQ 0 THEN plot_type = replicate(1,n_elements(plot_type))
timeseries  = plot_type[0]
correlation = plot_type[1]

aplot = data_type#plot_type

; Set up the plot layout (!p.multi)

p = total(aplot)
nx = ([0,1,2,2,2,2,2,2,2])[p]
ny = ([0,1,1,2,2,3,3,4,4])[p]

p_multi_old = !p.multi
!p.multi = [0,nx,ny]


InitVar, forecast, /key

; # tomographic time series (3rd dimension of ytomo)

ntomo = size(ytomo)
IF ntomo[0] EQ 2 THEN ntomo = 1 ELSE ntomo = ntomo[3]

x0t = 0.03	; 0.65	;0.4
x0c = 0.60
y0t = 0.90
y0c = 0.15
xw  = 0.12
yh  = 0.075
xs  = 0.20

InitVar, colored, !d.name NE !ThePrinter AND !d.name NE 'PS'
IF NOT colored THEN message, /info, 'black-and-white output'

; Colors used for adding tomographic time series (first entry not used)

more_color = [0, 22, 7, 36, 83, 127]
more_color = more_color*colored+!p.color*(1-colored)

; Color used for spacecraft insitu time series

sc_color = !p.color

IF colored THEN loadct, 12, /silent	; Load color table only if colors needed

here = where(aplot[*,0],nhere)

IF nhere GT 0 THEN BEGIN
	InitVar, fmin, -1
	InitVar, fmax, -1

	InitVar, fmin, replicate(fmin[0], nhere), count=1
	InitVar, fmax, replicate(fmax[0], nhere), count=1
ENDIF

uday = TimeUnit(/day)
umin = TimeUnit(/minute)

FOR ihere=0,timeseries*nhere-1 DO BEGIN

	p = here[ihere]

	ztomo = ytomo[*,p,0]
	zsitu = ysitu[*,p  ]

	ysymbol 	= vu_fnc(p, /symbol)
	yunit		= vu_fnc(p, /unit  )
	ytitle		= vu_fnc(p, /symbol)+' ('+yunit+')'
	tomo_color	= vu_fnc(p, /color )*colored+!p.color*(1-colored)

	; Set up time range for x-axis.

	CASE IsType(ut0, /defined) OF

	0: BEGIN

		title  = ''
		IF NOT IsTime(tx_range) THEN tx_range = TimeLimits(ttomo, /range)

	END

	1: BEGIN

		ut0 = Carrington(ut0, /get_time)
		InitVar, upto, umin
		title  = TimeGet(ut0,/ymd,upto=upto)+' UT'

		; If time range for x-axis is specified as a time structure, use it unmodified.
		; If it is a time difference structure, use it with ut0.
		; If it is defined, but is not a structure, interpret as a time
		; time difference in uday, and use it with ut0.

		CASE 1 OF
		IsTime(tx_range			):
		IsType(tx_range,/defined): tx_range = TimeOp(/add, ut0, TimeSet(/diff,tx_range,uday))
		ELSE: BEGIN

			; If tx_range not defined than set it such that ut0 is centered
			; in the plot, and tx_range brackets all times in ttomo.

			tx_range = TimeOp(/subtract, TimeLimits(ttomo, /range), ut0, uday)
			tx_range = max(abs(tx_range))*[-1,1]*( 2*(tx_range[0] LT 0)-1 )
			tx_range = TimeOp(/add, ut0, TimeSet(/diff, tx_range, uday))
		END
		ENDCASE

	END

	ENDCASE

	data_ok = where( finite(ztomo) )	; Check whether tomography data are available
	have_data = data_ok[0] NE -1

	; Plot the time axis, even when no data are available.
	; This will update !p.multi and defines !x.crange, and !y.crange (which are needed
	; even when no time series are plotted.

	TimeXAxis, /exact, tx_range, uday, /noxtitle, upto=uday, _extra=_extra

	; If there are no tomography data then erase the time axis again

	IF NOT have_data THEN BEGIN

		;message, /info, 'no data: '+ytitle

		n  = !p.multi[1]*!p.multi[2]-!p.multi[0]
		xn = (n-1) mod !p.multi[1]
		yn = !p.multi[2]-1-(n-1)/!p.multi[1]

		dx = !d.x_size/!p.multi[1]
		dy = !d.y_size/!p.multi[2]

		tv, replicate(!p.background,dx,dy),xn*dx, yn*dy ; Clear plot area

	ENDIF

	; Plot time label at top of plot
	; If in addition keyword 'forecast' is set, plot FORECASE in right half of plot

	IF IsTime(ut0) THEN BEGIN
		xyouts, mean(!x.crange), !y.crange[1]+0.04*(!y.crange[1]-!y.crange[0]), title, align=0.5, _extra=_extra
		IF forecast AND have_data THEN	$
			xyouts, total(([0,1]+0.25*[1,-1])*!x.crange),	$
					total(([0,1]+0.50*[1,-1])*!y.crange), 'FORECAST', align=0.5, _extra=_extra
	ENDIF

	CASE have_data OF

	0: BEGIN						; No tomography data available

		xn = (xn+0.5)/!p.multi[1]
		yn = (yn+0.5)/!p.multi[2]
		xyouts, xn, yn, /norm, align=0.5, 'no data: '+ytitle

	END

	1: BEGIN						; Tomography data available

		situ_ok  = where( finite(zsitu) )
		IF situ_ok [0] EQ -1 AND IsType(sc_label, /defined) THEN	$
			message, /info, 'no '+sc_label+' '+ytitle+' data'

		; Set up range of y-axis using input fmin, fmax (if available).

		yrange = [0,0]

		CASE fmin[ihere] EQ -1 OF
		0: yrange[0] = fmin[ihere]
		1: BEGIN
			yrange[0] = min(ztomo[data_ok],/nan)
			IF situ_ok[0] NE -1 THEN yrange[0] = min([yrange[0],zsitu[situ_ok]],/nan)
		END
		ENDCASE

		CASE fmax[ihere] EQ -1 OF
		0: yrange[1] = fmax[ihere]
		1: BEGIN
			yrange[1] = max(ztomo[data_ok],/nan)
			IF situ_ok[0] NE -1 THEN yrange[1] = max([yrange[1],zsitu[situ_ok]],/nan)
		END
		ENDCASE

		;PlotCurve, /oplot, /newyaxis, ttomo, ztomo, ytitle=ytitle,$
		PlotCurve, /oplot, /newyaxis, ttomo, ztomo, ytitle=ytitle,$
			yrange=yrange, color=tomo_color, _extra=_extra	; Plot 1st tomography series

		FOR i=1,ntomo-1 DO PlotCurve, /oplot, ttomo, ytomo[*,p,i],	$
			color=more_color[i], _extra=_extra			; Plot other tomography series

		;IF IsTime(ut0 THEN PlotCurve, /oplot, [ut0,ut0], !y.crange, linestyle=2
		IF IsTime(ut0) AND forecast THEN BEGIN			; Dashed line down center of plot
			n  = 20
			dx = [ut0,ut0]
			dy = gridgen(n, range=!y.crange)
			FOR i=0,n-1,2 DO PlotCurve, /oplot, dx, dy[i:i+1], /silent
		ENDIF

		FOR ii=0,n_elements(vert_ut)-1 DO BEGIN
			n  = 20
			dx = [vert_ut[ii],vert_ut[ii]]
			dy = gridgen(n, range=!y.crange)
			FOR i=0,n-1,2 DO PlotCurve, /oplot, dx, dy[i:i+1], /silent
		ENDFOR

		dx = !x.crange[1]-!x.crange[0]
		dy = !y.crange[1]-!y.crange[0]

		b  = !x.crange[0]+x0t*dx
		w  = xw*dx
		h  = !y.crange[0]+y0t*dy+yh*dy

		IF situ_ok[0] NE -1 AND IsType(sc_label, /defined) THEN BEGIN ; Overplot s/c data if present

			PlotCurve, /oplot, tsitu, zsitu, linestyle=2*(1-colored), color=sc_color, _extra=_extra

			; Correlation makes only sense if ttomo = tsitu

			h = h-yh*dy
			plots , b+[0.0,w],h*[1.0,1.0], linestyle=2*(1-colored), color=sc_color, _extra=_extra
			xyouts, b+(1+xs)*w,h, sc_label, _extra=_extra

			both_data_ok = where( finite(ztomo) AND finite(zsitu) )
			IF both_data_ok[0] NE -1 THEN BEGIN
				tmp = abs(ztomo[both_data_ok]-zsitu[both_data_ok])
				xyouts, !x.crange[1]-x0t*dx, h, align=1.0, _extra=_extra,	$
					'!4D!X'+ysymbol+' = '+	$
					anicenum( mean(tmp), digits=digits, /astring )+$
					' '+yunit
			ENDIF

		ENDIF

		h = h-yh*dy
		plots , b+[0.0,w],h*[1.0,1.0], color=tomo_color, _extra=_extra
		xyouts, b+(1+xs)*w,h, tm_label[0], _extra=_extra

		FOR i=1,ntomo-1 DO BEGIN
			IF total(finite(ytomo[*,p,i])) NE 0 THEN BEGIN
				h = h-yh*dy
				plots , b+[0.0,w],h*[1.0,1.0], color=more_color[i], _extra=_extra
				xyouts, b+(1+xs)*w,h, tm_label[i], _extra=_extra
			ENDIF
		ENDFOR

	END

	ENDCASE

ENDFOR

CASE IsType(sc_label, /defined) OF
0: nhere = 0
1: here = where(aplot[*,1],nhere)
ENDCASE

FOR ihere=0,correlation*nhere-1 DO BEGIN

	p = here[ihere]

	ytitle = vu_fnc(p,/symbol)+' ('+vu_fnc(p,/unit)+')'

	ztomo = ytomo[*,p,0]
	zsitu = ysitu[*,p]

	data_ok = where( finite(ztomo) and finite(zsitu) )

	CASE data_ok[0] NE -1 OF
	0: BEGIN

		i = !p.multi[1]*!p.multi[2]
		!p.multi[0] = (!p.multi[0]-1+i) mod i

		n  = !p.multi[1]*!p.multi[2]-!p.multi[0]
		xn = (n-1) mod !p.multi[1]
		yn = !p.multi[2]-1-(n-1)/!p.multi[1]

		dx = !d.x_size/!p.multi[1]
		dy = !d.y_size/!p.multi[2]

		tv, replicate(!p.background,dx,dy),xn*dx, yn*dy

		xn = (xn+0.5)/!p.multi[1]
		yn = (yn+0.5)/!p.multi[2]

		xyouts, xn, yn, /norm, align=0.5, 'no data: '+ytitle

	END

	1: BEGIN

		ztomo = ztomo[data_ok]
		zsitu = zsitu[data_ok]

		cc = correlate(ztomo, zsitu)
		range = [min([zsitu,ztomo],/nan),max([zsitu,ztomo],/nan)]
		title = 'correlation'+strcompress(string(cc,format='(F10.3)'))

		a = lsqNormalFit(ztomo, zsitu,	/plot,	$
			title	= title,	$
			xtitle	= tm_label[0]+' '+ytitle,	$
			ytitle	= sc_label   +' '+ytitle,	$
			xrange	= range,	$
			yrange	= range,	$
			_extra	= _extra)

		oplot, !x.crange,!y.crange, linestyle=2

	END

	ENDCASE

ENDFOR

!p.multi = p_multi_old

RETURN  &  END