ComboCode.cc.modeling.physics.EnergyBalance

1 # -*- coding: utf-8 -*- 2 3 """ 4 Module for calculating the energy balance. 5 6 Author: R. Lombaert, H. Olofsson & M. Maercker (Chalmers, Sweden) 7 8 For the use of the EnergyBalance module, please contact one of the three authors 9 listed here. 10 11 """ 12 13 import os, collections, functools, copy 14 import numpy as np 15 from scipy.interpolate import InterpolatedUnivariateSpline as spline1d 16 from scipy.interpolate import interp1d 17 from scipy.integrate import odeint, trapz, cumtrapz 18 from astropy import constants as cst 19 from astropy import units as u 20 21 import matplotlib.pyplot as plt 22 23 from cc.data import Data 24 from cc.plotting import Plotting2 25 from cc.tools.io import DataIO 26 from cc.tools.numerical import Operators as op 27 from cc.tools.numerical import Gridding 28 from cc.modeling.profilers import Velocity, Density, Profiler, Grainsize 29 from cc.modeling.profilers import Mdot, Radiance, Opacity, Temperature 30 from cc.tools.readers import CollisReader, PopReader, LamdaReader, MlineReader 31 from cc.tools.readers import RadiatReader 32 import cc.path 33 34 #-- Define a few global constants. Calling cst inside the functions slows them 35 # down too much 36 37 #-- The Boltzmann constant 38 global k_b 39 k_b = cst.k_B.cgs.value # erg/K 40 41 #-- Hydrogen mass 42 global mh 43 mh = cst.m_p.cgs.value # g 44 45 46

47 -def dTdr(T,r,v,gamma,rates=None,warn=1):

48 49 ''' 50 The differential equation for the kinetic temperature profile. 51 52 Currently only takes into account adiabatic expansion. 53 54 This function is used by EnergyBalance to iterate and solve the ODE. 55 56 @param T: The temperature at which to evaluate the differential equation 57 @type T: float 58 @param r: The radial point(s) (cm) 59 @type r: array/float 60 @param v: The velocity profile object 61 @type v: Velocity() 62 @param gamma: The adiabatic coefficient profile as function of T 63 @type gamma: Profiler() 64 65 @keyword rates: The heating and cooling rates, summed up (erg/s/cm3, H-C). 66 Includes the density factor, without the velocity component. 67 This speeds up iteration, limiting the radial grid 68 evaluation where possible. Default if only adiabatic cooling 69 is taken into account. 70 71 (default: None) 72 @type rates: Profiler() 73 @keyword warn: Warn when extrapolation occurs in an interpolation object. 74 Not applicable to functional evaluation, or gamma in which 75 case the extrapolation is deemed safe (constant value at 76 low and high end temperatures consistent with diatomic gases) 77 78 (default: 1) 79 @type warn: bool 80 81 @return: The derivative with respect to radius (K/cm) 82 @rtype: array/float 83 84 ''' 85 86 #-- Set the gamma-based factor. gamma is assumed to be a Profiler object. 87 tg = gamma.eval(T,warn=0) 88 factor1 = 2-2*tg 89 factor2 = tg-1 90 91 #-- Determine the velocity and its derivative from the Velocity() object 92 vi = v.eval(r,warn=warn) 93 dvi = v.diff(r,warn=warn) 94 95 #-- Calculate dT/dr 96 # 1) adiabatic cooling term 97 term1 = factor1 * (1./r + 0.5*1./vi*dvi) * T 98 99 # 2) other heating and cooling terms (can be none at all). This includes 100 # the density factor, except the velocity. 101 rr = rates.eval(r,warn=warn) if rates else 0. 102 term2 = factor2 * rr / vi 103 104 #-- And the final summation of the two terms 105 Tprime = term1 + term2 106 107 return Tprime

108 109 110

111 -class EnergyBalance(object):

112 113 ''' 114 The energy balance class. 115 116 Author: R. Lombaert, H. Olofsson & M. Maercker (Chalmers, Sweden) 117 118 For the use of the EnergyBalance module, please contact one of the three 119 authors listed here. 120 121 Calculates energy balance and provides tools for reading and writing I/O, 122 and setting input physics. 123 124 For now limited to the energy balance from the dust formation radius up to 125 the outer radius. 126 127 An example for calculating the energy balance (see the respective functions 128 for in-depth information): 129 >>> from cc.modeling.physics import EnergyBalance as EB 130 >>> eb = EB.EnergyBalance() 131 >>> eb.iterT() 132 >>> eb.plotT() 133 >>> eb.plotRates() 134 This calculates and plots the temperature profiles and heating/cooling 135 rates across all iterations using the default settings given in 136 cc.path.aux/inputEnergyBalance_standard.dat 137 138 Default in the inputEnergyBalance_standard file can be adapted through 139 keywords passed to EnergyBalance. A few examples: 140 >>> eb = EB.EnergyBalance(a=[0.005,0.25,300,1]) 141 which creates a grain-size grid from 0.005 to 0.25 micron with 300 142 points distributed in log scale. 143 144 >>> eb = EB.EnergyBalance(hterms=['dg','dt']) 145 adds dust-gas collisional heating and heat exchange to the adiabatic 146 cooling as heat/cool terms. 147 148 >>> eb = EB.EnergyBalance(template='gastronoom') 149 uses the input template to reproduce the GASTRoNOoM settings. 150 151 >>> eb = EB.EnergyBalance(gamma='h2') 152 calculates the adiabatic coefficient from the temperature-dependent measured 153 coefficient of molecular hydrogen. 154 155 >>> eb = EB.EnergyBalance(heatmode='gs2014') 156 calculates the heating and cooling terms from dust-gas interaction based on 157 the formalism of Gail & Sedlmayr 2014 as opposed to the classical 158 implementations of Decin et al. 2006 and Schoïer et al. 2001. 159 160 Furthermore, the iterative procedure can be adapted to user-specific 161 needs: 162 >>> eb.iterT(conv=0.005,imax=200,step_size=0.05) 163 sets the relative convergence criterion to 0.005 from default 0.01, the 164 maximum iterations to 200, and the step_size of the gradual maximum 165 allowed temperature change between iterations to 5%. 166 167 ''' 168

169 - def __init__(self,fn=None,template='standard',**kwargs):

170 171 ''' 172 Creating an EnergyBalance object used to calculate the energy balance. 173 174 Can work in conjunction with iteration with a RT code. 175 176 Takes into account several contributors to the heating and cooling of 177 an outflow: 178 - Adiabatic cooling 179 - dust-gas collisional indirect heating (elastic collisions) 180 - dust-gas direct heat exchange (accommodation) 181 182 Not yet implemented: 183 - Radiative line cooling and heating (for H2, CO, H2O, ...) 184 - Photoelectric heating 185 - Heating by cosmic rays 186 187 An instance of EnergyBalance contains three types of properties: 188 - basic input parameters given by the inputfile 189 - independent variables and basic profiles set before calculations 190 - profiles calculated from the basic parameters, variables, and 191 profiles, reset whenever any of the previous two properties are 192 changed. 193 194 First set is given in aux/inputEnergyBalance.dat. 195 196 Second set includes: radius, grain size, gas velocity, 197 gas/dust mass-loss rates, dust opacities, dust temperature, stellar 198 radiance. 199 200 Third set includes: gas temperature, gas/dust density, drift velocity. 201 202 Note that for now line cooling parameters are fixed. Hence after a new 203 RT model was calculated, a new energy balance must be calculated. In the 204 future, the level pops will be updateable. 205 206 @keyword fn: The parameter inputfile filename. If not given, a default 207 inputfile is read from aux/inputEnergyBalance.dat. 208 209 (default: aux/inputEnergyBalance.dat) 210 @type fn: str 211 @keyword template: The inputfile template used for initializing the 212 energy balance. Default is the standard, with as much 213 consistency as possible. Alternatives are 'mcp' and 214 'gastronoom' to reproduce the respective code's 215 results. Templates cannot be changed, but inputfiles 216 given by fn overwrite any input given there. 217 218 (default: 'standard') 219 @type template: str 220 221 @keyword kwargs: Input parameters listed in the default 222 inputEnergyBalance file can also be passed through the 223 initialisation of this instance. They will replace the 224 values given in the inputfile. 225 226 (default: {}) 227 @type kwargs: dict 228 229 ''' 230 231 #-- A welcome message, and credentials 232 m = "Welcome to the EnergyBalance module of ComboCode! \n"+\ 233 "This module was written at Chalmers University of Technology. "+\ 234 "The main author is R. Lombaert, with support from H. Olofsson "+\ 235 "and M. Maercker. Please contact one of these people if you wish "+\ 236 "to use the EnergyBalance with the purpose of publishing results." 237 print m 238 239 #-- Initialise all variables 240 self.__reset() 241 242 #-- Set default input template, and the custom parameter file 243 self.template = template.lower() 244 if self.template not in ['mcp','gastronoom']: self.template = 'standard' 245 self.fn = fn 246 self.include_lc = 0 247 248 #-- Read input parameters, set constants, coordinate grids, stellar 249 # properties, basic profiles, adiabatic coefficient, line cooling, 250 # velocity 251 self.readInputParameters(**copy.deepcopy(kwargs)) 252 self.setGrids() 253 self.setStar() 254 self.setProfiles() 255 self.setGamma() 256 self.setVelocity() 257 258 #-- Everything has been initialised. Set the current Temperature profile 259 self.__setT() 260 261 #-- Set the abundance for each molecule. The initial T has been set, so 262 # also calculate the line cooling terms if requested. 263 for m in self.molecules: 264 self.setAbundance(m) 265 if self.include_lc: self.setLineCooling(m)

266 267 268

269 - def formatInput(self,ilst):

270 271 ''' 272 Format an input keyword list for a variable. 273 274 Can be given as a string, split by spaces, or as the list output from 275 str.split(). 276 277 Assumes the first element is the requested function or a constant, 278 followed by key=value pairs that serve as the input keyword arguments 279 for the function if the first element is not a constant. The first 280 element thus does not contain a '=' sign. 281 282 The key=value pairs are formatted into a dictionary. 283 284 Example str: 285 'read_opacity species=AMCDHSPREI index=1 order=5' 286 Example list: 287 ['read_opacity','species=AMCDHSPREI','index=1','order=5'] 288 289 @param ilst: The input line, arguments split by spaces 290 @type ilst: list[str]/str 291 292 @return: 293 @rtype: list[str,dict] 294 295 ''' 296 297 #-- Check if a list or a str is given. Apply the ' ' split. 298 if not isinstance(ilst,collections.Iterable) or isinstance(ilst,str): 299 ilst = str(ilst).split() 300 301 #-- If the list is empty, return an empty entry 302 if not ilst: 303 return ['',{}] 304 305 #-- Read the first element and parse the rest a dictionary. Either the 306 # first element gives a constant, or it gives a read method for files 307 # The rest are the function input arguments in key=value pairs 308 ddk = {'convert_lists':1,'convert_floats':1,'convert_ints':1} 309 return [ilst[0],DataIO.readDict(lines=ilst[1:],**ddk)]

310 311 312

313 - def readInputParameters(self,**kwargs):

