Source Code for Module ivs.asteroseismology.granulation

  1  """ 
  2  Simulation of a granulation signal 
  3   
  4  Author: Joris De Ridder 
  5   
  6  Error messages are written to the logger "granulation". 
  7  """ 
  8   
  9  import numpy as np 
 10  from numpy.random import normal 
 11  import logging 
 12   
 13   
 14  # Setup the logger. 
 15  # Add at least one handler to avoid the message "No handlers could be found"  
 16  # on the console. The NullHandler is part of the standard logging module only  
 17  # from Python 2.7 on. 
 18   
 19 -class NullHandler(logging.Handler): 
 20 -    def emit(self, record): 
 21          pass 
 22           
 23  logger = logging.getLogger("granulation") 
 24  nullHandler = NullHandler() 
 25  logger.addHandler(nullHandler) 
 26   
 27   
 28   
 29 -def granulation(time, timescale, varscale): 
 30   
 31      """ 
 32      Simulates a time series showing granulation variations 
 33       
 34      A first-order autoregressive process is used, as this gives a Harvey 
 35      model in the frequency domain. See also: 
 36      De Ridder et al., 2006, MNRAS 365, pp. 595-605. 
 37       
 38      @param time: time points 
 39      @type time: ndarray 
 40      @param timescale: array of time scale "tau_i" of each granulation component 
 41                        of the granulation/magnetic activity. Same units as 'time'. 
 42      @type timescale: ndarray 
 43      @param varscale: array of variation scale "sigma_i" of each component of the  
 44                       granulation/magnetic activity in the appropriate passband. 
 45                       Same size as the timescale array. Unit: ppm 
 46      @type varscale: ndarray 
 47      @return: the granulation signal 
 48      @rtype: ndarray 
 49       
 50      Example: 
 51       
 52      >>> time = np.linspace(0,100,200)             # E.g. in days 
 53      >>> timescale = np.array([5.0, 20.])          # time scales in days 
 54      >>> varscale = np.array([10.0, 50.0])         # variation scale in ppm 
 55      >>> gransignal = granulation(time, timescale, varscale) 
 56      >>> flux = 100000.0                           # mean flux level 
 57      >>> signal = flux * (1.0 + gransignal)        # signal in flux 
 58       
 59      """ 
 60       
 61       
 62      Ntime = len(time) 
 63      Ncomp = len(timescale) 
 64   
 65      logger.info("Simulating %d granulation components\n" % Ncomp) 
 66           
 67      # Set the kick (= reexcitation) timestep to be one 100th of the 
 68      # shortest granulation time scale (i.e. kick often enough). 
 69   
 70      kicktimestep = min(timescale) / 100.0 
 71       
 72      logger.info("Kicktimestep = %f\n" % kicktimestep) 
 73   
 74      # Predefine some arrays 
 75   
 76      signal = np.zeros(Ntime) 
 77      granul = np.zeros(Ncomp) 
 78      mu = np.zeros(Ncomp) 
 79      sigma = np.sqrt(kicktimestep/timescale)*varscale 
 80   
 81      # Warm up the first-order autoregressive process 
 82   
 83      logger.info("Granulation process warming up...\n") 
 84           
 85      for i in range(2000): 
 86          granul = granul * (1.0 - kicktimestep / timescale) + normal(mu, sigma) 
 87   
 88      # Start simulating the granulation time series 
 89   
 90      logger.info("Simulating granulation signal.\n") 
 91   
 92      delta = 0.0 
 93      currenttime = time[0] - kicktimestep 
 94   
 95      for i in range(Ntime): 
 96   
 97          # Compute the contribution of each component separately. 
 98          # First advance the time series right *before* the time point i, 
 99   
100          while((currenttime + kicktimestep) < time[i]): 
101              granul = granul * (1.0 - kicktimestep / timescale) + normal(mu,sigma) 
102              currenttime = currenttime + kicktimestep 
103   
104          # Then advance the time series with a small time step right *on* time[i] 
105   
106          delta = time[i] - currenttime 
107          granul = granul * (1.0 - delta / timescale)    \ 
108                   + normal(mu, np.sqrt(delta/timescale)*varscale) 
109          currenttime = time[i] 
110   
111          # Add the different components to the signal.  
112   
113          signal[i] = sum(granul) 
114   
115   
116      # That's it! 
117   
118      return(signal) 
119