function hsi_find_spinaxis,cbe,det,VERB=verb,PLOT=plot
;+
; PURPOSE:
;  Finds the rhessi spin axis using a 4-parameter fit to the SAS
;  cbe.dx,cbe.dy, cbe.roll_angle arrays
;
; METHOD:
;  Finds the 4 parameters a,b,xc,yc that best (least square) fit
;   dx= -a*cos(phi) + b*sin(phi) + xc
;   dy = a*sin(phi) + b*cos(phi) + yc
;   where phi =cbe.roll_angle+orient
;
;  The (xc,yc) values are the coordinates of the "least-squares" spin axis 
;   and sqrt(a^2+b^2) = radius of spin circle

; INPUTS:
;  cbe = pointer to calibrated event list
;  det = scalar value of detector (0-8)
;
; OUTPUTS:
;   [a,b,xc,yc]
;   xc,yc = coordinates of the "least-squares" spin axis 
;   and sqrt(a^2+b^2) = radius of spin circle
;
; VERSION HISTORY:
;  ejs december 2002-jan 2003  schmahl@hessi.gsfc.nasa.gov
;-

det=det mod 9

orient=((hsi_grid_parameters()).orient)[det]
phi=orient+((*cbe[det,0]).roll_angle)
dx=((*cbe[det,0]).dx)
dy=((*cbe[det,0]).dy)


cphi=(cos(phi)) & mcphi=mean(cphi)
sphi=(sin(phi)) & msphi=mean(sphi)

; Use the method of least sqaures
; minimize the objective function
; S=sum((-a*cos(phi)+b*sin(phi)+xc-dx)^2 + (a*sin(phi)+b*cos(phi)+yc-dy)^2)
; with respect to a,b,xc,yc.

c1=-mean(dx*cphi)+mean(dy*sphi)
c2=mean(dx*sphi)+mean(dy*cphi)
c3=mean(dx)
c4=mean(dy)

m=[ [1,     0,    -mcphi,msphi], $
    [0,     1,     msphi,mcphi], $
    [-mcphi,msphi, 1,    0    ], $
    [ msphi,mcphi, 0,    1    ]  ]
abxy=invert(m)#[c1,c2,c3,c4]

return,abxy

end

;MAIN

restore,'cbe+p_2002-03-25-20:05:00.sav',/ver  ; xc,yc=[339.6,  0.05]
;restore,'cbe+p_2002-04-20-23:10:00.sav',/ver  ; xc,yc=[364.5,  -93.]
;restore,'cbe+p_2002-02-20-09:40:00.sav',/ver  ; xc,yc=]247., -95.]
;restore,'cbe+p_2002-06-02-11:44:00.sav',/ver  ; xc,yc=[174,146]


det=8
abxy=hsi_find_spinaxis(cbe,det)
print,'OUTPUT: a,b,xc,yc=',abxy






end