314 315 ''' 316 Read the input parameters for the energy balance. 317 318 If a filename is not given, the one set upon initialisation is used. If 319 none was given upon initialisation the default one from aux/ is read. 320 321 Resets any variables that depend on independent variables. 322 323 @keyword fn: The parameter inputfile filename. If None, a default 324 inputfile is read from aux/ 325 326 (default: None) 327 @type fn: str 328 329 @keyword kwargs: Input parameters listed in the default 330 inputEnergyBalance file can also be passed through the 331 initialisation of this instance. They will replace the 332 values given in the inputfile. 333 334 (default: {}) 335 @type kwargs: dict 336 337 ''' 338 339 #-- Only reset if self.T is already yet defined. 340 if self.T: self.__reset() 341 342 #-- Read the template 343 fnbase = os.path.join(cc.path.aux,'inputEnergyBalance_') 344 templates = {k: '{}{}.dat'.format(fnbase,k) 345 for k in ['standard','mcp','gastronoom']} 346 self.pars = DataIO.readDict(templates[self.template],convert_lists=1,\ 347 convert_floats=1,multi_keys=['molecule']) 348 if not self.pars.has_key('molecule') or self.pars['molecule'] == ['']: 349 self.pars['molecule'] = [] 350 351 #-- Read from the filename if it is given. 352 if not self.fn is None: 353 dd = DataIO.readDict(self.fn,convert_lists=1,convert_floats=1,\ 354 multi_keys=['molecule']) 355 m = dd.pop('molecule',[]) 356 self.pars.update(dd) 357 self.pars['molecule'] += m 358 359 #-- Insert any extra arguments given in function call 360 m = kwargs.pop('molecule',[]) 361 if isinstance(m, str): m = [m] 362 self.pars.update(kwargs) 363 self.pars['molecule'] += m 364 365 #-- Format the molecule information, and extract molecule names. 366 self.pars['molecule'] = [mlst.split() for mlst in self.pars['molecule']] 367 abuns = [] 368 self.molecules = [] 369 for mlst in self.pars['molecule']: 370 m = mlst.pop(0) 371 if not m in self.molecules: 372 self.molecules.append(m) 373 abuns.append(mlst) 374 self.pars['molecule'] = abuns 375 376 #-- Read left over information and parse it as a dictionary. This 377 # contains the abundance information. Only used when no pops are given 378 # and for Hpe in case of 12C16O 379 self.pars['molecule'] = [self.formatInput(mlst) 380 for mlst in self.pars['molecule']] 381 382 #-- Then initialise the heating and cooling terms requested 383 self.H = {term: {} for term in self.pars['hterms']} 384 self.C = {term: {} for term in self.pars['cterms'] if term != 'lc'} 385 386 #-- Adiabatic cooling is always calculated 387 self.C['ad'] = {} 388 389 #-- Prepare the settings for the molecular abundances. Note that these 390 # abundance are taken from the model output when possible in case 391 # level populations are used. If those are not available, the info 392 # given in the molecule input is used. 393 for m in self.molecules: 394 self.abun = {} 395 396 #-- Add the molecules for line cooling if applicable, and set the 397 # containers for collision rates, level pops, spectroscopy, Texc, 398 # abundances. Spectroscopy dict will refer to either collis or pop 399 # depending on the code (LamdaReader or MlineReader, respectively) 400 if 'lc' in self.pars['cterms']: 401 self.pars['cterms'] = [p for p in self.pars['cterms'] if p != 'lc'] 402 self.include_lc = 1 403 self.collis = {} 404 self.pop = {} 405 self.Texc = {} 406 self.ipop = {} 407 self.mol = {} 408 409 #-- Set the dict entries for each molecule. Also add a decorated 410 # function for each molecule that adds the molecule as an 411 # argument by default, so Clc can be called without arguments. 412 # (note that 'lc' does not occur in the self.C dict keys) 413 for m in self.molecules: 414 fname = 'lc_{}'.format(m) 415 self.C[fname] = {} 416 self.pars['cterms'].append(fname) 417 setattr(self,'C'+fname,functools.partial(self.Clc,m)) 418 419 #-- And also set the proper readers depending on the code. 420 if self.pars['rtcode'].lower() == 'gastronoom': 421 self.colread = CollisReader.CollisReader 422 self.popread = MlineReader.MlineReader 423 else: 424 self.colread = LamdaReader.LamdaReader 425 self.popread = PopReader.PopReader

426 427 428

429 - def setGrids(self,r=None,a=None,l=None):

430 431 ''' 432 Initialise the coordinate grids: radius (cm), grain size (cm), 433 wavelength (cm). 434 435 The grids can be given as lists, which are used as input for 436 Gridding.makeGrid. If not given, they are taken from the inputfile. 437 438 @keyword r: The radius grid (cm). [rmin,rmax,nsteps,log?] 439 440 (default: None) 441 @type r: list 442 @keyword a: The grain size grid (cm). [amin,amax,nsteps,log?] Can be 443 single value in case of one average grain size. 444 445 (default: None) 446 @type a: list/float 447 @keyword l: The wavelength grid (cm). [lmin,lmax,nsteps,log?] 448 449 (default: None) 450 @type l: list 451 452 ''' 453 454 #-- Only reset if self.T is already yet defined. 455 if self.T: self.__reset() 456 457 #-- Take the grids from the inputfile if not given here 458 if r is None: r = self.pars['r'] 459 if a is None: a = self.pars['a'] 460 if l is None: l = self.pars['l'] 461 462 #-- Check if a is a single value or a grid 463 a = Data.arrayify(a) 464 if a.size < 2: 465 self.a = a 466 else: 467 self.a = Gridding.makeGrid(*a) 468 469 #-- Set radius and wavelength 470 self.r = Gridding.makeGrid(*r) 471 self.l = Gridding.makeGrid(*l)

472 473 474

475 - def setStar(self,rstar=None,Tstar=None,Ltype=None):

476 477 ''' 478 Set the stellar properties. 479 480 If any of these are None, they are taken from the inputEnergyBalance.dat 481 file. 482 483 Defines the Radiance profile based on these properties. 484 485 @keyword rstar: The stellar radius (cm) 486 487 (default: None) 488 @type rstar: float 489 @keyword Tstar: The stellar effective temperature (K) 490 491 (default: None) 492 @type Tstar: float 493 @keyword Ltype: The type of stellar radiance profile (blackbody, ...) 494 495 (default: None) 496 @type Ltype: float 497 498 ''' 499 500 #-- Only reset if self.T is already yet defined. 501 if self.T: self.__reset() 502 503 #-- Take the profile definitions from the inputfile if not given here 504 if rstar is None: rstar = self.pars['rstar'] 505 if Tstar is None: Tstar = self.pars['Tstar'] 506 if Ltype is None: Ltype = self.pars['Ltype'] 507 508 self.rstar = rstar 509 self.Tstar = Tstar 510 511 #-- Set the stellar radiance profile and stellar surface area 512 self.rad = Radiance.Radiance(l=self.l,func=Ltype,T=Tstar) 513 self.rad.setSurface(rstar)

514 515 516

517 - def setProfiles(self,opac=None,mdot=None,mdot_dust=None,T=None,Td=None, 518 mu=None):

519 520 ''' 521 Initialise the independent profiles:opacity, mass-loss rate, 522 initial temperature, dust temperature 523 524 The grids can be given as a 1- or 2-item list, which are used as input 525 for the respective Profiler classes. If not given, they are taken from 526 the inputfile. 527 528 Format of the input lists: 529 - [constant value] 530 - [func/str,{kwargs}] 531 532 @keyword opac: The opacity profiles (l, cm^2/g). [func,{pars}] 533 534 (default: None) 535 @type opac: [func,dict] 536 @keyword mdot: The mass-loss rate profile (r, Msun/yr) [func,{pars}] or 537 constant. 538 539 (default: None) 540 @type mdot: [func,dict] 541 @keyword mdot_dust: The dust mass-loss rate profile (r, Msun/yr) 542 [func,{pars}] or constant. 543 544 (default: None) 545 @type mdot_dust: [func,dict] 546 @keyword T: The initial temperature profile (r, K). [func,{pars}] 547 The first arg is the function, followed by T0 and r0. Can be 548 string Td as well to set it to the dust temperature. 549 550 (default: None) 551 @type T: [func,dict] 552 @keyword Td: The dust temperature profile (r, K). [func,{pars}] 553 The first arg is the function. 554 555 (default: None) 556 @type Td: [func,dict] 557 @keyword mu: The mean molecular weight if nonstandard. Otherwise 558 calculated from fH and fHe. Could be specified in case, 559 e.g., some sort of dissociation is ongoing in the inner 560 wind. 561 562 (default: None) 563 @type mu: float 564 565 ''' 566 567 #-- Only reset if self.T is already yet defined. 568 if self.T: self.__reset() 569 570 #-- Take the profile definitions from the inputfile if not given here 571 if opac is None: opac = self.formatInput(self.pars['opac']) 572 if mdot is None: mdot = self.formatInput(self.pars['mdot']) 573 if mdot_dust is None: 574 mdot_dust = self.formatInput(self.pars['mdot_dust']) 575 if T is None: T = self.formatInput(self.pars['Tinit']) 576 if Td is None: Td = self.formatInput(self.pars['Td']) 577 if mu is None: mu = self.pars.get('mu',None) 578 579 #-- Set the opacity profile: l, func, pars 580 self.opac = Opacity.Opacity(self.l,opac[0],**opac[1]) 581 582 #-- Set the mass-loss-rate profiles. Profiler checks itself for constant 583 # If a constant, then mdot[1] is an empty dict 584 self.mdot = Mdot.Mdot(self.r,mdot[0],**mdot[1]) 585 586 #-- Set the dust temperature profile 587 self.Td = Temperature.Temperature(self.r,Td[0],**Td[1]) 588 589 #-- Use a power law for T in the inner wind? 590 self.inner = T[1]['inner'] 591 592 #-- Set the initial temperature profile 593 if T[0].lower() == 'td': 594 self.T_iter[0] = Temperature.Temperature(self.r,Td[0],\ 595 inner=self.inner,**Td[1]) 596 else: 597 self.T_iter[0] = Temperature.Temperature(self.r,T[0],**T[1]) 598 599 #-- Determine the T epsilon going from rstar to r0. This is used for the 600 # inner wind. This can be turned off by setting inner=0 in Tinit. 601 # Reset the inner wind power law for the T profiler as well. 602 self.T0 = self.T_iter[0].T0 603 self.r0 = self.T_iter[0].r0 604 self.inner_eps = -np.log(self.Tstar/self.T0)/np.log(self.rstar/self.r0) 605 self.T_iter[0].setInnerEps(self.inner_eps) 606 607 #-- Set the dust mass-loss rate, and add the inner radius if needed. 608 if not mdot_dust[1].has_key('r0'): mdot_dust[1]['r0'] = self.r0 609 self.mdot_dust = Mdot.Mdot(self.r,mdot_dust[0],**mdot_dust[1]) 610 611 #-- Set the mean molecular weight 612 if mu is None: 613 mu = Density.mean_molecular_weight(self.pars['fH'],self.pars['fHe']) 614 self.mu = mu

615 616 617

618 - def setGamma(self,gamma=None):

619 620 ''' 621 Set the adiabatic coefficient. 622 623 Three options: 624 - Constant value 625 - Step function at 350K (going 7/5 => 5/3 at lower T, MCP model) 626 - T-dependent gamma for H2 gas 627 628 @keyword gamma: The adiabatic coefficient profile (K, /). Is either a 629 constant (float), or 'h2' in which case the adiabatic 630 coefficient is read from a file as a function of T, or 631 h2_mcp in which case gamma is a step function: 7./5. for 632 T>350 K, 5./3. otherwise. 633 634 (default: None) 635 @type gamma: float/str 636 637 ''' 638 639 #-- When gamma == 'h2' or 'h2_mcp', vary gamma as function of T, 640 # appropriate for H2. 641 if isinstance(self.pars['gamma'],str): 642 gamma = self.pars['gamma'].lower() 643 #-- The MCP implementation 644 if gamma == 'h2_mcp': 645 #-- At low T the diatomic molecule H2 (gamma~1.4) loses internal 646 # degrees of freedom, and becomes like a mono-atomic gas 647 # (gamma~1.67). 648 kwargs = {'x': self.T_iter[0].eval(inner_eps=self.inner_eps,\ 649 warn=not self.inner), 650 'func': Profiler.step, 651 'ylow': 5./3., 'yhigh': 7./5., 'xstep': 350. } 652 gamma = Profiler.Profiler(**kwargs) 653 #-- Interpolate the adiabatic coefficient taken from NIST database 654 else: 655 fn = os.path.join(cc.path.aux,'h2_physical_properties.dat') 656 kwargs = {'x': self.T_iter[0].eval(inner_eps=self.inner_eps,\ 657 warn=not self.inner), 658 'func': Profiler.interp_file, 659 'ikwargs': {'ext': 3, 'k': 3}, 660 'filename': fn, 'ycol': -1} 661 gamma = Profiler.Profiler(**kwargs) 662 663 #-- Alternatively, gamma is constant, and gives the value of the 664 # coefficient 665 else: 666 kwargs = {'x': self.T_iter[0].eval(inner_eps=self.inner_eps,\ 667 warn=not self.inner), \ 668 'func': Profiler.constant, 'c': self.pars['gamma']} 669 gamma = Profiler.Profiler(**kwargs) 670 671 self.gamma = gamma

672 673 674

675 - def setVelocity(self,v=None):

676 677 ''' 678 Set the velocity profile, based on the requested function. 679 680 Done separately from other profiles, because this can depend on the 681 temperature at the inner radius (independent), and on the adiabatic 682 coefficient (for the sound velocity). 683 684 This is done before the line cooling term is set. 685 686 @keyword v: The velocity profile (r, cm/s). [func,{pars}] 687 688 (default: None) 689 @type v: [func,dict] or [cst] 690 691 ''' 692 693 #-- Take v from the inputfile if not given 694 if v is None: v, vdd = self.formatInput(self.pars['v']) 695 else: v, vdd = v[0], v[1] 696 697 #-- if sound velocity requested for v0, and not enough pars given, 698 # set automatically 699 if v == 'vbeta' and vdd.get('v0',0) == 'vs' and not vdd.has_key('T'): 700 vdd['T'], vdd['mu'] = self.T0, self.mu 701 vdd['gamma'] = self.gamma.eval(self.T0) 702 703 #-- Set the velocity profile: r, func, dfunc, order interp dfunc, pars 704 self.v = Velocity.Velocity(self.r,v,**vdd)

