;+
; NAME:
;	RemoteView_Init_Matrix
; PURPOSE:
;	Defines structure for storing 3D matrix and associated information
; CATEGORY:
;	Sky project: structure definition
; CALLING SEQUENCE:
	FUNCTION RemoteView_Init_Matrix,	$
		f3dmat		= F3Dmat,		$
		b3dmat		= B3Dmat,		$
		matrix_time	= matrix_time,	$
		matrix_rect	= matrix_rect,	$
		volsize		= volsize,		$
		start_lng	= start_lng,	$
		degrees		= degrees,		$

		; keyword view_time and rot_cme are accessed only if /cme is set.
		; They SHOULD NOT BE MODIFIED.

		cme			= cme,			$
		view_time	= view_time,	$
		rot_cme		= rot_cme,		$
		silent		= silent
; INPUTS:
;
;	F3Dmat	array[N,L,M,K]; type: float
;			array[N,L,M  ] (K  = 1)
;				Specifies a sequence of 3D heliospheric volume data to be rendered.
;				The array must have at 3 or 4 dimensions. The first three are the
;				spatial dimensions; the last usually is the time dimension. The matrix
;				can be specified on a Cartesian or a spherical grid, but always refers
;				to heliographic coordinates:
;
;				Each matrix F3Dmat[*,*,*,k] (k=0,K-1) represents a spatial
;				scalar solar wind quantity at time TT[k].

;				The matrix covers a whole Carrington rotation (360 degree
;				in heliographic longitude) starting at longitude start_lng (see below)
;
;				Rectangular: N=L=M; cube of data with the location of the Sun
;						as implied by 'volsize' keyword (see below).
;				Spherical: N=# longitudes covering 360 deg
;						   L=# latitudes  covering -90,90 deg
;						   M=# heliocentric distances covering [volsize[0],volsize[1]]
;				The Sun is located at radial distance 0.
;
;	B3Dmat	array[N,L,M,K]; type: float
;			array[N,L,M  ] (K  = 1)
;				Radial component of the magnetic field, corresponding to F3Dmat volume date.
;
; OPTIONAL INPUT PARAMETERS:
;	Parameters used to define the volume data matrix.
;
;	matrix_time=matrix_time
;			array[1] or array[K]; type: time structure
;				times for each of the K 3D matrices in F3Dmat
;				For K=1 no time needs to be specified.
;
;	/matrix_rect
;			if set then F3Dmat refers to a Cartesian (rectangular) coordinate system
;				(default is a spherical coordinate system).
;
;	volsize=volsize
;			scalar, array[2]; type: float; default: [0,1.5]
;				A scalar is interpreted as [0,volsize]
;
;				Cartesian: volsize[1] is half the size of the cube covered by F3Dmat
;					(volsize[0] is not used)
;				Spherical: heliocentric distance range covered by F3Dmat
;
;	volsize=volsize
;			Determines the spatial extent of the volume covered by F3Dmat (in AU).
;			If tomography files are read then this keyword can be used to limit the extent.
;			of the volume rendered.
;
;			If /matrix_rect NOT set (spherical data set)
;				scalar or array[2]; type: float; default: [0,1.5]
;				A scalar is interpreted as [0,volsize].
;				volsize[0] and volsize[1] are the inner and outer radial boundary of F3Dmat
;
;			If /matrix_rect set (rectangular data set)
;				scalar, array[2], array[3], array[2,3]
;				volsize determines the edges of the cube in x, y and z dimension, respectively.
;				-A scalar volsize is interpreted as
;				 [[-volsize,volsize], [-volsize,volsize], [-volsize,volsize]] i.e. the
;				 origin of the coordinate system is in the center of a cube with sides 2*volsize.
;				- An array[2] is interpreted as
;				 [[-volsize[0],volsize[1]], [-volsize[0],volsize[1]], [-volsize[0],volsize[1]]]
;				 i.e. the volume is a cube with sides volsize[1]-volsize[0], but the origin is
;				 not necessarily in the center.
;				-An array[3] is interpreted as
;				 [[-volsize[0],volsize[0]], [-volsize[1],volsize[1]], [-volsize[2],volsize[2]]]
;				 i.e. the size of the volume is different for all dimensions, but the
;				 origin is still in the center of the volume.
;
;	start_lng=start_lng
;			scalar or array[K]; type: float: default: 0.0
;				start heliographic longitude of the matrices in F3Dmat
;				If K > 1 then start_lng can be an array with different start
;				longitudes for each of the K matrices in F3Dmat.
;				(if scalar is used then the same value is used for all K matrices).
;
;	/degrees	if set then start_lng is assumed to be in degrees.
;
;	cme=cme		if set then href=RemoteView_CMEDensity= is called to add the density of a CME to
;				the input F3Dmat density matrix. This is an old trick which superposes
;				a 3D HELIOS CME density on the input volume data.
;				if /matrix_rect is set then this keyword is ignored
;
;				If /cme is set then:
;
;				The CME density matrix is assigned the same time as the matrix time.
;				The CME is expanded from matrix_time to view_time (usually the view times
;				will be a couple of days later than the matrix time).
;				rot_cme are two angles (in radians) used to rotate the CME matrix to
;					the desired part of the heliosphere
; OUTPUTS:
;	M3D				array[1]; type: Sky_Matrix3D structure
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; RESTRICTIONS:
;	For matrices referring to Cartesian coordinates the grid spacings in all three dimensions
;	must correspond to the same physical distance in AU. So it is not permitted e.g to use
;	a grid spacing of 0.2 AU for the x-dimension and 0.1 for y and z.
; CALLS:
;	InitVar, ToRadians, TimeSystem, IsType, BadValue, SuperArray, AngleRange
;	remoteview_CMEDensity
; PROCEDURE:
; >	The 3D matrix F3Dmat is specified in a heliographic coordinate system, with
;	x and y axis in the solar equator and z-axis toward solar north. The Sun is
;	located in the origin. The heliographic longitude of the x-axis is given by
;	'start_lng'.
;
; >	If only one spatial matrix is specified (i.e. K=1) then the time keyword need not be specified
;	(unless /cme is set)
;
; >	The structure used to carry around information about the 3D matrix.
;		time		array[K]		time of matrix		(pointer to heap variable)
;		F3Dmat		array[n,m,l,K]	the matrix itself	(pointer)
;		Lng0		array[K]		heliographic longitude at start of matrix (pointer)
;
;		dR			float			spatial resolution (AU/grid spacing)
;		rEdge		fltarr(2)		inner and outer edge (grid spacings)
;		Origin		float(3)		location of Sun (grid spacings) (0,0,0) for spherical matrix)
;		rectangular	integer			0: spherical coordinates; 1: Cartesian coordinates
;
; MODIFICATION HISTORY:
;	MAY-2000, Paul Hick (UCSD/CASS)
;	DEC-2000, Paul Hick (UCSD/CASS; pphick@ucsd.edu), adapted to handle a sequence of
;		matrices from time-dependent tomography
;-

