function correct_corr_vec_tails,calibdata,fraction=fraction,nored=nored,noblue=noblue
     if(n_elements(fraction) ne 1) then fraction=0.05
     v=calibdata.corr_vec
     wlmin=min(calibdata.wave,max=wlmax)
     wlr=wlmax-wlmin
     btail=where(calibdata.wave ge wlmin and calibdata.wave lt (wlmin+wlr*fraction))
     breg=where((calibdata.wave ge wlmin+wlr*fraction) and $
                (calibdata.wave lt (wlmin+wlr*2.0*fraction)) and $
                (finite(v) eq 1), cbreg)
     rtail=where(calibdata.wave le wlmax and calibdata.wave gt (wlmax-wlr*fraction))
     rreg=where((calibdata.wave le wlmax-wlr*fraction) and $
                (calibdata.wave gt (wlmax-wlr*2.0*fraction)) and $
                (finite(v) eq 1), crreg)
     if((cbreg gt 2) and (not keyword_set(noblue))) then begin
         klb=linfit(calibdata.wave[breg],calibdata.corr_vec[breg])
         ww=findgen(cbreg)/(double(cbreg-1))
         v[breg]=(klb[0]+klb[1]*calibdata.wave[breg])*(1.0-ww)+calibdata.corr_vec[breg]*ww
         v[btail]=klb[0]+klb[1]*calibdata.wave[btail]
     endif
     if((crreg gt 2) and (not keyword_set(nored))) then begin
         klr=linfit(calibdata.wave[rreg],calibdata.corr_vec[rreg])
         ww=findgen(crreg)/(double(crreg-1))
         v[rreg]=(klr[0]+klr[1]*calibdata.wave[rreg])*ww+calibdata.corr_vec[rreg]*(1.0-ww)
         v[rtail]=klr[0]+klr[1]*calibdata.wave[rtail]
     endif
     return,v
end


function create_weight_vec,datavec,win=win
    if(n_elements(win) ne 1) then win=fix(n_elements(datavec)*0.03)

    n_dv=n_elements(datavec)
    wvec=dblarr(n_dv)
    

    b_datavec=where((finite(datavec) ne 1) or (datavec eq 0), cb_datavec)
    if(cb_datavec gt 0 and win gt 0) then begin
        krnl=dblarr(win*2+1)
        krnl[0:win]=findgen(win+1)
        krnl[win+1:*]=reverse(krnl[0:win-1])
        krnl=krnl/max(krnl)

        for i=0,cb_datavec-1 do begin
            ndmin=(b_datavec[i]-win) > 0
            ndmax=(b_datavec[i]+win) < (n_dv-1)
            kdmin=(win-b_datavec[i]) > 0
            kdmax=kdmin+(ndmax-ndmin)
            wvec[ndmin:ndmax]=(krnl[kdmin:kdmax] > wvec[ndmin:ndmax])
        endfor

        wvec=wvec/max(wvec)
    endif

    return,1.0-wvec
end

function mask_exclreg,wl,exclreg
    s_e=size(exclreg)
    res=byte(wl*0.0)
    if(s_e[0] ne 2) then return,res

    for i=0,s_e[2]-1 do begin
        minwl=exclreg[0,i]-exclreg[1,i]
        maxwl=exclreg[0,i]+exclreg[1,i]
        badreg=where((wl ge minwl) and (wl le maxwl),cbad)
        if(cbad gt 0) then res[badreg]=1
    endfor

    return,res
end

pro correct_indous_spec_bintable,star_id,$
    dir_bintable=dir_bintable,dir_nightcorr=dir_nightcorr,tail_frac=tail_frac,$
    plot=plot,wl_out=wl_out,minsnr=minsnr,win_sm=win_sm,corr_ext=corr_ext,debug=debug,$
    corr_2o=corr_2o,basel_atm=basel_atm

nights_2o=read_asc('calib_5450_second_order.txt',/str)

iu_atm_par=get_atm_par()

wlcorr=read_asc('extnames_indous_all_wlcorr.txt',/str)
goodwlcorr=where(wlcorr[0,*] eq star_id,cgoodwlcorr)
wlcorr=(cgoodwlcorr gt 0)? wlcorr[*,goodwlcorr] : ['nothing','a0000009999','-1','0.0','0.0','0.0']
if(cgoodwlcorr eq 0) then message,/inf,'Wavelength shift information not found for '+star_id

