;+
; NAME:
;	vu_correlate
; PURPOSE:
;	Correlate time series for tomography and in situ data;
; CALLING SEQUENCE:
	PRO vu_correlate, t1, f1, t2, f2,	$

		density		= density		, $
		speed		= speed			, $
		bradial 	= bradial		, $
		btang		= btang			, $
		bnorm		= bnorm			, $

		bstrength	= bstrength		, $
		data_type	= data_type		, $

		sc_label	= sc_label		, $
		tm_label	= tm_label 		, $

		result_fit	= result_fit	, $
		result_str	= result_str
; INPUTS:
;	t1			array[N]; type: time structure
;					times of available data points for tomographic data
;	f1			array[N,M]; 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: radial magnetic field component
;					m=5: bt^2+bn^2
;					m=6: total field strength
;
;					Currently only m=0,1,2,3 are provided. The others are
;					calculated internally. The normal component is set to zero.
;
;	t2			array[N2]; type: time structure
;					times of available data point for in situ data
;					(always the same as t1???)
;	f2			array[N2,M]; type: float
;					spacecraft in situ time series with 'bad data' indicated
;					by !values.f_nan containing same quantities as f1.
; OPTIONAL INPUT PARAMETERS:
;	/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
;	/bstrength		select magnetic field strength
;
;	tm_label=tm_label
;				 	scalar; type: string; default: 'corotating'
;					label used for the tomography time series
;	sc_label=sc_label
;					scalar; type: string; default: 'in situ'
;					label used to indicate origin of in situ data
;					(usually this is obtained from href=Instrument= with
;					the keyword /label set).
; OPTIONAL OUTPUT PARAMETERS:
;	result_str=result_str
;					array[n]; type: scalar
;					identifying string of correlation results for all time
;					series for which correlations were calculated.
;	result_fit=result_fit
;					array[6,n]; type: float
;					result[  0,n]: correlation coefficient
;					result[1:5,n]: best fit result from function href=lsqNormalFit=
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, lsqNormalFit, boost
; 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:
;	JUL-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Extracted from PlotInsituCurve (now called vu_insitu_curve)
;-

InitVar, speed		, /key
InitVar, density	, /key
InitVar, bradial	, /key
InitVar, btang		, /key
InitVar, bnorm		, /key
InitVar, bstrength	, /key
InitVar, label		, 'corotating'

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

InitVar, data_type, [speed,density,bradial,btang,bnorm,bstrength]

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

ytitles = [	'velocity (km s!E-1!N)'	,	$
			'density (cm!E-3!N)'	,	$
			'B!Drad!N (nT)'			,	$
			'B!Dtang!N (nT)'		,	$
			'B!Dnorm!N (nT)'		,	$
			'!MS!X(B!Drad!N!E2!N+B!Dtang!N!E2!N) (nT)',	$
			'|B| (nT)']

ytomo = f1		; Contains V,n,Br,Bt

; Add zero normal component for magnetic field

IF (size(ytomo,/dim))[1] EQ 4 THEN boost, ytomo, pushdim=1, plus=1
boost, ytomo, pushdim=1, plus=2

ytomo[*,5] = sqrt( total(ytomo[*,2:3]*ytomo[*,2:3],2, /nan) )	; Tomography: Br^2+Bt^2 = B^2
ytomo[*,6] = sqrt( total(ytomo[*,2:4]*ytomo[*,2:4],2, /nan) )

ysitu = f2		; Contains V,n,Br,Bt,Bn,B

boost, ysitu, pushdim=1, plus=1

ysitu[*,5] = sqrt( total(ysitu[*,2:3]*ysitu[*,2:3],2, /nan) )		; In situ: Br^2+Bt^2 ( Bn is NOT zero)
ysitu[*,6] = f2[*,5]												; In situ: B

here = where(data_type,nhere)

FOR ihere=0,nhere-1 DO BEGIN

	p = here[ihere]

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

	both_data_ok = where( finite(ztomo) AND finite(zsitu) )
	have_data = both_data_ok[0] NE -1

	boost, result_str, ytitles[p]+' ('+label+' vs. '+sc_label+')'

	CASE have_data OF
	0: BEGIN
		message, /info, 'no data: '+ytitles[p]
		boost, result_fit, replicate(BadValue(ztomo),6)
	END
	1: BEGIN
		ztomo = ztomo[both_data_ok]
		zsitu = zsitu[both_data_ok]

		cc = correlate(ztomo, zsitu)
		aa = lsqNormalFit(ztomo, zsitu)

		boost, result_fit, [cc,aa[*,0]]
	END
	ENDCASE

ENDFOR

RETURN  &  END