pro robust_reg,xx,yy,a,b,yfit,Log,miss=miss

; Log is the lun for the log file; -1 means no file is connected and no data is logged
; performs robust regression between x and y according to
; Emerson, J. D., and D. C. Hoaglin. 1983.  Resistant lines 
;  for y versus x.  pp. 129-65. In D. C. Hoaglin, F. Mosteller,
;  and J. W. Tukey, eds., Understanding Robust and               
;  Exploratory Data Analysis.  John Wiley & Sons.
;
if n_elements(miss) eq 0 then miss=-9999
nan=!values.f_nan
x=float(xx) & y=float(yy)

; Tim Mitchell routine added to turn all Miss into NaN
allmiss = where(finite(x) eq 0 or finite(y) eq 0 or x eq miss or y eq miss)
x(allmiss)=nan
y(allmiss)=nan

a=nan & b=nan
yfit=replicate(nan,n_elements(xx))
; Mark New: m=where(finite(x) eq 0 and finite(y) eq 0,mm)
m=where(finite(xx) eq 0 and finite(yy) eq 0,mm)
if mm gt 0 then x(m)=miss
if mm gt 0 then y(m)=miss
; Mark New: n=where(x ne miss and y ne miss,nn)
n=where(xx ne miss and yy ne miss,nn)
if nn lt 9 then return

if (Log GE 0) then printf, Log, mm, nn, format='(2i6)'

srt=sort(x) & nsrt=n_elements(srt) & n3=nsrt/3
x=x(srt) & y=y(srt)

;first iteration
xl=median(x(0:n3-1)) & xm=median(x(n3:2*n3-1))
xr=median(x(n3*2:nsrt-1))
yl=median(y(0:n3-1)) & ym=median(y(n3:2*n3-1))
yr=median(y(n3*2:nsrt-1))
b0=(yr-yl)/(xr-xl)
a0=median(y-b0*(x-xm))
yf=b0*(x-xm)
r0=y-yf
drb0=median(r0(n3*2:nsrt-1))-median(r0(0:n3-1))

if (Log GE 0) then begin
	printf, Log, xl, xm, xr, format='(3e12.4)'
	printf, Log, yl, ym, yr, format='(3e12.4)'
	printf, Log, a0, b0, drb0, format='(3e12.4)'
endif

; calculate b1 and drb1 iteratively
for i=0,5 do begin
 db1=drb0/(xr-xl)
 if i eq 0 then b1=b0+db1
 r1=y-(b1)*(x-xm)
 drb1=median(r1(n3*2:nsrt-1))-median(r1(0:n3-1))

 if (Log GE 0) then begin
	printf, Log, i, format='(i4)'
	printf, Log, db1, b1, drb1, format='(3e12.4)'
 endif

 if drb1 gt -0.001 and drb1 lt 0.001 then begin
  a=median(y-b1*(x))
  b=b1
  yfit(n)=a+b*xx(n)

  if (Log GE 0) then printf, Log, a, b, format='(2e12.4)'
  return
 endif
 ;
 ;interpolate between first and second estimates
 b2=b1-drb1*((b1-b0)/(drb1-drb0))
 b0=b1
 b1=b2
 drb0=drb1

 if (Log GE 0) then printf, Log, b0, b1, b2, drb0, format='(4e12.4)'
endfor

end