;+
; NAME:
;	vu_radialcut
; PURPOSE:
;	Makes cut along raidal line through tomography output
;	from Sun to the position of any of bodies supported by the
;	ephemeris program big_eph.
; CATEGORY:
;	sat/idl/util/vupack
; CALLING SEQUENCE:
	PRO vu_radialcut, hdr, ff, ihdr=ihdr, ut0=ut0,	$
		radius		= radius		, $
		type		= type			, $
		printer		= printer		, $
		silent		= silent		, $
		module		= module		, $
		body		= body			, $
		degrees		= degrees		, $
		textfile	= textfile		, $
		_extra		= _extra
; INPUTS:
;	hdr		array[1]; type: structure
;			    header information about the tomographic reconstruction
;			    as returned from vu_select. If 'hdr' does not exist,
;			    vu_select is called.
;	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:
;	ihdr=ihdr	scalar; type: integer: default: 0
;			    If 'hdr' is an array (and 'ff' the corresponding data arrays)
;			    then a synoptic map for data at index 'ihdr' is displayed.
;
;	ut0=ut0 	scalar or array[1]; type: time structure, float or integer
;				(passed to vu_getdata)
;			    Time (converted to an Earth-based Carrington variable).
;				If not specified then the time in 'hdr' is used.
;			    If a sequence of files from a time-dependent tomography run is
;			    used, and ihdr is not specified, then the hdr closest to ut0 is used.
;
;	type=type	array[1]; type: structure
;			    output structure from vu_type_insitu
;	radius=radius
;				scalar; type: float: default: outer boundary of data sphere
;			    The planarcut covers a square with sides 2*radius with the
;				Sun in the center.
;	/circle		Shows the planar cut always as a circular area with radius
;				specified by keyword 'radius'
;				(For values of radius larger/equal to the radius of the
;				data sphere, the presence of keyword /circle does not
;				make a difference).	
;	zoffset=zoffset
;				scalar; type: float; default: 0.0
;					distance of plane of the cut to the origin (the Sun)
;
;	body=body	passed to PlotPlanarCut
;
;	abg_euler	array[3]; type: float; default: none
;					See PROCEDURE
;       textfile        scalar; type: string
;                           File name with or without path to write time and values to
;
; OUTPUTS:
;       Output is displayed on the screen or written to file.
;       See href=vu_get_page= for information on the keywords involved.
; OPTIONAL OUTPUT PARAMETERS:
;	hdr		array[1]; type: structure
;	ff		array[n,m,l,2]; type: float
;			    Returns headers and arrays read by vu_select and vu_read
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
	@vu_fnc_include.pro 	;
