pro mideshift_xy,wavelength,flux,xmin,ymin,sigmax,npoints=npoints,order=order,$
interp=interp,minabs=minabs,noisy=noisy,noplot=noplot,tpb=tpb

;+
; Measuring the wavelength of the center of a line.
;
;
; IN:  wavelength - dblarr - wavelengths (angstroms); absolute values
;      flux 	  - fltarr - fluxes 
;
; OUT: xmin       - double  - wavelength of the line bottom
;      ymin       - double  - flux at xmin
;      sigmax	  - double  - sigma associated with the measurement of xmin
;
;
; KEYWORDS: npoints- number of pixels around the minimum to enter the fit.
;			(default: 7) It has to be an even number.
;	    order  - order of the polynomial (defult: 3 = third order)
;	    interp - when on, we use a spline interpolation to improve sampling
;			(step= 0.005 A)
;	    minabs - min. central absorption of a line in order to be considered
;			(default: 0.98)
;	    noisy  - if set, we measure line shifts of lines with irregular
;			shapes close to the line center. Otherwise, we don't.
;	    noplot - when set, it does not produce any plot
;	    tpb	   - two-point-bisector tecnique (Hamilton & Lester 1999)
;
; NOTE: when the measurement cannot be performed, -1000 is return in xmin,ymin
;	and sigmax 
;
; C. Allende Prieto, UT,  1999
; C. Allende Prieto, UT,  2001 included spline interp. to improve sampling 
; C. Allende Prieto, UT,  2001 more rigorous determination of error in lambda
;				implementation of 'security' checks
; C. Allende Prieto, UT,  2006 fixed a bug in the determination of the error
;				and added a 2nd error estimate, the average of
;				the two is kept as tests indicate it is more robust
;				than either method alone.
;-
npar = n_params()
if (npar eq 0) then begin
	print,'mideshift_xy,x,y,xmin,ymin,lsig'
	return
endif 

step=wavelength(1)-wavelength(0)

keepwavelength=wavelength & keepflux=flux

if keyword_set(interp) then begin
; spline interp. to improve sampling
	;if(step gt 5e-3) then begin
		step=5e-3
		nns=floor((max(wavelength)-min(wavelength))/step)
		nwavelength=findgen(nns)*step+min(wavelength)
	;endif
	;nflux=spline(wavelength,flux,nwavelength)
	nflux=interpol(flux,wavelength,nwavelength,/spline)
	wavelength=nwavelength & flux=nflux
endif

if not keyword_set(npoints) then npoints=7
if not keyword_set(order) then order=3
if not keyword_set(minabs) then minabs=0.98
if not keyword_set(noisy) then begin
	noisy=0
endif else begin
	noisy=100
endelse
nhalf=(npoints-1)/2

data=wavelength(min(where(flux eq min(flux))))
big=size(flux)
pixelmin=fix(min(where(flux eq min(flux))))


; we  request at least nhalf points each side of the  minimum
;		a consistant slope for each group of nhalf points on each side 
;			(unless /noisy is set)
;	 	a significant absorption to exist (minabs of the continuum level)
cool=1
if (pixelmin+nhalf gt big(1)-1 or pixelmin-nhalf lt 0 or $
		min(flux) gt minabs) then begin
	cool=0
endif else begin
	if(max(deriv(flux(pixelmin-nhalf:pixelmin-1))) gt noisy or $
		min(deriv(flux(pixelmin+1:pixelmin+nhalf))) lt -noisy) then cool=0
endelse

if (cool eq 0) then begin			
		print,'% mideshift_xy: measurement failure'
		xmin=-1000
		ymin=-1000
		sigmax=-1000
		sigmay=-1000
		wavelength=keepwavelength & flux=keepflux
		if not keyword_set(noplot) then begin
			plot,wavelength,flux,psym=2,title=data(0),charsize=1.7
		endif	
		return
endif else begin

if  keyword_set(tpb) then begin ; 2pd if

	first_blue=0.5*(wavelength(pixelmin-1)+$
	interpol(wavelength(pixelmin:n_elements(wavelength)-1),$
	flux(pixelmin:n_elements(wavelength)-1),flux(pixelmin-1)))
	first_red=0.5*(wavelength(pixelmin+1)+$
	interpol(wavelength(0:pixelmin),$
	flux(0:pixelmin),flux(pixelmin+1)))
	xmin=mean([first_blue,first_red])
	sigmax=0. ; unknown
	ymin=flux(pixelmin); approximately
	