705 706 707

708 - def setLineCooling(self,m):

709 710 ''' 711 Set the line cooling parameters. 712 713 This includes reading the collision rates and level populations. 714 715 This is done per molecule, and upon initialisation of the EnergyBalance. 716 717 A distinction is made between the source of the spectroscopy/level pops: 718 - GASTRoNOoM: The info is taken from the MlineReader and CollisReader 719 - ALI: The info is taken from the LamdaReader, PopReader and .par 720 output file 721 722 @param m: The molecule name from the input molecules list. 723 @type m: str 724 725 ''' 726 727 #-- Pops were read when the abundances were set. Read the collision 728 # rates 729 imol = self.molecules.index(m) 730 self.collis[m] = self.colread(self.pars['collis'][imol]) 731 732 #-- ipop remembers the iteration number for which the level 733 # populations were set, so they can be updated later. 734 self.ipop[m] = self.i 735 736 #-- If no level populations are given, define some default ones. 737 if not self.pars['pop'][imol]: 738 self.setPopInitial(m) 739 740 #-- Pops are code-dependent. Spectroscopy as well. Set that here 741 elif self.pars['rtcode'].lower() == 'gastronoom': 742 #-- Note that the abundance is in the mline output for GASTRoNOoM 743 # Hence, the populations have already been read. 744 # Refer the spectroscopy to the mline file for later use 745 self.mol[m] = self.pop[m] 746 else: 747 self.pop[m] = self.popread(self.pars['pop'][imol]) 748 749 #-- Refer the spectroscopy to the lamda file for later use 750 self.mol[m] = self.collis[m] 751 752 #-- Set the coll rate and level pop interpolators 753 self.pop[m].setInterp(itype='spline',k=3,ext=3) 754 self.collis[m].setInterp(itype='spline',k=1,ext=0) 755 756 #-- Calculate the cooling term for the initial temperature 757 self.Clc(m,update_pop=0) 758 759 #-- Calculate the excitation temperature 760 self.setTexc(m)

761 762 763

764 - def setPopInitial(self,m):

765 766 ''' 767 Set an initial set of level populations. 768 769 Based on the Boltzmann distribution for a given excitation temperature. 770 The kinetic temperature is taken for now. Should likely be scalable in 771 the future. 772 773 Solves the set of equations: 774 - Sum(n_i) = 1 775 - n_i = n_0 * (g_l/g_0) exp((E_0 - E_l)/Tkin) 776 777 @param m: The molecule name from the input molecules list. 778 @type m: str 779 780 ''' 781 782 #-- In case of GASTRoNOoM, need the radiat file to be read 783 if self.pars['rtcode'] == 'gastronoom': 784 imol = self.molecules.index(m) 785 fn = self.pars['collis'][imol].replace('collis','radiat') 786 ny = max(self.collis[m]['coll_trans']['lup']) 787 self.mol[m] = RadiatReader.RadiatReader(fn,ny=ny) 788 789 #-- Get T profile and other information 790 T = self.T.eval(inner_eps=self.inner_eps,warn=not self.inner) 791 g0 = self.mol[m].getLWeight(1) 792 E0 = self.mol[m].getLEnergy(1,unit='erg') 793 ny = self.mol[m]['pars']['ny'] 794 g = self.mol[m].getLWeight(range(2,ny+1)) 795 E = self.mol[m].getLEnergy(range(2,ny+1),unit='erg') 796 797 #-- Force gups/Eups into a column vector for array multiplication 798 g.shape = (g.size,1) 799 E.shape = (E.size,1) 800 801 #-- The equation that calculates n_0 802 Ediff = E0 - E 803 nfactor = g/g0*np.exp(np.outer(Ediff,1./(k_b*0.999*T))) 804 denominator = 1 + sum(nfactor) 805 n0 = 1./denominator 806 807 #-- Calculate all the other levels 808 n = n0*nfactor 809 810 #-- Create a PopReader object with filename "None" 811 self.pop[m] = PopReader.PopReader(None) 812 self.pop[m].setP(self.r) 813 self.pop[m].setNY(ny) 814 self.pop[m].setPop(index=1,n=n0) 815 for i in range(2,ny+1): 816 self.pop[m].setPop(index=i,n=n[i-2,:])

817 818 819

820 - def setAbundance(self,m):

821 822 ''' 823 Set the abundance profile for a molecule. 824 825 Two possibilities: 826 - Level populations are given, which include molecular abundance 827 profiles as well, and are RT code dependent 828 - No level pops are given, so use the default abundance profile 829 830 The second possibility requires a read method and a filename to be given 831 in the molecule definition, e.g. 832 molecule=12C16O np.loadtxt fname=waql.par usecols=[1,4] skiprows=9 unpack=1 833 Alternatively, a constant value can be given as well (not advised). 834 835 @param m: The molecule name from the input molecules list. 836 @type m: str 837 838 ''' 839 840 #-- Which molecule? 841 imol = self.molecules.index(m) 842 843 #-- Check if level pops are given (those files contain an abundance 844 # profile as well). If not, use the molecule definition. 845 if not self.pars['pop'] or len(self.pars['pop']) <= imol \ 846 or not self.pars['pop'][imol]: 847 mlst = self.pars['molecule'][imol] 848 #-- Run a check if an abundance is actually given. 849 if not mlst[0]: 850 q = "No abundance input given in molecule={} ".format(m)+\ 851 "definition. Cannot calculate default level populations "+\ 852 "or set photoelectric heating." 853 raise IOError(q) 854 855 #-- A constant value might be requested, otherwise read the file 856 # If constant value, the kwargs dict will be empty. 857 if len(mlst[1]) == 0: 858 p, abun = self.r, mlst[0] 859 else: 860 dd = DataIO.read(func=mlst[0],**mlst[1]) 861 862 #-- Check whether the values are given in rstar or in cm 863 if dd[0][0]/self.pars['rstar']<0.1: 864 p = dd[0]*self.pars['rstar'] 865 else: 866 p = dd[0] 867 868 #-- In case three columns are read, nmol and nh2 are read 869 # separately 870 if len(dd) == 2: 871 abun = dd[1] 872 elif len(dd) == 3: 873 abun = dd[1]/dd[2] 874 875 #-- If pops are given, extract the abundance profile, depending on the 876 # rt code. Molecule information in inputfile is then not used. 877 elif self.pars['rtcode'].lower() == 'gastronoom': 878 #-- Abundances for mline in GASTRoNOoM are in the mline output, 879 # which contains the level populations. So read that here already 880 self.pop[m] = self.popread(self.pars['pop'][imol]) 881 882 #-- Extract the abundance profile as a function of impact parameter 883 p = self.pop[m].getP() 884 abun = self.pop[m].getProp('amol') 885 else: 886 #-- Abundance profile for MCP/ALI is in a separate parameter output 887 # file 888 fn = self.pars['pop'][imol].replace('pop','log') 889 890 #-- Find out number of cells. Note vols 0 and 1 are not used, so -2. 891 nvol = DataIO.getKeyData(filename=fn,incr=0,keyword='nvol')[0] 892 nvol = int([line.split('=')[1] 893 for line in nvol 894 if 'nvol' in line.lower()][0]) 895 896 #-- Find out where the data start + 2 lines for units and separator 897 data = DataIO.readFile(fn,' ') 898 ixmol = DataIO.findKey(i=0,key='xmol',data=data) + 2 899 p, abun = np.genfromtxt(fn,usecols=[1,4],skip_header=ixmol+1,\ 900 unpack=1,max_rows=nvol-2) 901 902 #-- Constant value is requested 903 if len(p) != len(abun): 904 abun_interp = mlst[0] 905 906 #-- We need flexible extrapolation, use interp1d and linear interp 907 # Linear because cubic is too unstable for the last radial step in 908 # a steep abundance decline. 909 else: 910 abun_interp = interp1d(x=p,y=abun,kind='linear',bounds_error=0,\ 911 assume_sorted=1,fill_value=(abun[0],0.)) 912 self.abun[m] = Profiler.Profiler(x=self.r,func=abun_interp)

913 914 915

916 - def setTexc(self,m):

917 918 ''' 919 Calculate the excitation temperature for the level populations given in 920 self.pop. This method is called when the pops are read for the first 921 time (or set by setInitialPop), and when the level populations are 922 updated. 923 924 The excitation temperature is set as an array of dimensions 925 (number of transition indices, number of impact parameter), hence for a 926 given transition index, you can retrieve the excitation temperature as 927 self.Texc[m][i-1,:] as a function of impact parameter. Alternatively, 928 you can access Texc by calling getTexc and passing either the indices or 929 the lup/llow. 930 931 The impact parameter grid can be retrieved from self.pop[m].getP(). 932 933 Note that the excitation temperature is calculated for all radiative 934 transitions included in the molecular spectroscopy; not for all 935 collisional transitions. 936 937 @param m: The molecule tag 938 @type m: str 939 940 ''' 941 942 #-- Retrieve the transition indices, and then extract relevant 943 # information for all transitions. 944 indices = self.mol[m].getTI(itype='trans') 945 lup = self.mol[m].getTUpper(indices) 946 llow = self.mol[m].getTLower(indices) 947 gl = self.mol[m].getLWeight(llow) 948 gu = self.mol[m].getLWeight(lup) 949 El = self.mol[m].getLEnergy(llow,unit='erg') 950 Eu = self.mol[m].getLEnergy(lup,unit='erg') 951 952 #-- Retrieve the associated populations 953 nl = self.pop[m].getPop(llow) 954 nu = self.pop[m].getPop(lup) 955 956 #-- Force g/E into column vectors for array multiplication 957 gl.shape = (gl.size,1) 958 gu.shape = (gu.size,1) 959 El.shape = (El.size,1) 960 Eu.shape = (Eu.size,1) 961 962 #-- The equation that calculates Texc 963 self.Texc[m] = -1./np.log(nu/nl*gl/gu)*(Eu-El)/k_b

964 965 966

967 - def getTexc(self,m,indices=None,lup=None,llow=None):

968 969 """ 970 Return the excitation temperature profile as a function of impact 971 parameter for given transition indices of a given molecule. 972 973 Can also take upper and/or lower level indices to determine the 974 transition indices. 975 976 The excitation temperature is returned as an array of dimensions 977 (number of transition indices, number of impact parameters). 978 979 @param m: The molecule tag 980 @type m: str 981 982 @keyword indices: The transition indices. Can be indirectly given 983 through upper and/or lower level indices. If default 984 and no llow/lup are given, all Texc are returned for 985 this molecule. 986 987 (default: None) 988 @type indices: list/array 989 @keyword lup: The upper level indices used to extract the transition 990 indices. If this or llow are defined, the keyword indices 991 is overwritten. 992 993 (default: None) 994 @type lup: list/array 995 @keyword llow: The upper level indices used to extract the transition 996 indices. If this or lup are defined, the keyword indices 997 is overwritten. 998 999 (default: None) 1000 @type llow: list/array 1001 1002 @return: The excitation temperature, with array shape (n_index,n_p) 1003 @rtype: array 1004 1005 """ 1006 1007 #-- If either lup or llow are given, get transition indices from them 1008 if not (lup is None and llow is None): 1009 if not lup is None: lup = Data.arrayify(lup) 1010 if not llow is None: llow = Data.arrayify(llow) 1011 indices = self.mol[m].getTI(lup=lup,llow=llow) 1012 1013 #-- If indices is still None, all transitions are returned 1014 if indices is None: return self.Texc[m] 1015 1016 #-- Make sure the indices are given as an array, and select. Note 1017 # python's 0-based indexing vs the fortran 1-based indexing. 1018 indices = Data.arrayify(indices) 1019 return self.Texc[m][indices-1,:]

1020 1021 1022

1023 - def __reset(self):

1024 1025 ''' 1026 Reset all variables that depend on the radius, velocity, etc. 1027 1028 ''' 1029 1030 #-- Density profiles 1031 self.gdens = None 1032 self.ddens = None 1033 self.nd = None 1034 1035 #-- Dust velocity profile 1036 self.vd = None 1037 self.w = None 1038 1039 #-- Temperature profile and iterations 1040 self.T = None 1041 self.T_iter = dict() 1042 self.i = 0 1043 1044 #-- And the heating and cooling terms 1045 self.H = {} 1046 self.C = {} 1047 self.rates = {}

1048 1049 1050

1051 - def __next_iter(self):

