;+
; NAME:
;	vu_extract
; PURPOSE:
;	Interpolate data in region of interest for points
;	specified in any one of the reference frames
;	supported by CvSky
; CATEGORY:
;	sat/idl/util/vupack
; CALLING SEQUENCE:
	FUNCTION vu_extract, hdr, ff	, $

		speed		= speed 		, $
		density 	= density		, $
		brad		= brad			, $
		btang		= btang 		, $

		rr			= rr			, $

		ndim		= ndim			, $
		range		= range 		, $

		grid		= grid			, $
		coord1		= coord1 		, $
		coord2		= coord2 		, $
		coord3		= coord3 		, $

		rectangular = rectangular	, $
		degrees 	= degrees		, $

		write		= write 		, $
		format		= format		, $
		destination = destination	, $
		outfiles	= outfiles 		, $
		silent		= silent		, $

		 _extra 	= _extra
; INPUTS:
;	hdr				array[1]; type: structure
;						information about the tomographic reconstruction.
;						if 'hdr' does not exist, then user is prompted for
;						a data file using the 'pickfile' dialog.
;	ff				array[n,m,l,2]; type: float
;						3D volume data; if not present then the files
;						specified in the 'hdr' argument are read.
; OPTIONAL INPUT PARAMETERS:
;	/speed			if set, processes the velocity data
;					This is the default if neither /speed or /density are set.
;	/density		if set, processes the density data
;	/bradial		if set, process radial magnetic field
;	/btang			if set, process tangential magnetic field
;
;	The reference frame is selected using one of the from_* keywords
;	of href=CvSky= (e.g. if /from_ihg is set then all coordinates are
;	assumed to be IHG coordinates.
;
;	/degrees		if set, then all angles are in degrees (default: radians)
;	/rectangular	if set, then input coordinates are assumed to be
;						in rectangular coordinates (default: spherical
;						coordinates).
;
;	The heliospheric locations where interpolated data are requested
;	can be specified using anyone of three ways:
;
;	(1) Explicitly (using keyword "rr"):
;
;	rr = rr 	array[3,*]; type: float; default: none
;				Points specified in spherical or rectangular
;				coordinates, depending on setting of /rectangular:
;				If /rectangular NOT set:
;					rr[0,*]: longitudes
;					rr]1,*]: latitudes
;					rr[2,*]; heliocentric distances (AU)
;				If /rectangular SET:
;					rr[0,*]: x-coordinates (AU)
;					rr[1,*]: y-coordinates (AU)
;					rr[2,*]: z-coordinates (AU)
;
;	If rr is defined then the keywords below (ndim,range,coord1
;	coord2,coord3) are ignored.
;
;	(2) Specify the range of the three spatial dimensions
;	    using the keywords coord1, coord2, coord3.
;		When used, all three keyword must be specified.
;
;	coord1 = coord1 scalar or array[n1]; type: float
;				x-coordinates or longitudes of points
;	coord2 = coord2 scalar or array[n2]; type: float
;				y-coordinate of latitudes of points
;	coord3 = coord3 scalar or array[n2]; type: float
;				z-coordinates or heliocentric distances of points
;
;				If /rectangular is SET then rectangular
;				(x,y,z) coordinates are assumed.
;				If /rectangular is NOT set, then spherical
;				(longitude, latitude, heliocentric distance)
;				coordinates are assumed.
;
;				If /grid is SET then coord1,coord2,coord3 are
;				assumed to be 1-dim grids (not necessarily
;				equally spaced) for the three spatial
;				coordinates. The 3-D grid used has dimension
;				[3,n1,n2,n3], combining all elements in
;				coord1, coord2 and coord3.
;	
;				In this case keyword "ndim" is ignored.
;
;				If /grid is NOT set, then coord1, coord2 and
;				coord3 are assumed to be 2-element array
;				representing ranges (start and end values),
;				and the 3-element array "ndim" specifies the
;				number of grid points in each dimension.
;				The 3-D grid has dimensions
;				[3,ndim[0],ndim[1],ndim[2] covering ranges
;				coord1, coord2 and coord3 in each dimension
;
;	(3) Use keywords range and ndim, to set up an equally
;		spaced spherical or rectangular grid.
;
;	ndim = ndim		scalar; array[3]; type: integer
;					grid resolution for three spatial dimensions
;					/rectangular NOT set: default: [72,36,31]
;					/rectangular SET    : default: [31,31,31]
;	range = range 	array[2,3]; ranges for the three spatial dimensions
;					Angles are in units set by keyword /degrees;
;					distances are in AU.
;					/rectangular NOT set:
;						longitude (default: [0,360])
;						latitude (default: [-90,90])
;						heliocentric distance (default: [0,3]
;					/rectangular SET:
;						rectangular x,y,z coordinates
;						(default: [-3,3] AU for all three coordinates)
;
;
;	/write			If set, then write function values into
;					ascii files.
;	destination = destination
;					scalar; type: string; default: same
;						directory as tomography file
;					Destination directory for output files.
;					Only used if /write is set
;
;	Additional keywords:
;
;	Keywords passed to href=vu_getdata= are explained there.
;	/from_* keywords are passed to href=CvSky=.
;
; OUTPUTS:
;	rr = rr 		array[3,*]; type: float
;					grid for which function values are returned.
;					This is returned in the coordinate system
;					implied by the selected from_* keyword to CvSky.
; OPTIONAL OUTPUTS:
;	outfiles=outfiles	array; type: string
;					list of names of all files created
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	vu_getdata, vu_get, AngleRange, boost, gridgen, Carrington
;	SetFileSpec, GetFileSpec, PutFileSpec, TimeSet; CvSky
;	destroyvar, InitVar, IsType, vu_type_insitu, ToDegrees
; EXAMPLE:
;	gg = vu_extract(hdr,ff,coord1=[0,360],coord2=[-90,90],	$
;		coord3=1.0,ndim=[360,180],/degrees,/from_ihg,/density)
;
;	or
;
;	gg = vu_extract(hdr,ff,range=[[0,360],[-90,90],[1,1]],	$
;		ndim=[360,180,1],/degrees,/from_ihg,/density)
;
;	returns a float array of size [360,180,1] for the density
;	at 1 AU for a 1-degree grid in IHG coordinates.
; PROCEDURE:
; >	The name of the output file is based on the name of the
;	tomography file, e.g. t3d_1962.345 will be output as
;	t3d_1962_345.out
;
;	The selected matrix is used to calculate function values
;	by linear interpolation.
; MODIFICATION HISTORY:
;	JUN-2008, Paul Hick (UCSD/CASS)
;-

InitVar, silent 	, 0
InitVar, write		, /key
InitVar, grid		, /key
InitVar, rectangular, /key

InitVar, format, ['F8.3','E12.5']

dpm = ToDegrees(degrees=degrees)

twopi 	= 360.0d0/dpm
pi		= twopi/2.0d0
halfpi	= twopi/4.0d0
zeropi	=   0.0d0/dpm

InitVar, speed		, /key
InitVar, density	, /key
InitVar, brad		, /key
InitVar, btang		, /key

index = 0*speed+1*density+2*brad+3*btang

utt = vu_getdata(hdr, ff, /get_roi, _extra=_extra, count=count)

IF count EQ 0 THEN RETURN, -1

type = vu_type_insitu(density=density, speed=speed, brad=brad, btang=btang)
id_type = (where(type.data))[0]

nhdr = n_elements(hdr)
destroyvar, outfiles

; Set up the output grid

CASE 1 OF

; Points explicitly defined.

IsType(rr,/defined): BEGIN

	szrr = size(rr,/dim)		; Needed to restore input array structure
	IF ssrr[0] NE 3 THEN message, '1st dimension of rr MUST be 3'

	ndim = szrr[1:*]
	rr = reform(rr,3,n_elements(rr)/3,/overwrite)

END

IsType(coord1,/defined): BEGIN

	CASE grid OF

	0: BEGIN

		InitVar, ndim,		  n_elements(coord1)
		InitVar, ndim, [ndim, n_elements(coord2)], count=1
		InitVar, ndim, [ndim, n_elements(coord3)], count=2

		InitVar, coord1, [coord1,coord1], count=1
		InitVar, coord2, [coord2,coord2], count=1
		InitVar, coord3, [coord3,coord3], count=1

		range = [[coord1],[coord2],[coord3]]

		CASE rectangular OF
		0: rr = gridgen(ndim,range=range,open=[1,1,0])
		1: rr = gridgen(ndim,range=range)
		ENDCASE

	END

	1: BEGIN

		ndim = [n_elements(coord1),n_elements(coord2),n_elements(coord3)]

		rr = fltarr([3,ndim[0]*ndim[1]*ndim[2]],/nozero)

		rr[0,*] = SuperArray(coord1,ndim[1]*ndim[2],/trail)
		rr[1,*] = SuperArray(SuperArray(coord2,ndim[0],/lead),ndim[2],/trail)
		rr[2,*] = SuperArray(coord3,ndim[0]*ndim[1],/lead )

	END

	ENDCASE

	szrr = [3,ndim]

END

ELSE: BEGIN

	CASE rectangular OF
	0: BEGIN
		InitVar, range, [[zeropi,twopi],[-halfpi,halfpi],[0,3]]
		InitVar, ndim, [360/5,180/5,1+round(3/0.1)]

		rr = gridgen(ndim,range=range,open=[1,1,0])
	END
	1: BEGIN
		InitVar, range, [[-3,3],[-3,3],[-3,3]]
		InitVar, ndim, (1+round(3/0.1))*[1,1,1]

		rr = gridgen(ndim,range=range)
	END
	ENDCASE

	szrr = [3,ndim]

END

ENDCASE

whatis, rr

formatn = '('+strcompress(ndim[0])+format+')'

FOR ihdr=0L,nhdr-1 DO BEGIN

	utt = vu_get(hdr[ihdr], /uttime)

	; Start lng (heliographic) of data array

	start_lng = (vu_get(hdr[ihdr], /array_range))[0]
	start_lng = carrington(start_lng,/get_longitude,degrees=degrees)

	nLng = vu_get(hdr[ihdr],/nlng)
	nLat = vu_get(hdr[ihdr],/nlat)
	nRad = vu_get(hdr[ihdr],/nrad)

	dR = vu_get(hdr[ihdr],/distance_step)
	R0 = vu_get(hdr[ihdr],/start_distance)			; Inner radius (AU)
	R1 = R0+(nRad-1)*dR								; Outer radius (AU)

	; Convert input grid to heliographic coordinates

	CASE rectangular OF
	0: vv = rr
	1: vv = cv_coord(from_rect=rr,/to_sphere,degrees=degrees)
	ENDCASE

	vv = CvSky(utt,vv,/to_heliographic,degrees=degrees,_extra=_extra)

	; The longitudes vv[0,*] need to be defined relative to
	; the start longitude of the matrix.

	vv[0,*] -= start_lng
	vv[0,*]  = AngleRange(vv[0,*], degrees=degrees)	; Map to [0,360]

	vv[0,*] =  vv[0,*]/twopi*(nLng-1)				; Convert to array indices
	vv[1,*] = (vv[1,*]+halfpi)/pi*(nLat-1)
	vv[2,*] = (vv[2,*]-R0)/(R1-R0)*(nRad-1)

	F = ff[*,*,*,index mod 2,ihdr]
	F = interpolate(F, vv[0,*], vv[1,*], vv[2,*], missing=0.0)
	F = reform(F, ndim, /overwrite)

	IF silent LE 1 THEN message, /info, $			; Display time
		vu_fnc(id_type,/name)+' @ '+TimeGet(utt,/ymd)

	IF write THEN BEGIN 		; Write the array F into an ascii file

		file = vu_get(hdr[ihdr], /file)
		SetFileSpec, file

		IF IsType(destination,/defined) THEN	$
			PutFileSpec, filepath(root=destination,''), upto='DIR'

		outfile =  GetFileSpec(upto='name')   +'_'+	$
			strmid(GetFileSpec(part='type'),1)+type.label+'.out'

		IF silent LE 1 THEN message, /info, hide_env(outfile)

		openw, iu, /get_lun, outfile
		printf, iu, ndim
		FOR i=0,2 DO printf, iu, rr[i,*], format=formatn[0]
		printf, iu, f, format=formatn[1]
		free_lun, iu

		boost, outfiles, outfile

	ENDIF

ENDFOR

rr = reform(rr,szrr,/overwrite)

RETURN, F  &  END