' subroutine to plot an ellipse represented as a quadratic form
'	(x-x0)'B(x-x0) = c
'
' where the inputs are
'
' 	  x0: 2x1 vector
' 	   B: 2x2 symmetric matrix
' 	   c: scalar
' %gname: name for graph object
'
' 1/22/99 h.
' revised 12/19/2000 for version 4 (h)

subroutine drawEllipse(vector x0, sym B, scalar c, string %gname)

	' input error check
	!rows = @rows(B)
	!cols = @columns(B)
	if (@rows(x0)<>2 or !rows<>2 or !cols<>2) then
		statusline input arguments must have dimension 2
		return
	endif

	' get eigenvalues/eigenvectors
	matrix ve = @eigenvectors(B)
	vector va = @eigenvalues(B)

	' construct the rotation matrix
	!ang = @atan(ve(2,1)/ve(1,1))			' angle
	matrix(2,2) rot
	rot.fill(b=r) @cos(!ang),-@sin(!ang),@sin(!ang),@cos(!ang)

	' draw the unrotated ellipse
	vector(101) theta	
	!pi = @acos(-1)
	for !i = 1 to 101
		theta(!i) = (!i-1)*2*!pi/100
	next

	!a = @sqrt(c/va(1))			' width of ellipse
	!b = @sqrt(c/va(2))
	matrix(2,101) xy
	for !i=1 to 101
		xy(1,!i) = !a*@cos(theta(!i))
		xy(2,!i) = !b*@sin(theta(!i))
	next

	' rotate 
	matrix rx = rot*xy			
	' offset
	matrix ones = @filledmatrix(1,101,1)
	rx = rx + x0*ones
	rx = @transpose(rx)

	' plot
	freeze({%gname}) rx.scat
	' connect scatter without symbols
	{%gname}.setelem(1) lpat(solid) symbol(none) 
endsub

 