;+
; NAME:
;	RemoteView_ZEclipticPlane
; PURPOSE:
;	Get coordinate arrays describing ecliptic plane in
;	a Cartesian heliographic coordinate system.
; CATEGORY:
;	Sky graphics
; CALLING SEQUENCE:
	FUNCTION RemoteView_ZEclipticPlane, TT, X, $
		corners = Corners	, $
		grid	= Grid		, $
		angle	= Angle 	, $
		degrees = Degrees	, $
		color	= color
; INPUTS:
;	TT				array[1]; type: standard time structure;
;						time
;	X				array[2,*]
;						X and Y coordinates of points in the solar
;						equatorial plane.
; OPTIONAL INPUT PARAMETERS:
;	/corners		if set then only the first 2 points X[*,0:1] are used. The output
;						arrays will be an array[3,5], in the sequence
;						(0,0),(1,0),(1,1),(0,1),(0,0), i.e. the Z array can be used to
;						draw a square representing the ecliptic plane.
;	/grid			if set, every value in X[0,*] is combined with every X[1,*]. This
;						defines a square grid of points covering the ecliptic plane
;	angle=Angle		scalar; default: 0
;						the input coordinates are rotated over Angle in the solar
;						equatorial plane before the Z-coordinates are calculated
;	/degrees		if set, Angle is in degrees (default: radians)
; OUTPUTS:
;	ZEcliptic		array[*], array[*,*]
;	X				array[*], array[*,*]
;						X,Y,Z coordinates of points in ecliptic plane
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	ToRadians, CvSky
; PROCEDURE:
;	The X-Y-Z coordinate system is associated with heliographic coordinates:
;	X-Y plane is the solar equator; X-axis has heliographic longitude zero;
;	Z-axis points to solar North pole.
;	Ecliptic north has has ecliptic latitude !pi/2 (the longitude doesn't matter,
;	use 0). Converting this to heliographic coordinates gives the heliographic
;	longitude, rLng, and latitude, rLat of ecliptic North.
;	If [A,B,C] is a unit vector pointing to ecliptic north, then
;		A = cos(rLat)*cos(rLng)
;		B = cos(rLat)*sin(rLng)
;		C = sin(rLat)
;	The equation of the ecliptic plane in heliographic coordinates is given
;	by the A*x+B*y+C*z=0. So for given x,y: z = -(A*x+B*y)/C
; MODIFICATION HISTORY:
;	28-Feb-1997, Paul Hick (UCSD, pphick@ucsd.edu)
;-

rpm = ToRadians(degrees=Degrees)

InitVar, Angle  , 0
InitVar, Corners, /key
InitVar, Grid   , /key

Grid = Grid OR Corners

bOne = (size(X))[0] EQ 1			; One point only

XX = reform(X[0,*])					; X,Y coordinates relative to Sun
YY = reform(X[1,*])

IF Grid THEN BEGIN					; Combine each X which each Y
	nDim = n_elements(XX)
	S = replicate(1.,nDim)
	XX = XX#S
	YY = S#YY
ENDIF

IF Angle NE 0 THEN BEGIN			; Rotate the X,Y vectors over angle A
	cosA = cos(Angle*rpm)
	sinA = sin(Angle*rpm)
	S  = XX*cosA-YY*sinA
	YY = XX*sinA+YY*cosA
	XX = S
ENDIF

					; Ecliptic North Pole in heliographic coordinates
Pole = CvSky(TT, from_ecliptic=[0.0,!pi/2], /to_heliographic, /silent)

ZZ = -(cos(Pole[0])*XX+sin(Pole[0])*YY)/tan(Pole[1])

CASE 1 OF

Corners: BEGIN
	ZZ =[ [XX[0,0],YY[0,0],ZZ[0,0]],	$
		  [XX[1,0],YY[1,0],ZZ[1,0]],	$
		  [XX[1,1],YY[1,1],ZZ[1,1]],	$
		  [XX[0,1],YY[0,1],ZZ[0,1]],	$
		  [XX[0,0],YY[0,0],ZZ[0,0]] ]

	IF IsType(Color,/defined) THEN plots, ZZ, /t3d, color=Color
END

bOne: ZZ = [XX[0],YY[0],ZZ[0]]

ELSE: BEGIN
	IF IsType(Color,/defined) THEN BEGIN
		Xmin = min(XX,max=Xmax)
		Ymin = min(YY,max=Ymax)
		Zmin = min(ZZ,max=Zmax)

		Position = Convert_Coord([Xmin,Xmax],[Ymin,Ymax],[Zmin,Zmax], /t3d, /data, /to_normal)
		Position = [Position[0:1],Position[3:4],Position[2],Position[5]]

		surface, ZZ, XX,YY, $
			color	= Color , $
			/t3d			, $
			xstyle	= 5 	, $
			ystyle	= 5 	, $
			zstyle	= 5 	, $
			;thick	= 0.1	, $
			position= Position
	ENDIF

	ZZ = transpose([[[XX]],[[YY]],[[ZZ]]],[2,0,1])
END

ENDCASE

RETURN, ZZ  &  END