FUNCTION emission_lines, wave, flux, fluxivar, plot=plot
;gal='NGC2410_restframe.txt' & readcol, gal, F='(F,F,F)', wave, flux,fluxivar
; USING fitcont.pro TO DETERMINE CONTINUUM
; To see how well emission is subtracted do
; out=emission_lines(wave,flux,fluxivar)
; plot,wave,flux,background='FFFFFF'x,col=255
; oplot,wave,out,col='000000'x

flux_out = flux
n = n_elements(flux)
; List of gas emission lines in vacuum to be fitted by Gaussian and subtracted
 lines = [3203.15 ,3345.81 ,3425.81 ,3726.03 ,3728.73 ,3868.69 ,3889.05 ,3967.40 ,3970.07 ,4101.73 , $ 
; Hg      [OIII]    HeII    [ArIV]   [ArIV]    Hb     [OIII]    [OIII]
;[NI]     [NI]    sky    HeI    NaI    NaI       
 4340.46, 4363.15 ,4685.74 ,4711.30 ,4740.10 ,4861.32 ,4958.83 ,5006.77 ,5197.90 ,5200.39 , 5577.0,  5875.60, 5890.0, 5896.0, $ 
; [OI]     [OI]     [NII]     Ha      [NII]    [SII]    [SII] [ArIII]  
 6300.20, 6363.67 ,6547.96 ,6562.80 ,6583.34 ,6716.31 ,6730.68, 7135.67, $ 
; from Appenzeller and Östreicher 1988 (AJ 95, 45)
; [OII]    [OII]  [FeXI]
  7319.9, 7330.2, 7891.9, $ 
; from Osterbrock (Table 3.12)
; [SIII] [SIII] [SIII]
 8829.9, 9068.9, 9531.0]

ilines = where(lines GE wave[0] AND lines LE wave[n-1], lcount)
lines = lines[ilines]
filter = 300
eps = 3e-3
; gas emission lines inside spectrum
for i=0,lcount-1 do begin
   ctrl = 1
   center = (where(wave GE lines[i]))[0]
;   print, 'center', wave[center]
   if (center LT filter/2) then cut = center-1 else $ 
   if (center GT n-filter/2) then cut = n-1-center $ 
   else cut = filter/2
;   llim = (center-filter GT 0)   ? center-filter : 0 
;   hlim = (center+filter LT n-1) ? center+filter : n-1 
   x = wave[center-cut:center+cut]
   y = flux[center-cut:center+cut]
   var = 1./fluxivar[center-cut:center+cut]
   fitcont, y, continuum=continuum, filter=cut, rms=rms
if (i EQ 0) then tmplines=[x[cut],lines[1:*]] else if (i EQ lcount-1) then $ 
                 tmplines=[lines[0:lcount-2],x[cut]] else $ 
                 tmplines=[lines[0:i-1],x[cut],lines[i+1:*]]
   ;print, '3 x rms:',3*sqrt(total((continuum-y)^2)/(2*cut+1))
   ;lower =  where([alog(y[0:cut])/alog(y[1:cut])]-[(((x[0:cut]-x[cut])*1.d)/((x[1:cut]-x[cut])*1.d))^2] LE eps AND (y[0:cut] GT continuum[0:cut]),lowcount)
   ;lower =  where([alog(y[0:cut]/y[1:cut])] LT 0 AND [alog(y[1:cut]/y[2:cut])] LT 0 AND abs([alog(y[0:cut]/y[1:cut])/alog(y[1:cut]/y[2:cut])]-[(((x[0:cut]-x[cut])*1.d)^2 - ((x[1:cut]-x[cut])*1.d)^2)/(((x[1:cut]-x[cut])*1.d)^2 - ((x[2:cut]-x[cut])*1.d)^2)]) LE eps,lowcount)