1052 1053 ''' 1054 Reset some of the properties that depend on temperature so they can be 1055 calculated anew. 1056 1057 ''' 1058 1059 #-- Dust and drift velocity profiles depend on T if w_thermal != none 1060 self.vd = None 1061 self.w = None 1062 1063 #-- Set the current temperature profile to the new calculation 1064 self.__setT() 1065 1066 #-- Recalculate the dust velocity. Drift will be done as well. 1067 self.setVdust()

1068 1069 1070

1071 - def setDrift(self):

1072 1073 ''' 1074 Calculate the drift velocity profile. 1075 1076 This assumes a stellar radiance profile is given, of which the 1077 wavelength grid is used to integrate Q_l * L_l, thereby interpolating 1078 the opacity. 1079 1080 The two available modes are standard and beta: 1081 - standard: Calculates the drift velocity based on the balance between 1082 the radiation pressure and the drag force. 1083 - vbeta: Follows MCP, where the terminal drift velocity is calculated 1084 from the terminal gas velocity, based on the balance between radiation 1085 pressure and drag force, with a beta law going to that max velocity. 1086 1087 Note that the vbeta mode follows the F mode for dust velocity in MCP. 1088 1089 1090 ''' 1091 1092 if self.w is None: 1093 #-- Set density for the mean molecular weight. 1094 self.setDensity('gas') 1095 1096 #-- Define shared input parameters for driftRPDF function 1097 kwargs = {'a': self.a, 'l': self.l, 'P': self.pars['P'], 1098 'sd': self.pars['sd'], 'opac': self.opac, 1099 'radiance': self.rad, 'alpha': self.pars['alpha'], 1100 'mu': self.gdens.getMeanMolecularWeight()} 1101 1102 #-- Check if the MCP method is requested, otherwise do 'standard' 1103 method = self.formatInput(self.pars['w_mode']) 1104 if method[0] == 'vbeta2D': 1105 #-- Calculates the terminal drift velocity from the terminal gas 1106 # velocity, and sets a beta law for the inner wind up to it 1107 #-- For this the terminal gas velocity is needed, and the mdot 1108 # at the inner radius (no variable mass loss or thermal comp) 1109 kwargs['v'] = max(self.v.eval()) 1110 kwargs['mdot'] = self.mdot.eval(self.r0) 1111 kwargs['w_thermal'] = 'none' 1112 wmax = Velocity.driftRPDF(r=self.r0,**kwargs)[0] 1113 1114 #-- Set the drift as a beta law with sensible values, unless 1115 # they were defined in the w_mode key. Even the max velocity 1116 # can be fixed in the input through vinf. 1117 v0 = method[1].get('v0',self.v.eval(self.r0)) 1118 r0 = method[1].get('r0',self.r0) 1119 beta = method[1].get('beta',1.0) 1120 wmax = method[1].get('vinf',wmax) 1121 self.w = Velocity.Drift(r=self.r,a=self.a,\ 1122 func=Velocity.vbeta2D,\ 1123 r0=r0,v0=v0,vinf=wmax,beta=beta) 1124 1125 else: 1126 #-- If standard, add the relevant driftRPDF keywords and set the 1127 # profile 1128 kwargs['w_thermal'] = self.pars['w_thermal'] 1129 kwargs['v'] = self.v 1130 kwargs['r'] = self.r 1131 kwargs['mdot'] = self.mdot 1132 kwargs['func'] = Velocity.driftRPDF 1133 1134 #-- Include T if a thermal velocity term is needed 1135 vtherm_types = ['kwok','mean','rms','prob','epstein'] 1136 if self.pars['w_thermal'].lower() in vtherm_types: 1137 #-- T is always set upon initialisation 1138 kwargs['T'] = self.T 1139 1140 self.w = Velocity.Drift(**kwargs)

1141 1142 1143

1144 - def setVdust(self):

1145 1146 ''' 1147 Set the dust velocity profile. The drift is averaged over grain size. 1148 1149 ''' 1150 1151 if self.vd is None: 1152 self.setDrift() 1153 kwargs = {'r': self.r, 'func': Velocity.vdust, 1154 'w': self.w, 'v': self.v, 'norm_type': 'standard'} 1155 self.vd = Velocity.Velocity(**kwargs)

1156 1157 1158

1159 - def setDensity(self,dtype):

1160 1161 ''' 1162 Calculate the density profile for gas or dust. 1163 1164 The dust density profile is not dependent on grain size, since the drift 1165 averaged over grain size. 1166 1167 @param dtype: The density profile type, being 'gas' or 'dust'. 1168 @type dtype: str 1169 1170 ''' 1171 1172 #-- Set the density profile for either dust or gas. Nothing is done if 1173 # dtype is not recognized. 1174 if dtype.lower() == 'gas': 1175 if self.gdens is None: 1176 kwargs = {'r': self.r,'func':Density.dens_mdot, 1177 'mdot': self.mdot, 'v': self.v} 1178 self.gdens = Density.Density(**kwargs) 1179 kwargs = {'fH': self.pars['fH'], 'fHe': self.pars['fHe']} 1180 self.gdens.calcNumberDensity(**kwargs) 1181 1182 #-- Replace mean molecular weight in case a nonstandard value is 1183 # requested 1184 self.gdens.mu = self.mu 1185 1186 if dtype.lower() == 'dust': 1187 #-- Calculate the mass density 1188 if self.ddens is None: 1189 #-- Calculate the dust velocity profile. The drift is averaged 1190 # over the grain size. 1191 self.setVdust() 1192 kwargs = {'r': self.r,'func':Density.dens_mdot, 'dust': 1, 1193 'mdot': self.mdot_dust, 'v': self.vd} 1194 self.ddens = Density.Density(**kwargs) 1195 1196 #-- Calculate the number density 1197 if self.nd is None and self.a.size > 1: 1198 #-- Set the grain size distribution. self.nd refers to a 1199 # Grainsize.Distribution object! 1200 kwargs = {'r': self.r, 'a': self.a, 'sd': self.pars['sd'], 1201 'func': getattr(Grainsize,self.pars['GSD']), 1202 'w': self.w, 1203 'v': self.v, 'mdot_dust': self.mdot_dust} 1204 self.nd = Grainsize.Distribution(**kwargs) 1205 1206 elif self.nd is None: 1207 #-- Working with average grain size. Calculate number density 1208 # from the mass density: rho/(4/3 pi rho_specific a^3) 1209 self.ddens.calcNumberDensity(sd=self.pars['sd'],a=self.a) 1210 #-- In this case, self.nd just refers to the Density() object 1211 self.nd = self.ddens

1212 1213 1214

1215 - def updatePop(self,m,fn='',updateLC=0):

1216 1217 ''' 1218 Update the level populations. 1219 1220 This reads the level populations from the same file as listed upon 1221 creation of the EnergyBalance and assumes those level populations are 1222 appropriate for the temperature profile of the iteration number self.i. 1223 1224 This is done per molecule. 1225 1226 @param m: The molecule name from the input molecules list. 1227 @type m: str 1228 1229 @keyword fn: The filename of a populations file, in case it changes or 1230 is added anew after default populations have been used. 1231 Ignored if the file does not exist. 1232 1233 (default: '') 1234 @type fn: str 1235 @keyword updateLC: The flag that moves the self.ipop to the current (or 1236 next) iteration, so that calcClc knows a new line 1237 cooling term must be updated. This is always set to 1 1238 if the populations are updated. Setting this to 1 1239 upon function call, tells calcClc the LC term must be 1240 calculated anew regardless of reading new pops (and 1241 thus uses the old pops with the new T profile of this 1242 iteration). This key can be set in the inputfile. 1243 1244 (default: 0) 1245 @type updateLC: bool 1246 1247 ''' 1248 1249 #-- First check if populations are given. If not, don't update them. 1250 imol = self.molecules.index(m) 1251 1252 #-- Maybe a new file is added for pops 1253 if os.path.isfile(str(fn)): 1254 self.pars['pop'][imol] = fn 1255 1256 #-- If a file name is available, check if updating is needed 1257 if self.pars['pop'][imol]: 1258 #-- Select the files and read them with the Reader objects. 1259 npop = self.popread(self.pars['pop'][imol]) 1260 if False in [np.array_equal(npop.getPop(i),self.pop[m].getPop(i)) 1261 for i in npop.getLI()]: 1262 self.pop[m] = npop 1263 1264 #-- Set the level pop interpolators and calc Texc 1265 self.setTexc(m) 1266 self.pop[m].setInterp(itype='spline',k=3,ext=3) 1267 1268 #-- Make sure to tell calcClc to recalculate the LC term 1269 updateLC = 1 1270 1271 #-- Tell calcClc to recalculate the cooling term for this molecule 1272 if updateLC or self.pars.get('updateLC',0): 1273 #-- If there's a key for this iteration index already, the cooling 1274 # term was already calculated, so set the index for the update to 1275 # the next iteration. 1276 if self.C['lc_{}'.format(m)].has_key(self.i): 1277 self.ipop[m] = self.i + 1 1278 else: 1279 self.ipop[m] = self.i

1280 1281 1282

1283 - def __setT(self):

1284 1285 ''' 1286 Set the current temperature profile for this iteration. When i == 0, 1287 this is the initial guess. 1288 1289 This is done explicitly, because upon initialisation the state of the 1290 object depends on self.T being None or not. 1291 1292 ''' 1293 1294 self.T = self.T_iter[self.i]

1295 1296 1297

1298 - def iterT(self,conv=0.01,imax=50,step_size=0.,dTmax=0.10,warn=1,\ 1299 *args,**kwargs):

1300 1301 ''' 1302 Iterate the temperature profile until convergence criterion is reached. 1303 1304 Extra arguments are passed on to the spline1d interpolation of the 1305 total cooling and heating terms through the calcT() call. 1306 1307 @keyword conv: The maximum relative allowed change in T for convergence. 1308 1309 (default: 0.01) 1310 @type conv: float 1311 @keyword imax: Maximum number of allowed iterations. Code stops when 1312 this number is reached, but not before dTmax = 1. 1313 1314 (default: 50) 1315 @type imax: int 1316 @keyword step_size: The minimum step size in maximum allowed relative 1317 temperature change. The code decides dynamically 1318 based on how close the iteration is to convergence 1319 which multiple of this step size the next iteration 1320 allows. If T change percentage grows to fast, 1321 decrease this number. Default is 0, in which case 1322 the maximum allowed T change is kept constant at 1323 dTmax. 1324 1325 (default: 0) 1326 @type step_size: float 1327 @keyword dTmax: The starting value for the maximum allowed relative 1328 temperature change between iterations. This maximum 1329 increases as the code reaches convergences, through a 1330 set multiple of step_size. If step_size is 0, dTmax 1331 stays constant 1332 1333 (default: 0.10) 1334 @type dTmax: float 1335 @keyword warn: Warn when extrapolation occurs. 1336 1337 (default: 1) 1338 @type warn: bool 1339 1340 ''' 1341 1342 print('-- Iterating T(r) now.') 1343 dTnsteps = (1.-dTmax)/step_size if step_size else 0. 1344 steps = 0. 1345 1346 #-- If first iteration in this call of iterT, calculate it in any case 1347 # This is not necessarily iteration 0! The pops may have been updated. 1348 print('Iteration {} for T(r)...'.format(self.i+1)) 1349 self.calcT(dTmax,warn=warn,*args,**kwargs) 1350 1351 #-- Then continue iterating until the relative difference between this 1352 # and the previous result is smaller than convergence criterion in all 1353 # radial points, or when the maximum number of iterations is reached 1354 while True: 1355 dTdiff = np.abs(1.-self.T.eval(warn=not self.inner)\ 1356 /self.T_iter[self.i-1].eval(warn=not self.inner)) 1357 if (np.all(dTdiff<conv) and steps == dTnsteps) or self.i == imax: 1358 break 1359 if np.all(dTdiff<(dTnsteps-steps)*conv): 1360 if dTnsteps-steps > 8.: steps += 4. 1361 elif dTnsteps-steps > 4.: steps += 2. 1362 else: steps += 1. 1363 print('Iteration {} for T(r)...'.format(self.i+1)) 1364 self.calcT(dTmax+step_size*steps,warn=warn,*args,**kwargs) 1365 if not np.all(np.isfinite(self.T.eval(inner_eps=self.inner_eps,\ 1366 warn=not self.inner))): 1367 print('NaNs found in T-profile. Breaking off iteration.') 1368 break

1369 1370 1371

1372 - def calcT(self,dTmax=1,warn=1,ode_kwargs={},*args,**kwargs):

