Source Code for Module ComboCode.cc.ivs.tools.spectra_combine

 1  # -*- coding: utf-8 -*- 
 2   
 3  """ 
 4  The combine method from the ivs.spectra.tools module. 
 5   
 6  """ 
 7   
 8  import logging 
 9  import numpy as np 
10  import scipy.stats 
11  from cc.ivs.units import conversions 
12   
13  logger = logging.getLogger("SPEC.TOOLS") 
14   
15 -def combine(list_of_spectra,R=200.,lambda0=(950.,'AA'),lambdan=(3350.,'AA')): 
16      """ 
17      Combine and weight-average spectra on a common wavelength grid. 
18       
19      C{list_of_spectra} should be a list of lists/arrays. Each element in the 
20      main list should be (wavelength,flux,error). 
21       
22      If you have FUSE fits files, use L{cc.ivs.fits.read_fuse}. 
23      If you have IUE FITS files, use L{cc.ivs.fits.read_iue}. 
24       
25      After Peter Woitke. 
26       
27      @param R: resolution 
28      @type R: float 
29      @param lambda0: start wavelength, unit 
30      @type lambda0: tuple (float,str) 
31      @param lambdan: end wavelength, unit 
32      @type lambdan: tuple (float,str) 
33      @return: binned spectrum (wavelengths,flux, error) 
34      @rtype: array, array, array 
35      """ 
36      l0 = conversions.convert(lambda0[1],'AA',lambda0[0]) 
37      ln = conversions.convert(lambdan[1],'AA',lambdan[0]) 
38      #-- STEP 1: define wavelength bins 
39      Delta = np.log10(1.+1./R) 
40      x = np.arange(np.log10(l0),np.log10(ln)+Delta,Delta) 
41      x = 10**x 
42      lamc_j = 0.5*(np.roll(x,1)+x) 
43   
44      #-- STEP 2: rebinning of data onto newly defined wavelength bins 
45      Ns = len(list_of_spectra) 
46      Nw = len(lamc_j)-1 
47      binned_fluxes = np.zeros((Ns,Nw)) 
48      binned_errors = np.inf*np.ones((Ns,Nw)) 
49   
50      for snr,(wave,flux,err) in enumerate(list_of_spectra): 
51          wave0 = np.roll(wave,1) 
52          wave1 = np.roll(wave,-1) 
53          lam_i0_dc = 0.5*(wave0+wave) 
54          lam_i1_dc = 0.5*(wave1+wave) 
55          dlam_i = lam_i1_dc-lam_i0_dc 
56           
57          for j in range(Nw): 
58              A = np.min(np.vstack([lamc_j[j+1]*np.ones(len(wave)),lam_i1_dc]),axis=0) 
59              B = np.max(np.vstack([lamc_j[j]*np.ones(len(wave)),lam_i0_dc]),axis=0) 
60              overlaps = scipy.stats.threshold(A-B,threshmin=0) 
61              norm = np.sum(overlaps) 
62              binned_fluxes[snr,j] = np.sum(flux*overlaps)/norm 
63              binned_errors[snr,j] = np.sqrt(np.sum((err*overlaps)**2))/norm 
64       
65      #-- STEP 3: all available spectra sets are co-added, using the inverse 
66      #   square of the bin uncertainty as weight 
67      binned_fluxes[np.isnan(binned_fluxes)] = 0 
68      binned_errors[np.isnan(binned_errors)] = 1e300 
69      weights = 1./binned_errors**2 
70       
71      totalflux = np.sum(weights*binned_fluxes,axis=0)/np.sum(weights,axis=0) 
72      totalerr = np.sqrt(np.sum((weights*binned_errors)**2,axis=0))/np.sum(weights,axis=0) 
73      totalspec = np.sum(binned_fluxes>0,axis=0) 
74       
75      #-- that's it! 
76      return x[:-1],totalflux,totalerr,totalspec 
77