; CALLS:
;	InitVar, ToRadians, vu_getdata, Carrington, TimeGet, TimeUnit
;	vu_gettime, vu_get, vu_type_insitu, gridgen, BadValue
;	AngleRange, vu_atlocation, PlotPlanarCut
;	set_page, vu_get_page
; PROCEDURE:
;	The cut plane is along the equator of the selected spherical
;	coordinate system using one of the from_* keywords
;	to href=CvSky=. E.g. setting
;		/from_centered_ecliptic, /heliocentric
;	will produce a cut along the ecliptic plane with the Earth
;	on the positive x-axis.
;
;	The keyword zoffset is used to displace the cut plane
;	parallel to the equator.
;
;	The Euler angles abg_euler are used to make cuts along
;	planes other the equator. To set up a user-defined cut,
;	define a spherical coordinate system with the user-defined
;	cut plane as equator (i.e., containing x-axis and y-axis
;	(with z-axis completing a right-handed coordinate system).
;	The thus-defined x-axis will be pointing to the right in
;	the plot; the y-xaxis will be up.
;	Then abg_euler are the Euler angles needed to rotate the
;	user-defined coordinate system to the selected CvSky
;	coordinate system.
;
;	For meridional cuts perpendicular to the ecliptic with the
;	Earth on the x-axis use:
;		/from_centered_ecliptic, /heliocentric, $
;		abg_euler=[-90,-90,+90], /degrees
;	The 2nd Euler angle (-90) sets the tilt of the cut plane
;	relative to the ecliptic. E.g. to make a cut at a 30-deg
;	tilt with the ecliptic with the west side of the plane
;	above the ecliptic use -30 deg for the 2nd angle.
;	The 3rd Euler angle (+90 deg) sets the longitude of the
;	nodes of the meridianal plane on the ecliptic.
;	E.g. to make a cut at 30 deg west of the Sun use +120 deg
;	for the 3rd angle.
; MODIFICATION HISTORY:
;	SEP-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, silent 	, 0
InitVar, degrees	, /key
InitVar, body		, 'Earth'

InitVar, printer	, /key

ihdr_set = IsType(ihdr, /defined)

utt = vu_getdata(hdr, ff, ut0=ut0, ihdr=ihdr, _extra=_extra, count=count, silent=silent)
IF count EQ 0 THEN RETURN

; utt determines the data matrix and is used to set the relation between
; the coordinate systems.

utt = Carrington(utt, /get_time)

say, threshold=1, silent=silent, TimeGet(utt,/ymd,upto=TimeUnit(/minute))+strcompress(ihdr)

CASE ihdr_set OF
0: hdri = vu_gettime(hdr, ff, ut=utt, ff_ut=ffi)	; Interpolate at time utt
1: hdri = vu_gettime(hdr, ff, ut=vu_get(hdr[ihdr], /uttime), ff_ut=ffi, /fixgrid); ihdr set explicitly
ENDCASE


; Return spherical heliographic coordinates at time utt for selected body
; (z to heliographic north; x to heliographic longitude zero)

rrbody = big_eph( utt,	$
	center	= 'Sun'	, $
	body	= body	, $
	/to_heliographic, $
	/to_sphere	, $
	/precess	, $
	;degrees	=degrees, $
	/onebody	)

IF NOT finite(rrbody[0]) THEN BEGIN
	say, threshold=2, silent=silent, ['no position for "'+body+'" at '+TimeGet(/ymd,utt),'... ABORTING']
	RETURN
ENDIF

type = vu_type_insitu('planarcut', _extra=_extra, hdr=hdri, type=type)
destroyvar, module

FOR itype=0,n_elements(type)-1 DO BEGIN

	id = (where( type[itype].data ))[0]

	; Extract fields from 'type' structure

	ctable	= type[itype].ctable
	Fmin	= type[itype].Fmin
	Fmax	= type[itype].Fmax
	display	= type[itype].display
	label	= type[itype].label
	format	= type[itype].format
	charsize= type[itype].charsize
	xysize  = type[itype].xysize

	; Done with 'type' structure
	;===========================

	InitVar, radius, vu_get(hdri,/stop_distance)

	; Pick up all radial distances in grid
	; Only keep grid points inside 'radius'

	rx = vu_get(hdri, /all_radii)
	inside = where(rx LE radius)
	IF inside[0] NE -1 THEN rx = rx[inside]

	rr = SuperArray(rrbody,n_elements(rx),/trail)
	rr[2,*] = rx		; Heliographical spherical coordinates

	xci = vu_get(hdri, /array_range	  )
	rri = vu_get(hdri, /start_distance)
	dri = vu_get(hdri, /distance_step )
	tti = vu_get(hdri, /uttime		  )

	zz = vu_atlocation(xci,rri,dri, vu_fnc(id,ffi,/t3d_array),	$
		sequence=tti, location=rr, times=tti, /cv2carrington, /periodic)
	zz = reform(zz)

	old_device = !D.name					; Store current plot device

	;set_page, $
	;	zbuffer		= 1-printer		, $
	;	printer		= printer		, $
	;	/silent						, $
	;	xysize		= xysize		, $
	;	fullsize	= printer		, $
	;	old_device	= old_device	, $
	;	ctable		= ctable		, $
	;	_extra		= _extra

	title  = vu_fnc(id,/symbol,/norm)
	title += ' ('+vu_fnc(id,/unit)+')'
	IF IsType(textfile,/defined) THEN BEGIN
		ix = where(finite(zz))
		IF ix[0] NE -1 THEN BEGIN
			openw, lu, textfile, /get_lun
			printf, lu, transpose(reform([rx[ix],zz[ix]], n_elements(rx[ix]) , 2))
			close, lu
		ENDIF ELSE $
			say, threshold=0, silent=silent, 'No data to write to file'
	ENDIF

	PlotCurve, rx, zz, _extra=_extra
	PlotCurve, rrbody[2]*[1,1], !y.crange, linestyle=0, /oplot

	boost, module, display+label
	;vu_get_page, hdri, title='', module[itype], ut=utt, device=old_device,	$
	;	silent=silent, new=itype ne 0, _extra=_extra

ENDFOR

RETURN  &  END