ComboCode.cc.modeling.profilers.Velocity

1 # -*- coding: utf-8 -*- 2 3 """ 4 Module for calculating the velocity profile as a function of radius. 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 Mdot 17 18 19

20 -def driftRPDF(r,a,l,v,mdot,opac,sd,radiance,T=None,P=0,alpha=0.,mu=2.,\ 21 w_thermal='none'):

22 23 ''' 24 Calculate the drift velocity as a function of radius and grain size by 25 equating the radiation pressure force and drag force. 26 27 If the grain size is given as a constant value, the profile is calculated 28 only for an average grain size, and is essentially 1d, in which case the 29 drift does not depend on grain size. 30 31 For now in the approximation that the thermal velocities are small compared 32 to the drift. This is only valid in the outer, cool regions of the wind. 33 34 The function integrates the luminosity and the opacity over the wavelength 35 grid to estimate the component for the radiation pressure. Note that the 36 opacities are assumed to be independent of grain size for this to work, 37 hence the Rayleigh regime should be applicable. 38 39 Inclusion of the thermal velocity requires T profile. Mean molecular weight 40 assumed to be one in that case. 41 42 @param r: The radial points (cm) 43 @type r: array/float 44 @param a: The grain size grid (cm) 45 @type a: array/float 46 @param l: The wavelength grid (cm) 47 @type l: array/float 48 @param v: The velocity profile (cm/s) 49 @type v: float/Velocity() 50 @param mdot: The mass-loss profile (Msun/yr) 51 @type mdot: float/Mdot() 52 @param opac: The opacity profile (cm2/g), must include l-dependence. Can be 53 extinction, absorption, or scattering, but note that the 54 calculation assumes this is an absorption opacity. 55 @type opac: Opacity() 56 @param sd: The specific density of the dust species 57 @type sd: float 58 @param radiance: The luminosity profile (ergs/s), must include l-dependence 59 @type radiance: Radiance() 60 61 @keyword T: The temperature profile. Only relevant if w_thermal is not 62 none. 63 64 (default: None) 65 @type T: Temperature() 66 @keyword P: The porosity of the grain (or equivalent, to represent 67 effective grain surface). Default is the spherical case. Is 68 given as the ratio of vacuum per unit volume of a grain. 69 70 (default: 0) 71 @type P: float 72 @keyword alpha: Sticking coefficient (Kwok 1975, Van Marle 2011). Default is 73 current value in MCP/GASTRoNOoM. Should be 0.25 according to 74 Kwok and van Marle. 1-alpha gives the fraction of elastic 75 collisions. 76 77 (default: 0.) 78 @type alpha: float 79 @keyword mu: The mean molecular weight. Default is a wind with purely H2. 80 To get the GASTRoNOoM value of 1.4, set fH to 1.5. 81 82 (default: 2.0) 83 @type mu: float 84 @keyword w_thermal: Type of thermal velocity. 85 - sqrt(2kT/m): most probably speed df(v)/dr = 0 86 - sqrt(8/pi kT/m): mean speed (f(v) * v dv)/f(v) dv 87 - sqrt(3kT/m): rms speed sqrt((f(v) * v2 dv)/f(v) dv) 88 - Kwok version: 3/4 sqrt(3kT/m) 89 - Epstein: 4/3 sqrt(8kT/mpi) 90 Choose from: epstein, kwok, prob, mean, rms, none 91 92 (default: none) 93 @type w_thermal: str 94 95 @return: The drift velocity (cm/s) as a function of r and a 96 @rtype: array/float 97 98 ''' 99 100 #-- Define some constants 101 c = cst.c.cgs.value 102 m_p = cst.m_p.cgs.value 103 k_b = cst.k_B.cgs.value 104 w_thermal = w_thermal.lower() 105 106 #-- Check whether v and mdot are constants 107 vi = v.eval(r) if isinstance(v,Velocity) else v 108 mdoti = mdot.eval(r) if isinstance(mdot,Mdot.Mdot) else mdot 109 mdoti = (mdoti*u.Msun/u.yr).to(u.g/u.s).value 110 111 #-- Integrate the opacity and luminosity profiles over wavelength 112 # For this: calculate the emitting surface of the central source 113 L = radiance.getLuminosity(l=l,ftype='flambda') 114 #Ltot = trapz(x=l,y=L) 115 #Lsun = 3.846e+33 116 #print 'The total stellar luminosity is %f Lsun.'%(Ltot/Lsun) 117 OpacL = trapz(x=l,y=opac.eval(l)*L)#+scat.eval(l)*(1-g)*L) 118 119 #-- Calculate the drift for each grain size 120 # 1) create the 2d array with r-dependent and a-dependent 1d arrays 121 # 2) Then add in anything that's constant 122 # Note that Q(a) = kappa*4/3*a*sd*(1-P)^(2/3) 123 arr = np.outer(vi/mdoti,a) 124 vK = np.sqrt(arr/(c*(1-alpha))*OpacL*4./3.*sd*(1.-P)**(2./3.)) 125 126 #-- Calculate the thermal term (see Decin 2006) or return vK (default) 127 if not w_thermal in ['kwok','mean','rms','prob','epstein']: 128 return vK 129 130 #-- RMS, ie sqrt(Int(f(v) v^2 dv)/Int(f(v) dv)) 131 if w_thermal == 'rms': 132 vT = (3.*k_b*T.eval(r)/(mu*m_p))**0.5 133 134 #-- Most probable speed: d(f(v)/dv = 0) 135 elif w_thermal == 'prob': 136 vT = (2.*k_b*T.eval(r)/(mu*m_p))**0.5 137 138 #-- Mean, ie Int(f(v) v dv)/Int(f(v) dv) ==> likely the most appropriate 139 elif w_thermal == 'mean': 140 vT = (8.*k_b*T.eval(r)/(mu*m_p*np.pi))**0.5 141 142 #-- Implementation by Epstein 1924, with the correct factor 4./3. 143 elif w_thermal == 'epstein': 144 vT = 4/3.*(8.*k_b*T.eval(r)/(mu*m_p*np.pi))**0.5 145 146 #-- kwok: Unknown factor 3/4 and rms velocity 147 else: 148 vT = 0.75*(3.*k_b*T.eval(r)/(mu*m_p))**0.5 149 150 #-- Make sure the right axis of vT is multiplied with vK^-1 151 vT = np.transpose([vT]) 152 xT = 0.5*(vT/vK)**2. 153 factor = ((1.+xT**2)**0.5-xT)**0.5 154 155 return factor*vK