endif else begin


coef=poly_fit(wavelength(pixelmin-nhalf:pixelmin+nhalf)-data(0),$
flux(pixelmin-nhalf:pixelmin+nhalf),order,p1,p2,sigmay,p4)

; old style for getting the minimum and sigma
;ran=3*step;
;xgrid=dindgen(ran/0.000001)*0.000001+wavelength(pixelmin-1)-data(0)
;ygrid=coef(0)+coef(1)*xgrid+coef(2)*xgrid^2+coef(3)*xgrid^3;+coef(4)*xgrid^4
;xmin=xgrid(where(ygrid eq min(ygrid)))+data(0)
;ymin=min(ygrid(where(ygrid eq min(ygrid))))
;New error estimate; Carlos UT Feb 2001 
;P(x)=sum(i=1,2,3) A(i) X**i
; sigma^2(x)=sigma^2(P) / (sum[i=1,2,3] A^2(i) i^2 x^(2*(i-1)))
;ii=transpose(findgen(n_elements(coef)))
;sigmax=sigmay/sqrt(total((coef*ii*(xmin(0)-data(0))^((ii-1)))^2))

;New estimate of minimum (derivative of the polynomial=0) and error
; Carlos UT Sep 2001
derivativecoef=shift(coef,-1)
derivativecoef=derivativecoef(0:n_elements(coef)-2)*$
	(findgen(n_elements(coef)-1)+1)
	
if order eq 3 then begin
	roots=[(-coef(2)+sqrt(coef(2)^2-3.*coef(3)*coef(1)))/3./coef(3),$
	(-coef(2)-sqrt(coef(2)^2-3.*coef(3)*coef(1)))/3./coef(3)]; I'll doit analitically for order=2
	if (max(where(imaginary(roots) eq 0)) eq -1) then begin
		print,'% mideshift_xy: Complex roots! I quit!'
			print,'% mideshift_xy: measurement failure'
			xmin=-1000
			ymin=-1000
			sigmax=-1000
			sigmay=-1000
			wavelength=keepwavelength & flux=keepflux
			if not keyword_set(noplot) then begin
				plot,wavelength,flux,psym=2,title=data(0),charsize=1.7
			endif
		return
	endif
	if (max(where(finite(roots))) eq -1) then begin
		print,'% mideshift_xy: Not finite roots! I quit!'
			print,'% mideshift_xy: measurement failure'
			xmin=-1000
			ymin=-1000
			sigmax=-1000
			sigmay=-1000
			wavelength=keepwavelength & flux=keepflux
			if not keyword_set(noplot) then begin
				plot,wavelength,flux,psym=2,title=data(0),charsize=1.7
			endif
		return
	endif	
endif else begin
	roots=fz_roots(derivativecoef,/double) ; this one get's them numerically
	if (max(where(imaginary(roots) eq 0)) eq -1) then begin
		print,'% mideshift_xy: Complex roots! I quit!'
			print,'% mideshift_xy: measurement failure'
			xmin=-1000
			ymin=-1000
			sigmax=-1000
			sigmay=-1000
			wavelength=keepwavelength & flux=keepflux
			if not keyword_set(noplot) then begin
				plot,wavelength,flux,psym=2,title=data(0),charsize=1.7
			endif
		return
	endif
	if (max(where(finite(roots))) eq -1) then begin
		print,'% mideshift_xy: Not finite roots! I quit!'
			print,'% mideshift_xy: measurement failure'
			xmin=-1000
			ymin=-1000
			sigmax=-1000
			sigmay=-1000
			wavelength=keepwavelength & flux=keepflux
			if not keyword_set(noplot) then begin
				plot,wavelength,flux,psym=2,title=data(0),charsize=1.7
			endif
		return
	endif	
	roots=double(roots)
endelse

xmin=roots(where(min(abs(roots)) eq abs(roots)))+data(0)
ymin=poly(xmin-data(0),coef)

;now error bars, two ways:
;1) sigma(x) directly from the rms scatter in the x axis     -> sigmax
;2) 1./sigma(x_i)^2 = sum_i (1./sigma(y_i)^2) (@y_i/@x_i)^2  -> sigmax2
; then we keep the mean of the two in sigmax