1373 1374 ''' 1375 Calculate the temperature profile based on the differential equation 1376 given by Goldreich & Scoville (1976). The function is defined in dTdr. 1377 1378 The initial temperature is taken to be the second argument of the 1379 initial temperature profile given by T in inputEnergyBalance.dat. 1380 This is typically the condensation temperature. 1381 1382 Additionals arguments are passed on to the spline1d interpolation of the 1383 total cooling and heating terms, e.g. k=3, ext=0 are defaults. 1384 1385 @keyword dTmax: The maximum allowed relative temperature change for this 1386 T calculation. Set to 100% by default. 1387 1388 (default: 1) 1389 @type dTmax: float 1390 @keyword warn: Warn when extrapolation occurs. 1391 1392 (default: 1) 1393 @type warn: bool 1394 @keyword ode_kwargs: Extra arguments for the ODE solver. They are added 1395 to the dict ode_args made by the function, and 1396 hence overwrites any defaults. In principle, the 1397 defaults are fine, but can be overwritten if needed 1398 1399 (default: {}) 1400 @type ode_kwargs: dict 1401 1402 ''' 1403 1404 #-- Note that T is always set upon initialisation and at the end of a 1405 # new T calculation. 1406 1407 #-- Move to the next iteration. Cooling and heating terms belong to this 1408 # next iteration (even tho they are based on initial guess T profile, 1409 # except for the adiabatic cooling term) 1410 self.i += 1 1411 1412 #-- Calculate the gas density if it hasn't been calculated yet 1413 self.setDensity('gas') 1414 1415 #-- Calculate the heating and cooling terms 1416 # Note that the adiabatic term is used explicitly in the dTdr 1417 # function. The rate is calculated with the new T profile at the end 1418 # for plotting. 1419 for term in self.pars['hterms']: 1420 getattr(self,'H'+term)() 1421 for term in self.pars['cterms']: 1422 getattr(self,'C'+term)() 1423 1424 #-- The abundance term can be included here as well to reduce the number 1425 # of profile evaluations in dTdr. 1426 nh2 = self.gdens.eval(dtype='nh2') 1427 1428 #-- The ascale factor is the correction for non-zero He and atomic H 1429 # abundances (Decin et al. 2006). Reduces to original if fH = fHe = 0. 1430 ascale = self.gdens.fH + 1. + self.gdens.fHe * (self.gdens.fH + 2.) 1431 1432 #-- The density term does not yet include the velocity factor. This is 1433 # included when dTdr is calculated. 1434 dens_term = nh2 * k_b * ascale 1435 1436 #-- Then set the heating and cooling terms as profiler objects 1437 yrates = sum([self.H[k][self.i] for k in self.H.keys()]) \ 1438 - sum([self.C[k][self.i] for k in self.C.keys() if k != 'ad']) 1439 1440 #-- Remember the rates profiler for each iteration 1441 self.rates[self.i] = yrates/dens_term 1442 1443 #-- Note that if no rates are requested, this will give an empty array 1444 if self.r.size == Data.arrayify(yrates).size: 1445 rp = Profiler.Profiler(self.r,spline1d(self.r,self.rates[self.i],\ 1446 *args,**kwargs)) 1447 else: 1448 rp = None 1449 1450 #-- Calculate the next iteration of the temperature profile 1451 ode_args = {'func': dTdr, 'y0': self.T0, 't': self.r, 1452 'args': (self.v,self.gamma,rp,warn)} 1453 ode_args.update(ode_kwargs) 1454 Tr = odeint(**ode_args)[:,0] 1455 1456 #-- Check the temperature change. Limit it to the requested maximum 1457 # allowed change. Only do this from the second iteration, to allow 1458 # a big jump from the initial condition. 1459 if self.i < 0: 1460 dTmax = self.pars['dTmax'] 1461 print('Changing T by {:.1f}%.'.format(dTmax*100.)) 1462 Ti = self.T.eval(inner_eps=self.inner_eps,warn=not self.inner) 1463 Tr = np.where(Tr<=0.,np.ones_like(Tr),Tr) 1464 Tr = Ti + dTmax*(Tr-Ti) 1465 1466 #-- Create a profiler for the new temperature structure. Extrapolation 1467 # done by returning boundary values, ie T0 and the temp at outer 1468 # boundary 1469 Tinterp = spline1d(self.r,Tr,k=3,ext=3) 1470 keys = {'inner':self.inner,'inner_eps':self.inner_eps,'r0':self.r0, 1471 'T0':self.T0} 1472 self.T_iter[self.i] = Temperature.Temperature(self.r,Tinterp,**keys) 1473 1474 #-- Initialize some T-dependent properties for the next iteration. Also 1475 # sets the new temperature profile as the current one. 1476 self.__next_iter() 1477 1478 #-- Calculate the adiabatic term for plotting. Also sets the current 1479 # temperature profile as the new calculation. 1480 self.Cad()

1481 1482 1483

1484 - def Cad(self):

1485 1486 ''' 1487 Calculate the adiabatic cooling rate. This is not used by the solver of 1488 the dTdr differential equation, and is only for plotting purposes. 1489 1490 Uses the latest calculation of the T-profile. 1491 1492 Equation derived from Decin et al. 2006 and Danilovich 2016. 1493 1494 ''' 1495 1496 #-- if cooling rate has already been calculated: don't do anything 1497 if self.C['ad'].has_key(self.i): return 1498 1499 #-- Set the gamma-based factor evaluated at the new T profile 1500 # calculated. 1501 gamma = self.gamma.eval(self.T.eval(inner_eps=self.inner_eps,\ 1502 warn=not self.inner),warn=0) 1503 1504 #-- Makes the rate negative, in accordance with -C in the second term of 1505 # the diff eq. 1506 factor1 = 2-2*gamma 1507 factor2 = gamma-1 1508 1509 #-- Determine the velocity and its derivative from the Velocity() object 1510 # Note that these will never complain bout extrapolation unless they 1511 # are profiles read beforehand instead of functions such as vbeta 1512 vi = self.v.eval() 1513 dvi = self.v.diff() 1514 1515 #-- Calculate the conversion term K/cm => erg/s/cm3: allows to compare 1516 # C terms with adiabatic term. 1517 self.setDensity('gas') 1518 nh2 = self.gdens.eval(dtype='nh2') 1519 ascale = self.gdens.fH + 1. + self.gdens.fHe * (self.gdens.fH + 2.) 1520 conversion = nh2 * k_b * vi * ascale / factor2 1521 1522 #-- Calculate the temperature dependent term 1523 main_term = (1./self.r + 0.5*1./vi*dvi) \ 1524 * self.T.eval(inner_eps=self.inner_eps,\ 1525 warn=not self.inner) 1526 1527 #-- Calculate the rate. Multiply by -1 to have positive cooling rate. 1528 # [erg/s/cm3] 1529 self.C['ad'][self.i] = -1. * conversion * factor1 * main_term

1530 1531 1532

1533 - def Hdg(self):

1534 1535 ''' 1536 Calculate the heating rate by dust-gas collisions. 1537 1538 Note that this term evaluates to zero for the inner wind at r<r0, if the 1539 dust density is 0 there. This can be enforced by choosing an mdot_step 1540 function. 1541 1542 Sputtering is applied if a grain size distribution is used. For a 1543 constant grain size, this is not done, since it's not realistic to 1544 remove all dust if the drift becomes too large. 1545 1546 The equation is derived from Goldreich & Scoville 1976, based on Schoier 1547 et al 2001, and Decin et al. 2006: 1548 1549 General case: 1550 Hdg = n_d sigma_d w x 1/2 rho_g w^2 1551 => 1552 Hdg = pi/2 n_H2 m_H (fH+2) (1+4fHe) (1-P)^(2/3) Int(n_d w^3 a^2 da) 1553 1554 MCP/ALI case (average grain size a): 1555 Hdg = pi n_H2 m_H n_d w^3 a^2 1556 1557 GASTRoNOoM case (MRN): 1558 Hdg = pi/2 A n_H2^2 m_H (fH+2)^2 (1+4fHe) Int(w^3 a^-1.5 da) 1559 1560 ''' 1561 1562 #-- if heating rate has already been calculated: don't do anything 1563 if self.H['dg'].has_key(self.i): return 1564 1565 #-- Initialize the current temperature profile as the initial guess 1566 # This is not done if i > 0, since the last solution is set as the 1567 # current T profile after a calculation 1568 if not self.T: 1569 self.__setT() 1570 1571 #-- Initialise density and a few parameters 1572 self.setDrift() 1573 self.setDensity('gas') 1574 self.setDensity('dust') 1575 nh2 = self.gdens.eval(dtype='nh2') 1576 fH = self.pars['fH'] 1577 fHe = self.pars['fHe'] 1578 P = self.pars['P'] 1579 mode = self.pars['heatmode'].lower() 1580 1581 #-- Set the constant factor. Same for both modes 1582 cfac = np.pi*nh2*mh*(fH+2.)*(1.+4.*fHe)*(1.-P)**(-2./3.) 1583 1584 #-- Hdg according to Gail & Sedlmayr 1585 if mode == 'gs2014': 1586 #-- The thermal velocity and factor 4/3. for specular reflection: 1587 vT = 4./3.*np.sqrt(8.*k_b*self.T.eval(inner_eps=self.inner_eps,\ 1588 warn=not self.inner)\ 1589 /(self.gdens.mu*np.pi*mh)) 1590 1591 #-- Hdg according to Decin et al. 2006: extra factor 0.5 for Ekin, no 1592 # thermal component. 1593 else: 1594 cfac *= 0.5 1595 vT = 0. 1596 1597 #-- Differ between having a grain size distribution or not. Note that 1598 # this factor remains constant for the same gas density, regardless of 1599 # fH, fHe. It is a different matter for gsfac where drift enters. 1600 if self.a.size > 1: 1601 #-- GSD: integrate drift, dust number density, a^2 over grain size 1602 # self.nd is a Grainsize.Distribution object! Apply sputtering. 1603 vT = np.transpose([vT]) 1604 nd = self.nd.eval(w=self.w.eval(),w_sputter=self.pars['w_sputter']) 1605 w = self.w.eval() 1606 afac = nd*self.a**2.*w**2.*(vT**2.+w**2.)**0.5 1607 gsfac = trapz(x=self.a,y=afac,axis=1) 1608 else: 1609 #-- Avg grain size: no integration 1610 # self.nd refers to the self.ddust, or a dust Density() object! 1611 # no sputtering for a single grain size 1612 nd = self.nd.eval(dtype='ndust') 1613 w = np.squeeze(self.w.eval()) 1614 gsfac = nd*self.a**2*w**2.*(vT**2.+w**2.)**0.5 1615 1616 #-- Set dust-gas collisional heating [erg/s/cm3] 1617 self.H['dg'][self.i] = cfac*gsfac

1618 1619 1620

1621 - def Hdt(self):

