;Modified from P.P.Hick -- Hsiu-Shan Yu/2013/06/12
;The equation for the fit is A*(x-Tx)+B*(y-Ty)=0
;==================================================================
FUNCTION LSQNormalFit_hsyu,Xin,Yin, $
         PLOTX=PLOTX,OPLOTX=OPLOTX,WORST=WORST,XFIT=XFIT,YFIT=YFIT, $
         FCOLOR=FCOLOR,THICK=THICK,IPSY=IPSY,SCOLOR=SCOLOR
;==================================================================
IF KEYWORD_SET(IPSY)   EQ 0 THEN IPSY=2
IF KEYWORD_SET(SCOLOR) EQ 0 THEN SCOLOR=0
IF KEYWORD_SET(FCOLOR) EQ 0 THEN FCOLOR=0
IF KEYWORD_SET(THICK)  EQ 0 THEN THICK=5
;--------------------
N=N_ELEMENTS(XIN)
IF N_ELEMENTS(YIN) NE N THEN MESSAGE,'Both x and y arrays must have the same size'
;--------------------
XARR=XIN
YARR=YIN
;---
AVGX=TOTAL(XARR)/N
AVGY=TOTAL(YARR)/N
;---
XX=XARR-AVGX
YY=YARR-AVGY
DISXX=SQRT(TOTAL(XX*XX)/N)
DISYY=SQRT(TOTAL(YY*YY)/N)
;--------------------
XARR=XARR/DISXX
YARR=YARR/DISYY
;---
TX=TOTAL(XARR)/N
TY=TOTAL(YARR)/N
;---
XX=XARR-TX
YY=YARR-TY
;---
RX=TOTAL(XX*XX)/N
RY=TOTAL(YY*YY)/N
RXY=TOTAL(XX*YY)/N
;---
RP=(RX+RY)/2.
RM=(RX-RY)/2.
ZXY=(RM*RM+RXY*RXY) > 0
ZXY=SQRT(ZXY)
;--------------------
A1=-TX*RXY+TY*(RM+ZXY)          ; Worst fit
B1= TY*RXY+TX*(RM-ZXY)
F1= SQRT(RP+ZXY)                ; Root of average of dist^2
;---
A2=-TX*RXY+TY*(RM-ZXY)          ; Best fit
B2= TY*RXY+TX*(RM+ZXY)
F2= SQRT((RP-ZXY) > 0)          ; Root of average of dist^2
;---
D    = SORT(XARR)               ; Do not sort the arrays prior to the calculation.
XARR = XARR[D]
YARR = YARR[D]                  ; This destroys the symmetry between x and y.
;------------------------------------------------------------------
IF KEYWORD_SET(PLOTX) THEN BEGIN
  OPLOT,XIN[D],YIN[D],PSYM=IPSY,COLOR=SCOLOR
ENDIF
;---
IF KEYWORD_SET(PLOTX) OR KEYWORD_SET(OPLOTX) THEN BEGIN
  D = 3.*SQRT(RP)
  S = SQRT(A2*A2+B2*B2)
  DX= D*(B2/S)
  DY=-D*(A2/S)
  OPLOT,DISXX*(TX+DX*[-1,1]),DISYY*(TY+DY*[-1,1]),COLOR=FCOLOR,THICK=THICK
  ;---
  IF KEYWORD_SET(worst) AND (A1 NE 0 OR B1 NE 0) THEN BEGIN
    s  = sqrt(A1*A1+B1*B1)
    dx =  d*(b1/s)
    dy = -d*(a1/s)
    oplot, Dxx*(Tx+dx*[-1,1]), Dyy*(Ty+dy*[-1,1])
  ENDIF
  CATCH,IERR
  IF IERR EQ 0 THEN S=POLY_FIT(DISXX*XARR,DISYY*YARR,1,D)
  CATCH,/CANCEL
ENDIF
;------------------------------------------------------------------
YFIT= AVGY-(A2/B2)*(DISYY/DISXX)*(XARR-AVGX)
XFIT= AVGX-(B2/A2)*(DISXX/DISYY)*(YARR-AVGY)
;---
RETURN, [[A2/DISXX,AVGX,B2/DISYY,AVGY,F2],[A1/DISXX,AVGX,B1/DISYY,AVGY,F1]]  &  END