ComboCode.cc.modeling.profilers.Opacity

1 # -*- coding: utf-8 -*- 2 3 """ 4 Module for calculating the opacity. 5 6 Author: R. Lombaert 7 8 """ 9 10 import collections 11 import numpy as np 12 from astropy import constants as cst 13 from scipy.interpolate import InterpolatedUnivariateSpline as spline1d 14 15 from cc.modeling.profilers import Profiler 16 from cc.modeling.profilers import Density 17 from cc.tools.readers import KappaReader as KR 18 from cc.tools.numerical import Operators as op 19 from cc.data import Data 20 21

22 -def read_opacity(x,species,index=1,unit='cm',*args,**kwargs):

23 24 ''' 25 Read the opacities with the KappaReader and interpolate. 26 27 Returns a tuple with the original grid, and the interpolation object. Used 28 by the Opacity class to set an interpolation object for a .opac/.particle 29 file through species given in usr/Dust.dat. 30 31 Additional args and kwargs can be passed to the interpolate method. 32 33 @param x: The x grid requested for the interpolation. Note that this is a 34 dummy variable to allow Profiler to work with this function. 35 @type x: array 36 @param species: The dust species (from Dust.dat) 37 @type species: string 38 39 @keyword index: The index of the kappas in the .opacity/.particle file. 40 0: extinction, 1: absorption, 2: scattering 41 42 (default: 1) 43 @type index: int 44 @keyword unit: The unit of the wavelength. Can be given as u.Unit() 45 object or as a string representation of those objects. 46 Can range from length, to frequency, and energy 47 48 (default: cm) 49 @type unit: str/u.Unit() 50 51 @return: The interpolator for the mass extinction/absorption/scattering 52 coefficients. (in cgs!) 53 @rtype: spline1d 54 55 ''' 56 57 kr = KR.KappaReader() 58 return (kr.getWavelength(species,unit=unit),kr.getKappas(species,index), 59 kr.interpolate(species,index,unit=unit,*args,**kwargs)) 60 61 62

63 -class Opacity(Profiler.Profiler):

64 65 ''' 66 An interface for an opacity profile and any type of extinction quantity. 67 68 ''' 69

70 - def __init__(self,l,func,dfunc=None,order=3,*args,**kwargs):