156 157 158

159 -def vsound(T,gamma,mu):

160 161 ''' 162 The sound velocity at T for a given mean molecular weight and gamma. 163 164 @param T: The temperature (K) 165 @type T: array/float 166 @param gamma: The adiabatic constant 167 @type gamma: float 168 @param mu: the mean molecular weight in units of hydrogen mass (g) 169 @type mu: float 170 171 @return: The sound velocity 172 @rtype: array/float 173 174 ''' 175 176 #-- Define constants explicitly for quick calculation 177 k_B = 1.3806488e-16 # astropy.constants.k_B.cgs.value 178 m_p = 1.672621777e-24 # astropy.constants.m_p.cgs.value 179 180 return np.sqrt(gamma*k_B*T/mu/m_p)

181 182 183

184 -def vbeta2D(r,a,*args,**kwargs):

185 186 ''' 187 Convenience function for vbeta in 2D. Used by Drift() to calculate a beta 188 law for a single grain size. 189 190 Ignores the grain size array, and passes arguments on to vbeta. 191 192 Additional args and kwargs are passed on to vbeta. 193 194 @param r: The radial points (cm) 195 @type r: array 196 @param a: Grain size (cm), dummy variable that is ignored. 197 @type a: array 198 199 @return: The beta law velocity profile (cm/s) 200 @rtype: array 201 202 ''' 203 204 return vbeta(r,*args,**kwargs)

205 206 207

208 -def vbeta(r,r0,v0,vinf,beta,vi_mode='constant',rstar=0.0,beta_inner=0.5,\ 209 *args,**kwargs):

