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,verbose=verbose,$ corr_2o=corr_2o,$ basel_atm=basel_atm,calib_model=calib_model,ndeg_mod=ndeg_mod,$ noc3450=noc3450,thr3450=thr3450,cdeg=cdeg,$ crd_corr_bvr=crd_corr_bvr,$ teff_min=teff_min,dir_output=dir_output,ignore_blacklist=ignore_blacklist,png=png if(n_elements(thr3450) ne 1) then thr3450=0.02 if(n_elements(cdeg) ne 1) then cdeg=3 if(n_elements(dir_output) ne 1) then dir_output='I_corr_by_seg/' if(n_elements(ndeg_mod) ne 1) then ndeg_mod=3 if(n_elements(teff_min) ne 1) then teff_min=4100.0 nights_2o=read_asc('calib_5450_second_order.txt',/str) if(keyword_set(corr_2o) and (n_elements(basel_atm) lt 1)) then basel_atm=read_basel_atm() iu_atm_par=get_atm_par() ;;;; coefficients for global corrections determined using "comp_phot" kcbv=[-0.021596466d,7.0420563d-05,-1.2461569d-05,6.8544417d-08] k2dcvr=[-0.017889214d,-0.0015763134d,2.1666991d-05,-1.6652234d-07,0.00066270271d,$ 3.2158038d-06,-7.5264603d-08,-1.3061167d-06,-2.2686052d-09,-1.1031633d-09] crd_corr_bv=poly(iu_atm_par._dej2000,kcbv) crd_corr_vr=poly2d(iu_atm_par._raj2000,iu_atm_par._dej2000,k2dcvr,deg1=3,deg2=3,/irr) ;iu_atm_par[where(iu_atm_par.name eq '29763')].teff=8000.0 ;iu_atm_par[where(iu_atm_par.name eq '29763')]._fe_h_=-1.0 ;iu_atm_par[where(iu_atm_par.name eq '29763')].log_g_=4.0 if(keyword_set(cont_corr)) then wbvr_data=get_wbvr() 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 if((iu_atm_par[atm_idx].TEFF) eq 0) then iu_atm_par[atm_idx].TEFF=3000.0 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 if(keyword_set(calib_model)) then begin if(iu_atm_par[atm_idx].TEFF gt teff_min) then begin best_model=get_best_model(iu_atm_par[atm_idx],$ basel_atm,k_teff=3.0,$ ;;; corr_ext=corr_ext,ext_k_norm=ext_k_norm,airmass=airmass,$ ;;; spectra are already ext.corrected /redden,bvobs=bvobs) best_model.flux=best_model.flux/mean(best_model.flux[where(abs(best_model.wave-5500.0) le 50.0)]) endif else begin calib_model=0 message,/inf,'Too cool star, not performing the correction' endelse endif endif else begin calib_model=0 message,/inf,'Atmosphere parameters not found endelse 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(verbose)) 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) bpix=[1411,2076,2077,2926] idxvec=findgen(4096) nbpix=n_elements(bpix) for i=0,num_ext-1 do begin r=mrdfits(b_filename,exttouse[i],h,/silent) r.spectrum[bpix]=!values.f_nan flagarr=bytarr(n_elements(r.spectrum)) flagarr[bpix]=1 gpix=where(flagarr eq 0) r.spectrum[bpix]=interpol(r.spectrum[gpix],idxvec[gpix],idxvec[bpix],/quad) 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) or (strmid(extname,0,7) eq 'a951219')) and $ (not keyword_set(ignore_blacklist))) then begin if(keyword_set(verbose)) 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(verbose)) 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(verbose)) then print,'Processing: ',extname,' with ',calibname,' nsp=',calibdata.n_spec g_calvec=where((finite(calibdata.corr_vec) eq 1) and $ (mask_exclreg(calibdata.wave,$ [[5876.,10.],[6563.0,40.],[6678.,10.],[6721.,5.]]) eq 0),cg_calvec) if(cg_calvec lt 0.1*n_elements(calibdata.corr_vec)) then begin if(keyword_set(verbose)) 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(verbose)) 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 ; if(keyword_set(calib_model)) then continue 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(verbose)) 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 flux_out_corr=wl_out*0.0 errvec_out_corr=wl_out*0.0 weight_out_corr=wl_out*0.0 corrvec_out_corr=wl_out*0.0 ref_r=mrdfits(b_filename,exttouse[0],h,/silent) ref_r.spectrum[[1411,2076,2077,2926]]=!values.f_nan flagarr=bytarr(n_elements(r.spectrum)) flagarr[bpix]=1 gpix=where(flagarr eq 0) ref_r.spectrum[bpix]=interpol(ref_r.spectrum[gpix],idxvec[gpix],idxvec[bpix],/quad) 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(verbose)) then print,' ' seg_info=replicate({star_id:star_id,extname:'',extnum:0,$ wave_min:0.0,wave_max:0.0,wave_offset:0.0,$ snr:0.0,corr_model:0,model_id:0L,rms_fit:0.0,rms_poly:0.0},num_ext) for i=0,num_ext-1 do begin r=mrdfits(b_filename,exttouse[i],h,/silent) flagarr=bytarr(n_elements(r.spectrum)) flagarr[bpix]=1 gpix=where(flagarr eq 0) r.spectrum[bpix]=interpol(r.spectrum[gpix],idxvec[gpix],idxvec[bpix],/quad) 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) seg_info[i].extname=extname seg_info[i].wave_min=min(r.wavelength,max=wwmax) seg_info[i].wave_max=wwmax seg_info[i].wave_offset=-wloffset[i] seg_info[i].snr=snr if(i eq 0) then min_used_wl=seg_info[i].wave_min 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(verbose)) 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 corr_err=r.sigma/cal_vec_res ; if(keyword_set(calib_model)) then begin ; corr_flux_med=1.0 ; endif else begin if(keyword_set(verbose)) then print,'Number of overlapping pixels:', cwl_ol 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 ; endelse corr_err_corr=corr_err corr_flux=corr_flux/corr_flux_med corr_err=corr_err/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(verbose)) 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 corr_flux_corr=(keyword_set(calib_model))? corr_flux*corr_flux_med : corr_flux ref_corr_flux=corr_flux if(keyword_set(calib_model)) then begin exclreg1=[[4861.0,12.0],[5893.0,8.0],[6564.0,12.0],[7625.0,60.0],[8620.0,60.0],[9000,55]] min_wl_seg=min(r.wavelength,max=max_wl_seg) wave_seg=where(best_model.wave ge min_wl_seg-10.0 and best_model.wave le max_wl_seg+10.0) corr_str=corrvec_indous_lowres(wl_iu=r.wavelength,fl_iu=corr_flux_corr,$ wl_a=best_model.wave,fl_a=best_model.flux,$ ndeg=ndeg_mod,spr_lowres=0.0,exclreg=exclreg1,tail_frac=0.07,verbose=verbose,/const,$ rms_poly=rms_poly,rms_fit=rms_fit) if(keyword_set(verbose)) then print,'start_wl, rms_poly, rms_fit: ',min_wl_seg,rms_poly,rms_fit if(((rms_fit gt thr3450) or keyword_set(noc3450)) and $ (abs(min_wl_seg-3450.0) lt 150.0)) then begin ;;; bad blue segment mwlbluetr=3300.0 ;;;3750.0 wltruncred=where(r.wavelength gt mwlbluetr,compl=wltruncblue,ncompl=cwltruncblue) min_wltruncred=min(r.wavelength[wltruncred],min_wltruncred_idx,max=max_wltruncred,subscript_max=max_wltruncred_idx) if(min_wltruncred gt min_used_wl) then min_used_wl=min_wltruncred corr_str_a1=corrvec_indous_lowres(wl_iu=r.wavelength[wltruncred],$ fl_iu=(corr_flux_corr)[wltruncred],wl_a=best_model.wave,fl_a=best_model.flux,$ ndeg=ndeg_mod,spr_lowres=0.0,exclreg=exclreg1,tail_frac=0.07,verbose=verbose,/const,$ rms_poly=rms_poly_a1,rms_fit=rms_fit_a1) if((rms_fit_a1 le thr3450) and (rms_fit_a1 lt rms_fit) and $ (not keyword_set(noc3450)) and (cwltruncblue gt 0)) then begin corr_str.c_vec_int[wltruncred]=corr_str_a1.c_vec_int corr_str.c_vec_int[wltruncblue]=corr_str_a1.c_vec_int[min_wltruncred_idx] rms_poly=rms_poly_a1 rms_fit=rms_fit_a1 if(keyword_set(verbose)) then print,'Using correction computed in the truncated wavelength range' if(keyword_set(verbose)) then print,'start_wl, rms_poly, rms_fit: ',min_wl_seg,rms_poly,rms_fit endif else begin ;;;; using red end rms_poly=9.999 rms_fit =9.999 corr_str.c_vec_int[*]=corr_str_a1.c_vec_int[max_wltruncred_idx] if(max_wltruncred gt min_used_wl) then min_used_wl=max_wltruncred if(keyword_set(verbose)) then print,'Using max(c_vec_int)' endelse endif seg_info[i].corr_model=1 seg_info[i].model_id=atm_idx seg_info[i].rms_fit=rms_fit seg_info[i].rms_poly=rms_poly corr_flux_corr=corr_flux_corr/corr_str.c_vec_int corr_err_corr=corr_err_corr/corr_str.c_vec_int endif 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(calib_model)) then $ oplot,r.wavelength,corr_str.c_vec_int/median(corr_str.c_vec_int),linest=0 ;;;col=50+(205.0/num_ext)*i,linest=0 if(keyword_set(debug)) then begin aaa='' read,aaa endif endif fl_tmp=interpol(corr_flux,r.wavelength,wl_out) err_tmp=interpol(corr_err,r.wavelength,wl_out) fl_tmp_corr=interpol(corr_flux_corr,r.wavelength,wl_out) err_tmp_corr=interpol(corr_err_corr,r.wavelength,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 weight_tmp_corr=fl_tmp_corr*0.0+1.0 if(cbwl gt 0) then begin fl_tmp[bwl]=0.0 err_tmp[bwl]=0.0 weight_tmp[bwl]=0.0 fl_tmp_corr[bwl]=0.0 err_tmp_corr[bwl]=0.0 weight_tmp_corr[bwl]=0.0 endif wvec=create_weight_vec(fl_tmp,win=win_sm) flux_out+=fl_tmp*wvec*snr errvec_out=sqrt(errvec_out^2+(err_tmp*wvec*snr)^2) weight_out+=weight_tmp*wvec*snr flux_out_corr+=fl_tmp_corr*wvec*snr errvec_out_corr=sqrt(errvec_out_corr^2+(err_tmp_corr*wvec*snr)^2) weight_out_corr+=weight_tmp_corr*wvec*snr if(keyword_set(plot)) then $ oplot,wl_out,wvec,col=50+(205.0/num_ext)*i endfor flux_out=flux_out/weight_out errvec_out=errvec_out/weight_out flux_out_corr=flux_out_corr/weight_out_corr errvec_out_corr=errvec_out_corr/weight_out_corr if(keyword_set(plot)) then begin oplot,wl_out,flux_out oplot,wl_out,errvec_out,col=40 if(keyword_set(calib_model)) then $ oplot,best_model.wave,best_model.flux,col=254 endif 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 begin flux_norm=total(flux_out[fl_norm])/cfl_norm flux_out=flux_out/flux_norm errvec_out=errvec_out/flux_norm flux_norm_corr=total(flux_out_corr[fl_norm])/cfl_norm flux_out_corr=flux_out_corr/flux_norm_corr errvec_out_corr=errvec_out_corr/flux_norm_corr endif else begin if(keyword_set(calib_model)) then begin flux_norm=mean(best_model.flux[where(abs(best_model.wave-5500.0) le 50.0)],/nan) flux_norm_corr=flux_norm endif else begin flux_norm=1.0 flux_norm_corr=1.0 endelse endelse if(keyword_set(calib_model)) then begin flux_diff=flux_out_corr/flux_out wl_r=where(finite(flux_diff) eq 1 and wl_out ge min_used_wl,cwl_r) if(cwl_r gt 10) then begin cfit_cont=poly_fit(wl_out[wl_r]/60.0,flux_diff[wl_r],cdeg) corr_cont=poly(wl_out/60.0,cfit_cont) wl_c_b=where(wl_out lt min_used_wl,cwl_c_b) if(cwl_c_b gt 0) then corr_cont[wl_c_b]=corr_cont[wl_r[0]] flux_out_corr=flux_out_corr/corr_cont errvec_out_corr=errvec_out_corr/corr_cont mod_cont=poly(best_model.wave/60.0,cfit_cont) mod_cont[where(best_model.wave le min_used_wl)]=corr_cont[wl_r[0]] mod_cont[where(best_model.wave gt wl_out[wl_r[cwl_r-1]])]=corr_cont[wl_r[cwl_r-1]] best_model.flux=best_model.flux/mod_cont endif endif if(keyword_set(crd_corr_bvr) and (catm_idx eq 1)) then begin leff_B=4361.0 leff_V=5448.0 leff_R=6407.0 if(keyword_set(verbose)) then print,'Correcting continuum' xcrd1=[leff_B,leff_V] ycrd1=[-crd_corr_bv[atm_idx],0]*0.9 kcv1=poly_fit(xcrd1,ycrd1,1) xcrd2=[leff_V,leff_R] ycrd2=[0,crd_corr_vr[atm_idx]]*0.6 kcv2=poly_fit(xcrd2,ycrd2,1) sgn_h=(kcv1[1] ge kcv2[1])? -1.0 : 1.0 a=(kcv1[1]+kcv2[1])/2.0 b=-(kcv1[1]-kcv2[1])/2.0 corrvec=10^(-0.4*(sgn_h*sqrt((wl_out-leff_V)^2*b^2+2e-5)+a*(wl_out-leff_V))) wl5500=where(abs(wl_out-5500) lt 50.0) corrvec=corrvec/median(corrvec[wl5500]) flux_out=flux_out*corrvec errvec_out=errvec_out*corrvec if(keyword_set(calib_model)) then begin flux_out_corr=flux_out_corr*corrvec errvec_out_corr=errvec_out_corr*corrvec best_model.flux=best_model.flux*interpol(corrvec,wl_out,best_model.wave) endif endif if(not keyword_set(plot)) then begin plot,wl_out,flux_out,xr=[3400,9500],xs=1,ys=1,yr=[0.0,2.0],title=star_id,xtitle='Wavelength' oplot,wl_out,flux_out_corr,col=80 ;oplot,wl_out,flux_out/flux_out_corr,thick=3,col=192 if(keyword_set(calib_model)) then oplot,best_model.wave,best_model.flux/flux_norm_corr,col=254 oplot,wl_out,errvec_out,col=80 endif if(keyword_set(png)) then write_png,dir_output+'/png/'+star_id+'.png',tvrd(/true) if(keyword_set(calib_model)) then begin flux_out=flux_out_corr errvec_out=errvec_out_corr endif 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,dir_output+star_id+'.fits',float(flux_out),h_out h_ext1=h_out sxdelpar,h_ext1,'SIMPLE' sxaddpar,h_ext1,'XTENSION','IMAGE' mwrfits,float(errvec_out),dir_output+star_id+'.fits',h_ext1 mwrfits,seg_info,dir_output+star_id+'.fits' end