;+
; NAME:
;	RemoteView_FOV_loc
; PURPOSE:
; CATEGORY:
;	3D viewing
; CALLING SEQUENCE:
	FUNCTION RemoteView_FOV_loc, vector, FOVInfo, $
		rectangular = rectangular	, $
		dloc		= dloc
; INPUTS:
;	vector			array[3,*]; type: float
;						heliographic coordinates (spherical or, if
;						/rectangular, Cartesian) of points to be
;						converted to display coordinates.
;						Distances are in data units (usually AU).
;	FOVInfo			array[1]; type: structure
;						structure set up by href=RemoteView_FOV=
; OPTIONAL INPUT:
;	/rectangular	if set then 'vector' is in rectangular heliographic
;						coordinates (default is spherical coordinates).
; OUTPUTS:
;	Result			arrays[3,*]
;						input vector points converted to sky coordinates.
;						These are normal coordinates, i.e. points inside
;						the fov range from 0 to 1.
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	cvt3d, RemoteView_EqTrans, RemoteView_FOV_Edge, SyncDims
; PROCEDURE:
; >	Read the RemoteView document (remoteview.doc) for
;	a detailed description.
; >	The 3D data set is specified in an arbitrary Cartesian
;	coordinate system x-y-z.
; >	The data set is contained within a sphere of radius DataRadius
;	centered on the origin.
; >	Put a rectangle with sides 2*tan(Fov) in a plane
;	perpendicular to the viewer direction at a distance of one unit.
;	Define a rectangular grid of dimension nDimX x nDimY covering
;	the rectangle. Connect the origin with each of the nDimX x nDimY
;	grid points to define an angular grid.
; >	Define a radial grid of dimension nDimZ covering the range
;	FOVInfo.view_extent
; >	The radial and angular grid define a 'sky grid' covering the
;	portion of space along the viewing direction which may
;	contain data.
; STRUCTURES:
;	FOVInfo.view_loc, FOVInfo.view_dir,  FOVInfo.view_fov,
;	FOVInfo.view_tilt, FOVInfo.view_eqtrans, FOVInfo.view_extent
; MODIFICATION HISTORY:
;	FEB-1998, Paul Hick (UCSD/CASS)
;	DEC-2000, Paul Hick (UCSD/CASS)
;		Rewrite
;	SEP-2008, Paul Hick (UCSD/CASS, pphick@ucsd.edu)
;		Bug fix?. The calculation of the z-component
;		(depth) for the FOV coordinates did not
;		preserve the sign of the z-component of the
;		native camera coordinates system. As a result
;		point "behind" the camera where put "in front"
;		of the camera.
;		Coordinates for points "behind" the camera are
;		now explicitly set to NaN in the return array.
;-

InitVar, rectangular, /key

x = reform(vector, 3, n_elements(vector)/3)
IF NOT rectangular THEN x = cv_coord(from_sphere=x, /to_rect)

view_eqtrans = FOVInfo.view_eqtrans
view_loc	 = cv_coord(from_sphere=FOVInfo.view_loc,/to_rect)
view_dir	 = FOVInfo.view_dir
view_tilt	 = FOVInfo.view_tilt
view_extent  = FOVInfo.view_extent
view_fovedge = RemoteView_FOV_Edge(FOVInfo.view_fov, view_eqtrans)

; Rotate to native camera coordinate system

InitVar, dloc, 0
x = cvt3d(translate = view_loc-dloc	, $
		  rotate	= [0.0,0.0,!pi+view_dir[0], 0.0,-(!pi/2.0-view_dir[1]),view_tilt], $
		  /xyexch	, $
		  vector 	= x)

; z is the distance from the observed
; Preserve the sign of x[2,*] to be able to
; distinguish "in front" and "behind" the camera.

z = sqrt(total(x*x,1))*(2*(x[2,*] GT 0)-1)

CASE view_eqtrans OF

RemoteView_EqTrans(/eqdist): BEGIN

	xy = x[2,*]

	;i = where(xy eq 0.0 and finite(xy))
	;IF i[0] ne -1 THEN xy = -1.0				; Avoid divide by zero

	x[0:1,*] /= view_fovedge#xy					; Range inside fov: [-1,1]

	; The above transformation.projects points (X,Y,Z) and -(X,Y,Z) to the
	; same location. To avoid problems we exclude all points with negative z.

	i = where(xy LE 0.0)
	IF i[0] NE -1 THEN x[0:1,i] = BadValue(x)	;2*x[0:1,i]/abs(x[0:1,i])

END

RemoteView_EqTrans(/eqangle): BEGIN

	xy = sqrt(total(x[0:1,*]*x[0:1,*],1))

	i = where(xy NE 0.0)
	IF i[0] NE -1 THEN xy[i] = atan(xy[i],x[2,i])/xy[i]
	i = where(xy eq 0.0)
	IF i[0] NE -1 THEN xy[i] = 1.0

	x[0:1,*] *= (1.0/view_fovedge)#[xy]			 ; Range inside fov: [-1,1]

END

RemoteView_EqTrans(/eqarea): BEGIN

	xy = sqrt(total(x[0:1,*]*x[0:1,*],1))

	i = where(xy NE 0.0)
	IF i[0] NE -1 THEN xy[i] = sqrt(2*(1.0-x[2,i]/z[i]))/xy[i]
	i = where(xy EQ 0.0)
	IF i[0] NE -1 THEN xy[i] = 1.0

	x[0:1,*] *= (1.0/view_fovedge)#[xy]			; Range inside fov: [-1,1]

END

ENDCASE

; Map to normal coordinates (range [0,1])
; z=0 corresponds to view_extent[1]
; i.e. lies farthest away from the viewer.
; z=1 corresponds to view_extent[0] (usually zero),
; i.e. z>1 is "behind" the camera.

x[0:1,*] = (x[0:1,*]+1)*0.5
x[2  ,*] = (z-view_extent[1])/(view_extent[0]-view_extent[1])

; Not sure this is a good idea.
; The plots function with /normal set will still plot
; points with z > 1 (behind the camera).
; We don't want that, so we set points with z > 1
; to bad here.

i = where(x[2,*] GT 1.0)
IF i[0] NE -1 THEN x[*,i] = BadValue(x)

SyncDims, x, size=size(vector)

RETURN, x  &  END