71 72 ''' 73 Create an instance of the Opacity() class. Requires a wavelength grid. 74 75 The function can also be given as an interpolation object. 76 77 The optional args and kwargs give the additional arguments for the 78 opacity function, which are ignored in case func is an interpolation 79 object. 80 81 Eval and diff work with the mass extinction coefficient in cm2/g. 82 83 In case func refers to an interpolation object in KappaReader, the 84 args/kwargs should always contain an index indicating whether extinction 85 (default), absorption or scattering is requested: 86 - 0: extinction, 87 - 1: absorption, 88 - 2: scattering 89 90 Note that if func is an interpolator object, the original input x and y 91 grids can be passed as additional keywords xin/lin and yin, which would 92 then be arrays. Otherwise, the x and the interpolator(x) are set as 93 xin/lin and yin. xin/lin and yin are ignored if func is a function, even 94 if it returns an interpolator (in which case the original grids are 95 known) 96 97 Extrapolation is done outside of the "original" (depending on if xin/yin 98 are given) wavelength ranges according to a power law ~l^-alpha, with 99 alpha given upon eval() calls. 100 101 @param l: The wavelength points (cm) 102 @type l: array 103 @param func: The function that describes the opacity profile. Can be 104 given as an interpolation object, in which case lin and yin 105 keywords can be passed as arrays for the original grids. 106 @type func: function 107 108 @keyword dfunc: Function that describes the derivative of the profile 109 with respect to x. Can be given as an interpolation 110 object. If None, a generic central difference is taken & 111 interpolated with a spline of which the order can be 112 chosen. 113 114 (default: None) 115 @type dfunc: function/interpolation object 116 @keyword order: Order of the spline interpolation. Default is cubic. 117 118 (default: 3) 119 @type order: int 120 121 @keyword args: Additional parameters passed to the functions when eval 122 or diff are called. 123 124 (default: []) 125 @type args: tuple 126 @keyword kwargs: Additional keywords passed to the functions when eval 127 or diff are called. 128 129 (default: {}) 130 @type kwargs: dict 131 132 ''' 133 134 #-- Make sure lin is recognized 135 if kwargs.has_key('lin'): 136 kwargs['xin'] = kwargs['lin'] 137 del kwargs['lin'] 138 139 #-- Do not name func, dfunc, etc in function call, or *args won't work 140 super(Opacity, self).__init__(l,func,dfunc,order,*args,**kwargs) 141 142 143 #-- Remember the index of the opacities file. If relevant. 144 if kwargs.has_key('index'): 145 self.index = kwargs['index'] 146 else: 147 self.index = 0 148 149 #-- Set the default alpha. Not used if func is not an interpolator. 150 self.alpha = 2. 151 152 #-- If the requested function is not an interpolator of a file, do not 153 # change anything. 154 if not self.interp_func: 155 self.l = self.x 156 self.kappa = self.y 157 return 158 159 #-- For Opacity(), in case of an interpolator object as function, 160 # the extrapolation must be re-done for it to be safe. 161 # The eval() method is rewritten here to function with a 162 # power law extrapolation, for which the settings can be chosen. 163 # The diff method must not be rewritten since it is based on the y 164 # profile. 165 self.l = None 166 self.kappa = None 167 self.y = None 168 169 #-- Run the class' own eval() method, with l=self.x, which will not be 170 # equal to the class' self.l attribute, and thus it will evaluated 171 # properly. Evaluated with alpha==2. Reevaluated if alpha is different 172 self.kappa = self.eval(self.x,alpha=2.) 173 self.y = self.kappa 174 175 #-- Now also redo the derivative method, in case it is an interpolator 176 if self.interp_dfunc: 177 self.dfunc = spline1d(x=self.x,y=op.diff_central(self.y,self.x),\ 178 k=self.order) 179 180 #-- Evaluate the derivative with the default grid 181 self.dydx = self.dfunc(self.x,*self._dargs,**self._dkwargs) 182 183 #-- Now set l, dydl. 184 self.l = self.x 185 self.dydl = self.dydx

186 187 188

189 - def __call__(self,l=None,*args,**kwargs):

190 191 ''' 192 Evaluate the profile function at a coordinate point. 193 194 l can be any value or array. No warning is printed in case of 195 extrapolation, which is done with a power law at small or large l values 196 by the eval() function. 197 198 The default y-grid is returned if l is None. 199 200 Passes the call to eval. Additional keywords for the power law 201 extrapolation can be passed this way as well. 202 203 @keyword l: The primary coordinate point(s). If None, the default 204 coordinate grid is used. 205 206 (default: None) 207 @type l: array/float 208 209 @return: The profile evaluated at l 210 @rtype: array/float 211 212 ''' 213 214 return self.eval(l,*args,**kwargs)

215 216 217

218 - def eval(self,l=None,alpha=2.):

219 220 ''' 221 Evaluate the profile function at a coordinate point. 222 223 l can be any value or array. No warning is printed in case of 224 extrapolation, which is done with a power law at small or large l values 225 226 The default y-grid is returned if l is None and alpha is 2 (the default) 227 228 @keyword l: The coordinate point(s). If None, the default 229 coordinate grid is used. 230 231 (default: None) 232 @type l: array/float 233 @keyword alpha: The exponent of the wavelength-dependent power law 234 extrapolation, such that kappa ~ lambda^-alpha 235 236 (default: 2.) 237 @type alpha: float 238 239 @return: The profile evaluated at l 240 @rtype: array/float 241 242 ''' 243 244 #-- Return self.y since l was given as None 245 if l is None and alpha == self.alpha: 246 return self.y 247 248 #-- l can still be None. Apparently different alpha requested. So calc. 249 if l is None: 250 l = self.l 251 252 #-- First retrieve the original profile. Don't warn since we re-do the 253 # extrapolation anyway. If this is the first call from the __init__ 254 # method of Opacity(), this will be the standard grid. 255 y = super(Opacity,self).eval(l,warn=0) 256 257 #-- If self.func is not an interpolator, no warnings needed, and no 258 # extrapolation needed either. 259 if not self.interp_func: return y 260 261 #-- Determine the regions where extrapolation is done, i.e. outside the 262 # original l-grid's range. 263 lmin, lmax = self.xin[0], self.xin[-1] 264 ymin, ymax = self.yin[0], self.yin[-1] 265 266 #-- Replace the extrapolated values with the new power law. Make sure l 267 # and y are an array for this. 268 larr, y = Data.arrayify(l), Data.arrayify(y) 269 y[larr<lmin] = ymin*(larr[larr<lmin]/lmin)**(-alpha) 270 y[larr>lmax] = ymax*(larr[larr>lmax]/lmax)**(-alpha) 271 272 return y if isinstance(l,collections.Iterable) else y[0]