lower = where(y[1:cut]-y[0:cut] GT 0 OR (y[0:cut] - continuum[0:cut] GT 0),lowcount)
higher =  where([y[cut:*]-y[cut+1:*]] GT 0 OR (y[cut:*] GT continuum[cut:*]),highcount)
ltest = (where(y[0:cut] LE continuum[0:cut],ccc))[ccc-1]
htest = (where(y[cut:*] LE continuum[cut:*]))[0]
;print, '============================================================'
;print, 'left',where(tmplines GT x[ltest] AND tmplines LT x[cut],nlines), x[0], x[ltest]
;print, 'right',where(tmplines GT x[cut] AND tmplines LT x[cut+htest],nlines),x[cut],x[cut+htest]
lprevious = -1 & found = 0 & hprevious = -1 & ll = -1 & hh = -1 
isThereLineL = where(tmplines GT x[ltest] AND tmplines LT x[cut],llines)
isThereLineR = where(tmplines GT x[cut] AND tmplines LT x[cut+htest],rlines)
if (llines NE 0) then leftlim=x[(where(x GE (tmplines[isThereLineL])[llines-1]))[0]]
if (rlines NE 0) then rightlim=x[(where(x GE (tmplines[isThereLineR])[0]))[0]]
for tt=0,cut-1 do begin
; if there is not another line until this one reaches the continuum 
; then expand until the line crosses continuum 
if (llines EQ 0) then $ ; expand until continuum 
ltmp =  (y[tt] - continuum[tt] GT 0) ? tt : -1 $
else $ ; stop since there is a line and the two are blended
ltmp =  ((y[1+tt] GT y[tt]) AND (y[tt] GT continuum[tt]) AND (x[tt] GE leftlim)) ? tt : -1
if (lprevious EQ -1 AND ltmp NE -1) then ll = ltmp
;print, tt, y[1+tt] - y[tt], y[tt]-continuum[tt],lprevious, ltmp, ll, nlines
lprevious = ltmp
if (rlines EQ 0) then $ ; expand until continuum 
htmp = (found EQ 0 AND (y[cut+tt] - continuum[cut+tt] GT 0)) ? tt : -1 $
else $ ; stop since there is a line and the two are blended
htmp =  ((found EQ 0) AND (y[cut+tt] - y[cut+1+tt] GT 0) AND (y[cut+tt] - continuum[cut+tt] GT 0) AND (x[cut+tt] LE rightlim)) ? tt : -1
if (htmp EQ -1 AND found EQ 0) then begin 
   found = 1 
   if (tt EQ 0) then hh = htmp else hh = hprevious
endif
;print, tt, y[cut+tt] - y[cut+1+tt], y[cut+tt] - continuum[cut+tt],hprevious, htmp, hh, found, nlines
hprevious = htmp
endfor
if (hh EQ 0) then hh = 1 ; if the line is strongly blended
;print, '============================================================'
;print,'ll',ll
;print,'hh',hh
   ;higher =  where([alog(y[cut:*])/alog(y[cut+1:*])]-[(((x[cut:*]-x[cut])*1.d)/((x[cut+1:*]-x[cut])*1.d))^2] GE eps AND (y[cut:*] GT continuum[cut:*]),highcount)
;higher =  where([alog(y[cut:*]/y[cut+1:*])] GT 0 AND [alog(y[cut+1:*]/y[cut+2:*])] GT 0 AND abs([alog(y[cut:*]/y[cut+1:*])/alog(y[cut+1:*]/y[cut+2:*])]-[(((x[cut:*]-x[cut])*1.d)^2 - ((x[cut+1:*]-x[cut])*1.d)^2)/(((x[cut+1:*]-x[cut])*1.d)^2 - ((x[cut+2:*]-x[cut])*1.d)^2)]) LE eps,highcount)


