PRO line_normal_form,px,py,rho,theta
;+
;line_normal_form,px,py,rho,theta returns the normal
;rho,theta representation of the regression line
;through the points represented by px and py.
;
;rho is in the same units as px and py and is
;always non-negative.
;
;theta is in radians and ranges from (-pi,pi)
;
;H. Rhody
;July, 2004
;-

n=n_elements(px)<n_elements(py)
IF n LE 1 THEN MESSAGE,'Insufficient number of points'

qx=float(px[0:n-1])
qy=float(py[0:n-1])
x0=min(qx)
x1=max(qx)
y0=min(qy)
y1=max(qy)

CASE Abs(x0-x1) LE abs(y0-y1) OF
0: BEGIN ;Line more horizontal
 Sx=total(qx)
 Sy=total(qy)
 Sxy=Total(qx*qy)
 Sxx=Total(qx^2)
 den=n*Sxx-Sx^2
 b=(Sy*Sxx-Sx*Sxy)/den
 m=(n*Sxy-Sx*Sy)/den
 IF b GE 0 THEN BEGIN
 	rho=b/sqrt(1+m^2)
 	theta=Atan(1,-m)
 ENDIF ELSE BEGIN
 	rho=-b/sqrt(1+m^2)
 	theta=Atan(-1,m)
 ENDELSE
 END ;Horizontal case

1: BEGIN ;Line more vertical
Sx=Total(qx)
Sy=Total(qy)
Sxy=Total(qx*qy)
Syy=Total(qy^2)
den=n*Syy-Sy^2
m=(n*Sxy-Sx*Sy)/den
b=(Sx*Syy-Sy*Sxy)/den
IF b GE 0 THEN BEGIN
	rho=b/sqrt(1+m^2)
	theta=atan(-m,1)
ENDIF ELSE BEGIN
	rho=-b/sqrt(1+m^2)
	theta=atan(m,-1)
ENDELSE
END
ENDCASE

END