273 274 275

276 - def diff(self,l=None,alpha=2.):

277 278 ''' 279 Evaluate the derivative of the profile function at a coordinate point. 280 281 l can be any value or array. No warning is printed in case of 282 extrapolation, which is done with a power law at small or large l values 283 284 The default y-grid is returned if l is None and alpha is 2 (the default) 285 286 @keyword l: The coordinate point(s). If None, the default 287 coordinate grid is used. 288 289 (default: None) 290 @type l: array/float 291 @keyword alpha: The exponent of the wavelength-dependent power law 292 extrapolation, such that kappa ~ lambda^-alpha 293 294 (default: 2.) 295 @type alpha: float 296 297 @return: The derivative evaluated at l 298 @rtype: array/float 299 300 ''' 301 302 #-- Return self.y since l was given as None 303 if l is None and alpha == self.alpha: 304 return self.dydx 305 306 #-- l can still be None. Apparently different alpha requested. So calc. 307 if l is None: 308 l = self.l 309 310 #-- First retrieve the original profile. Don't warn since we re-do the 311 # extrapolation anyway. This will have the wrong alpha if alpha is not 312 # 2, but we re-do the extrapolation anyway. The original profile will 313 # be untouched. 314 dydl = super(Opacity,self).diff(l,warn=0) 315 316 #-- If self.func or self.dfunc is not an interpolator, no warnings 317 # needed, and no extrapolation needed either. 318 if not (self.interp_dfunc and self.interp_func): return dydl 319 320 #-- Determine the regions where extrapolation is done, i.e. outside the 321 # original l-grid's range. 322 lmin, lmax = self.xin[0], self.xin[-1] 323 ymin, ymax = self.yin[0], self.yin[-1] 324 325 #-- Replace the extrapolated values with the new power law. Make sure l 326 # and y are an array for this. 327 larr, dydl = Data.arrayify(l), Data.arrayify(dydl) 328 dydl[larr<lmin] = ymin*(larr[larr<lmin]/lmin)**(-alpha-1.)*-alpha/lmin 329 dydl[larr>lmax] = ymax*(larr[larr>lmax]/lmax)**(-alpha-1.)*-alpha/lmax 330 331 return dydl if isinstance(l,collections.Iterable) else dydl[0]

332 333 334

335 - def getEfficiency(self,l,a,sd,P=0.,alpha=2.):

336 337 ''' 338 Calculate the extinction/absorption/scattering (depending on index) 339 efficiency for given grain size and porosity on a wavelength grid. 340 341 Follows: Q = (extinction cross section) / (geometric cross section) 342 where geometric cross section is the effective surface, hence 343 pi*a^2/(1-P)^(2/3), taking into account the porosity of the grain. 344 345 @param l: The wavelength points (cm) 346 @type l: array 347 @param a: The grain size (cm) 348 @type a: float 349 @param sd: The specific density of the dust species (g/cm3) 350 @type sd: float 351 @keyword P: The porosity of the grain (or equivalent, to represent 352 effective grain surface). Default is the spherical case. Is 353 given as the ratio of vacuum per unit volume of a grain. 354 355 (default: 0) 356 @type P: float 357 @keyword alpha: The exponent of the wavelength-dependent power law 358 extrapolation, such that kappa ~ lambda^-alpha 359 360 (default: 2.) 361 @type alpha: float 362 363 @return: The dimensionless extinction efficiencies 364 @rtype: array 365 366 ''' 367 368 kappas = self.eval(l,alpha=alpha) 369 return kappas * 4./3. * sd * a * (1-P)**(2./3.)

370

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