;+
; NAME:
;	vu_solardisk
; PURPOSE:
;	Projects tomographic data at fixed distance from the Sun onto the
;	plane of the sky.
; CATEGORY:
;	Tomography
; CALLING SEQUENCE:
	PRO vu_solardisk, hdr, ff, ihdr=ihdr, ut0=ut0, radius=radius, type=type,	$
		hardcopy=hardcopy, focus_center=focus_center, degrees=degrees,			$
		silent=silent, module=module, _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
;			    Time (converted to an Earth-based Carrington variable)on which
;			    the map is centered. 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: 1 AU
;			    heliocentric distance for synoptic map (AU)
;			    The radius is rounded to the nearest grid distance in the 3D
;			    heliospheric matrix.
; OUTPUTS:
;	(plot to screen)
; 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, TimeGet, TimeUnit
;	vu_gettime, vu_get, vu_type_insitu, gridgen, BadValue, big_eph, jpl_body
;	CvSky_Heliographic, CvT3D, AngleRange, vu_atlocation, PlotSolarDisk
;	set_page, vu_get_page
; PROCEDURE:
; MODIFICATION HISTORY:
;	AUG-2002, Paul Hick (UCSD/CASS)
;	JAN-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Minor bug fix. The diameter variable passed to PlotSolarDisk was misspelled.
;-

InitVar, silent 	, 0

InitVar, hardcopy	, /key
InitVar, fill		, /key
InitVar, earth		, /key
InitVar, north		, /key
InitVar, south		, /key

IF IsType(focus_center, /defined) THEN	$
	focus = focus_center*ToRadians(degrees=degrees)

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 = Carrington(utt, /get_time)

IF silent LE 1 THEN message, /info, 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
												; Pick up a default 'type' if none is input
type = vu_type_insitu('earth_solardisk', _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
	;===========================

	diameter = 0.8*xysize[0]

	; x,y coordinates over a square from -radius to + radius.

	InitVar, radius, vu_get(hdri, /start_distance)
	IF radius EQ 0.0 THEN message, 'no source surface distance found'

	rr = gridgen(diameter*[1,1], range=radius*[-1,1,-1,1])
	r_xy = total(rr*rr,1)				; Distance to disk center

	rr[*,where(r_xy gt radius*radius)] = BadValue(rr)	; Set pixels outside radius bad

	inside = where( finite(rr[0,*]) )	; Pixels on disk

	zz = reform(rr[0,*])				; Z-values across disk
	zz[inside] = sqrt( (radius*radius-r_xy[inside]) > 0 )

	boost, rr, zz						; x,y,z, coordinates

	rr = rr[*,inside]

	; In the disk image the x-y axes define the screen
	; In the real world this corresponds to a heliocentric spherical
	; coordinate system:

	CASE IsType(focus, /defined) OF

	0: BEGIN

		; The default is the disk as seen from Earth, so the +z axis
		; is the vector pointing from Sun to Earth.

		; xyz_ecl = Euler angles from x-y-z to ecliptic coordinates
		; ecl_hel = Euler angles from ecliptic to heliographic
		; rot = rotation sequence from x-y-z to heliographic

		; z-axis is Sun-Earth direction.
		; x-axis points west
		; y-axis is to ecliptic north.

		xyz_ecl = big_eph(utt,		$
			body = jpl_body(/earth,/string)	, $
			/to_ecliptic, $
			/to_sphere	, $
			/precess	, $
			/onebody	, $
			/silent 	)

		xyz_ecl = -[!pi/2, !pi/2, xyz_ecl[0]]
		ecl_hel = CvSky_Heliographic(utt)
		rot = [0,0,xyz_ecl[0],0,xyz_ecl[1:2],0,0,ecl_hel[0],0,ecl_hel[1:2]]

		earth = 1

	END

	1: BEGIN

		; 'focus' defines the direction of the +z axis in heliographic
		; coordinates (heliographic longitude and latitude)
		; The rotation 'rot' rotates from these x,y,z coordinates to
		; heliographic coordinates

		; z = out of the screen = direction to viewer (default: Sun-Earth direction)
		; x = right = points west from Sun center
		; y = up    = toward heliographic north

		rot = [0,0,-!pi/2, 0, focus[1]-!pi/2, -focus[0]]

	END

	ENDCASE

	; Rotate vector 'rr' from x,y,z to heliographic coordinates, and
	; convert from rectangular to spherical coordinates

	rr = CvT3D(rot=rot, vector=rr)
	rr = cv_coord(from_rect=rr, /to_sphere, /double)

	rr[0,*] = AngleRange(rr[0,*])		; Keep inside [0,2*!pi]
	rr[2,*] = rr[2,*] > radius			; Needed at source surface

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

	zz[inside] = 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)

	old_device = !D.name						; Store current plot device

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

	reset_colors, /reverse						; White foreground, black background
	stretch, !d.n_colors-1, 0					; Invert color table

	PlotSolarDisk, zz, ut=utt, diameter=diameter, radius=radius, rot=rot, earth=earth,	$
		breakval=BreakVal, xysize=xysize, charsize=charsize, _extra=_extra, $
		title=vu_fnc(id,/symbol,/norm)+' ('+vu_fnc(id,/unit)+')'

	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