if(n_elements(win_sm) ne 1) then win_sm=300.0

exclreg=[[5893.0,8.0],[10000.0,10.0]]

if(n_elements(minsnr) ne 1) then minsnr=10.0

wl_out=3465.0+0.4*findgen(15011)

bl6800=read_asc('blacklist_6800.txt',/str)

wlmin_ref=[3450.0,4100.0,4700.0,5450.0,5950.0,6800.0,7300.0,8100.0]

if(n_elements(dir_bintable) ne 1) then dir_bintable='BINTABLE/'
if(n_elements(dir_iraf) ne 1) then dir_iraf='IRAF/'
if(n_elements(dir_nightcorr) ne 1) then dir_nightcorr='night_seg_corr/'
b_filename=dir_bintable+star_id+'.fits'
i_filename=dir_iraf+star_id+'.fits'
dir_published='I/'
o_filename=dir_published+star_id+'.fits'

r_pub=readfits(o_filename,h_pub,/silent)
bvobs=sxpar(h_pub,'OBJBMAG')-sxpar(h_pub,'OBJVMAG')
parse_spechdr,h_pub,wl=wl_pub

atm_idx=where(iu_atm_par.name eq star_id, catm_idx)
if(catm_idx eq 1) then begin
    sxaddpar,h_pub,'TEFF_WU',iu_atm_par[atm_idx].TEFF
    sxaddpar,h_pub,'ETEFF_WU',iu_atm_par[atm_idx].E_TEFF
    sxaddpar,h_pub,'LOGG_WU',iu_atm_par[atm_idx].LOG_G_
    sxaddpar,h_pub,'ELOGG_WU',iu_atm_par[atm_idx].E_LOG_G_
    sxaddpar,h_pub,'FEH_WU',iu_atm_par[atm_idx]._FE_H_
    sxaddpar,h_pub,'EFEH_WU',iu_atm_par[atm_idx].E__FE_H_
    sxaddpar,h_pub,'RFEH_WU',iu_atm_par[atm_idx].R__FE_H_
    sxaddpar,h_pub,'SIMBADID',iu_atm_par[atm_idx].SIMBADNAME
endif



fits_open,b_filename,fcb
n_ext=fcb.nextend
fits_close,fcb
exttouse=indgen((n_ext-1)/2)+(n_ext-1)/2+2
if(keyword_set(debug)) then print,'Extensions:',exttouse
num_ext=n_elements(exttouse)

if(keyword_set(plot)) then $
    plot,[0],[0],/nodata,xs=1,ys=1,xr=[3400,9500],yr=[0,2.5],title=star_id,xtitle='Wavelength'

wlmin_arr=dblarr(num_ext)
flag_arr=bytarr(num_ext)
calibname_arr=strarr(num_ext)
date_wl_arr=strarr(num_ext)
wloffset=dblarr(num_ext)

for i=0,num_ext-1 do begin
    r=mrdfits(b_filename,exttouse[i],h,/silent)

    snr=median(r.spectrum/r.sigma)
    extname=sxpar(h,'EXTNAME')
    wlcchk=where(wlcorr[1,*] eq extname, cwlcchk)
    wloffset[i]=(cwlcchk eq 1)? double(wlcorr[5,wlcchk]) : 0d
    r.wavelength=r.wavelength-wloffset[i]