210 211 ''' 212 The beta-power law description of the velocity profile. 213 214 Below the cut-off r0, the velocity is taken constant and equal to v0. This 215 behaviour can be changed with the vi_mode keyword. 216 217 @param r: The radial points (cm) 218 @type r: array 219 @param r0: The inner radius (typically at dust condensation radius, cm) 220 @type r0: float 221 @param v0: The initial velocity at r0 (cm/s). Can be 'vs' for the sound 222 velocity, in which case additional args and kwargs must be passed 223 to the function call (see method vsound()). 224 @type v0: float 225 @param vinf: The terminal wind velocity (cm/s) 226 @type vinf: float 227 @param beta: The exponent of the beta law 228 @type beta: float 229 230 @keyword vi_mode: The type of velocity law used in the inner wind at r<r0. 231 Options: 232 - constant: A constant value at v0 233 - beta: an increasing velocity up to v0 according to a 234 beta law with beta=beta_inner, r0~rstar. 235 Lower case letters required. 236 237 (default: 'constant') 238 @type vi_mode: str 239 @keyword rstar: The stellar radius (cm), only relevant when 240 vi_mode==beta 241 242 (default: 0) 243 @type rstar: float 244 @keyword beta_inner: Only relevant for vi_mode==beta, the beta exponent for 245 the inner wind at r<r0. Default is value for GASTRonOoM 246 247 (default: 0.5) 248 @type beta_inner: float 249 250 @return: The beta law velocity profile (cm/s) 251 @rtype: array 252 253 ''' 254 255 #-- In case sound velocity is requested, calculate it 256 if v0 == 'vs': 257 v0 = vsound(*args,**kwargs) 258 259 #-- Set the velocity law 260 v = v0 + (vinf-v0)*(1-r0/r)**beta 261 262 #-- Set the inner wind velocity 263 if vi_mode == 'constant': 264 vi = v0 265 266 #-- To reproduce the velocity law of GASTRoNOoM exactly, use v0 = sound 267 # velocity at r0. See source_cooling/vellaw.f: 268 # R0 = R_STAR*(1d0 - (0.01d0**(1d0/0.5d0))) 269 # vr = VEL_SOUND * (1d0-(R0/R_INNER))**(-0.5d0) * 270 # $ (1d0-(R0/rad1))**(0.5d0) 271 #-- Can be used with any v0, for a smooth transition r<r0 => r>r0 272 elif vi_mode == 'beta': 273 #-- Because starting from 0 km/s at rstar can lead to numerical issues, 274 # instead scale the rstar such that we start from 0.01*v0, maintaining 275 # 0 km/s at rstar itself. 276 rstar_scaled = rstar*(1 - 0.01**(1./beta_inner)) 277 278 #-- To have a properly scaled inner-wind beta law, with a v0 at r0 for 279 # beta == 0.5, one must parametrise vinf_inner and r0_inner for this 280 # law (whereby r0_inner is actually rstar_scaled), assuming 281 # v0_inner == 0, and filling in values in the beta law. 282 # v0 = vinf_inner (1 - r0_inner/r0)**0.5 283 # leading to: vinf_inner = v0 * (1-r0_inner/r0)**-0.5 284 v0_factor = (1.-rstar_scaled/r0)**(-beta_inner) 285 vi = v0*v0_factor*(1-rstar_scaled/r)**beta_inner 286 287 #-- Return the full velocity law, with the inner wind velocity of choice. 288 return np.where(r<=r0,vi,v)

289 290 291

292 -def dvbetadr(r,r0,v0,vinf,beta,vi_mode='constant',rstar=0.0,beta_inner=0.5,\ 293 *args,**kwargs):

