function ZEclipticPlane, TT, X,		$
		sun=Sun, corners=Corners, grid=Grid, angle=Angle, degrees=Degrees, $
		color = color
;+
; NAME:
;	ZEclipticPlane
; PURPOSE:
;	Get coordinate arrays describing ecliptic plane in a
;	Cartesian heliographic coordinate system.
; CATEGORY:
;	Sky graphics
; CALLING SEQUENCE:
;	ZEcliptic = ZEclipticPlane (TT, X, sun=Sun, /corners, angle=Angle, /degrees] )
; INPUTS:
;	TT		scalar; 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[2,0:1]
;						are used. The output arrays will be array[5],
;						in the sequence (0,0),(1,0),(1,1),(0,1),(0,0),
;						i.e. the X,Y,Z arrays 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
;	sun=Sun		scalar, array[2], array[3]; default: [0,0,0])
;						x,y,z coordinates of Sun
;						If Sun is scalar then x=y=z is assumed.
;						If Sun is array[2], z=0 is assumed
;	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[*,*]
;	Y				array[*], array[*,*]
;									X,Y,Z coordinates of points in ecliptic plane
; OPTIONAL OUTPUT PARAMETERS:
; CALLS:
;	CvSky
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
;	28-Feb-1997, Paul Hick (UCSD, pphick@ucsd.edu)
;-
on_error, 2

rpd = ToRadians(degrees=Degrees)

if n_elements(Sun  ) eq 0 then Sun   = [0,0,0]
if n_elements(Angle) eq 0 then Angle = 0

Corners = keyword_set(Corners)
Grid    = keyword_set(Grid) or Corners

bOne = (size(X))[0] eq 1

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

if Grid then begin
	nDim = n_elements(XX)
	S = replicate(1.,nDim)		; Combine each X which each Y
	XX = XX#S
	YY = S#YY
endif

if Angle ne 0 then begin		; Rotate the X,Y vectors over angle A
	cosA = cos(A*rpd)
	sinA = sin(A*rpd)
	S  = XX*cosA-YY*sinA
	YY = XX*sinA+YY*cosA
	XX = S
endif

; 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 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

					; Ecliptic North Pole in heliographic coordinates
Pole = CvSky(TT, fromecliptic=[0.,!pi/2], /toheliographic)

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

if Corners then 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 n_elements(Color) ne 0 then plots, ZZ, /t3d, color=Color

endif else if bOne then begin
	ZZ = [XX[0],YY[0],ZZ[0]]

endif else begin
	Xmin = min(X,max=Xmax)
	Ymin = min(Y,max=Ymax)
	Zmin = min(Z,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]]
	if n_elements(Color) ne 0 then $
		 surface, ZZ, XX,YY, color=Color, /t3d,	$
			xstyle=5,ystyle=5,zstyle=5,			$
;			thick=0.1,	$
			position=Position
	ZZ = transpose ( [ [[XX]], [[YY]], [[ZZ]] ], [2,0,1] )

endelse

return, ZZ  &  end