;    print,'Extname,wloffset:',sxpar(h,'EXTNAME'),wloffset[i]

    blchk=where(bl6800 eq extname, cblchk)
    if(cblchk gt 0  or  snr le minsnr) then begin
        if(keyword_set(debug)) then print,'extname: '+extname+' is blacklisted'
        continue
    endif
    date=strmid(extname,0,7)
    wlmtmp=min(abs(wlmin_ref-min(r.wavelength)),wlidx)
    wlmin_arr[i]=min(r.wavelength) ;;; wlmin_ref[wlidx]
    date_wl=date+'_'+string(wlmin_ref[wlidx],format='(i4)')
    calibname=dir_nightcorr+'/'+date_wl+'.fits'
    date_wl_arr[i]=date_wl
    calibname_arr[i]=calibname
    if(file_test(calibname) ne 1) then begin
        if(keyword_set(debug)) then print,'The calibration file '+calibname+' does not exist'
        continue ;;;; calibration does not exist for that night
    endif
    calibdata=mrdfits(calibname,1,/silent)
    if(keyword_set(debug)) then print,'Processing: ',extname,' with ',calibname,' nsp=',calibdata.n_spec
    g_calvec=where(finite(calibdata.corr_vec) eq 1,cg_calvec)
    if(cg_calvec lt 0.1*n_elements(calibdata.corr_vec)) then begin
        if(keyword_set(debug)) then print,'Too few good points in: ',calibname
        continue ;;;; too few good points
    endif
    nored=(wlmin_arr[i] lt 6700)?0 : 1
    noblue=(wlmin_arr[i] lt 8000)?0 : 1
    corr_vec_tails=correct_corr_vec_tails(calibdata,fraction=tail_frac,nored=nored,noblue=noblue)
    ;;if(wlmin_arr[i] lt 6700) then 
    calibdata.corr_vec=corr_vec_tails

    if(keyword_set(corr_ext)) then begin
        r_ref=mrdfits(i_filename,exttouse[i]-1,hr_raw,/silent)
        airmass=sxpar(hr_raw,'AIRMASS')
        ext_vec=get_ext_vec(r.wavelength,k_norm=calibdata.ext_k_norm)*airmass
        if(keyword_set(debug)) then print,'Disperser: ',sxpar(hr_raw,'DISPERSE'),'  airmass=',airmass
    endif else ext_vec=r.wavelength*0.0

    r.spectrum=r.spectrum*10^(0.4*ext_vec)

    flag_arr[i]=1
    cal_vec_res=interpol(calibdata.corr_vec[g_calvec],calibdata.wave[g_calvec],r.wavelength)
    corr_flux=r.spectrum/cal_vec_res
    corr_flux_med=median(corr_flux)

    cont_2o_model=0.0*corr_flux
    if(keyword_set(corr_2o) and (catm_idx gt 0)) then begin
        bnight=where(date_wl_arr[i] eq nights_2o[0,*] and nights_2o[1,*] eq 'yes',cbnight)
        if(cbnight gt 0) then begin
            if(keyword_set(debug)) then print,'Correcting the 2nd order contamination'
            cont_2o_model=get_second_order_cont(r.wavelength,corr_flux/corr_flux_med,$
                basel_atm=basel_atm,star_atm_par=iu_atm_par[atm_idx],/redden,bvobs=bvobs,dqemult=dqemult,$
                corr_ext=corr_ext,ext_k_norm=calibdata.ext_k_norm,airmass=airmass)
            cont_2o_model=cont_2o_model*dqemult
        endif
    endif
    if(keyword_set(plot)) then $
        oplot,r.wavelength,corr_flux/corr_flux_med-cont_2o_model,col=50+(205.0/num_ext)*i
endfor
;;print,wloffset
good_seg=where(flag_arr eq 1, cgs)
if(cgs eq 0) then begin
    message,'No good segments in '+b_filename,/inf
    return
endif

;print,wloffset
;print,wlmin_arr
;print,exttouse

exttouse=exttouse[good_seg]
wlmin_arr=wlmin_arr[good_seg]
calibname_arr=calibname_arr[good_seg]
date_wl_arr=date_wl_arr[good_seg]
wloffset=wloffset[good_seg]
num_ext=cgs

s_wl=sort(wlmin_arr)

wlmin_arr=wlmin_arr[s_wl]
exttouse=exttouse[s_wl]
calibname_arr=calibname_arr[s_wl]
date_wl_arr=date_wl_arr[s_wl]
wloffset=wloffset[s_wl]

;print,wloffset
;print,wlmin_arr
;print,exttouse

if(keyword_set(plot)) then $
    plot,[0],[0],/nodata,xs=1,ys=1,xr=[3400,9500],yr=[0,2.5],title=star_id,xtitle='Wavelength'

flux_out=wl_out*0.0
errvec_out=wl_out*0.0
weight_out=wl_out*0.0
corrvec_out=wl_out*0.0