InitVar, matrix_rect, /key
InitVar, cme		, /key
InitVar, silent 	, 0

cme = cme AND NOT matrix_rect

InitVar, matrix_time, TimeSystem(/utc,/silent)  ; If no time specified then use system time
InitVar, start_lng	, 0.0

rpm = ToRadians(degrees=degrees)

nsize = size(F3Dmat, /dim)

; nsize  : the spatial dimensions of the density matrix
; nMatrix: the number of density matrices

CASE n_elements(nsize) OF
3: nMatrix = 1
4: BEGIN
	nMatrix = nsize[3  ]
	nsize   = nsize[0:2]
	IF nMatrix EQ 1 THEN F3Dmat = reform(F3Dmat, /overwrite); Remove degenerate trailing dimension
END
ELSE: message, 'invalid matrix for volume data'
ENDCASE

IF IsType(volsize,/undefined ) THEN rEdge = 1.5 ELSE rEdge = volsize
IF IsType(rEdge, /generic_int) THEN rEdge = float(rEdge)

CASE matrix_rect OF
0: BEGIN
	CASE n_elements(rEdge) OF					; Edges of data cube (AU)
	1   : rEdge = [0,rEdge[0]]
	ELSE: rEdge = rEdge[0:1]
	ENDCASE

	dR		= (rEdge[1]-rEdge[0])/(nsize[2]-1)	; # AU per radial grid unit

	Origin	= replicate(BadValue(0.0),3)		; Not used