1622 1623 ''' 1624 Calculate the heating rate by dust-gas heat exchange. 1625 1626 The equation is derived from Burke & Hollenbach 1983, based on 1627 Groenewegen 1994, Schoier et al. 2001, and Decin et al. 2006: 1628 1629 Note that this term evaluates to zero for the inner wind at r<r0, if the 1630 dust density is 0 there. This can be enforced by choosing an mdot_step 1631 function. 1632 1633 Sputtering is applied if a grain size distribution is used. For a 1634 constant grain size, this is not done, since it's not realistic to 1635 remove all dust if the drift becomes too large. 1636 1637 The older implementations used by MCP/ALI and GASTRoNOoM follow Burke & 1638 Hollenbach 1983 and Groenwegen 1994. The general-case implementation 1639 follows the more recent work of Gail & Sedlmayr, adding in proper vT and 1640 drift terms. The GS2014 is not yet fully implemented, but the Hdt term 1641 is already available here. 1642 1643 General case (following Gail & Sedlmayr 2014, see Eq 15.19): 1644 Hdt = alpha pi k_b n_h2 (fH+2.)(1.-P)^(-2./3.) Int(n_d a^2 da) (T_d - T) 1645 sqrt(8*vT^2+drift^2) (1/(gamma-1)) 1646 1647 MCP/ALI case (n_d: dust number density, average grain size a): 1648 Hdt = 2 alpha pi k_b n_h2 (n_d a^2) (T_d - T) vT 1649 1650 GASTRoNOoM case (n_d: distribution following MRN): 1651 Hdt = 2 alpha pi k_b n_h2 (fH+2.) Int(n_d a^2 da) (T_d - T) vT 1652 1653 In all cases, alpha is the accommodation coefficient (from Groenewegen 1654 1994): 1655 alpha = 0.35 exp(-sqrt(0.002*(T_d + T)))+0.1 1656 1657 and vT is the thermal velocity: 1658 vT = sqrt(8 k_b t/(mu pi m_H)) 1659 1660 ''' 1661 1662 #-- if heating rate has already been calculated: don't do anything 1663 if self.H['dt'].has_key(self.i): return 1664 1665 #-- Initialise density and a few parameters 1666 self.setDrift() 1667 self.setDensity('gas') 1668 self.setDensity('dust') 1669 nh2 = self.gdens.eval(dtype='nh2') 1670 fH = self.pars['fH'] 1671 P = self.pars['P'] 1672 mode = self.pars['heatmode'] 1673 1674 #-- Extract the dust and gas temperatures 1675 td = self.Td.eval() 1676 t = self.T.eval(inner_eps=self.inner_eps,warn=not self.inner) 1677 1678 #-- Calculate the average accommodation coefficient 1679 accom = 0.35*np.exp(-np.sqrt((td+t)*0.002))+0.1 1680 1681 #-- Calculate the constant factor as a function of r 1682 cfac = np.pi*k_b*nh2*(fH+2.)*(1.-P)**(-2./3.) 1683 1684 #-- The thermal velocity: 1685 vT = np.sqrt(8*k_b*t/(self.gdens.mu*np.pi*mh)) 1686 1687 #-- Calculate the temp factor 1688 tdiff = (td-t) 1689 1690 #-- Differ between having a grain size distribution or not. 1691 if self.a.size > 1: 1692 #-- GSD: integrate dust number density, a^2 over grain size 1693 # self.nd is a Grainsize.Distribution object! Apply sputtering. 1694 nd = self.nd.eval(w=self.w.eval(),w_sputter=self.pars['w_sputter']) 1695 gsfac = trapz(x=self.a,y=nd*self.a**2,axis=1) 1696 else: 1697 #-- Avg grain size: no integration 1698 # self.nd refers to the self.ddust, or a dust Density() object! 1699 # A single grain size: no sputtering. 1700 nd = self.nd.eval(dtype='ndust') 1701 gsfac = nd*self.a**2 1702 1703 #-- Set the dust-gas collisional heating term. First mode is Gail & 1704 # Sedlmayr's equation from 2014, based on Draine 1980. Default is the 1705 # accommodation heat exchange by Burke & Hollenbach 1983 1706 if mode.lower() == 'gs2014': 1707 #-- Avg over a^2. Normalisation disappears when combined with gsfac 1708 # Sputtering is not included for the average of w, since it would 1709 # skew the distribution of drift velocities to a lower value. We 1710 # cannot approximate what happens to destroyed grains at this 1711 # point, hence we don't know how this would affect the drift 1712 drift = self.w.avgDrift(norm_type='collisional',nd=self.nd) 1713 1714 #-- fg/2. based on gamma, with fg = 2/(gamma-1) 1715 fg = 1./(self.gamma.eval(self.T.eval(inner_eps=self.inner_eps,\ 1716 warn=not self.inner),\ 1717 warn=0)-1.) 1718 1719 #-- Velocity factor 1720 vfac = (8.*vT**2+drift**2)**0.5 1721 self.H['dt'][self.i] = accom * cfac * fg * tdiff * gsfac * vfac 1722 1723 else: 1724 #-- Extra fac 2 enters through 2kT (BH1983) 1725 self.H['dt'][self.i] = 2. * accom * cfac * tdiff * vT * gsfac

1726 1727 1728

1729 - def Hpe(self):

1730 1731 ''' 1732 Calculate the heating rate by the photoelectric effect. 1733 1734 Two methods are available at present: 1735 1) Following Draine 1978 and Huggins et al. 1988, as used by MCP 1736 (see also Crosas & Menten 1997) 1737 2) Following Bakes & Tielens 1994., as implemented by Decin et al. 1738 2006 in GASTRoNOoM. 1739 1740 No options to tweak these methods has been implemented yet, but a lot of 1741 consistency checks should be done, and some of the assumptions can be 1742 improved on with a consistent calculation (e.g. Av in method 2). 1743 1744 Photoelectric heating as derived by Bakes & Tielens is maybe a good 1745 basis for further development of the heating term, but see also Woitke 1746 2015 (2015EPJWC.10200011W) for recent developments. 1747 1748 The derivation in method 2 makes a lot of assumptions and are applied 1749 specifically to the case of GASTRoNOoM (e.g. fixed grain size 1750 distribution). Hence, the implementation is considered appropriate for 1751 the default settings in the inputEnergyBalance_gastronoom.dat file. Any 1752 deviations from that must be treated carefully. 1753 1754 No sputtering is applied. Small grains are the most relevant for Hpe, 1755 and any sputtering would reduce the size of existing dust grains, thus 1756 increasing the amount of small grains. Since we cannot have Hpe directly 1757 depend on the grain size distribution (and instead work with a scaling 1758 factor amin_scale), we do not apply sputtering to the photoelectric 1759 heating. 1760 1761 ''' 1762 1763 #-- if heating rate has already been calculated: don't do anything 1764 if self.H['pe'].has_key(self.i): return 1765 1766 #-- Initialise density 1767 self.setDensity('gas') 1768 self.setDensity('dust') 1769 1770 #-- Method 1: Scale the heating rate with the gas density and a constant 1771 # factor (Draine 1978) -- Maybe ngas is better here 1772 method = self.formatInput(self.pars['pe_method']) 1773 if method[0].lower() == 'draine': 1774 nh2 = self.gdens.eval(dtype='nh2') 1775 #-- Default Kpe 1e-26 erg/s, following Huggins et al. 1988 and 1776 # Schoier & Olofsson 2001 1777 Kpe = method[1].get('Kpe',1e-26) 1778 self.H['pe'][self.i] = Kpe * nh2 1779 1780 #-- Method 2: More elaborate calculation of the heating rate, though 1781 # caution is advised due to the uncertainties in quite a few of the 1782 # assumptions. 1783 elif method[0].lower() == 'bakes': 1784 T = self.T.eval(inner_eps=self.inner_eps,warn=not self.inner) 1785 G0 = method[1].get('G0',1.) 1786 amin_scale = method[1].get('amin_scale',0.2) 1787 d2g = self.ddens.eval()/self.gdens.eval() 1788 nhtot = self.gdens.eval(dtype='nhtot') 1789 1790 #-- Calculate H_pe-hat from Decin 2006 EQ. 9 (see Sect 2.2.2) 1791 # First: electron density. Assumed to come from photodissociation 1792 # of CO into C and O, whereby C is almost immediately ionised. 1793 # Hence, difference between C/O > 1 and C/O < 1, since the former 1794 # has free C available for ionisation. 1795 # Note that we want the relative CO abundance compared to the max 1796 # at R_STAR or R_INNER. 1797 aco = self.abun['12C16O'].eval() 1798 xco = aco/aco[0] 1799 abun_c = method[1].get('abun_c') 1800 abun_o = method[1].get('abun_o') 1801 if abun_c > abun_o: 1802 ne = nhtot * abun_c * (1-xco*abun_o/abun_c) 1803 1804 #-- If C/O < 1, all C is in CO: ne = 0, in the inner wind. More e in 1805 # outer wind where CO gets photodissociated. 1806 # Note that a small bump is formed possibly in the inner wind, but 1807 # this is due to the GASTRoNOoM interpolation, and having values 1808 # very close to zero. It's not important 1809 # considering how small the rate is in the inner wind. 1810 else: 1811 ne = nhtot * abun_c * (1-xco) 1812 1813 #-- Calculate the factor itself, noting the a_min scaling and G0 1814 Hpehat = amin_scale*1e-24*nhtot 1815 Hpehat[ne>0] *= 0.03/(1+2e-4*G0*np.sqrt(T[ne>0])/ne[ne>0]) 1816 Hpehat[ne==0] = 0. 1817 1818 #-- Calculate the correction on H_pe-hat based on the optical depth 1819 # Note the use of G0 again. tau_uv = 1.8 * A_v = 1.8 * 1.6e-22 N_H 1820 # H Column density from R_outer integrated inwards, as a function 1821 # of self.r. Note r_outer-r to have a positive column density. 1822 Nh = cumtrapz(nhtot[::-1],x=(self.r[-1]-self.r)[::-1],\ 1823 initial=0.)[::-1] 1824 1825 #-- Extra scaling factor for the dust-to-gas ratio in the extinction 1826 Av = 1.6e-22*Nh/0.01*d2g 1827 self.H['pe'][self.i] = Hpehat*G0*np.exp(-1.8*Av)

1828 1829 1830

1831 - def Hcr(self):

1832 1833 ''' 1834 Calculate the heating rate by cosmic rays. 1835 1836 The rate is taken from Goldsmith & Langer 1978, and is also used by 1837 Groenwegen 1994 and Decin et al 2006. This is a very approximate formula 1838 and will need updating in the future. 1839 1840 Two modes are available: The standard one, using n(H2), and the one used 1841 by Groenewegen 1994 and Decin et al. 2006 (where all gas mass is placed 1842 into H2, even if fH or fHe are non-zero). Controlled through the keyword 1843 cr_method == 'groenewegen' or cr_method == 'standard' in the inputfile. 1844 1845 ''' 1846 1847 #-- if heating rate has already been calculated: don't do anything 1848 if self.H['cr'].has_key(self.i): return 1849 1850 #-- Initialise density 1851 self.setDensity('gas') 1852 1853 #-- Goldsmith & Langer 1978 give the rate for ionisation of H2 by cosmic 1854 # rays. However Groenewegen 1994 and Decin 2006 use this rate by 1855 # placing all gas mass into H2 (even if fH and fHe are not zero), ie 1856 # by dividing the density rho by 2*m_H 1857 # Technically this should only be nh2, but we follow Groenewegen/Decin 1858 # for now. Why?! Both options are implemented. 1859 nh2 = self.gdens.eval(dtype='nh2') 1860 fH = self.gdens.fH 1861 fHe = self.gdens.fHe 1862 1863 #-- Decin et al 2006/Groenewegen 1994 1864 if self.pars['cr_method'].lower() == 'groenewegen': 1865 self.H['cr'][self.i] = 6.4e-28 * nh2 * (fH/2.+1.)*(1.+4.*fHe) 1866 else: 1867 self.H['cr'][self.i] = 6.4e-28 * nh2

1868 1869 1870

1871 - def Ch2(self):

1872 1873 ''' 1874 Calculate the line cooling rate by vibrational excitation of H_2. 1875 1876 Two modes are available (used by MCP/ALI and GASTRoNOoM, respectively): 1877 1) Following Groenewegen 1994, based on Hartquist et al 1980. Based 1878 on fitting of tabulated data of H2 cooling under LTE conditions 1879 2) Following Decin et al 2006, based on GS1976, Hollenbach & McKee 1880 1979 and Hollenbach & McKee 1989. 1881 1882 The keyword h2_method (groenewegen, or decin) determines which of the 1883 two is used. 1884 1885 ''' 1886 1887 #-- if cooling rate has already been calculated: don't do anything 1888 if self.C['h2'].has_key(self.i): return 1889 1890 #-- Initialise density 1891 self.setDensity('gas') 1892 nh2 = self.gdens.eval(dtype='nh2') 1893 T = self.T.eval(inner_eps=self.inner_eps,warn=not self.inner) 1894 1895 #-- Method by Groenewegen 1994 1896 if self.pars['h2_method'].lower() == 'groenewegen': 1897 self.C['h2'][self.i] = 2.6111e-21*nh2*(T/1000.)**4.74 1898 1899 #-- Method by Decin et al 2006 1900 else: 1901 #-- Based on deexcitation collisions between H2-H2, and H2-H, so 1902 # need H-number density 1903 nh = self.gdens.eval(dtype='nh') 1904 1905 #-- Spontaneous emission rate from first vibrational state of H2 1906 # (in s^-)1 1907 A10 = 3e-7 1908 1909 #-- Energy (in erg) of the emitted photon (0.6 eV converted) 1910 hnu10 = 9.61305939e-13 1911 1912 #-- The collisional deexcitation rate factor alpha 1913 coll_deex_h = 1.0e-12*np.sqrt(T)*np.exp(-1000./T) 1914 coll_deex_h2 = 1.4e-12*np.sqrt(T)*np.exp(-18100./(T+1200.)) 1915 alpha = nh*coll_deex_h + nh2*coll_deex_h2 1916 1917 #-- Nominator and denominator for n1 1918 nomin = alpha*np.exp(-hnu10/k_b/T) 1919 denomin = alpha*(1+np.exp(-hnu10/k_b/T))+A10 1920 1921 #-- Number density of vibrationally excited molecules: 1922 n1 = nh2*nomin/denomin 1923 self.C['h2'][self.i] = A10 * hnu10 * n1

1924 1925 1926

1927 - def Clc(self,m,update_pop=1):

