;+
; NAME:
;	vu_planarcut
; PURPOSE:
;	Makes planar cut through tomography output
; CATEGORY:
;	sat/idl/util/vupack
; CALLING SEQUENCE:
	PRO vu_planarcut, hdr, ff, ihdr=ihdr, ut0=ut0,	$
		radius		= radius		, $
		zoffset 	= zoffset		, $
		type		= type			, $
		printer		= printer		, $
		silent		= silent		, $
		module		= module		, $
		body		= body			, $
		abg_euler	= abg_euler		, $
		degrees		= degrees		, $
		edge		= edge			, $
		circle		= circle		, $
		_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
; 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, printer	, /key
InitVar, fill		, /key
InitVar, earth		, /key
InitVar, circle 	, /key

ihdr_set = IsType(ihdr, /defined)

utt = vu_getdata(hdr, ff, ut0=ut0, ihdr=ihdr, _extra=_extra, count=count, silent=silent+1)
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, TimeGet(utt,/ymd,upto=TimeUnit(/minute))+strcompress(ihdr), threshold=0, silent=silent

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

IF IsType(edge,/defined) THEN BEGIN
	;iedge = round( (edge-hdri.distance_start)/hdri.distance_step )
	;ffi[*,*,iedge+1:*,*] = BadValue(ffi)

        iedge = round( (edge-vu_get(hdri,/start_distance))/vu_get(hdri,/distance_step) )+1
        IF 0 LE iedge AND iedge LT vu_get(hdri,/nrad) THEN ffi[*,*,iedge:*,*] = BadValue(ffi)

ENDIF


; 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 abg_euler exists then apply that rotation first

IF IsType(abg_euler,/string) THEN BEGIN

	; This is a fudge to allow making a planar cut through the orbital
	; plane of a body whose orbital elements are stored in
	; $EPHEM/orbits/orbital_elements.txt (mostly comets, I think).

	tmp_body = abg_euler
	tmp = strpos(tmp_body,'&')
	IF tmp NE -1 THEN tmp_body = strmid(tmp_body,tmp)

	tmp = filepath(root=getenv('EPHEM'),subdir='orbits','orbital_elements.txt')
	IF flt_read(tmp,mask=strjoin(replicate('1',20)),data,crumbs=crumbs,/double) THEN BEGIN
		crumbs = strtrim(reform(crumbs[0,*]))		; Body name
		tmp = (where( strlowcase(tmp_body) EQ strlowcase(crumbs) ))[0]
		IF tmp NE -1 THEN BEGIN
			boost, body, abg_euler
		    print, data[5:6,tmp]			; Asc node and inclination
		    abg_euler = [-90,-data[6,tmp],90-data[5,tmp]]/ToDegrees(degrees=degrees)
		ENDIF
	ENDIF

	IF IsType(abg_euler,/string) THEN destroyvar, abg_euler

ENDIF

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('planarcut', _extra=_extra, hdr=hdri, type=type)

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)
	InitVar, zoffset, 0.0

	diameter = round(0.8*xysize[0])			; Diameter of cross section in pixels
	diameter = (diameter-1)/2*2+1
											; Rectangular 2D-vectors in plane of cross section
	rr = gridgen(diameter*[1,1], range=[-1.0,1.0,-1.0,1.0])
											; Rectangular 3D-vectors in coordinate frame of data array
	rr = (zoffset*runit[*,2])#replicate(1,diameter*diameter)+radius*runit[*,0:1]#rr

	IF circle THEN outside = where(total(rr*rr,1) GT radius*radius)

	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, diameter, diameter, /overwrite)
	IF circle THEN IF outside[0] NE -1 THEN zz[outside] = BadValue(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)+')'

	PlotPlanarCut, zz, ut=utt		, $
		radius		= radius		, $
		euler_angles= euler_angles	, $
		degrees		= degrees		, $	
		euler_info	= euler_info	, $
		earth		= earth 		, $
		breakval	= breakval		, $
		charsize	= charsize		, $
		title		= title 		, $
		body		= body			, $
		_extra		= _extra

	IF IsType(module, /undefined) OR n_elements(type) EQ 2 THEN BEGIN
		IF IsType(module, /defined) AND n_elements(type) EQ 2 AND itype EQ 0 THEN destroyvar, module
		boost, module, display+label
	ENDIF
        
	vu_get_page, hdri, title='', module[itype], ut=utt, device=old_device,	$
		silent=silent, new=itype ne 0, _extra=_extra

ENDFOR

RETURN  &  END