Source Code for Module ComboCode.cc.modeling.profilers.Grainsize

  1  # -*- coding: utf-8 -*- 
  2   
  3  """ 
  4  Module for calculating the grain size distribution.  
  5   
  6  Author: R. Lombaert 
  7   
  8  """ 
  9   
 10  import numpy as np 
 11  from astropy import constants as cst 
 12  from astropy import units as u 
 13  from scipy.integrate import trapz 
 14   
 15  from cc.modeling.profilers import Profiler 
 16  from cc.modeling.profilers import Velocity, Mdot 
 17   
 18   
 19 -def MRN(r,a,w,v,mdot_dust,sd,a_min=0.005e-4,a_max=0.25e-4,returnA=0): 
 20   
 21      ''' 
 22      Calculate the MRN grain size distribution for an AGB wind.  
 23       
 24      This assumes an average drift velocity independent of grain size, and a  
 25      constant dust-to-gas ratio throughout the wind. 
 26       
 27      The MRN distribution is taken from Mathis, Rumpl and Nordsieck (1977), and  
 28      reads:  
 29      n_d(a,r) = A(r) n_Htot(r) a^-3.5 
 30       
 31      A(r) is a scaling factor that is derived from the dust mass-loss rate. This  
 32      causes the H-number density to drop out of the equation, so this information 
 33      is not needed: 
 34      A(r) = md / [nHtot 16/3 pi^2 sd r^2 Int((vg + w) a^-0.5 da)] 
 35       
 36      The MRN distribution assumes minimum grain size of .005 micron and maximum 
 37      grain size of 0.25 micron. While the method allows going outside this range 
 38      (e.g. such as in the case of calculating a local value just outside the  
 39      range when solving a differential equation), it is recommended changing them 
 40      if the grain size distribution extends beyond it. These are the values used 
 41      to derive A(r). 
 42       
 43      @param r: The radial grid (cm) 
 44      @type r: array 
 45      @param a: The grain size grid (cm) 
 46      @type a: array 
 47      @param w: The drift velocity profile (cm/s) as a function of r and a. 
 48      @type w: Drift() 
 49      @param v: The gas velocity (cm/s). Either as a profile, or as a cst 
 50      @type v: float/Velocity() 
 51      @param mdot_dust: The dust mass-loss rate (msun/yr). Either as a profile, or 
 52                        as a cst 
 53      @type mdot_dust: float/Mdot() 
 54      @param sd: The average specific density of the dust grains in g/cm3. 
 55      @type sd: float 
 56       
 57      @keyword a_min: The minimum grain size in cm for MRN 
 58       
 59                      (default: 0.005e-4) 
 60      @type a_min: float 
 61      @keyword a_max: The maximum grain size in cm for MRN 
 62       
 63                      (default: 0.25e-4) 
 64      @type a_max: float 
 65      @keyword returnA: Return the scale factor separately as well. This is the  
 66                        second element of the tuple. Note that this does not  
 67                        contain Htot, ie A = A_local/nhtot 
 68                         
 69                        (default: 0) 
 70      @type returnA: bool 
 71       
 72      @return: The grain size distribution as a function of a 
 73      @rtype: array 
 74       
 75      ''' 
 76       
 77      #-- Check if the mass-loss rate is given as a constant or a profile 
 78      mddi = mdot_dust.eval(r) if isinstance(mdot_dust,Mdot.Mdot) else mdot_dust 
 79      mddi = (mddi*u.Msun/u.yr).to(u.g/u.s).value 
 80      vi = v.eval(r) if isinstance(v,Velocity.Velocity) else v 
 81       
 82      #-- Calculate the integration over a for A. Use the drift's default a grid, 
 83      #   because A(r) does not depend on the specific a requested here. 
 84      term1 = 2*(np.sqrt(a_max)-np.sqrt(a_min))*vi 
 85      term2 = trapz(x=w.a,y=w.eval(x=r)*w.a**-0.5) 
 86   
 87      #-- Calculate the factors for A 
 88      denominator = 16*np.pi**2*sd*r**2*(term1+term2) 
 89      numerator = mddi*3. 
 90      A = numerator/denominator 
 91       
 92      #-- Calculate the distribution 
 93      #   Note that we do not need the TOTAL hydrogen density, since it drops out  
 94      #   through A where nHtot ends up in the denominator 
 95      nd = np.outer(A,a**-3.5) 
 96       
 97      if returnA: return (nd,A) 
 98       
 99      return nd 
