;+
; NAME:
;	vu_vox_write
; PURPOSE:
;	Convert ascii t3d file to vox file for Volume Pro board
; CATEGORY:
;	sat/idl/util/vupack
; CALLING SEQUENCE:
	PRO vu_vox_write, hdr, ff,		$
		destination = destination	, $
		speed		= speed 		, $
		density 	= density		, $
		brad		= brad			, $
		btang		= btang 		, $
		ndim		= ndim			, $
		range		= range 		, $
		lngrange	= lngrange		, $
		latrange	= latrange		, $
		fix_sun 	= fix_sun		, $
		fix_earth	= fix_earth 	, $
		fix_stars	= fix_stars 	, $
		fix_ihg 	= fix_ihg		, $
		 _extra 	= _extra		, $
		show_sun	= show_sun		, $
		show_earth	= show_earth	, $
		show_orbit	= show_orbit	, $
		show_heliogrid = show_heliogrid, $
		show_eclipgrid = show_eclipgrid, $
		show_hcs	= show_hcs		, $
		obj_val 	= obj_val		, $
		nodata		= nodata		, $
		hdr_bb		= hdr_bb		, $
		bb			= bb			, $
		nearest 	= nearest		, $
		vox_files	= vox_files 	, $
		silent		= silent		, $
		bb_filter	= bb_filter
