Viewing contents of file '../idllib/astron/contrib/freudenreich/xyfit.pro'
FUNCTION XYFIT,X,Y,W,CCX,CCY,DIST,FX,FY
;+
; NAME:
; XYFIT
; PURPOSE:
; Returns the bisecting line of the regressions: Y on X and X on Y.
; OR the mean of the two lines, with FX and FY the weights.
; CALLED BY FITXYERRS.
; CALLING SEQUENCE:
; YFIT = XYFIT( X, Y, W, CCX, CCY, DIST, FX, FY )
; INPUTS:
; X,Y,W = points and weights
; FX, FY = weight given to X|Y and Y|X regressions in their average
; OUTPUT:
; CCX = coefficients of Y on X regression
; CCY = same for X on Y
; DIST = perpendicular distance of points to line
; AUTHOR:
; H.T. Freudenreich
;-
on_error,2
EPS=1.0E-20
SX=TOTAL(W*X) & SY=TOTAL(W*Y) & SXY=TOTAL(W*X*Y) & SXX=TOTAL(W*X*X)
SYY=TOTAL(W*Y*Y)
USE_X=1
USE_Y=1
IF N_PARAMS() EQ 8 THEN BEGIN
IF FX LT .001 THEN USE_X=0
IF FY LT .001 THEN USE_Y=0
ENDIF
IF USE_Y EQ 1 THEN BEGIN
; The Y vs X line.
D=SXX-SX*SX
IF ABS(D) LT EPS THEN BEGIN
PRINT,'XYFIT: Determ=0. No fit possible.'
RETURN,0.
ENDIF
YSLOP=(SXY-SX*SY)/D & YYINT=(SXX*SY-SX*SXY)/D
SLOP=YSLOP & YINT=YYINT
IF USE_X EQ 0 THEN BEGIN
CCX=0.
CCY=[YINT,SLOP]
DIST=Y-(YINT+SLOP*X)
RETURN,CCY
ENDIF
ENDIF
IF USE_X EQ 1 THEN BEGIN
; Get the X vs Y line.
D=SYY-SY*SY
IF ABS(D) LT EPS THEN BEGIN
PRINT,'XYFIT: Determ=0. No fit possible.'
RETURN,0.
ENDIF
TSLOP=(SXY-SY*SX)/D & TYINT=(SYY*SX-SY*SXY)/D
; Now invert it to get the form Y=a+bX:
IF ABS(TSLOP) LT EPS THEN BEGIN
PRINT,'XYFIT: X vs Y line uninvertable. No fit possible.'
RETURN,0.
ENDIF
XSLOP=1./TSLOP & XYINT=-TYINT/TSLOP
IF USE_Y EQ 0 THEN BEGIN
CCY=0.
CCX=[XYINT,XSLOP]
DIST=X-(TYINT+TSLOP*Y)
RETURN,CCX
ENDIF
ENDIF
IF N_PARAMS() EQ 8 THEN BEGIN
; Calculate a weighted mean line:
YINT=FX*XYINT+FY*YYINT
SLOP=FX*XSLOP+FY*YSLOP
ENDIF ELSE BEGIN
; Calculate the equation of the bisector of the 2 lines:
IF YSLOP GT XSLOP THEN BEGIN
A1=YYINT & B1=YSLOP & R1=SQRT(1.+YSLOP^2)
A2=XYINT & B2=XSLOP & R2=SQRT(1.+XSLOP^2)
ENDIF ELSE BEGIN
A2=YYINT & B2=YSLOP & R2=SQRT(1.+YSLOP^2)
A1=XYINT & B1=XSLOP & R1=SQRT(1.+XSLOP^2)
ENDELSE
YINT=(R1*A2+R2*A1)/(R1+R2)
SLOP=(R1*B2+R2*B1)/(R1+R2)
ENDELSE
R=SQRT(1.+SLOP^2) & IF YINT GT 0. THEN R=-R
U1=SLOP/R & U2=-1./R & U3=YINT/R
DIST=U1*X+U2*Y+U3 ; = orthog. distance to line
CC=[YINT,SLOP]
CCX=[XYINT,XSLOP] & CCY=[YINT,YSLOP]
RETURN,CC
END