Source Code for Module ivs.asteroseismology.solarosc

  1  """ 
  2  Simulation of a solarlike oscillations 
  3   
  4  Author: Joris De Ridder 
  5   
  6  Error messages are written to the logger "solarosc". 
  7   
  8  Notes: 
  9   
 10  Suppose Delta is the peak FWHM linewidth (muHz) and eta is the damping rate (muHz), then:: 
 11   
 12      Delta = eta/pi 
 13  """ 
 14   
 15  import numpy as np 
 16  from numpy.random import uniform, normal 
 17  from math import sin,cos, floor, pi 
 18  import logging 
 19   
 20   
 21  # Setup the logger. 
 22  # Add at least one handler to avoid the message "No handlers could be found"  
 23  # on the console. The NullHandler is part of the standard logging module only  
 24  # from Python 2.7 on. 
 25   
 26 -class NullHandler(logging.Handler): 
 27 -    def emit(self, record): 
 28          pass 
 29           
 30  logger = logging.getLogger("solarosc") 
 31  nullHandler = NullHandler() 
 32  logger.addHandler(nullHandler) 
 33   
 34   
 35   
 36   
 37 -def solarosc(time, freq, ampl, eta): 
 38       
 39      """ 
 40      Compute time series of stochastically excited damped modes 
 41       
 42      See also De Ridder et al., 2006, MNRAS 365, pp. 595-605. 
 43   
 44      Example: 
 45       
 46      >>> time = np.linspace(0, 40, 100)      # in Ms 
 47      >>> freq = np.array([23.0, 23.5])       # in microHz 
 48      >>> ampl = np.array([100.0, 110.0])     # in ppm 
 49      >>> eta = np.array([1.e-6, 3.e-6])      # in 1/Ms 
 50      >>> oscsignal = solarosc(time, freq, ampl, eta) 
 51      >>> flux = 1000000.0                      # average flux level 
 52      >>> signal = flux * (1.0 + oscsignal) 
 53      >>> # The same with a logger 
 54      >>> import sys, logging, logging.handlers 
 55      >>> myLogger = logging.getLogger("solarosc") 
 56      >>> myLogger.addHandler(logging.StreamHandler(sys.stdout)) 
 57      >>> myLogger.setLevel(logging.INFO) 
 58      >>> oscsignal = solarosc(time, freq, ampl, eta, myLogger) 
 59      Simulating 2 modes 
 60      Oscillation kicktimestep: 3333.333333 
 61      300 kicks for warm up for oscillation signal 
 62      Simulating stochastic oscillations 
 63       
 64      @param time: time points [0..Ntime-1] (unit: e.g. Ms) 
 65      @type time: ndarray 
 66      @param freq: oscillation freqs [0..Nmodes-1] (unit: e.g. microHz) 
 67      @type freq: ndarray 
 68      @param ampl: amplitude of each oscillation mode 
 69                   rms amplitude = ampl / sqrt(2.) 
 70      @type ampl: ndarray 
 71      @param eta: damping rates (unit: e.g. (Ms)^{-1}) 
 72      @type eta: ndarray 
 73      @return: signal[0..Ntime-1] 
 74      @rtype: ndarray 
 75      """ 
 76     
 77      Ntime = len(time) 
 78      Nmode = len(freq) 
 79   
 80      logger.info("Simulating %d modes" % Nmode) 
 81   
 82      # Set the kick (= reexcitation) timestep to be one 100th of the 
 83      # shortest damping time. (i.e. kick often enough). 
 84   
 85      kicktimestep = (1.0 / max(eta)) / 100.0 
 86       
 87      logger.info("Oscillation kicktimestep: %f" % kicktimestep) 
 88   
 89      # Init start values of amplitudes, and the kicking amplitude 
 90      # so that the amplitude of the oscillator will be on average be 
 91      # constant and equal to the user given amplitude 
 92   
 93      amplcos = 0.0 
 94      amplsin = 0.0 
 95      kick_amplitude = ampl * np.sqrt(kicktimestep * eta) 
 96   
 97      # Warm up the stochastic excitation simulator to forget the 
 98      # initial conditions. Do this during the longest damping time. 
 99      # But put a maximum on the number of kicks, as there might 
100      # be almost-stable modes with damping time = infinity 
101     
102      damp = np.exp(-eta * kicktimestep) 
103      Nwarmup = min(20000, int(floor(1.0 / min(eta) / kicktimestep))) 
104   
105      logger.info("%d kicks for warm up for oscillation signal" % Nwarmup) 
106           
107      for i in range(Nwarmup): 
108          amplsin = damp * amplsin + normal(np.zeros(Nmode), kick_amplitude) 
109          amplcos = damp * amplcos + normal(np.zeros(Nmode), kick_amplitude) 
110   
111   
112      # Initialize the last kick times for each mode to be randomly chosen 
113      # a little before the first user time point. This is to avoid that 
114      # the kicking time is always exactly the same for all of the modes. 
115   
116      last_kicktime = uniform(time[0] - kicktimestep, time[0], Nmode) 
117      next_kicktime = last_kicktime + kicktimestep 
118   
119   
120      # Start simulating the time series. 
121   
122      logger.info("Simulating stochastic oscillations") 
123           
124      signal = np.zeros(Ntime) 
125   
126      for j in range(Ntime): 
127   
128          # Compute the contribution of each mode separatly 
129   
130          for i in range(Nmode): 
131   
132              # Let the oscillator evolve until right before 'time[j]' 
133   
134              while (next_kicktime[i] <= time[j]): 
135   
136                  deltatime = next_kicktime[i] - last_kicktime[i] 
137                  damp = np.exp(-eta[i] * deltatime) 
138                  amplsin[i] = damp * amplsin[i] + kick_amplitude[i] * normal(0.,1.) 
139                  amplcos[i] = damp * amplcos[i] + kick_amplitude[i] * normal(0.,1.) 
140                  last_kicktime[i] = next_kicktime[i] 
141                  next_kicktime[i] = next_kicktime[i] + kicktimestep 
142   
143              # Now make the last small step until 'time[j]' 
144   
145              deltatime = time[j] - last_kicktime[i] 
146              damp = np.exp(-eta[i] * deltatime) 
147              signal[j] = signal[j] + damp * (amplsin[i] * sin(2*pi*freq[i]*time[j])    \ 
148                                    + amplcos[i] * cos(2*pi*freq[i]*time[j])) 
149   
150      # Return the resulting signal 
151   
152      return(signal) 
153