294 295 ''' 296 The derivative of the vbeta function: analytic. 297 298 Below the cut-off r0, the velocity is taken constant and equal to v0. This 299 behaviour can be changed with the vi_mode keyword. (NYI) 300 301 @param r: The radial points (cm) 302 @type r: array 303 @param r0: The inner radius (typically at dust condensation radius, cm) 304 @type r0: float 305 @param v0: The initial velocity at r0 (cm/s). Can be 'vs' for the sound 306 velocity, in which case additional args and kwargs must be passed 307 to the function call (see method vsound()). 308 @type v0: float 309 @param vinf: The terminal wind velocity (cm/s) 310 @type vinf: float 311 @param beta: The exponent of the beta law 312 @type beta: float 313 314 @keyword vi_mode: The type of velocity law used in the inner wind at r<r0. 315 Options: 316 - constant: A constant value at v0 317 - beta: an increasing velocity up to v0 according to a 318 beta law with beta=beta_inner, r0~rstar. 319 Lower case letters required. 320 321 (default: 'constant') 322 @type vi_mode: str 323 @keyword rstar: The stellar radius (cm), only relevant when 324 vi_mode==beta 325 326 (default: 0) 327 @type rstar: float 328 @keyword beta_inner: Only relevant for vi_mode==beta, the beta exponent for 329 the inner wind at r<r0. Default is value for GASTRonOoM 330 331 (default: 0.5) 332 @type beta_inner: float 333 334 @return: The analytic derivative of the beta law velocity profile (cm^2/s) 335 @rtype: array 336 337 ''' 338 339 #-- In case sound velocity is requested, calculate it 340 if v0 == 'vs': 341 v0 = vsound(*args,**kwargs) 342 343 #-- Set the derived velocity law 344 dvdr = beta*(vinf-v0)*(1.-r0/r)**(beta-1)*r0/(r*r) 345 346 #-- Set the inner wind velocity 347 if vi_mode == 'constant': 348 dvi = 0.0 349 350 #-- To reproduce the velocity law of GASTRoNOoM exactly, use v0 = sound 351 # velocity at r0. See source_cooling/vellaw.f and vbeta() 352 elif vi_mode == 'beta': 353 #-- Because starting from 0 km/s at rstar can lead to numerical issues, 354 # instead scale the rstar such that we start from 0.01*v0, maintaining 355 # 0 km/s at rstar itself. 356 rstar_scaled = rstar*(1 - 0.01**(1./beta_inner)) 357 358 #-- To have a properly scaled inner-wind beta law, with a v0 at r0 for 359 # beta == 0.5, one must parametrise vinf_inner and r0_inner for this 360 # law (whereby r0_inner is actually rstar_scaled), assuming 361 # v0_inner == 0, and filling in values in the beta law. 362 # v0 = vinf_inner (1 - r0_inner/r0)**0.5 363 # leading to: vinf_inner = v0 * (1-r0_inner/r0)**-0.5 364 #-- Then calc derivative. 365 v0_factor = (1.-rstar_scaled/r0)**(-beta_inner) 366 deriv_factor = rstar_scaled/(r*r) 367 dvi = v0*v0_factor*0.5*(1-rstar_scaled/r)**(beta_inner-1.)*deriv_factor 368 369 #-- Return the derivative of the full velocity law, with the inner wind 370 # velocity of choice 371 return np.where(r<=r0,dvi,dvdr)

372 373 374

375 -def vdust(r,v,w,*args,**kwargs):

376 377 ''' 378 Calculate the dust velocity profile from a gas velocity and a drift velocity 379 profile. 380 381 Requires averaging the drift velocity over grain size. Several normalisation 382 possibilities available. See Drift().avgDrift. 383 384 @param r: The radial points (cm) 385 @type r: array 386 @param v: The velocity profile (cm/s) 387 @type v: float/Velocity() 388 @param w: The drift velocity profile (cm/s) as a function of r and a. 389 @type w: float/Drift() 390 391 @keyword args: Additional parameters passed to the avgDrift method 392 393 (default: []) 394 @type args: tuple 395 @keyword kwargs: Additional keywords passed to the avgDrift method 396 397 (default: {}) 398 @type kwargs: dict 399 400 @return: The dust velocity profile (cm/s) 401 @rtype: array 402 403 ''' 404 405 vi = v.eval(r) if isinstance(v,Velocity) else v 406 wi = w.avgDrift(r,*args,**kwargs) if isinstance(w,Drift) else w 407 408 return vi + wi

409 410 411

412 -class Velocity(Profiler.Profiler):

413 414 ''' 415 An interface for a velocity profile and its derivative 416 417 ''' 418

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

420 421 ''' 422 Create an instance of the Velocity() class. Requires a radial grid and 423 a velocity function. 424 425 The function can also be given as an interp1d object. 426 427 The optional args and kwargs give the additional arguments for the 428 velocity function, which are ignored in case func is an interp1d object. 429 430 @param r: The radial points (cm) 431 @type r: array 432 @param func: The function that describes the velocity profile. Can be 433 given as an interp1d object. 434 @type func: function 435 436 @keyword dfunc: The function that describes the derivative of the 437 profile with respect to r. Can be given as an interp1d 438 object. If None, a generic central difference is taken 439 and interpolated. 440 441 (default: None) 442 @type dfunc: function/interp1d object 443 @keyword order: Order of the spline interpolation. Default is cubic. 444 445 (default: 3) 446 @type order: int 447 448 @keyword args: Additional parameters passed to the functions when eval 449 or diff are called. 450 451 (default: []) 452 @type args: tuple 453 @keyword kwargs: Additional keywords passed to the functions when eval 454 or diff are called. 455 456 (default: {}) 457 @type kwargs: dict 458 ''' 459 460 #-- In case of vbeta, use the analytic function for the derivative 461 if func == vbeta or func == 'vbeta': 462 dfunc = dvbetadr 463 464 #-- Do not name func, dfunc, etc in function call, or *args won't work 465 super(Velocity, self).__init__(r,func,dfunc,order,*args,**kwargs) 466 467 self.r = self.x 468 self.v = self.y 469 self.dvdr = self.dydx

