pro krEllipsoidFit, x,y,z, weighting=weighting, maxK=maxK, evec=evec, evals=evals, $ kArray=kArray, Xm=Xm, P=P ;+ ; :Description: ; finds the best fit ellipse of ellipsoid of the input data ; ; :Params: ; x,y,[z] = vertices of the data to be fitted ; ; :Keywords: ; maxK = The max K value that will form an ellipse/ellipsoid that will encompass all ; the data ; kArray = array of all the K values as calculated by the program ; evec = eigen vector of the ellipse/ellipsoid ; evals = eigen vales of the ellipse/ellipsoid ; P = second weighted moment ; weighting = optional weighting for each vertices. i.e. the intensity of the ; pixel at each point ; Xm = center of mass (if weighting is all 1). Otherwise this is the first weighted moment ; ; :Author: Ronn Kling, based upon an idea by Jerry Lefever ; ;- numDims = n_params() if numDims le 1 then message,'You must have at least two dimensions' numPoints = n_elements(x) ;calculate the first moment if numDims eq 2 then begin xVec = transpose([[x],[y]]) endif else begin xVec = transpose([[x],[y],[z]]) endelse if n_elements(weighting) ne 0 then begin if numDims eq 2 then begin W = transpose([[weighting],[weighting]]) endif else begin W = transpose([[weighting],[weighting],[weighting]]) endelse Xm = total(xVec*W,2)/total(weighting) endif else begin weighting = fltarr(numPoints)+1.0 Xm = total(xVec,2)/total(weighting) endelse ;calculate the second moment P dXVec = xVec - (Xm # replicate(1.0,numPoints)) P = 0.0 for i=0L,numPoints-1 do P = P + weighting[i]*(dXVec[*,i]#dXVec[*,i]) P = P/total(weighting) ;find the max K kArray = fltarr(numPoints,/nozero) invP = invert(P) for i=0L,numPoints-1 do kArray[i] = (xVec[*,i] - Xm)#invP#(xVec[*,i] - Xm) maxK = sqrt(max(kArray)) evals = eigenql(P, eigenvectors = evec, residual = residual) return & end