xfit=dblarr(npoints)  
for i=0,npoints-1 do begin
	obs=flux(pixelmin-nhalf+i)
	if (flux(pixelmin-nhalf+i) lt ymin(0)) then begin
		; if the observation is lower than the minimum of the line
		obs=obs+2.*(ymin-flux(pixelmin-nhalf+i))
	endif  
	coef2=coef
	coef2[0]=coef2[0]-obs[0]
	roots=fz_roots(transpose(coef2),/double)
	;plot,wavelength(pixelmin-nhalf:pixelmin+nhalf)-data(0),$
	;flux(pixelmin-nhalf:pixelmin+nhalf)
	;oplot,wavelength(pixelmin-nhalf:pixelmin+nhalf)-data(0),$
	;	poly(wavelength(pixelmin-nhalf:pixelmin+nhalf)-data(0),coef),$
	;	col=140,thick=2
	;stop
	if (max(where(imaginary(roots) eq 0)) eq -1) then begin
		print,'% mideshift_xy: Complex roots! I quit!'
			print,'% mideshift_xy: measurement failure'
			xmin=-1000
			ymin=-1000
			sigmax=-1000
			sigmay=-1000
			wavelength=keepwavelength & flux=keepflux
			if not keyword_set(noplot) then begin
				plot,wavelength,flux,psym=2,title=data(0),charsize=1.7
			endif
		return
	endif
	roots=double(roots(where(imaginary(roots) eq 0)))
	xfit(i)=roots(where(abs(roots-wavelength(pixelmin-nhalf+i)+data(0)) eq $
	min(abs(roots-wavelength(pixelmin-nhalf+i)+data(0)))))
	;print,roots
	;print,'pick:',xfit(i),'   ref=',wavelength(pixelmin-nhalf+i)-data(0)
endfor

sigmax=stdev(xfit-wavelength(pixelmin-nhalf:pixelmin+nhalf)-data(0))
ydiff=poly(wavelength(pixelmin-nhalf:pixelmin+nhalf)-data(0),coef)-$
	    flux(pixelmin-nhalf:pixelmin+nhalf)
xmin=xmin(0)
ymin=ymin(0)
sigmax=sigmax(0)/sqrt(npoints*1.)
sigmax2=poly(wavelength(pixelmin-nhalf:pixelmin+nhalf)-data(0),derivativecoef)
sigmax2=total(1.d0/ydiff^2*sigmax2^2)
sigmax2=1./sqrt(sigmax2)

;print,sigmax,sigmax2,(sigmax+sigmax2)/2.0
sigmax=(sigmax+sigmax2)/2.0d0

endelse; 2pd if

if not keyword_set(noplot) then begin
	plot,wavelength(pixelmin-nhalf:pixelmin+nhalf)-data(0),$
	flux(pixelmin-nhalf:pixelmin+nhalf),psy=8,$
	yr=[min(flux(pixelmin-nhalf:pixelmin+nhalf))*0.98,$
	max(flux(pixelmin-nhalf:pixelmin+nhalf))*1.02],$
	xr=[min(wavelength(pixelmin-nhalf:pixelmin+nhalf))-data(0)-0.01,$
	max(wavelength(pixelmin-nhalf:pixelmin+nhalf))-data(0)+0.01],$
	xstyle=1,ystyle=1,thick=2,charsize=2.,/nodata,$
	title=string(xmin,format='(f10.4)')
	
	plotsym,0,2,/fill
	oplot,wavelength(pixelmin-nhalf:pixelmin+nhalf)-data(0),$
	flux(pixelmin-nhalf:pixelmin+nhalf),psy=8
	ran=10*step
	xgrid=dindgen(ran/0.0001)*0.001+wavelength(pixelmin-nhalf)-data(0)
	if not keyword_set(tpb) then begin
		oplot,xgrid,poly(xgrid,coef),thick=2
		oplot,[xmin-data(0),xmin-data(0)],[ymin-sigmay,ymin+sigmay]
	endif	
	xx=[xmin-data(0),xmin-data(0)]
	yy=[0,1]
	oplot,xx,yy
	oplot,[xmin-data(0)-sigmax,xmin-data(0)+sigmax],[ymin,ymin],thick=2
endif

xx=[xmin,xmin]
yy=[0,1]

; get the input wave.-flux vectors back
wavelength=keepwavelength & flux=keepflux

endelse

end