470 471 472

473 -class Drift(Profiler.Profiler2D):

474 475 ''' 476 An interface for the drift velocity as a function of radius and grain size. 477 478 ''' 479

480 - def __init__(self,r,a,func,*args,**kwargs):

481 482 ''' 483 Create an instance of the Drift() class. Requires a radial grid, a grain 484 size grid and a velocity function. 485 486 The function can also be given as an interpolator object. 487 488 The optional args and kwargs give the additional arguments for the 489 drift function, which are ignored in case func is an interpolator object 490 491 @param r: The radial points (cm) 492 @type r: array 493 @param a: The grain size grid (cm) 494 @type a: array 495 @param func: The function that describes the drift profile. Can be 496 given as an interpolator object. 497 @type func: function 498 499 @keyword args: Additional parameters passed to the functions when eval 500 or diff are called. 501 502 (default: []) 503 @type args: tuple 504 @keyword kwargs: Additional keywords passed to the functions when eval 505 or diff are called. 506 507 (default: {}) 508 @type kwargs: dict 509 510 ''' 511 512 super(Drift, self).__init__(r,a,func,*args,**kwargs) 513 self.r = self.x 514 self.a = self.y 515 self.w = self.z

516 517 518

519 - def avgDrift(self,r=None,norm_type='standard',nd=None,warn=1):

520 521 ''' 522 Calculate the average drift velocity as a function of radius, 523 normalising with respect to the average grain size (requires the 524 grain size distribution). 525 526 @keyword r: The radial points (cm). If None, the default r grid is used. 527 528 (default: None) 529 @type r: float/array 530 @keyword norm_type: The weighting of normalisation for the average drift 531 velocity: 532 - collisional: Int(w*nd*a^2*da)/Int(nd*a^2*da) 533 - dens: Int(w*nd*a^3*da)/Int(nd*a^3*da) 534 - nd: Int(w*nd*da)/Int(nd*da) 535 - a: Int(w*a*da)/Int(a*da) 536 - standard: Int(w*da)/Int(da) 537 538 (default: 'standard') 539 @type norm_type: str 540 @keyword nd: The grain size distribution. Only needed when norm_type is 541 dens or nd. 542 543 (default: None) 544 @type nd: Distribution() 545 @keyword warn: Warn when extrapolation occurs. 546 547 (default: 1) 548 @type warn: bool 549 550 @return: The average drift velocity evaluated at r (cm/s). Array is 551 reduced by one dimension since a drops out as variable. 552 @rtype: array/float 553 554 ''' 555 556 if self.a.size == 1: 557 return np.squeeze(self.eval(r,warn=warn)) 558 559 norm_type = norm_type.lower() 560 #-- Normalise over grain size 561 if norm_type == 'a': 562 wsum = trapz(y=self.eval(x=r,warn=warn)*self.a,x=self.a,axis=1) 563 norm = trapz(y=self.a*self.a,x=self.a) 564 565 #-- Normalize over grain surface * nd (collisional drift) 566 elif norm_type == 'collisional': 567 #-- Constant pi drops out due to normalisation 568 cs_tot = nd.eval(x=r,y=self.a,warn=warn)*self.a**2 569 wsum = trapz(y=self.eval(x=r,warn=warn)*cs_tot,x=self.a,axis=1) 570 norm = trapz(y=cs_tot,x=self.a,axis=1) 571 572 #-- Normalise over number density 573 elif norm_type == 'nd': 574 ndens = nd.eval(r,self.a,warn=warn) 575 wsum = trapz(y=self.eval(x=r,warn=warn)*ndens,x=self.a,axis=1) 576 norm = trapz(y=ndens,x=self.a,axis=1) 577 578 #-- Normalise over dust density 579 elif norm_type == 'dens': 580 #-- Constant 4/3*pi*spec_dens drops out due to normalisation. 581 rho = nd.eval(r,self.a,warn=warn)*self.a**3 582 wsum = trapz(y=self.eval(x=r,warn=warn)*rho,x=self.a,axis=1) 583 norm = trapz(y=rho,x=self.a,axis=1) 584 585 #-- Default is 'standard', so if neither a or nd or dens, do standard. 586 else: 587 wsum = trapz(y=self.eval(x=r,warn=warn),x=self.a,axis=1) 588 norm = trapz(y=np.ones_like(self.a),x=self.a) 589 590 return wsum/norm

591

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