;********************************************************************
;NAME       : linecenter (procedure)
;FUNCTION   : get minimum of absorption line with pol fit method
;						K.I. '91/08/26
;====================================================================

;-----------------------------------------------------------------------
function nearestindex,array,value	
; return i at which the value of array(i) is the nearest to the 'value'
diff=abs(array-value)
index=where(diff eq min(diff))
return,index(0)
end 
;-----------------------------------------------------------------------

pro linecenter,a,xy,width,xc,yc,zc
; a     : spectrum image array
; xy	: direction of dispersion
; width : width of absorption line
; xc,yc	: position of line center
; z  	: intensity at line center
kk=4

; -------------   get absorption line centers with mouse   ------------
on_error,2              ;Return to caller if an error occurs
s = size(a)
if s(0) ne 2 then message, 'Image parameter not 2d.'
nx = s(1)-1		;To n-1
ny = s(2)-1
!err=0
!order=0
device,cursor_standard=132	; cursor   to arrow
print,'click line position'
tvrdc,x,y,3,/dev
;	if (!err and 3) ne 0 then begin	;New line?
;	   print,form="($,a)",string("12b)
;	   while (!err ne 0) do begin wait,.1 & tvrdc,x,y,0,/dev & end
;	  endif
if (x ge nx) or (y ge ny)  or (x lt 0) or (y lt 0) then goto,labe
if ( !err eq 4 ) then begin
	xc=-1
	return
endif

if( xy eq 'x' ) then begin
	z=a(x-width/2:x+width/2,y)
	xmin=where(z eq min(z))
	x=x+xmin(0)-width/2
	z=a(x-width/2:x+width/2,y)
	yc=y
endif else begin
	z=transpose(a(x,y-width/2:y+width/2))
	ymin=where(z eq min(z))
	y=y+ymin(0)-width/2
	z=transpose(a(x,y-width/2:y+width/2))
	xc=x
endelse
ii=indgen(width/2+width/2+1)

coeff=transpose(poly_fit(ii,z,kk))
aa=fltarr(kk)
for i=0,kk-1 do aa(i)=coeff(i+1)*(i+1)
zroots,aa,root
reindex=where(imaginary(root) eq 0)	; index where roots are real
reroot=float(root(reindex))
value=where(z eq min(z))
if( xy eq 'x' ) then begin
	xc0=reroot(nearestindex(reroot,value(0)))
		; xc0 is counted from left side of the fitting range
	xc=float(x-width/2)+xc0
		; xc is counted from left side of the image
	zc=total(coeff*xc0^indgen(n_elements(coeff)))
endif else begin
	yc=reroot(nearestindex(reroot,value(0)))
		; yc0 is counted from bottom of the fitting range
	yc=float(y-width/2)+yc0
		; yc is counted from bottom of the image
	zc=total(coeff*yc0^indgen(n_elements(coeff)))
endelse
xyouts,xc,yc,'!I!9n',size=1,alignment=0.5,color=255,/dev

labe:
device,cursor_standard=34	; cursor to crosshair

end