pro testEllipseFit device, decomposed=0 loadct, 39 ;create a guassian ellipse that we know is correct nx = 1000L & ny = 1000L ; Define input function parameters: A = [0.0, 10., nx/6., ny/10., nx/2., .6*ny] ; Create X and Y arrays: X = FINDGEN(nx) # REPLICATE(1.0, ny) Y = REPLICATE(1.0, nx) # FINDGEN(ny) ; Create an ellipse: U = ((X-A[4])/A[2])^2 + ((Y-A[5])/A[3])^2 ; Create gaussian Z: Z = A[0] + A[1] * EXP(-U/2) ; Fit the function, no rotation: yfit = GAUSS2DFIT(Z, B) ; Report results: PRINT, 'Should be: ', STRING(A, FORMAT='(6f10.4)') PRINT, 'Is: ', STRING(B(0:5), FORMAT='(6f10.4)') xCenter = nx/2 yCenter = .6*ny ;we need 1D vectors x = reform(x,nx*ny) y = reform(y,nx*ny) ;weight by the intensity weighting = reform(z,nx*ny) ;color by the weighting color = bytscl(weighting) plot,x,y,psym=1,/isotropic,/nodata for i = 0L,n_elements(x)-1 do begin plots,x[i],y[i],color=color[i],psym=3 ;plot the points by color endfor krEllipsoidFit, x,y, weighting=weighting, maxK=maxK, evec=evec, evals=evals, $ kArray=kArray, Xm=Xm, P=P ;k = sqrt(sigma(kArray)) ;k = sqrt(max(kArray)) k = 1 ; will show the 1 sigma ellipse print,'k = ',k semiMajor = sqrt(evals[0]) semiMinor = sqrt(evals[1]) angle = atan(evec[1,0],evec[0,0]) print,'angle 1 = ',angle/!dtor majorVecX = [xCenter, xCenter + semiMajor*cos(angle)*k] majorVecY = [yCenter, yCenter + semiMajor*sin(angle)*k] ;oplot, majorVecX, majorVecY,color=255 angle = atan(evec[1,1],evec[0,1]) print,'angle 2 = ',angle/!dtor minorVecX = [xCenter, xCenter + semiMinor*cos(angle)*k] minorVecY = [yCenter, yCenter + semiMinor*sin(angle)*k] ;oplot, minorVecX, minorVecY,color=255 count = 0 for theta=0.0,2*!pi,0.1 do count++ x = fltarr(count) y = x i=0 ;plot the fitted ellipse points for theta=0.0,2*!pi,0.1 do begin u = ([cos(theta),sin(theta)]) denom = sqrt((u # invert(p)#u)) a = k/denom x[i] = a*cos(theta) + xCenter y[i] = a*sin(theta) + yCenter plots,x[i],y[i],psym=1,color=255 i++ endfor ;overplot the 1 sigma axis from the gaussian fit plots,B[4],B[5],psym=1,symsize=3 oplot,[B[4],B[4]+B[2]],[B[5],B[5]],thick=3 oplot,[B[4],B[4]],[B[5],B[5]+B[3]],thick=3 return & end