100       
101       
102       
103 -class Distribution(Profiler.Profiler2D):  
104       
105      ''' 
106      An interface for a grain size distribution as a function of grain size and 
107      radius, and for given minimum and maximum grain sizes. 
108       
109      ''' 
110       
111 -    def __init__(self,r,a,func,*args,**kwargs): 
112       
113          ''' 
114          Create an instance of the Distribution() class. Requires a radial 
115          and a grain size grid. 
116           
117          The function can also be given as an interpolation object. 
118           
119          The optional args and kwargs give the additional arguments for the  
120          temperature function, which are ignored in case func is an interpolation  
121          object. 
122           
123          @param r: The radial points (cm) 
124          @type r: array 
125          @param a: The grain size points (cm). Must include minimum and maximum 
126                    grain sizes, unless manually declared in kwargs. 
127          @type a: array 
128          @param func: The function that describes the grain size distribution.  
129                       Can be given as an interpolation object. 
130          @type func: function/interpolation object 
131                   
132          @keyword args: Additional parameters passed to the functions when eval 
133                         or diff are called.  
134                          
135                         (default: []) 
136          @type args: tuple 
137          @keyword kwargs: Additional keywords passed to the functions when eval 
138                           or diff are called.  
139                          
140                           (default: {}) 
141          @type kwargs: dict 
142           
143          ''' 
144   
145          #-- Set min and max grain size if needed 
146          if not kwargs.has_key('a_min'): 
147              kwargs['a_min'] = a[0] 
148          if not kwargs.has_key('a_max'): 
149              kwargs['a_max'] = a[-1] 
150               
151          #-- Do not name func, dfunc, etc in function call, or *args won't work             
152          super(Distribution, self).__init__(r,a,func,*args,**kwargs) 
153           
154          self.r = self.x 
155          self.a = self.y 
156          self.nd = self.z 
157           
158           
159           
160 -    def __call__(self,*args,**kwargs): 
161       
162          ''' 
163          Evaluate the profile function at a coordinate point. 
164                   
165          Passes all args and kwargs on the call to self.eval(). 
166           
167          @return: The profile evaluated at x and y 
168          @rtype: array/float 
169           
170          ''' 
171       
172          return self.eval(*args,**kwargs) 
173       
174       
175       
176 -    def eval(self,x=None,y=None,warn=1,w_sputter=0.,w=np.empty(0)): 
177           
178          ''' 
179          Evaluate the grain-size distribution function at a coordinate point. 
180           
181          x/y can be any value or array. If func is an interpolation object, it is 
182          in principle limited by the x/y-range of the interpolator. It is advised 
183          not to extend much beyond the given x/y-range.  
184           
185          If one of the two variables is None, it is replaced by the default grid. 
186           
187          Sputtering can be applied here by passing the drift array and the max 
188          drift velocity to the call. Any grain sizes with a drift above w_sputter 
189          will have their number density set to 0. This depends on both the grain  
190          size and the radius. 
191                   
192          @keyword x: The primary coordinate point(s). If None, the default  
193                      coordinate grid is used. 
194           
195                      (default: None) 
196          @type x: array/float 
197          @keyword y: The secondary coordinate point(s). If None, the default  
198                      coordinate grid is used. 
199           
200                      (default: None) 
201          @type y: array/float 
202          @keyword warn: Warn when extrapolation occurs. 
203           
204                         (default: 1) 
205          @type warn: bool 
206          @keyword w_sputter: The sputtering drift velocity above which the number 
207                              density is set to zero. Default in case no  
208                              sputtering is to be applied.  
209                               
210                              (default: 0.) 
211          @type w_sputter: float 
212          @keyword w: The Drift() object. Only relevant when w_sputter is nonzero. 
213           
214                      (default: np.empty(0)) 
215          @type w: array 
216           
217          @return: The profile evaluated at x and y 
218          @rtype: array/float 
219           
220          ''' 
221           
222          #-- Grab the original grain size distribution. 
223          z = super(Distribution,self).eval(x=x,y=y,warn=warn) 
224           
225          #-- w_sputter is zero, don't change anything 
226          if w_sputter == 0.: return z 
227           
228          #-- Set the distribution to zero wherever the drift velocity is too  
229          #   large, and return 
230          z = np.where(w>w_sputter,np.zeros_like(z),z) 
231          return z 
232