;+
; NAME:
;	RemoteView
; PURPOSE:
;	Calculates image of sky from 3D density matrix for a given viewer location
;	and viewing direction.
; CATEGORY:
;	Volume rendering
; CALLING SEQUENCE:
	PRO RemoteView,					$
		getcube		= getcube		, $
		fovcube		= FOVCube		, $		; Output only when /getcube is set

		cube		= Cube			, $
		drawecliptic= DrawEcliptic	, $
		meshecliptic= MeshEcliptic	, $
		meshhcs		= MeshHCS		, $
		;mergehcs	= MergeHCS		, $

		voxel		= Voxel			, $
		base_name	= base_name		, $
		noplot		= NoPlot		, $

		upto		= upto			, $
		_extra		= _extra		, $

		images		= images		, $

		M3D			= M3D			, $
		V3D			= V3D			, $
		D3D			= D3D			, $

		silent		= silent
; OPTIONAL INPUT PARAMETERS:
;	Most keywords are those described in href=RemoteView_Init=.
;
;	/cube			add bounding cube
;	/drawecliptic	add ecliptic boundary
;	/meshecliptic	add mesh in ecliptic
;	/noplot			suppresses all plotting activity
;	/meshhcs		add heliospheric current sheet as mesh surface
;
;	voxel		array, type: integer; default: 1
;				up to a 4 element array of values 0,1,2 or 3.
;				Values 0,1 and 2 result in a 3D remote view (see
;				href=RemoteView_Display3D=:
;				0: skip rgbo calculation (bodies and orbits only)
;				1: use IDL Voxel_Proj function
;				2: home-brew rgbo rendering (needs work)
;				Value 3 results in a 2D image (see href=RemoteView_Display2D=.
;				The keyword f3dfnc determines the content of the 2D image
;				(the default is the Thomson scattering brightness)
;	base_dir	scalar; type: string; default: getenv('TUB')
;				destination directory for output image files
;	base_name	array; type: string
;				base name used for output image files. A 3-digit counter for
;				each of the views rendered is appended to the base name.
;				Especially if keyword 'voxel' is specified as an array (i.e.
;				multiple images are calculated for each view) then a separate
;				base_name should be specified for each of the 'voxel' entries.
;
;	/getcube	?
;
;	upto = upto	passed to remoteview_display2d and remoteview_display3d
; OUTPUTS:
;	(to screen or image file)
;
;	images=images
;				scalar or array; type: string
;				fully-qualified names of all output images created
;	fovcube=FOVCube
;				Output only when /getcube is set
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	InitVar, boost, destroyvar, IsType, CvSky, gridgen, SuperArray, big_orbit
;	RemoteView_Init, RemoteView_FOV, RemoteView_FOV_loc, RemoteView_Display2D
;	RemoteView_Display3D, RemoteView_CubeFrame, RemoteView_EarthOrbit
;	RemoteView_BodyLoc, RemoteView_BodySize, RemoteView_Destroy
;	RemoteView_Destroy_FOV, RemoteView_ZEclipticPlane, RemoteView_VertexMap
;	RemoteView_CurrentSheet
; PROCEDURE:
; MODIFICATION HISTORY:
;	APR-2001, Paul Hick (UCSD/CASS)
;	JUN-2004, Paul Hick (UCSD/CASS)
;		Added option to plot orbits for planets and Moon (depends on JPL
;		ephemeris). A negative Voxel value now suppresses labeling in
;		RemoteView_Display3D.
;	JUN-2010, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Explicitly added keyword upto=upto.
;-
InitVar, silent 	 , 0
InitVar, cube		 , /key
InitVar, DrawEcliptic, /key
InitVar, MeshEcliptic, /key
InitVar, MeshHCS	 , /key
;InitVar, MergeHCS	 , /key
InitVar, NoPlot		 , /key

DoPlot = 1-NoPlot
InitVar, Voxel, 1

z_buffering = 0 LE abs(Voxel) AND abs(Voxel) LE 2
init_zbuffer = (where(z_buffering))[0] NE -1


RemoteView_Init, M3D, V3D, D3D, silent=silent, _extra=_extra

IF NOT IsType(M3D, /struct) THEN RETURN


t3d, /reset			; Make sure !p.t is set to identity matrix

; Define Z-buffer; make sure to switch back to terminal device before continuing

;IF init_zbuffer*DoPlot THEN BEGIN
;	set_plot,'Z'
;	device, /z_buffering, set_resolution=D3D.xysize; ,set_colors=n_elements(D3D.ct_red)	;!d.n_colors
	set_plot, !TheTerminal
;ENDIF

rOuter = M3D.rEdge[1]

destroyvar, images
nvox  = n_elements(Voxel)
nbase = n_elements(base_name)

IF V3D.nView GT 1 AND nbase EQ nvox THEN base_name = SuperArray(base_name, V3D.nView, /trail)

FOR iView=0,V3D.nView-1 DO BEGIN					; Loop over all views

	; Process view (V3D). Output are the FOVCube and the FOVInfo structure.

	IF RemoteView_FOV(V3D,M3D,FOVCube,FOVInfo,FOVLoc, $
		view		= iview 	, $
		FOVCubeB	= FOVCubeB	, $
		silent		= silent	, $
		_extra		= _extra) THEN BEGIN

		FOR ivox=0,nvox-1 DO BEGIN

			ibase = nvox*iview+ivox
			IF nbase GT ibase THEN vox_name = base_name[ibase] ELSE destroyvar, vox_name

			IF DoPlot THEN BEGIN
				set_plot,'Z'
				device, /z_buffering, set_resolution=D3D.xysize; ,set_colors=n_elements(D3D.ct_red)	;!d.n_colors

				tvlct, D3D.colortable[*,0], D3D.colortable[*,1], D3D.colortable[*,2]

				!p.background	= D3D.White[0]	; Background white
				!p.color		= D3D.Black[0]	; Draw in black

				erase							; Erase with background color
			ENDIF


			TT = FOVInfo.view_time

			; All plotting is done in normalized coordinates withouth data
			; coordinates, i.e. !x.s, !y.s, !x.s are not needed yet, and
			; are established only when dealing with the volume data in
			; RemoteView_Display.

			IF cube THEN BEGIN
				z = RemoteView_CubeFrame(2*rOuter)-rOuter
				z = RemoteView_FOV_loc(z, FOVInfo, /rect)

				IF DoPlot THEN BEGIN
					plots, z[*,*,0], /t3d, color=D3D.Green[0], /norm
					plots, z[*,*,1], /t3d, color=D3D.Green[0], /norm
				ENDIF
			ENDIF

			IF MeshEcliptic THEN BEGIN
				nMesh = 16

				;z = gridgen(nMesh*[1,1], range=rOuter*[[-1,1],[-1,1]])
				;z = RemoteView_ZEclipticPlane(TT, z)
				;z = RemoteView_FOV_loc(z, FOVInfo, /rect)

				x = gridgen(nMesh, range=rOuter*[-1,1])
				y = replicate(rOuter, nMesh)
				z = transpose([[x,-y, x,y], [-y,x, y,x]])

				z = RemoteView_ZEclipticPlane(TT, z)
				z = RemoteView_FOV_loc(z, FOVInfo, /rect)

				IF DoPlot THEN FOR i=0,2*nMesh-1 DO plots, z[*,i+[0,2*nMesh]], /t3d, color=D3D.Black[0], /norm
			ENDIF

			IF DrawEcliptic THEN BEGIN
				z = rOuter*[ [-1,-1], [1,1] ]

				z = RemoteView_ZEclipticPlane(TT, z, /corners)
				z = RemoteView_FOV_loc(z, FOVInfo, /rect)

				IF DoPlot THEN plots, z, /t3d, color=D3D.Black[0], /norm
			ENDIF

			IF MeshHCS AND n_elements(FOVInfo.matrixb_index) NE -1 THEN BEGIN
				z = RemoteView_CurrentSheet(FOVInfo, P)

				IF DoPlot THEN BEGIN
					n = n_elements(P)
					IF n GT 0 THEN z = RemoteView_FOV_loc(z, FOVInfo, /rect)

					i = -1

					; Plot all the polygons separately.
					; P[i] = # vertices in polygon
					; P[i+1:i+P[i]],P[i+1]]: Positions of vertices in vertex array
					; z[*,[P[i+1:i+P[i]],P[i+1]]]: x,y,z coordinates of all vertices in polygon

					WHILE i LT n-1 DO BEGIN
						i += 1
						plots, z[*,[P[i+1:i+P[i]],P[i+1]]], /t3d, color=D3D.Green[0], /norm
						i += P[i]
					ENDWHILE
				ENDIF

			ENDIF

			;z = RemoteView_EarthOrbit(TT)
			;z = RemoteView_FOV_loc(z, FOVInfo, /rect)
			;IF DoPlot THEN plots, z, /t3d, color=D3D.Black[0], /norm


			Loc = cv_coord(from_sphere=FOVInfo.view_loc, /to_rect)

			m = n_elements(D3D.Vertex)/3

			FOR n=0,n_elements(D3D.MarkBody)-1 DO BEGIN; Markers for celestial bodies: Sun, Earth, etc.

				IF D3D.MarkBody[n] THEN BEGIN
														; Heliocentric position in Cartesian heliographic coordinates
					zbody = RemoteView_BodyLoc(TT, body=D3D.MarkName[n], silent=silent)

					; Only insert body and orbit if a finite distance is returned

					IF finite(zbody[2]) THEN BEGIN 

						IF D3D.MarkOrbit[n] THEN BEGIN
							z = big_orbit(TT, body=D3D.MarkName[n], /to_heliographic, /precess, silent=silent)
							z = RemoteView_FOV_loc(z, FOVInfo, /rect)
							IF DoPlot THEN plots, z, /t3d, color=D3D.Black[0], /norm, _extra=_extra
						ENDIF

						; dz exists only for bodies close to Earth (Moon and SGP satellites)

						z = zbody[*,0]
						CASE size(zbody,/n_dim) EQ 2 OF
						0: BEGIN
							destroyvar, dz
							S = z-Loc				; Position relative to viewer location
						END
						1: BEGIN
							dz = zbody[*,1] 		; Heliocentric position Earth
							S = dz-Loc+z			; Position relative to viewer location
						END
						ENDCASE

						S = RemoteView_BodySize(sqrt(total(S*S)), D3D.RealBodyRadius[n],	$
							D3D.FakeBodyRadius[n], rOuter, FOVInfo.view_fov)

						IF S NE 0.0 THEN BEGIN

							CASE 1 OF

							1-finite((RemoteView_FOV_loc(z, FOVInfo, /rect, dloc=dz))[0]): BEGIN
								message, /info, 'bad location returned by RemoteView_FOV_loc for '+D3D.MarkName[n]
							END

							S LT 0.1: BEGIN						; Very small angular size

								z = RemoteView_FOV_loc(z, FOVInfo, /rect, dloc=dz)
								z = SuperArray(z,m,/trail)+S/2*D3D.Vertex

								IF DoPlot THEN BEGIN
									set_shading, values=D3D.BodyColor[*,n]
									z = polyshade(z, D3D.Polygon, /t3d, /norm)
								ENDIF

							END

							D3D.AddBodyMap[n]: BEGIN			; Wrap 2D map around sphere

								CASE D3D.MarkName[n] OF
								jpl_body(/sun,/string): BEGIN	; Wrap 2D map around Sun
									z = D3D.RealBodyRadius[n]*D3D.Vertex
									z = RemoteView_FOV_loc(z, FOVInfo, /rect, dloc=dz)
									zColor = RemoteView_VertexMap(D3D.Vertex,*D3D.BodyMap[n],reverse(D3D.Greyish))
									IF DoPlot THEN z = polyshade(z, D3D.Polygon, /t3d, shades=zColor, /norm)
								END
								jpl_body(/earth,/string): BEGIN	; Wrap 2D map around Earth
									z = SuperArray(z,m,/trail)+D3D.RealBodyRadius[n]*		$
										CvSky(TT,from_geographic=D3D.Vertex, /to_heliographic, /rect,silent=silent)
									z = RemoteView_FOV_loc(z, FOVInfo, /rect)
									zColor = RemoteView_VertexMap(D3D.Vertex,*D3D.BodyMap[n])
									IF DoPlot THEN z = polyshade(z, D3D.Polygon, /t3d, shades=zColor, /norm)
									;zColor = RemoteView_VertexMap(D3D.Vertex,*D3D.BodyMap[1],polygon=D3D.Polygon)
									;IF DoPlot THEN z = polyshade(z, D3D.Polygon, /t3d, poly_shades=zColor, /norm)
								END
								ENDCASE

							END

							ELSE: BEGIN					; Mark with shaded sphere

								; This doesn't work for some reason

								z = SuperArray(z,m,/trail)+D3D.RealBodyRadius[n]*D3D.Vertex
								z = RemoteView_FOV_loc(z, FOVInfo, /rect, dloc=dz)

								IF DoPlot THEN BEGIN
									set_shading, values=D3D.BodyColor[*,n]
									z = polyshade(z, D3D.Polygon, /t3d, /norm)
								ENDIF

							END

							ENDCASE

						ENDIF

					ENDIF

				ENDIF

			ENDFOR

			IF DoPlot THEN BEGIN

				CASE z_buffering[ivox] OF
				0: RemoteView_Display2D, D3D, FOVInfo, FOVCube, FOVLoc,	$
						iView, voxel=Voxel[ivox], base_name=vox_name,	$
						silent=silent, _extra=_extra, images=tmp_images, upto=upto
				1: RemoteView_Display3D, D3D, FOVInfo, FOVCube, FOVLoc,	$
						iView, voxel=Voxel[ivox], base_name=vox_name,	$ 
						silent=silent, _extra=_extra, images=tmp_images,$
						FOVCubeB=FOVCubeB, upto=upto
				ENDCASE

			ENDIF

			boost, images, tmp_images

		ENDFOR

	ENDIF

	RemoteView_Destroy_FOV, FOVInfo

ENDFOR


IF NOT arg_present(M3D) THEN RemoteView_Destroy, M3D=M3D
IF NOT arg_present(V3D) THEN RemoteView_Destroy, V3D=V3D
IF NOT arg_present(D3D) THEN RemoteView_Destroy, D3D=D3D

nimage = n_elements(images)

CASE nimage OF
0:
1: images = images[0]
ELSE: BEGIN
    nimage = nimage/(V3D.nView*nvox)
    images = reform( reform( images, nimage, nvox, V3D.nView ) )
END
ENDCASE

RETURN  &  END