pro mask_skyline,specname,path=path,dir=dir,skywave,noplot=noplot,outdir=outdir

; locate the spectrum
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if keyword_set(path) and keyword_set(dir) then begin
       filename=strtrim(path,2)+strtrim(dir,2)+strtrim(specname,2)
 endif else begin
       filename=specname
endelse 

if not file_test(filename) then begin
	print,'The spectrum does not exist, please check again...'
    return
endif 
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

;read the spectrum
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
spec = MRDFITS(filename, 0, hdr)
flux = spec[*,0]
;invertvar = spec[1,*]
;andmask = spec[3,*]
;ormask = spec[4,*]

wave0 = sxpar(hdr, 'COEFF0')
dwave = sxpar(hdr, 'COEFF1')
Ra=sxpar(hdr,'RA')
Dec=sxpar(hdr,'DEC')
nx=sxpar(hdr,'NAXIS1')
logwaverange = lindgen(nx)*dwave + wave0
waverange = 10.0^(logwaverange)
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


; mark skyline in wavelength
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
skymask = intarr(nx)*0
skynum = n_elements(skywave)
for i = 0, skynum-1 do begin
	sel=where(abs(waverange-skywave[i]) lt 0.8,Nsel)
	if Nsel gt 0 then skymask[sel] = 1
endfor
mask=where(skymask eq 1,Nm)
print,nm,' sky emission line mask'
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


;%%% fit continuum and finding peaks
flux_fit = robust_poly_fit(waverange, flux, 6, flux_mean)
flux_residual = abs(flux - flux_mean)
flux_sigma = robust_sigma(flux_residual)
flux_peakfind = where((flux gt (flux_mean + flux_sigma*3.0) or flux lt flux_mean-flux_sigma*3.0) and flux ne 0,Npeak)
sky_peak= where((flux gt (flux_mean + flux_sigma*3.0) or flux lt flux_mean-flux_sigma*3.0) and flux ne 0 and skymask eq 1,Nskypeak)
;%%% subtract skyflux peaks caused peaks
sel=where((flux gt flux_mean+flux_sigma*3 or flux lt flux_mean-flux_sigma*3) and skymask eq 1 and flux ne 0.,fcount)
flux_mask=flux
if fcount gt 0 then flux_mask[sel] = flux_mean[sel]
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


; plot
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if keyword_set(noplot) then goto,next

psfile=replace_string(specname,'.fits','_msky.ps')
ps_start,psfile,/metric
;  define yrange
  if abs(max(flux_mean) - flux_sigma) gt 100. then begin
    ymax = max(flux_mean) + 1.2*max(flux_mean)
    ymin = min(flux_mean) - 1.2*max(flux_mean)
  endif else begin
    ymax = max(flux_mean) + 25.*flux_sigma
    ymin = min(flux_mean) - 25.*flux_sigma
  endelse

; plot origin spectrum, continuum and peaks
  cgplot, waverange, flux, xrange=[min(waverange)-150, max(waverange)+150], xstyle=1, $
  yrange=[ymin, ymax], ystyle=1, xthick=2, ythick=2, $
  charsize=1.5, charthick=2, xtitle=textoidl('Wavelength (\AA)'), ytitle='Flux',title=specname,noclip=0
  cgplots, waverange, flux_mean, color='green'
    
  if Nskypeak gt 0 then $
  cgplots, waverange[sky_peak], flux[sky_peak], psym=9, symsize=0.2, color='red',noclip=0
    
; plot sky peaks
  skywavey = fltarr(n_elements(skywave)) + (max(flux_mean) + 10.*flux_sigma)
  cgplots, skywave, skywavey, psym=11, symsize=0.4, color='blue'
  
; smooth and plot smooth spectrum
  flux_s = smooth(flux_mask, 5)
  cgplots, waverange, flux_s-0.55*ymax,noclip=0 
  cgtext, 0.75, 0.88, charthick=2, 'RA= '+strtrim(string(Ra),2), /normal
  cgtext, 0.75, 0.84, charthick=2, 'Dec= '+strtrim(string(Dec),2), /normal
ps_end   
; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

; no plot jump to here: output the spectrum after mask
next:

;%%% write reduced fits file
if not keyword_set(outdir) then outdir = 'redu_spec'
if File_test(outdir, /dir) eq 0 then file_mkdir, outdir

outfile=outdir+'/'+replace_string(specname,'.fits','_msky.fits')
sxdelpar, hdr, 'NAXIS2'
sxaddpar, hdr, 'NAXIS2',2

outtab=dblarr(nx,2)
outtab[*,0]=flux_mask
outtab[*,1]=skymask
mwrfits,outtab,outfile, hdr, /create

END