pro testEllipsoidFit

device, decomposed=1

;create some random vertices
vertices = randomn(seed,3,2000,/normal)

;define a center
xCenter = 20
yCenter = 10
zCenter = 15
x = reform(vertices[0,*]) + xCenter
y = reform(vertices[1,*]) + yCenter
z = reform(vertices[2,*]) + zCenter

numPoints = n_elements(x)

krEllipsoidFit, x,y,z, weighting=weighting, maxK=maxK, evec=evec, evals=evals, $
                            kArray=kArray, Xm=Xm, P=P

;position should match the center points from above
print, 'mean position = ',Xm
print,'Second moment = ', P
print, 'evec = ',evec
print, 'evals = ',evals

oPoly1 = obj_new('IDLgrPolygon',x,y,z,color=[255,255,255],style=0)

;k = maxK
;curious to see the min and am
print, max(kArray), MIN(KARRAY)
k = 1 ;one sigma

;put the axes on
oSmall = obj_new('IDLgrPolyline', $
  [Xm[0],Xm[0]+evec[0,0]*k*evals[0]],[Xm[1],Xm[1]+evec[1,0]*k*evals[0]],[Xm[2],Xm[2]+evec[2,0]*k*evals[0]],$
  color=[255,0,0])

oMiddle = obj_new('IDLgrPolyline', $
 [Xm[0],Xm[0]+evec[0,1]*k*evals[1]],[Xm[1],Xm[1]+evec[1,1]*k*evals[1]],[Xm[2],Xm[2]+evec[2,1]*k*evals[1]],$
 color=[0,0255,0])

oBig = obj_new('IDLgrPolyline', $
 [Xm[0],Xm[0]+evec[0,2]*k*evals[2]],[Xm[1],Xm[1]+evec[1,2]*k*evals[2]],[Xm[2],Xm[2]+evec[2,2]*k*evals[2]],$
 color=[0,0,255])

background = [0,0,0]
    
delta = 5*!dtor
count = 0L
x = fltarr(5329)
y = x
z = y

;create the ellipsoid
for theta=0.0,2*!pi,delta do begin
  for phi=0.0,2*!pi,delta do begin
    ;formula from the readme.pdf
    u = ([cos(theta)*cos(phi),sin(theta)*cos(phi), sin(phi)])
    denom = sqrt((u # invert(p)#u))
    a = k/denom
    x[count] = a*cos(theta)*cos(phi) + Xm[0]
    y[count] = a*sin(theta)*cos(phi) + Xm[1]
    z[count] = a*sin(phi) + Xm[2]
    count++
  endfor
endfor   

oPoly = obj_new('IDLgrPolygon',x,y,z,color=[255,255,0],style=0)

xobjview, [oBig, oMiddle, oSmall, oPoly1, oPoly], background=background

return & end