;+
; NAME:
;	RemoteView_CurrentSheet
; PURPOSE:
;	Find the polygons of the currentsheet surface in cartesian coordinates
;	for a given helographic-sperical B3D matrix.
; CATEGORY:
;  tool
; CALLING SEQUENCE:
	FUNCTION RemoteView_CurrentSheet, FOVInfo, P
; INPUTS:
;	FOVInfo			array[1]; type: structure
;						structure set up by href=RemoteView_Init_FOV=
;	B3D				array[lng,lat,rad]; type: float
;						3D magnetic data;
; OUTPUTS:
;	Result			array[3,*]; type: float
;						polygon coordinates in rectangular coordinates
;						Units are data units (usually AU).
;	P				array[*]; type: integer
;						information about polygons (# points for each polygon,
;						and index lists for all points in each polygon)
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	gridgen, AngleRange, BadValue
; PROCEDURE:
; >	Taking a B3D array in heliographic spherical  convert to cartesian.
; >	Create a currentsheet based on B3D information
; >	Build a contoured surface at that currentsheet out of polygons
; >	Return the list of vertices and polygon-connectivity
; MODIFICATION HISTORY:
;	OCT-2001, T. Dunn (tdunn@ucsd.edu)
;	JUL-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Integrated into RemoteView package
;-

; Set up the rectangular heliographic grid for the mesh.
; Convert to spherical coordinates

nMesh = 16

; Cartesian coordinates in AU; convert to spherical coordinates

x = gridgen(nMesh*[1,1,1], range=(FOVInfo.matrix_rEdge[1]*[-1,1])#[1,1,1])
x = cv_coord(from_rect=x, /to_sphere)

; For the interpolation on the magnetic field matrix the longitudes need
; to be defined relative to the start longitude of the matrix.

x[0,*] = AngleRange(x[0,*]-FOVInfo.matrix_lng0)	; Map to [0,360]


; Convert to array indices needed for the interpolation

sz = size( *FOVInfo.matrixb )

x[0,*] =  x[0,*]/(2*!pi)*(sz[1]-1)
x[1,*] = (x[1,*]+!pi/2)/!pi*(sz[2]-1)
x[2,*] = (x[2,*]-FOVInfo.matrix_rEdge[0])/FOVInfo.matrix_dR

; Interpolate on data array

bb = (*FOVInfo.matrixb)[*,*,*,FOVInfo.matrixb_index]
bb = interpolate(bb,x[0,*], x[1,*], x[2,*], missing=BadValue(bb))
bb = reform(bb, nMesh, nMesh, nMesh, /overwrite)

shade_volume, bb, 0.0, V, P

;whatis, V,P
; V is an array index between 0 and nMesh-1; convert to AU
; V and P may not exist at this point

CASE IsType(V, /defined) OF
0: V = BadValue(0.0)
1: V = (2.0*V/(nMesh-1)-1.0)*FOVInfo.matrix_rEdge[1]
ENDCASE

RETURN, V  &  END