END
1: BEGIN
	CASE n_elements(rEdge) OF
	1	: rEdge = [[-rEdge   ,rEdge   ], [-rEdge   ,rEdge   ], [-rEdge   ,rEdge   ]]
	2   : rEdge = [[-rEdge[0],rEdge[1]], [-rEdge[0],rEdge[1]], [-rEdge[0],rEdge[1]]]
	3	: rEdge = [[-rEdge[0],rEdge[0]], [-rEdge[1],rEdge[1]], [-rEdge[2],rEdge[2]]]
	ELSE: rEdge = reform(rEdge[0:5],2,3)
	ENDCASE

	; Here we assume that all three dimensions use the same scale (AU/grid spacing)

	dR	   = (rEdge[1,0]-rEdge[0,0])/(nsize[0]-1); # AU per grid unit

	Origin = reform(-rEdge[0,*]*(nsize-1)/(rEdge[1,*]-rEdge[0,*])) ; Origin in grid spacings

	; Get the radius of a sphere containing all the data.

	rEdge = [[0,0,0],[0,1,0],[1,0,0],[1,1,0],[0,0,1],[0,1,1],[1,0,1],[1,1,1]]
	rEdge = rEdge*SuperArray(nsize-1,8)-SuperArray(Origin,8)
	rEdge = sqrt(total(rEdge*rEdge,1))
	rEdge = [0, max(rEdge)]*dR					; Edges of data cube (AU)
END
ENDCASE

IF silent LE 0 THEN $
	message, /info, (['spherical','rectangular'])[matrix_rect]+ $
		' data; radial range ['+strjoin(strcompress(rEdge,/rem),',')+'] AU'

; If keyword /cme is set then more than one time may have been specified in view_time, while
; only a single 3D matrix is specified in F3Dmat.

IF cme THEN BEGIN
	IF nMatrix EQ 1 AND n_elements(view_time) GT 1 THEN BEGIN	; More than one view time specified
		nMatrix = n_elements(view_time)
		F3Dmat = SuperArray(F3Dmat, nMatrix, /trail)
		matrix_time = replicate(matrix_time, nMatrix)
	ENDIF

	cme_average = mean(F3Dmat, /nan)
	cme_smooth  = 9

	FOR i=0,nMatrix-1 DO $
		F3Dmat[*,*,*,i] += remoteview_CMEDensity(matrix_time[i], view_time[i], nsize, rEdge, $
			rot_cme=rot_cme, degrees=degrees, cme_average=cme_average, cme_smooth=cme_smooth)

	matrix_time = view_time
ENDIF


; A magnetic field must have the same spatial dimensions as the volume data.
; If more than one 3D magnetic field matrix is specified, than it must be the same
; number as the number of time steps for the volume data. If this is not the case
; than only the first magnetic field matrix is used. The others are tossed.

; nsize  : the spatial dimensions of the magnetic field matrix
; nMatrix: the number of density matrices

nsizeB = size(B3Dmat, /dim)

CASE size(B3Dmat, /n_dim) OF
3: BEGIN
	nMatrixB = 1
	IF (where(nsize NE nsizeB))[0] EQ -1 THEN Bmat = B3Dmat
END
4: BEGIN
	nMatrixB = nsizeB[3  ]
	nsizeB   = nsizeB[0:2]
	IF (where(nsize NE nsizeB))[0] EQ -1 THEN BEGIN

		CASE nMatrixB OF
		1      : Bmat = reform(B3Dmat, /overwrite); Remove degenerate trailing dimension
		nMatrix: Bmat = B3Dmat
		ELSE   : BEGIN
			nMatrixB = 1
			Bmat = B3Dmat[*,*,*,0]
			message, /info, 'inconsistent magnetic field data; only retaining first 3D-matrix'
		END
		ENDCASE

	ENDIF
END
ELSE: nMatrixB = 0
ENDCASE

CASE nMatrixB OF
0   : B3D = ptr_new()
ELSE: B3D = ptr_new(Bmat, /no_copy)
ENDCASE


lng0 = AngleRange(start_lng*rpm)			; Convert to radians
IF n_elements(lng0) EQ 1 AND nMatrix GT 1 THEN lng0 = replicate(lng0[0], nMatrix)

; Set up the structure

rtn = {SMEIremoteMatrix}

rtn.nmatrix		= nMatrix
rtn.nmatrixb	= nMatrixB

rtn.time		= ptr_new(matrix_time, /no_copy)	; Pointer to time (time structure)
rtn.F3D			= ptr_new(F3Dmat, /no_copy)			; Pointer to 3D matrix for volume data
rtn.B3D			= B3D								; Pointer to magnetic field data, or null pointer
rtn.lng0		= ptr_new(lng0  , /no_copy)			; Heliographic longitude at start of matrix (radians)

rtn.dR			= dR								; Spatial resolution (AU)

rtn.rEdge		= rEdge								; Spatial boundaries (inner and outer edge in AU)
rtn.Origin		= Origin							; Origin = Location of Sun (grid spacings)

rtn.rectangular	= matrix_rect						; 0: spherical coordinates; 1: Cartesian coordinates

RETURN, rtn  &  END