; 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:
;	ndim=ndim		scalar; type: integer; default: 128
;						resolution of the binary file stored in the vox file
;						(128 or 256)
;	range=range		scalar or array[2]; type: float; default=[0.0,1.5]
;						range of radial distances (AU) written to binary file.
;						Only the part of the array inside range[0],range[1] is
;						extracted and stored in the vox file.
;	lngrange=lngrange scalar or array[2];type;integer;default = [-180,180]
;					  range of longitudes will be passed as a parameter to
;					  functions of the show_heliogrid and show_eclipgrid
;	latrange=latrange scarlar or array[2];type;integer;default = [-90,90]
;					  range of latitudes will be passed as a parameter to
;					  functions of the show_heliogrid and show_eclipgrid
;
;	By default the start of the region of interest of a matrix is put along
;	the x-axis in the vox volume. For an animated sequence of vox files it is
;	important to set one of the following three keywords.
;
;	/fix_sun		makes the x-y-z coordinates of the vox file coincide with
;					heliographic coordinates (i.e. heliographic longitude zero
;					is put along the x-axis)
;	/fix_earth		puts the Sun-Earth direction along the x-axis
;	/fix_stars		puts the heliographic longitude of the equinox (ecliptic
;					longitude zero, ecliptic latitude zero) on the x-axis
;					(i.e. uses a sidereal coordinate frame)
;	/fix_ihg		puts the ascending node of the solar equator on the
;					ecliptic on the x-axis. This also is essentially
;					a sidereal coordinate frame (except for the slow motion
;					of the ascending node over time (1 deg per 72 year).
;
;	/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
;
;	destination=destination
;					scalar; type: string; default: same directory as t3d file
;					Destination directory
;
;	Keywords for merging objects with volume data:
;
;	/show_sun		show Sun at center of vox volume
;	/show_earth		show Earth
;	/show_orbit		show Earth's orbit
;	/show_heliogrid show longitudes and latitudes in heliographic coordinates
;	/show_eclipgrid show longitudes and latitudes in ecliptic coordinates
;	/show_hcs		try to add current sheet
;
;	hdr_bb=hdr_bb, bb=bb
;					as an alternative to setting /show_hcs headers
;					and volume data for the magnetic field can be specified
;					directly using these two keywords.
;
;	obj_val=obj_val	scalar; type: real; default: 1.1 x max value in volume data
;					Function values used to mark objects added to vox volume
;					with show_* keywords. The units are the same as used for
;					the volume data (e.g. km/s for velocity, cm^-3 for number
;					densities, etc.)
;
;	/nodata			if set, the volume data are suppressed (by muliplying with
;					zero). (useful for showing merged objects only)
;
;	Additional keywords:
;
;	Keywords passed to href=vu_getdata= are explained there.
;
;	Keywords passed to href=vox_write= are explained there.
;	The most important ones are.
;
;	minf=minf		these three keywords are used in vox_write to convert
;	maxf=maxf		the volume data scale to byte (they are used as arguments
;	topb=topb		to the idl 'bytscl' function).
;
; OUTPUTS:
;	(*.vox file)
; OPTIONAL OUTPUTS:
;	vox_files=vox_files	array; type: string
;				list of file names of all vox files created
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	vu_getdata, vu_select, AngleRange, SetFileSpec, GetFileSpec, PutFileSpec
;	vox_write, gridgen, TimeSet, TimeUnit, TimeOp, big_eph, jpl_body
;	CvSky, ArrayLocation, Carrington, BadValue, destroyvar, vu_get, vu_prefix
;	vu_vox_drawhcs, vu_vox_draworbit, vu_vox_showgrid,vu_vox_drawsphere
;	InitVar, IsType, vu_type_insitu, boost
; PROCEDURE:
; >	When processing a time sequence of tomography files one of the keywords
;	/fix_sun, /fix_earth, /fix_stars or /fix_ihg should be set, to control
;	the	heliographic longitude that is put along the x-axis in the vox volume.
;	Also the minf, maxf and topb keywords should be used to make sure that the
;	same color scale is used for each of the images.
;
; >	The name of the vox file is based on the name of the t3d file
;	e.g.  t3d_1962.345 will be output as t3d_1962_345.vox
;
;	The selected density or velocity matrix is used to calculate
;	function values on a Cartesian grid by linear interpolation.
; MODIFICATION HISTORY:
;	JAN-2001, Paul Hick (UCSD/CASS)
;	SEP-2001, Kevin Nguyen, Paul Hick (UCSD/CASS)
;		Added options to merge objects (Sun, Earth, Earth's orbit)
;		Added option to add current sheet
;	JUL-2002, Cindy Wang
;		Added options to show grid in heliographic coordination
;		Added options to show grid in ecliptic coordination
;	JUN-2003, Paul Hick (UCSD/CASS)
;		Added postfix to name of output vox file, depending on type
;		of data (density, speed, etc.).
;	OCT-2006, Paul Hick (UCSD/CASS)
;		Replaced NewcombSun by big_eph
;	JUN-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added /fix_ihg
;-

InitVar, silent 	, 0

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

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

InitVar, ndim, 128

InitVar, fix_sun	, /key
InitVar, fix_earth	, /key
InitVar, fix_stars	, /key
InitVar, fix_ihg	, /key

InitVar, show_sun	, /key
InitVar, show_earth	, /key
InitVar, show_orbit	, /key
InitVar, show_hcs	, /key
InitVar, show_heliogrid	, /key ;show grid in heliographic coord
InitVar, show_eclipgrid	, /key ;show grid in ecliptic coord
InitVar, nodata		, /key
InitVar, nearest	, /key

InitVar, range, 1.5
InitVar, range, [0,range], count=1

InitVar, lngrange, 180*[-1,1]
InitVar, lngrange, [0,lngrange], count=1

InitVar, latrange,  90*[-1,1]
InitVar, latrange, [0,latrange], count=1

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

IF count EQ 0 THEN RETURN

IF nodata THEN ff = 0*ff

IF show_hcs THEN BEGIN

	; /show_hcs is set then try to pick up magnetic field data.
	; Find the magnetic field files by replacing 'nv' by 'bb' in hdr.file

	SetFileSpec, vu_get(hdr, /file)

	InitVar, bb_filter, vu_prefix(/magnetic,/silent)

	CASE nearest OF
	0: file = GetFileSpec(upto='directory')+bb_filter+strmid(GetFileSpec(from='name'),4)
	1: BEGIN
		file = (GetFileSpec(upto='directory'))[0]+bb_filter+'*.*'
		file = file_search(file, count=count_bb)
		show_hcs = count_bb GT 0
	END
	ENDCASE

	IF show_hcs THEN BEGIN
		hdr_bb = vu_select(file, /check, /read, count=count_bb, ff=bb)
		show_hcs = count_bb EQ count
	ENDIF

ENDIF

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

IF silent LE 0 THEN BEGIN
	IF fix_sun    THEN message, /info, 'using heliographic reference frame'
	IF fix_earth  THEN message, /info, 'using reference frame with Earth stationary'
	IF fix_stars  THEN message, /info, 'using sidereal reference frame'
	IF fix_ihg	  THEN message, /info, 'using IHG reference frame'

	IF show_sun   THEN message, /info, 'adding sun'
	IF show_earth THEN message, /info, 'adding earth'
	IF show_orbit THEN message, /info, 'adding orbit of earth'
	IF show_hcs   THEN message, /info, 'adding heliospheric current sheet'
	IF show_heliogrid THEN message, /info, 'adding latitude and longitudes in helio coord'
	IF show_eclipgrid THEN message, /info, 'adding latitude and longitudes in eclip coord'
ENDIF

nhdr = n_elements(hdr)
destroyvar, vox_files

FOR ihdr=0L,nhdr-1 DO BEGIN

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

	x = vu_get(hdr[ihdr], /array_range)
	start_lng = 360*(1-(x[0]-floor(x[0])))

	; Zero_lng is the heliographic longitude that is put along the x-axis in the vox volume.

	zero_lng = start_lng
	IF fix_sun   THEN zero_lng = 0
	IF fix_earth THEN zero_lng = 360*(1-(utc-floor(utc)))
	IF fix_stars THEN zero_lng = (CvSky(utt, from_ecliptic=[0.0d0,0.0d0], /to_heliographic, /degrees, /silent))[0]
	IF fix_ihg   THEN zero_lng = (CvSky(utt, from_ihg=[0.0d0,0.0d0], /to_heliographic, /degrees, /silent))[0]

	;off_lng = start_lng-zero_lng

	x = vu_get(hdr[ihdr], /dimension)
	nLng = x[0]
	nLat = x[1]
	nRad = x[2]

	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)

	iRad = round( ((range > R0) < R1)/dR )
	nRad = iRad[1]-iRad[0]+1							; Updated # radial grid points

	R0   = R0+iRad[0]*dR								; Updated inner radius (AU)
	R1   = R0+(nRad-1)*dR								; Updated outer radius (AU)

	; Set up the rectangular coordinates for all the voxels; then convert to spherical
	cntr = 0.5*(ndim-1.0)

	x    = gridgen(ndim*[1,1,1], range=(R1*[-1,1])#[1,1,1]); Cartesian coordinates
	x    = cv_coord(from_rect=x, /to_sphere, /degrees)	; Convert to spherical

	; The longitudes x[0,*] need to be defined relative to the start longitude
	; of the matrix. The longitude zero_lng (the x-axis of the vox volume) is
	; located zero_lng-start_lng from the beginning of the matrix.

	x[0,*] = x[0,*]+(zero_lng-start_lng)

	x[0,*] = AngleRange(x[0,*], /degrees)				; Map to [0,360]

	x[0,*] =  x[0,*]/360*(nLng-1)						; Convert to array indices
	x[1,*] = (x[1,*]+90)/180*(nLat-1)
	x[2,*] = (x[2,*]-R0)/(R1-R0)*(nRad-1)

	F = ff[*,*,iRad[0]:iRad[1],index mod 2,ihdr]
	F = interpolate(F, x[0,*], x[1,*], x[2,*], missing=0.)
	F = reform(F, ndim, ndim, ndim, /overwrite)			; Interpolate on data array

	IF IsType(obj_val, /undefined) THEN obj_val = 1.1*max(F, /nan)
	dvox = (ndim-1.0)/(2*R1)							; # Voxels per AU

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

	EP = big_eph(utt	, $								; Ecliptic location of Earth
		body = jpl_body(/earth,/string) , $
		/to_ecliptic	, $
		/to_sphere		, $
		/degrees		, $
		/precess		, $
		/onebody		, $
		/silent 		)

	EPlng_eclip = EP[0]
	EP = CvSky(utt, from_ecliptic=EP, /to_heliographic, /degrees, /silent)
	EP[0] = EP[0]-zero_lng								; Longitude offset with vox volume
	EPlng_helio = EP[0]

	EP = cv_coord(from_sphere=EP, /to_rect, /degrees)	; Convert to rectangular coordinates

	; The center of the vox volume is at cntr; dvox is the # vox elements per AU.
	; Convert to rectangular coordinates in vox volume.

	EP = cntr+EP*dvox

	IF show_sun THEN	$
		F = vu_vox_drawsphere(F, cntr*[1,1,1], obj_val, /sun)

	IF show_earth THEN	$
		F = vu_vox_drawsphere(F, EP, obj_val)			; Draw sphere at position Earth

	IF show_orbit THEN	$
		F = vu_vox_draworbit(F, utt, hdr[ihdr], zero_lng, cntr, dvox, obj_val)

	IF show_heliogrid THEN	$
		F = vu_vox_showheliogrid(F, EPlng_helio, cntr, dvox, obj_val,lngrange,latrange)

	IF show_eclipgrid THEN	$
		F = vu_vox_showeclipgrid(F, utt, zero_lng, EPlng_eclip,cntr, dvox, obj_val,lngrange,latrange)

	; this part is used to write 3 random points to header file.
	; copied from function draw_orbit

	OP = 366
	OP = gridgen(OP, range=-OP/2+[0,OP-1])
	OP = TimeSet(/diff, OP, TimeUnit(/days))
	OP = TimeOp(/add, utt, OP)						; One year centered on ut in steps of 1 day

	OP = big_eph(OP 	, $							; Ecliptic positions of Earth over one year
		body = jpl_body(/earth,/string) , $
		/to_ecliptic	, $
		/to_sphere		, $
		/degrees		, $
		/precess		, $
		/onebody		, $
		/silent 		)

	; Convert to heliographic locations at the matrix time ut

	OP = CvSky(utt, from_ecliptic=OP, /to_heliographic, /degrees, /silent)

	;OP = OP-[zero_lng,0,vu_get(hdr,/distance_start)]#replicate(1,n)
	OP[0,*] = OP[0,*]-zero_lng

	OP = cv_coord(from_sphere=OP, /to_rect, /degrees)
	OP = OP*dvox+cntr
	OP = OP[*,[0,9,39]]
	OP = reform(OP,n_elements(OP),/overwrite)

	; Merge magnetic data with density in vox volume

	IF show_hcs THEN BEGIN

		CASE nearest OF
		0: itmp = ihdr
		1: tmp  = min( abs(TimeOp(/subtract,vu_get(hdr_bb,/uttime),utt,TimeUnit(/day))), itmp)
		ENDCASE

		B = bb[*,*,iRad[0]:iRad[1],0,ihdr]
		B = interpolate(B, x[0,*], x[1,*], x[2,*], missing=BadValue(B))
		B = reform(B, ndim, ndim, ndim, /overwrite)		; Interpolate on data array

		F = vu_vox_drawhcs(F, B, obj_val)

	ENDIF

	; Write the array F into a vox file

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

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

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

	; The current vx.exe looks for 'Velocity' or 'Density'
	; This needs to be made a little more robust sometime, so this can be removed.

	title = vu_fnc(id_type,/name)
	;title = strmid(title,0,strpos(title,' '))
	title = strupcase(strmid(title,0,1))+strmid(title,1)

	vox_write, vox_file, F, silent=silent+1, _extra=_extra,	$
		volumename		 = title,		$
		general_comments = file ,		$
		volume_comments  = [			$
			'UT '+TimeGet(utt, /ymd),	$
			'Radial range (AU)'		+strjoin(strcompress([R0,R1])),	$
			'Center(CNTR)'			+strcompress(cntr),				$
			'Earth Position(EP)'	+strjoin(strcompress(EP)),		$
			'Oribit 3 random points'+strjoin(strcompress(OP))		]

	boost, vox_files, vox_file

ENDFOR

RETURN  &  END