ref_r=mrdfits(b_filename,exttouse[0],h,/silent)
ref_r.wavelength=ref_r.wavelength-wloffset[0]
;print,'         !!!!!!!!!!!!!!!!          '
;print,'Extname,wloffset:',sxpar(h,'EXTNAME'),wloffset[0]

calibdata=mrdfits(calibname_arr[0],1,/silent)

if(keyword_set(corr_ext)) then begin
    r_ref=mrdfits(i_filename,exttouse[0]-1,hr_raw,/silent)
    airmass=sxpar(hr_raw,'AIRMASS')
    ext_vec=get_ext_vec(ref_r.wavelength,k_norm=ext_k_norm)*airmass
endif else ext_vec=ref_r.wavelength*0.0

ref_r.spectrum=ref_r.spectrum*10^(0.4*ext_vec)

;;print,'Processing: ',sxpar(h,'EXTNAME'),' with ',calibname_arr[0],' nsp=',calibdata.n_spec
nored=(wlmin_arr[0] lt 6700)?0 : 1
noblue=(wlmin_arr[0] lt 8000)?0 : 1
corr_vec_tails=correct_corr_vec_tails(calibdata,fraction=tail_frac,nored=nored,noblue=noblue)
;;;if(wlmin_arr[0] lt 6700) then $
calibdata.corr_vec=corr_vec_tails
g_calvec=where(finite(calibdata.corr_vec) eq 1,cg_calvec)


cal_vec_res=interpol(calibdata.corr_vec[g_calvec],calibdata.wave[g_calvec],ref_r.wavelength)
ref_wl=ref_r.wavelength
ref_corr_flux=ref_r.spectrum/cal_vec_res
ref_corr_flux_med=median(ref_corr_flux)
ref_corr_flux=ref_corr_flux/ref_corr_flux_med
if(keyword_set(plot)) then $
    oplot,ref_wl,ref_corr_flux,col=50+(205.0/num_ext)*0

if(keyword_set(debug)) then print,' '

for i=0,num_ext-1 do begin
    r=mrdfits(b_filename,exttouse[i],h,/silent)
    r.wavelength=r.wavelength-wloffset[i]
;    print,'Extname,wloffset:',sxpar(h,'EXTNAME'),wloffset[i]

    snr=median(r.spectrum/r.sigma)
    extname=sxpar(h,'EXTNAME')
    date=strmid(extname,0,7)
    calibname=calibname_arr[i]
    calibdata=mrdfits(calibname,1,/silent)

    if(keyword_set(corr_ext)) then begin
        r_ref=mrdfits(i_filename,exttouse[i]-1,hr_raw,/silent)
        airmass=sxpar(hr_raw,'AIRMASS')
        ext_vec=get_ext_vec(r.wavelength,k_norm=ext_k_norm)*airmass
    endif else ext_vec=r.wavelength*0.0

    r.spectrum=r.spectrum*10^(0.4*ext_vec)

    if(keyword_set(debug)) then print,'Processing: ',extname,' with ',calibname,' nsp=',calibdata.n_spec
    nored=(wlmin_arr[i] lt 6700)?0 : 1
    noblue=(wlmin_arr[i] lt 8000)?0 : 1
    corr_vec_tails=correct_corr_vec_tails(calibdata,fraction=tail_frac,nored=nored,noblue=noblue)
    ;;;if(wlmin_arr[i] lt 6700) then $
    calibdata.corr_vec=corr_vec_tails
    g_calvec=where(finite(calibdata.corr_vec) eq 1,cg_calvec)

    cal_vec_res=interpol(calibdata.corr_vec[g_calvec],calibdata.wave[g_calvec],r.wavelength)
    wl_ol=where((r.wavelength le max(ref_wl)) and $
                (mask_exclreg(r.wavelength,exclreg) eq 0),cwl_ol)
    corr_flux=r.spectrum/cal_vec_res
    if(cwl_ol gt 11) then begin
        wl_olref=where((ref_wl gt min(r.wavelength)) and $
                       (mask_exclreg(ref_wl,exclreg) eq 0),cwl_olref)
        corr_flux_med=1.0/(median(ref_corr_flux[wl_olref])/median(corr_flux[wl_ol]))
    endif else begin
        owl_cur_min=min(ref_wl,max=owl_cur_max)
        owl_olmin=owl_cur_min+0.90*(owl_cur_max-owl_cur_min)
        owl_ol=where(ref_wl ge owl_olmin)
        owl_ol_pub=where(wl_pub ge owl_olmin and wl_pub le owl_cur_max)

        wl_cur_min=min(r.wavelength,max=wl_cur_max)
        wl_olmax=wl_cur_min+0.10*(wl_cur_max-wl_cur_min)
        wl_ol=where(r.wavelength le wl_olmax)
        wl_ol_pub=where(wl_pub ge wl_cur_min and wl_pub le wl_olmax)

        corr_flux_med=median(corr_flux[wl_ol])/(median(ref_corr_flux[owl_ol])*median(r_pub[wl_ol_pub])/median(r_pub[owl_ol_pub]))
    endelse
    corr_flux=corr_flux/corr_flux_med
    ref_wl=r.wavelength

    cont_2o_model=0.0*corr_flux
    if(keyword_set(corr_2o) and (catm_idx gt 0)) then begin
        bnight=where(date_wl_arr[i] eq nights_2o[0,*] and nights_2o[1,*] eq 'yes',cbnight)
        if(cbnight gt 0) then begin
            if(keyword_set(debug)) then print,'Correcting the 2nd order contamination'
            cont_2o_model=get_second_order_cont(r.wavelength,corr_flux,$
                basel_atm=basel_atm,star_atm_par=iu_atm_par[atm_idx],/redden,bvobs=bvobs,dqemult=dqemult,$
                corr_ext=corr_ext,ext_k_norm=calibdata.ext_k_norm,airmass=airmass)
            cont_2o_model=cont_2o_model*dqemult
        endif
    endif

    corr_flux=corr_flux-cont_2o_model
    ref_corr_flux=corr_flux

    if(keyword_set(plot)) then begin
        oplot,r.wavelength,corr_flux,col=50+(205.0/num_ext)*i
        oplot,r.wavelength,cal_vec_res,col=50+(205.0/num_ext)*i,linest=2
        if(keyword_set(debug)) then begin
            aaa=''
            read,aaa
        endif
    endif

    fl_tmp=interpol(corr_flux,r.wavelength,wl_out)
