;+
; NAME:
;	vu_spherecut
; PURPOSE:
;	Makes cut on sun-centered sphere in tomography output
; CATEGORY:
;	sat/idl/util/vupack
; CALLING SEQUENCE:
	PRO vu_spherecut, hdr, ff, ihdr=ihdr, ut0=ut0,	$
		radius		= radius		, $
		type		= type			, $
		printer		= printer		, $
		silent		= silent		, $
		module		= module		, $
		abg_euler	= abg_euler		, $
		degrees		= degrees		, $
		_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 sphere covers a sphere at heliocentric distance
;					of 'radius' AU.
;
;	abg_euler	array[3]; type: float; default: none
;					See PROCEDURE
; 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, PlotSphereCut
;	set_page, vu_get_page
; PROCEDURE:
;	The sphere is at the distance specified in the 'radius' keyword.
;	By default the horizontal is along the ecliptic plane.
;	This can be changed using one of the from_* keywords
;	to href=CvSky=.
; MODIFICATION HISTORY:
;	MAY-2010, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, silent , 0
InitVar, degrees, /key

pi  = !dpi/ToRadians(degrees=degrees)
pi2 = pi/2

InitVar, printer, /key
InitVar, fill	, /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


; Set up three unit vectors along x,y,z axes in the output frame
; (with z the normal to the plane, and x,y axes in the plane)
; The default output frame is the ecliptic seen from the north
; Then rotate these to the heliographic coordinate frame
; (z to heliographic north; x to heliographic longitude zero)

runit = [[1.0,0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]]


IF IsType(abg_euler,/defined) THEN	$
	runit = EulerRotate(abg_euler, runit, degrees=degrees, /rectangular)

; Rotate to heliographic coordinates

runit = CvSky(utt, runit, _extra=_extra, /to_heliographic,	$
	euler_angles=euler_angles, /rectangular, degrees=degrees, euler_info=euler_info, silent=silent)

; If abg_euler exists combine it with the rotation to heliographic coordinates

IF IsType(abg_euler,/defined) THEN	$
	euler_angles = CombineRotations(abg_euler, euler_angles, degrees=degrees)

type = vu_type_insitu('synoptic', _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

	breakval= type[itype].breakval

	i = where(finite(breakval), nbreak)
	IF nbreak GT 0 THEN BEGIN
		breakval = type[itype].breakval[i]
		IF format EQ '' THEN breakval = round(breakval)
	ENDIF

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

	InitVar, radius, vu_get(hdr,/stop_distance)

	height = 0.5*round(0.8*xysize[0])			; Vertical size of map in pixels
	height = (height-1)/2*2+1
	width  = 2*height

	; 3D grid covering sphere in longitude and latitude of selected coordinate system
	; at fixed distance from the sun (radius in AU)

	rr = dblarr(3,width*height,/nozero)
	rr[0:1,*] = gridgen([width,height], range=[0,2*pi,-pi2,pi2])
	rr[2  ,*] = radius

	; Convert to rectangular coordinates in selected coordinate system

	rr = cv_coord(from_sphere=rr,/to_rect,degrees=degrees)

	; Rectangular 3D-vectors in coordinate frame of data array (AU)

	rr = runit#rr

	rr = cv_coord(from_rect=rr, /to_sphere) ; Convert to spherical coordinates
	rr[0,*] = AngleRange(rr[0,*])			; Keep inside [0,2*!pi]

	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,width,height,/overwrite)

	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)+' ('+vu_fnc(id,/unit)+')'

	PlotSphereCut, zz, ut=utt		, $
		radius		= radius		, $
		euler_angles= euler_angles	, $
		degrees		= degrees		, $	
		euler_info	= euler_info	, $
		breakval	= breakval		, $
		charsize	= charsize		, $
		title		= title 		, $
		_extra		= _extra

	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