;+
; NAME:
;	RemoteView_FOV_xyz
; PURPOSE:
;	Internal use. Should be called by href=RemoteView_FOV_Cube= only.
; CATEGORY:
;	RemoteView project
; CALLING SEQUENCE:
	FUNCTION RemoteView_FOV_xyz, x, FOVInfo, bigmemory=bigmemory,	$
		fov=fov, eqtrans=eqtrans
; INPUTS:
;	x			float array[4,nDimX,nDimY,nDimZ]
;					Initialized in calling program RemoteView_FOV_Cube with:
;						x = fltarr(4,nDimX,nDimY,nDimZ, /nozero)
;					Only x[0:2,*,*,*] is filled here.
;	FOVInfo		array[1]; type: structure
;					structure set up by href=RemoteView_Init_FOV=
; OPTIONAL INPUT PARAMETERS:
;	/bigmemory	if set then no explicit do-loops are used in the calculations.
;				This will be significantly faster if lots of memory is available
;				(i.e. as long as no disk swapping is needed)
; OUTPUTS:
;	Result		float array[4,nDimX*nDimY*nDimZ]
;					Cartesian coordinates of locations where interpolated
;					densities are required. The coordinate system is the X,Y,Z
;					coordinate system defined in the RemoteView document.
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	RemoteView_FOV_Edge, RemoteView_EqTrans, gridgen, IsType, InitVar
; PROCEDURE:
;	If the output array is reformed to Result[4,nDimX,nDimY,nDimZ] then the nDimX and
;	nDimY dimensions cover field of view. The nDimZ dimension is the 'depth' dimension.
;	It stores the distance to the viewer location starting at the maximum distance from
;	the viewer, i.e. R[*,*,*,0] is the depth layer farthest away from the viewer along
;	the line of sight; R[*,*,*,nDimZ-1] is closest to the viewer.
; MODIFICATION HISTORY:
;	FEB-1998, Paul Hick (UCSD/CASS)
;	DEC-2000, Paul Hick (UCSD/CASS, pphick@ucsd.edu), adapted for handling of
;		time-dependent tomography
;-

InitVar, bigmemory, /key

IF IsType(FOVInfo,/defined) THEN BEGIN
	fov		= FOVInfo.view_fov
	eqtrans = FOVInfo.view_eqtrans
	extent	= FOVInfo.view_extent[[1,0]]
ENDIF

fovedge = RemoteView_FOV_Edge(fov, eqtrans)

r = size(x)

nDimR = r[1]
nDimX = r[2]
nDimY = r[3]
IF r[0] EQ 4 THEN nDimZ = r[4] ELSE nDimZ = 1

; The field of view is covered by a regular grid of dimension nDimX x nDimY
; Create r[2,nDimX,nDimY] array. r[0,*,*] and r[1,*,*] contain the
; x- and y-coordinates, respectively, across the field of view.

r = fltarr(2,nDimX*nDimY, /nozero)
r[0,*] = replicate(1.,nDimX)#gridgen(nDimY, range=fovedge[0]*[-1,1])
r[1,*] = gridgen(nDimX, range=fovedge[1]*[-1,1])#replicate(1.,nDimY)

;r = fovedge*r							; 2 x (nDimX * nDimY) grid

CASE eqtrans OF
RemoteView_EqTrans(/eqdist): begin			; Preserve distance
	ss = 1.0/sqrt(1+total(r*r,1))
	s  = ss
end
RemoteView_EqTrans(/eqangle): begin
	s = sqrt(total(r*r,1))				; [nDimX*nDimY] array

	ss = s								; [nDimX*nDimY] array
	i = where(s NE 0)
	IF i[0] NE -1 THEN ss[i] = sin(s[i])/s[i]
	i = where(s eq 0)
	IF i[0] NE -1 THEN ss[i] = 1.0

	s = cos(s)							; [nDimX*nDimY] array
END
RemoteView_EqTrans(/eqarea): BEGIN
	s  = total(r*r,1)
	ss = sqrt(1.-s/4.)
	s  = 1.-s/2.
END
ENDCASE


IF nDimZ EQ 1 THEN BEGIN
	x = reform(x,nDimR,nDimX*nDimY, /overwrite)

	x[0,*] = ss*r[0,*]
	x[1,*] = ss*r[1,*]
	x[2,*] = s

ENDIF ELSE BEGIN

	; The radial distance range FOVInfo.view_extent is covered by an equally
	; spaced grid of nDim elements. Note that rr[0] corresponds to the
	; largest distance from the viewer location.

	; Each grid point in the fov combines with each radial element to define
	; a 3-D grid.

	; Convert to the Cartesian X-Y-Z coordinate system (this amounts to
	; reversing a perspective transformation).

	; The Z-component is calculated first.
	; The X- and Y-component scale proportionally with the Z-component
	; The unit of length is the unit used for FOVInfo.view_extent.

	z = gridgen(nDimZ, range=extent)

	IF bigmemory THEN BEGIN

		x = reform(x,4,nDimX*nDimY*nDimZ, /overwrite)
		OneZ = replicate(1.,nDimZ)

		x[2,*] = s#z

		ss = ss#z
		x[0,*] = ss*(reform(r[0,*])#OneZ)	; All [nDimX*nDimY,nDimZ] arrays
		x[1,*] = ss*(reform(r[1,*])#OneZ)

	ENDIF ELSE BEGIN

		x = reform(x,4,nDimX*nDimY,nDimZ, /overwrite)

		;tmp = [SuperArray(ss,2, /lead)*r[0:1,*]. tramspose(s)]
		;for i=0,nDimZ-1 do x[*,*,i] = tmp*z[i]

		FOR i=0,nDimZ-1 DO BEGIN
			x[2,*,i] = s*z[i]

			tmp = ss*z[i]
			x[0,*,i] = tmp*r[0,*]
			x[1,*,i] = tmp*r[1,*]
		ENDFOR

		x = reform(x,4,nDimX*nDimY*nDimZ, /overwrite)

	ENDELSE

ENDELSE

RETURN, X  &  END