Source Code for Module ivs.spectra.moments

  1  """ 
  2  Compute the moments of a line profile. 
  3  """ 
  4  import numpy as np 
  5  import scipy.integrate 
  6  import pylab as pl 
  7   
  8 -def moments(velo,flux,SNR=200.,max_mom=3): 
  9      """ 
 10      Compute the moments from a line profile. 
 11      """ 
 12      mymoms,e_mymoms = np.zeros(max_mom+1),np.zeros(max_mom+1) 
 13       
 14      #-- compute zeroth moment (equivalent width) 
 15      m0 = scipy.integrate.simps( (1-flux),x=velo) 
 16       
 17      #-- to calculate the uncertainties, we need the error in each velocity 
 18      #   bin, and the error on the zeroth moment 
 19      sigma_i = 1./SNR / np.sqrt(flux) 
 20      Delta_v0 = np.sqrt(scipy.integrate.simps( sigma_i,x=velo)**2) 
 21      e_m0 = np.sqrt(2)*(Delta_v0/m0) 
 22       
 23      mymoms[0] = m0 
 24      e_mymoms[0] = e_m0 
 25           
 26      #-- calculate other moments 
 27      for n in range(1,max_mom+1): 
 28          mymoms[n] = scipy.integrate.simps((1-flux)*velo**n,x=velo)/ m0 
 29          #-- and the uncertainty on it 
 30          Delta_vn = np.sqrt(scipy.integrate.simps(sigma_i*velo**n)**2) 
 31          e_mymoms[n] = np.sqrt(  (Delta_vn/m0)**2 + (Delta_v0/m0 * mymoms[n])**2 ) 
 32       
 33      return mymoms,e_mymoms 
 34   
 35   
 36 -def moments_fromfiles(filelist,read_func,max_mom=3,**kwargs): 
 37      """ 
 38      Compute the moments from a list of files containg line profiles. 
 39       
 40      The C{read_func} should return 2 arrays and a float, representing velocities, 
 41      normalised fluxes and the SNR of the line profile. 
 42       
 43      m0: equivalent width 
 44      m1: radial velocity 
 45      m2: variance 
 46      m3: skewness 
 47       
 48      @param filelist: list of filenames 
 49      @type filelist: list of strings 
 50      @param read_func: function which reads in a file and returns velocities (array), 
 51      the line profile (array) and the SNR of the whole profile (float) 
 52      @type read_func: Python function 
 53      @param max_mom: maximum moment to compute 
 54      @type max_mom: integer 
 55      @return: a list containg the moments and a list containing the errors 
 56      @rtype: [max_mom x array],[max_mom x array] 
 57      """ 
 58       
 59      output = [np.zeros(len(filelist)) for i in range(max_mom+1)] 
 60      errors = [np.zeros(len(filelist)) for i in range(max_mom+1)] 
 61      for i,f in enumerate(filelist): 
 62          #-- read in the profile 
 63          velo,flux,SNR = read_func(f,**kwargs) 
 64          mymoms,e_mymoms = moments(velo,flux,SNR) 
 65          for n in range(len(mymoms)): 
 66              output[n][i] = mymoms[n] 
 67              errors[n][i] = e_mymoms[n] 
 68          pl.plot(velo,flux+i*0.01) 
 69      return output,errors 
 70   
 71 -def profiles_fromfiles(filelist,read_func,max_mom=3,**kwargs): 
 72      """ 
 73      Compute an average profile from a file list of profiles. 
 74       
 75      The C{read_func} should return 2 arrays and a float, representing velocities, 
 76      normalised fluxes and the SNR of the line profile. 
 77       
 78      m0: equivalent width 
 79      m1: radial velocity 
 80      m2: variance 
 81      m3: skewness 
 82       
 83      @param filelist: list of filenames 
 84      @type filelist: list of strings 
 85      @param read_func: function which reads in a file and returns velocities (array), 
 86      the line profile (array) and the SNR of the whole profile (float) 
 87      @type read_func: Python function 
 88      @param max_mom: maximum moment to compute 
 89      @type max_mom: integer 
 90      @return: a list containg the moments and a list containing the errors 
 91      @rtype: velo,av_prof,fluxes 
 92      """ 
 93       
 94      output = [np.zeros(len(filelist)) for i in range(max_mom+1)] 
 95      errors = [np.zeros(len(filelist)) for i in range(max_mom+1)] 
 96      for i,f in enumerate(filelist): 
 97          #-- read in the profile 
 98          velo,flux,SNR = read_func(f,**kwargs) 
 99          if i==0: 
100              template_velo = velo 
101              fluxes = np.zeros((len(filelist),len(flux))) 
102              fluxes[0] = flux 
103          else: 
104              fluxes[i] = np.interp(template_velo,velo,flux) 
105       
106      av_prof = fluxes.mean(axis=0) 
107      fluxes -= av_prof 
108       
109      return template_velo,av_prof,fluxes 
110