virusw_special_reduction.py 文件源码

def subtract_sky_from_sci(sci, sky_list, sky_model_list, wave, use_sci=False):
    sci.log.info('Subtracting sky from %s' %sci.basename)
    if use_sci:
        fn = op.join(sci.path, 'sky_model.txt')
        if True:#not op.exists(fn):
            F = fits.open('panacea/leoI_20131209.fits')
            G = get_fits(sci)
            get_master_sky(G['wavelength'].data, 
                           G['spectrum'].data - F[0].data*sci.exptime,
                           G['fiber_to_fiber'].data, sci.path, sci.exptime)
        wave, sky_model = np.loadtxt(fn, unpack=True)
    else:   
        weight, sorted_ind = get_interp_weights(sci, sky_list)
        arr_list = []
        for i,w in enumerate(weight):
            if w > 0.0:
                arr_list.append(w * sky_model_list[sorted_ind[i]])
        sky_model = np.sum(arr_list, axis=(0,))
    SCI = get_fits(sci)
    sky_subtracted = np.zeros(SCI['spectrum'].data.shape)
    for i, spec in enumerate(SCI['spectrum'].data):
        sky_subtracted[i,:] = (spec - np.interp(SCI['wavelength'].data[i,:],
                                                wave, sky_model*sci.exptime) 
                                      * SCI['fiber_to_fiber'].data[i,:])
    s = fits.ImageHDU(sky_subtracted)
    erase = []
    for i,S in enumerate(SCI):
        if S.header['EXTNAME'] == 'sky_subtracted':
            erase.append(i)
    for i in sorted(erase,reverse=True):
        del SCI[i]
    SCI.append(s)
    SCI[-1].header['EXTNAME'] = 'sky_subtracted'
    write_fits(SCI)            
    return SCI['wavelength'].data, sky_subtracted, wave, sky_model