;if (lower[0] NE -1 AND lowcount GT 1) then begin
; test = where(lower[1:*]-lower NE 1,cc)
; ll = (test[0] NE -1) ? (lower[test+1])[cc-1] : lower[0]
;print,'lower[test',lower[test+1]
;endif else ctrl=0
if (lower[0] EQ -1 OR lowcount LE 1) then ctrl = 0
;print,'ll',ll
;if (higher[0] NE -1 AND highcount GT 1) then begin
; test = where(higher[1:*]-higher NE 1)
; hh = (test[0] NE -1) ? (higher[test-1])[0] : higher[n_elements(higher)-1]
;endif else ctrl=0
if (higher[0] EQ -1 AND highcount LE 1) then ctrl = 0
;print,'higher',higher
;print,'ll',ll
;print,'hh',cut+hh
;  ----- CTRL NE 0 ------
; if there are enough points on both sides and the line is in emission
if (ctrl NE 0 AND y[cut] GT continuum[cut]) then begin 
   if (keyword_set(plot)) then begin
   print,i,' ---',lines[i];,'     --',scount,' lines ----------------------' 
   plot,x,y,background='FFFFFF'x,col='000000'x,charsize=1.5, psym=1,xrange=[x[cut-30],x[cut+30]]
   oplot,x,continuum,col='00FF00'x, thick=2
   oplot, [[x[cut]],[x[cut]]],[[0],[y[cut]]],thick=2, line=3, col='000000'x
   oplot, x[lower], y[lower], col='FF0000'x, psym=2 
   oplot, [x[ll],x[ll]], [0,y[ll]], col='FF0000'x, psym=-2, thick=2
   oplot, x[cut+higher], y[cut+higher], col='0000FF'x, psym=2
   oplot, [x[cut+hh],x[cut+hh]], [0,y[cut+hh]], col='0000FF'x, psym=-2, thick=2
   for kk=0, lcount-1 do $
   oplot,[lines[kk],lines[kk]],[0,200], col='000000'x, line=4
endif

; Fit emission lines with Gaussians and subtract best-fit from spectrum
xx = x[ll:hh+cut]
yy = y[ll:hh+cut]
err2 = var[ll:hh+cut]
nn = n_elements(xx)
nterms=4 & A=-1
tmpx = xx & tmpy = yy
quasiRMS = continuum[cut]+3*median(y/continuum)
; Conditions:
if (n_elements(tmpy) GT nterms AND tmpx[n_elements(tmpx)-1] NE wave[n-1]) AND $    
; 1) On the left and right from the center, flux should drop
;(y[cut]-y[cut-1] GE 0 AND y[cut] - y[cut+1] GE 0) AND $
; 2) Center flux should be significantlly above the local continuum
(y[cut] GT quasiRMS) AND $
; 3) Center line must be inside ll AND hh+cut
(x[cut] GT x[ll] AND x[cut] LT x[cut+hh]) then begin
; ------------  fit the line ----------------
expr = 'P[0] + GAUSS1(X, P[1:3])'
if (i EQ 0) then begin  
mmin = min(tmpx) & mmax = (tmpx[where(tmpy EQ max(tmpy))]-mmin) + mmin
endif else begin 
 mmax = max(tmpx) & mmin = tmpx[where(tmpy EQ max(tmpy))] - (mmax - tmpx[where(tmpy EQ max(tmpy))]) 
endelse
start = [median(continuum[ll:hh+cut])*1.D, .5D*(mmin+mmax), .5D*(mmax - mmin), .5D*(mmax-mmin)*max(tmpy)]
;gfit = gaussfit(tmpx,tmpy,coeffs,measure_errors=1/sqrt(err2),estimates=start,nterms=4,chisq=chisq)
gfit = MPFITEXPR(expr,tmpx,tmpy,err2,start,yfit=yfit,perror=perror,BESTNORM=bestnorm,dof=dof,/quiet)
;
if (keyword_set(plot)) then begin
print, 'fitting the line ...', gfit
oplot, tmpx, tmpy, col='FF00FF'x,thick=2
;oplot, tmpx, gfit(0)+gauss1(tmpx, gfit(1:3)),col='FF00FF'x,thick=2
oplot, tmpx, tmpy - (-continuum[ll:hh+cut]+gfit(0)+gauss1(tmpx,gfit(1:3))),col='000000'x,thick=2
pause
endif  ;plot

yy = tmpy - (-continuum[ll:hh+cut]+gfit(0)+gauss1(tmpx, gfit(1:3)))
endif else yy = tmpy ;dont fit the line - too few params


;endfor ; END blending lines ----------------------

index = (indgen(2*cut+1)+center-cut)[ll:hh+cut]
flux_out[min(index):max(index)] = yy 

endif ; ----- CTRL NE 0 -----
;oplot, wave[index],flux[index], psym=-3, thick=2, col='000000'x
endfor ; gas emission lines inside spectrum

return, flux_out

END