1928 1929 ''' 1930 Calculate the radiative line cooling rate for a molecule. 1931 1932 For GS2014, no correction term yet for the excitation temperature per 1933 level. 1934 1935 Follows Sahai 1990. 1936 1937 Cubic spline interpolation for the level populations. Linear 1938 interpolation and extrapolation for the collision rates (as for 1939 GASTRoNOoM). 1940 1941 Currently not yet implemented to use a sqrt(T/T0) extrapolation at lower 1942 boundary as is done by MCP/ALI. 1943 1944 Note that EnergyBalance will internally call this function with the 1945 molecule tag tacked onto the Clc function name. 1946 1947 @param m: The molecule name from the input molecules list. 1948 @type m: str 1949 @keyword update_pop: Check if the level populations have been updated 1950 1951 (default: 1) 1952 @type update_pop: bool 1953 1954 ''' 1955 1956 def calcLevelLC(llow,r,T): 1957 1958 ''' 1959 Calculate the line cooling contribution for a single transition 1960 index, with a fixed lower level i and for which j > i, for the 1961 entire radial grid. 1962 1963 Keep in mind, the goal is to include all transitions from every 1964 level to every level. It doesn't actually matter if the energy is 1965 lower or higher: The Sahai + Einstein equations change the sign if 1966 the lower level is really the upper level in terms of energy. 1967 1968 We can use the structure of the collision rate files, where all 1969 possible collisional transitions are included. By returning all 1970 transitions to a given lower level, we already have j > i. 1971 1972 A note must be made here. Normally one would want to work with 1973 all levels that have higher energy than Elow. However, for CO 1974 this leads to issues because some v=1 levels have lower energy 1975 than some v=0 levels, while they are still sorted going v=0 to 1976 jmax, then v=1 to jmax, ie not sorted by energy. The collision 1977 rates however assume that they are sorted by energy. Meaning 1978 that for some collisional transitions Eup-Elow becomes < 0 1979 because of how the CO spectroscopy is sorted. This is not 1980 necessarily a problem, hence why we assume Eup-Elow must be > 0 1981 in what follows, and force it to be through abs(Eup-Elow). 1982 We assume the collision rate files are sorted properly, thus 1983 retrieve all "higher energy levels" by simply passing the llow 1984 index to the getTI method. This leads to results that are 1985 identical with GASTRoNOoM CO cooling rates. 1986 1987 1988 @param llow: index of the transition lower level. 1989 @type llow: int 1990 @param r: The radial grid in cm 1991 @type r: array 1992 @param T: The temperature in K for which to calculate the cooling 1993 @type T: array 1994 1995 @return: The contribution to the line cooling for a given lower 1996 level i (in ergs * cm^3 / s). 1997 @rtype: float 1998 1999 ''' 2000 2001 #-- Energy, weight and population of the lower level 2002 Elow = self.mol[m].getLEnergy(index=llow,unit='erg') 2003 glow = self.mol[m].getLWeight(index=llow) 2004 poplow = self.pop[m].getInterp(llow)(r) 2005 2006 #-- Get the transition indices that go to the level with index llow 2007 # Enforce indices to be an array. 2008 indices = self.collis[m].getTI(itype='coll_trans',llow=(llow,)) 2009 2010 #-- In case no upper levels were found, llow is the highest level 2011 # available, and there are no collisions to be taken into account 2012 # return 0 2013 if not indices.size: 2014 return 0. 2015 2016 #-- Retrieve the level indices, energies, weights and populations of 2017 # upper levels 2018 lups = self.collis[m].getTUpper(index=indices,itype='coll_trans') 2019 Eups = self.mol[m].getLEnergy(lups,unit='erg') 2020 popups = np.array([self.pop[m].getInterp(lup)(r) for lup in lups]) 2021 gups = self.mol[m].getLWeight(index=lups) 2022 2023 #-- Force gups/Eups into a column vector for array multiplication 2024 gups.shape = (gups.size,1) 2025 Eups.shape = (Eups.size,1) 2026 2027 #-- Retrieve the rates between llow and all lups 2028 # Also get the respective weights and energies 2029 Culs = np.array([self.collis[m].getInterp(i)(T) for i in indices]) 2030 2031 #-- Calculate the reversed rate for lower to upper level. 2032 # Based on the Einstein relation: Cul/Clu = gl/gu exp(Eul/kT) 2033 expfac = np.exp(np.outer(-1*abs(Elow-Eups),1./(k_b*T))) 2034 Clus = Culs*np.multiply(expfac,gups/glow) 2035 2036 #-- Calculate the sum of the cooling contribution across all j>i 2037 # transitions 2038 return np.sum((Clus*poplow-Culs*popups)*abs(Eups-Elow),axis=0)

2039 2040 2041 #-- if cooling rate has already been calculated: don't do anything 2042 if self.C['lc_{}'.format(m)].has_key(self.i): 2043 return 2044 2045 #-- Initialise the gas density 2046 self.setDensity('gas') 2047 2048 #-- Check if the level populations have to be updated 2049 if update_pop: self.updatePop(m) 2050 2051 #-- Else if the level populations have not yet been updated, use the 2052 # cooling term calculated for the temperature for which the level pops 2053 # were calculated. Note that self.pop[m][0] always exists. It is 2054 # calculated upon initialisation. 2055 print "Passing through Clc at iteration %i"%self.i 2056 if self.i > self.ipop[m]: 2057 print "Taking previous Clc for iteration %i"%self.ipop[m] 2058 j = self.ipop[m] 2059 self.C['lc_{}'.format(m)][self.i] = self.C['lc_{}'.format(m)][j] 2060 return 2061 2062 print "Calculating new Clc for iteration %i"%self.i 2063 #-- Extract the relevant density profiles 2064 nh2 = self.gdens.eval(dtype='nh2') 2065 amol = self.abun[m].eval(warn=0) 2066 2067 #-- Calculate the line cooling term for this molecule. 2068 llows = self.pop[m].getLI() 2069 LCterm = np.empty(shape=(len(self.r),len(llows))) 2070 for i in llows: 2071 LCterm[:,i-1] = calcLevelLC(i,self.r,\ 2072 self.T.eval(inner_eps=self.inner_eps,\ 2073 warn=not self.inner)) 2074 LCtotal = nh2*nh2*amol*np.sum(LCterm,axis=1) 2075 2076 #-- Do NOT Multiply by -1. This already gives the net energy lost 2077 self.C['lc_{}'.format(m)][self.i] = LCtotal 2078 2079 2080

2081 - def plotRateIterations(self,iterations=[],dTsign='C',mechanism='ad',\ 2082 scale=1,fn=None,cfg=None,**kwargs):

2083 2084 ''' 2085 Plot the heating or cooling rates for several iterations of a given 2086 mechanism and heat exchange sign (cooling or heating). 2087 2088 @keyword iterations: The iteration indices to be plotted. If default, 2089 all iterations are plotted 2090 2091 (default: []) 2092 @type iterations: list 2093 @keyword dTsign: The dT type: heating (H) or cooling (C). 2094 2095 (default: 'C') 2096 @type dTsign: str 2097 @keyword mechanism: The cooling/heating mechanism to be plotted. 2098 2099 (default: 'ad') 2100 @type mechanism: type 2101 @keyword scale: Scale the heating and cooling rates with r^4/1e14cm. 2102 2103 (default: 1) 2104 @type scale: bool 2105 @keyword fn: The filename and path of the plot (without extension). If 2106 default, a filename can be given in a cfg file, and if that 2107 is not the case, the plot is simply shown but not saved. 2108 2109 (default: None) 2110 @type fn: str 2111 @keyword cfg: The filename to the cfg file for this plot with extra 2112 settings. 2113 2114 (default: None) 2115 @type cfg: str 2116 2117 @keyword kwargs: Additional keywords passed to the plotting method. 2118 Overwrites any keys given in the cfg file. 2119 2120 (default: {}) 2121 @type kwargs: dict 2122 2123 ''' 2124 2125 #-- Read the dict if available, and extract the filename. 2126 cfg_dict = Plotting2.readCfg(cfg) 2127 if cfg_dict.has_key('filename'): 2128 fn = cfg_dict.pop('filename') 2129 2130 #-- If not iterations are specified, plot all of them. 2131 if not iterations: 2132 iterations = range(1,self.i) 2133 2134 #-- Make sure we have a list, and proper input format. Adapt filename 2135 iterations = Data.arrayify(iterations) 2136 #if 0 in iterations: iterations = [i for i in iterations if i!=0] 2137 dTsign = dTsign[0].upper() 2138 if not fn is None: 2139 fn += '_{}{}'.format(dTsign,mechanism) 2140 fn += '_iter{}-{}'.format(iterations[0],iterations[-1]) 2141 ddict['filename'] = fn 2142 2143 #-- Set parameters. kwargs > cfg > locally defined 2144 pars = {'xlogscale': 1, 'ylogscale': 1, 'xaxis': 'r (cm)', 2145 'extension': 'pdf'} 2146 pars.update(cfg_dict) 2147 pars.update(kwargs) 2148 2149 #-- Additional settings for plot 2150 pmech = mechanism.replace('_','\_') 2151 ddict = {'yaxis': dTsign+pmech+'-rate (ergs s$^{-1}$ cm$^{-3}$)', 2152 'keytags': [str(i) for i in iterations], 'x': [], 'y': []} 2153 if scale: 2154 ddict['yaxis'] += r' $\times$ (r/10$^{14}$ cm)$^4$' 2155 2156 #-- Where are the mechanisms? Only select first two characters in mech 2157 mech_types = {'ad': 'C', 'dg': 'H', 'dt': 'H', 'lc': 'C', 2158 'pe': 'H', 'h2': 'C', 'cr': 'H' } 2159 mtype = mech_types[mechanism[:2]] 2160 2161 #-- Set the sign. Only select first two characters in mech 2162 Csigns = {'ad': 1, 'dg': -1, 'dt': -1, 'lc': 1, 2163 'pe': -1, 'h2': 1, 'cr': -1 } 2164 sign = Csigns[mechanism[:2]] 2165 if dTsign == 'H': 2166 sign = -1 * sign 2167 2168 #-- The heating/cooling terms for several iterations 2169 for i in iterations: 2170 isel = getattr(self,mtype)[mechanism][i]*sign > 0 2171 ix = self.r[isel] 2172 ddict['x'].append(ix) 2173 iy = getattr(self,mtype)[mechanism][i][isel]*sign 2174 if scale: iy = iy*(ix*1e-14)**4 2175 ddict['y'].append(iy) 2176 2177 #-- Add in other parameters 2178 if abs(min([min(yi) for yi in ddict['y']])) < 1e-20: 2179 ddict['ymin'] = 1e-20 2180 ddict.update(pars) 2181 fn = Plotting2.plotCols(**ddict) 2182 2183 if not fn is None: 2184 print('Heating and cooling rates plotted at:') 2185 print(fn+'.pdf')

2186 2187 2188

2189 - def plotRates(self,scale=1,fn=None,iteration=None,cfg=None,join=0,**kwargs):