;    err_tmp=interpol(err_iu/c_str.c_vec_int,wl_iu,wl_out)
;    cv_tmp=interpol(c_str.c_vec_int,wl_iu,wl_out)
    bwl=where((wl_out le min(r.wavelength,/nan)) or $
              (wl_out ge max(r.wavelength,/nan)),cbwl,comp=gwl,ncompl=cgwl)
    weight_tmp=fl_tmp*0.0+1.0
    if(cbwl gt 0) then begin
        fl_tmp[bwl]=0.0
;        err_tmp[bwl]=0.0
;        cv_tmp[bwl]=0.0
        weight_tmp[bwl]=0.0
    endif
    wvec=create_weight_vec(fl_tmp,win=win_sm)
    flux_out+=fl_tmp*wvec*snr
;    errvec_out+=err_tmp
    weight_out+=weight_tmp*wvec*snr
;    corrvec_out+=cv_tmp
    if(keyword_set(plot)) then $
        oplot,wl_out,wvec,col=50+(205.0/num_ext)*i

endfor

flux_out=flux_out/weight_out
if(keyword_set(plot)) then oplot,wl_out,flux_out


fl_norm=where(finite(flux_out) eq 1 and wl_out ge 5450.0 and wl_out le 5550.0, cfl_norm)
if(cfl_norm gt 0) then flux_out=flux_out/(total(flux_out[fl_norm])/cfl_norm)

if(not keyword_set(plot)) then $
    plot,wl_out,flux_out,xs=1,ys=1,yr=[0,1.8],title=star_id,xtitle='Wavelength'

sxaddpar,h_pub,'NCOMBINE',num_ext
h_out=h_pub
sxaddpar,h_pub,'CTYPE1','WAVE-AWAV'
sxaddpar,h_pub,'CRVAL1',wl_out[0]
sxaddpar,h_pub,'CDELT1',wl_out[1]-wl_out[0]
sxaddpar,h_pub,'CD1_1',wl_out[1]-wl_out[0]

writefits,'I_corr_by_seg/'+star_id+'.fits',float(flux_out),h_out

end