2190 2191 ''' 2192 Plot the heating and cooling rates for a single iteration. 2193 2194 @keyword scale: Scale the heating and cooling rates with r^4/1e14cm. 2195 2196 (default: 1) 2197 @type scale: bool 2198 @keyword fn: The filename and path of the plot (without extension). If 2199 default, a filename can be given in a cfg file, and if that 2200 is not the case, the plot is simply shown but not saved. 2201 2202 (default: None) 2203 @type fn: str 2204 @keyword iteration: The iteration for which to plot the rates. Default 2205 is the last iteration. 2206 2207 (default: None) 2208 @type iteration: int 2209 @keyword cfg: The filename to the cfg file for this plot with extra 2210 settings. 2211 2212 (default: None) 2213 @type cfg: str 2214 @keyword join: Join together the plot for heating and cooling terms. 2215 2216 (default: 0) 2217 @type join: bool 2218 2219 @keyword kwargs: Additional keywords passed to the plotting method. 2220 Overwrites any keys given in the cfg file. 2221 2222 (default: {}) 2223 @type kwargs: dict 2224 2225 ''' 2226 2227 #-- Read the dict if available, and extract the filename. 2228 cfg_dict = Plotting2.readCfg(cfg) 2229 if cfg_dict.has_key('filename'): 2230 fn = cfg_dict.pop('filename') 2231 2232 #-- Check which iteration to plot, use default if not give or invalid. 2233 if iteration is None or iteration > self.i: iteration = self.i 2234 2235 #-- Add the iteration if any but the last one is requested 2236 if fn and iteration != self.i: fn += '_iter{}'.format(self.i) 2237 2238 #-- Default line_types 2239 lts = {'ad': '-.g', 'dg': '-.b', 'dt': '--r', 2240 'pe': '--g', 'h2': '-.k', 'cr': '-k' } 2241 mkeys = [k for k in self.C.keys() if k[:2] == 'lc'] 2242 lts_extra = ['-r', '--b', '--k', '-m', '--m', '-y', '--y', '-g'] 2243 for mk in mkeys: 2244 lts[mk] = lts_extra.pop(0) 2245 2246 #-- Set parameters. kwargs > cfg > locally defined 2247 pars = {'xlogscale': 1, 'ylogscale': 1, 'xaxis': 'r (cm)', 2248 'extension': 'pdf'} 2249 pars.update(cfg_dict) 2250 pars.update(kwargs) 2251 2252 #-- Iterate over heating and cooling type, two plots total 2253 fns = [] 2254 for sign in ['H','C']: 2255 ddict = {'yaxis': sign+'-rate (ergs s$^{-1}$ cm$^{-3}$)', 2256 'keytags': [], 'x': [], 'y': [], 'line_types': []} 2257 2258 if scale: 2259 ddict['yaxis'] += r' $\times$ (r/10$^{14}$ cm)$^4$' 2260 2261 #-- Set the filename if one is given with a suffix 2262 if fn: 2263 ddict['filename'] = '{}_{}'.format(fn,sign) 2264 2265 #-- The heating (alt. cooling) terms 2266 for term in sorted(getattr(self,sign).keys()): 2267 sterm = '$_\\mathrm{{{s}}}$'.format(s=term.replace('_','\_')) 2268 ddict['keytags'].append(sign+sterm) 2269 isel = getattr(self,sign)[term][iteration] > 0 2270 ix = self.r[isel] 2271 ddict['x'].append(ix) 2272 iy = getattr(self,sign)[term][iteration][isel] 2273 if scale: iy = iy*(ix*1e-14)**4 2274 ddict['y'].append(iy) 2275 ddict['line_types'].append(lts[term]) 2276 2277 #-- The negative cooling (alt. heating) terms 2278 negsign = 'C' if sign == 'H' else 'H' 2279 for term in sorted(getattr(self,negsign).keys()): 2280 isel = getattr(self,negsign)[term][iteration] < 0 2281 if not isel.any(): continue 2282 sterm = '$_\\mathrm{{{s}}}$'.format(s=term.replace('_','\_')) 2283 ddict['keytags'].append(sign+sterm) 2284 ix = self.r[isel] 2285 ddict['x'].append(ix) 2286 iy = -1.*getattr(self,negsign)[term][iteration][isel] 2287 if scale: iy = iy*(ix*1e-14)**4 2288 ddict['y'].append(iy) 2289 ddict['line_types'].append(lts[term]) 2290 2291 #-- Add in other parameters 2292 if abs(min([min(yi) for yi in ddict['y'] if yi.size])) < 1e-20: 2293 ddict['ymin'] = 1e-20 2294 ddict.update(pars) 2295 ifn = Plotting2.plotCols(**ddict) 2296 fns.append(ifn) 2297 2298 #-- Print the plotted filenames, if any 2299 if join and fns[-1] and fns[-1][-3:] == 'pdf': 2300 DataIO.joinPdf(old=fns,new=fn+'.pdf',del_old=1) 2301 print('Heating and cooling rates plotted at:') 2302 print(fn+'.pdf') 2303 elif fns[-1]: 2304 print('Heating and cooling rates plotted at:') 2305 print('\n'.join(fns))

2306 2307 2308

2309 - def plotT(self,fn=None,iterations=[],cfg=None,**kwargs):

2310 2311 ''' 2312 Plot the temperature profiles of the different iterations. 2313 2314 @keyword fn: The filename and path of the plot (without extension). If 2315 default, a filename can be given in a cfg file, and if that 2316 is not the case, the plot is simply shown but not saved. 2317 2318 (default: None) 2319 @type fn: str 2320 @keyword iterations: The iteration indices to be plotted. If default, 2321 all iterations are plotted. The first and last 2322 iterations are always plotted. 2323 2324 (default: []) 2325 @type iterations: list 2326 @keyword cfg: The filename to the cfg file for this plot with extra 2327 settings. 2328 2329 (default: None) 2330 @type cfg: str 2331 @keyword kwargs: Additional keywords passed to the plotting method. 2332 Overwrites any keys given in the cfg file. 2333 2334 (default: {}) 2335 @type kwargs: dict 2336 2337 ''' 2338 2339 #-- Read the dict if available, and extract the filename. 2340 cfg_dict = Plotting2.readCfg(cfg) 2341 if cfg_dict.has_key('filename'): 2342 fn = cfg_dict.pop('filename') 2343 2344 #-- If no iterations are specified, plot all of them. 2345 if iterations: 2346 iterations = [i for i in iterations if i < self.i and i > 0] 2347 elif not iterations: 2348 iterations = range(1,self.i) 2349 2350 #-- Add first and last iteration 2351 iterations = [0] + iterations + [self.i] 2352 2353 #-- Set parameters. kwargs > cfg > locally defined 2354 pars = {'xlogscale': 1, 'ylogscale': 1, 'yaxis': 'T (K)', 2355 'xaxis': 'r (cm)', 'filename': fn, 2356 'keytags': [r'T$_\mathrm{d}$'] 2357 +['Iteration {}'.format(i) for i in iterations], 2358 'extension': 'pdf'} 2359 pars.update(cfg_dict) 2360 pars.update(kwargs) 2361 2362 #-- Set the T-profiles and plot 2363 y = [self.Td.eval()]+[self.T_iter[i].eval(inner_eps=self.inner_eps,\ 2364 warn=not self.inner) 2365 for i in iterations] 2366 pfn = Plotting2.plotCols(x=self.r,y=y,**pars) 2367 2368 #-- Print the plotted filename 2369 if pfn: 2370 print('The temperature profiles are plotted at:') 2371 print(pfn)

2372 2373 2374

2375 - def plotHCTerm(self,fn=None,iterations=[],dTsign='C',cfg=None,**kwargs):

2376 2377 ''' 2378 Plot the total heating and cooling term (excluding adiabatic cooling) 2379 calculated by calcT before the differential equation is solved. 2380 2381 This includes the density factor that enters, and so is essentially 2382 (H - C)/rho. The velocity does not enter here, since that is calculated 2383 explicitly in dTdr. Hence the y-axis units K/s. 2384 2385 @keyword fn: The filename and path of the plot (without extension). If 2386 default, a filename can be given in a cfg file, and if that 2387 is not the case, the plot is simply shown but not saved. 2388 2389 (default: None) 2390 @type fn: str 2391 @keyword iterations: The iteration indices to be plotted. If default, 2392 all iterations are plotted. The first and last 2393 iteration is always plotted. 2394 2395 (default: []) 2396 @type iterations: list 2397 @keyword dTsign: The dT type: heating (H) or cooling (C). Either the 2398 positive (C) or the negative (H) sum is plotted. 2399 2400 (default: 'C') 2401 @type dTsign: str 2402 @keyword cfg: The filename to the cfg file for this plot with extra 2403 settings. 2404 2405 (default: None) 2406 @type cfg: str 2407 @keyword kwargs: Additional keywords passed to the plotting method. 2408 Overwrites any keys given in the cfg file. 2409 2410 (default: {}) 2411 @type kwargs: dict 2412 2413 ''' 2414 2415 #-- Read the dict if available, and extract the filename. 2416 cfg_dict = Plotting2.readCfg(cfg) 2417 if cfg_dict.has_key('filename'): 2418 fn = cfg_dict.pop('filename') 2419 2420 #-- If no iterations are specified, plot all of them. 2421 if iterations: 2422 iterations = [i for i in iterations if i < self.i and i > 1] 2423 elif not iterations: 2424 iterations = range(2,self.i) 2425 2426 #-- Add first and last iteration 2427 iterations = [1] + iterations + [self.i] 2428 2429 #-- Set parameters. kwargs > cfg > locally defined 2430 pars = {'xlogscale': 1, 'ylogscale': 1, 2431 'yaxis': 'H - C (K s$^{-1}$)', 2432 'xaxis': 'r (cm)', 'filename': fn, 2433 'keytags': ['Iteration {}'.format(i) for i in iterations], 2434 'extension': 'pdf'} 2435 pars.update(cfg_dict) 2436 pars.update(kwargs) 2437 2438 #-- Set the H - C terms and plot 2439 dTsign = -1 if dTsign.upper() == 'H' else 1 2440 y = [self.rates[i]*dTsign for i in iterations] 2441 pfn = Plotting2.plotCols(x=self.r,y=y,**pars) 2442 2443 #-- Print the plotted filename 2444 if pfn: 2445 print('The H-C terms are plotted at:') 2446 print(pfn)

2447 2448 2449

2450 - def plotTexc(self,m,indices=None,llow=None,lup=None,fn=None,cfg=None,\ 2451 **kwargs):

2452 2453 ''' 2454 Plot the excitation temperature as a function of impact parameter for a 2455 selection of radiative transitions, given by either the transition 2456 indices or the upper and/or lower level indices. 2457 2458 Retrieves the excitation temperature from self.Texc. They are 2459 appropriate for the level populations used by the current iteration, and 2460 are updated when the level populations are updated. 2461 2462 At this time, plotting Texc for other iterations is not possible. 2463 2464 @param m: The molecule tag 2465 @type m: str 2466 2467 @keyword indices: The transition indices. Can be indirectly given 2468 through upper and/or lower level indices. If default 2469 and no llow/lup are given, all Texc are returned for 2470 this molecule. 2471 2472 (default: None) 2473 @type indices: list/array 2474 @keyword lup: The upper level indices used to extract the transition 2475 indices. If this or llow are defined, the keyword indices 2476 is overwritten. 2477 2478 (default: None) 2479 @type lup: list/array 2480 @keyword llow: The upper level indices used to extract the transition 2481 indices. If this or lup are defined, the keyword indices 2482 is overwritten. 2483 2484 (default: None) 2485 @type llow: list/array 2486 @keyword fn: The filename and path of the plot (without extension). If 2487 default, a filename can be given in a cfg file, and if that 2488 is not the case, the plot is simply shown but not saved. 2489 2490 (default: None) 2491 @type fn: str 2492 @keyword cfg: The filename to the cfg file for this plot with extra 2493 settings. 2494 2495 (default: None) 2496 @type cfg: str 2497 @keyword kwargs: Additional keywords passed to the plotting method. 2498 Overwrites any keys given in the cfg file. 2499 2500 (default: {}) 2501 @type kwargs: dict 2502 2503 2504 2505 ''' 2506 2507 #-- Read the dict if available, and extract the filename. 2508 cfg_dict = Plotting2.readCfg(cfg) 2509 if cfg_dict.has_key('filename'): 2510 fn = cfg_dict.pop('filename') 2511 2512 #-- Check which iteration to plot, use default if not give or invalid. 2513 if iteration is None or iteration > self.i: iteration = self.i 2514 2515 #-- Add the iteration if any but the last one is requested 2516 if fn and iteration != self.i: fn += '_iter{}'.format(self.i) 2517 2518 #-- If either lup or llow are given, get transition indices from them 2519 if not (lup is None and llow is None): 2520 if not lup is None: lup = Data.arrayify(lup) 2521 if not llow is None: llow = Data.arrayify(llow) 2522 indices = self.mol[m].getTI(lup=lup,llow=llow) 2523 2524 #-- Redefine lup and llow to match the given transitions. 2525 # Do it here because lup/llow are needed for plot keytags. 2526 lup = self.mol[m].getTUpper(indices) 2527 llow = self.mol[m].getTLower(indices) 2528 2529 #-- Get the impact parameter grid 2530 p = self.pop[m].getP() 2531 2532 #-- Get kinetic T profile and the excitation temperatures 2533 T = self.T.eval(p,inner_eps=self.inner_eps,warn=not self.inner) 2534 Texc = self.getTexc(m=m,indices=indices) 2535 2536 #-- Set parameters. kwargs > cfg > locally defined 2537 pars = {'xlogscale': 1, 'ylogscale': 0, 2538 'yaxis': 'T$_\mathrm{exc}$ (K)', 2539 'xaxis': 'p (cm)', 'filename': fn, 2540 'extension': 'pdf', 2541 'keytags': ['T$_\mathrm{kin}$']+ 2542 ['llup = {}, llow = {}'.format(nui,nli) 2543 for nui,nli in zip(lup,llow)]} 2544 pars.update(cfg_dict) 2545 pars.update(kwargs) 2546 2547 #-- Make the plot 2548 y = [T]+[Texc[i,:] for i in range(Texc.shape[0])] 2549 pfn = Plotting2.plotCols(x=p,y=y,**pars)

2550

Source Code for Module ComboCode.cc.modeling.physics.EnergyBalance