ComboCode.cc.modeling.objects.Star

1 # -*- coding: utf-8 -*- 2 3 """ 4 Module including functions for stellar parameters, and the STAR class and 5 its methods and attributes. 6 7 Author: R. Lombaert 8 9 """ 10 11 import types 12 from glob import glob 13 import os 14 from scipy import pi, log, sqrt 15 from scipy import array, exp, zeros 16 from scipy import integrate, linspace 17 from scipy import argmin,argmax, empty 18 from scipy.interpolate import interp1d 19 import operator 20 from numpy import savetxt 21 from astropy import units as u 22 from astropy import constants as cst 23 24 import cc.path 25 from cc.data import Data 26 from cc.tools.units import Equivalency as eq 27 from cc.tools.io import Database 28 from cc.tools.io import DataIO, Atmosphere 29 from cc.tools.numerical import Interpol 30 from cc.modeling.objects import Molecule 31 from cc.modeling.objects import Transition 32 from cc.modeling.tools import ColumnDensity 33 from cc.modeling.codes import MCMax 34 35

36 -def getStar(star_grid,modelid,idtype='GASTRONOOM'):

37 38 ''' 39 Grab a Star() object from a list of such objects, given a model id. 40 41 If no modelid is found, an empty list is returned. If Star() objects are 42 found (even only one), a list of them is returned. 43 44 Based on the cooling modelid. 45 46 @param star_grid: the Star() objects 47 @type star_grid: list[Star()] 48 @param modelid: the given modelid for which the selection is made. 49 @type modelid: string 50 51 @keyword idtype: The type of model id 52 53 (default: GASTRONOOM) 54 @type idtype: string 55 56 @return: The models matching the modelid 57 @rtype: list[Star()] 58 59 ''' 60 61 modelid, idtype = str(modelid), str(idtype) 62 return [s for s in star_grid if s['LAST_%s_MODEL'%idtype] == modelid]

63 64 65

66 -def makeStars(models,id_type,path,code):

67 68 ''' 69 Make a list of dummy Star() objects. 70 71 @param models: model_ids for the new models 72 @type models: list[string] 73 @param id_type: The type of id (PACS, GASTRONOOM, MCMAX) 74 @type id_type: string 75 @param path: Output folder in the code's home folder 76 @type path: string 77 @param code: The code (which is not necessarily equal to id_type, such as 78 for id_type == PACS) 79 @type code: string 80 81 @return: The parameter sets, mostly still empty! 82 @rtype: list[Star()] 83 84 ''' 85 86 extra_pars = dict([('path_'+code.lower(),path)]) 87 star_grid = [Star(example_star={'LAST_%s_MODEL'%id_type.upper():model},\ 88 **extra_pars) 89 for model in models] 90 return star_grid

91 92 93

94 -class Star(dict):

95 96 """ 97 Star class maintains information about a stellar model and its properties. 98 99 Inherits from dict. 100 101 """ 102 103 104

105 - def __init__(self,path_gastronoom='',path_mcmax='',extra_input=None,\ 106 example_star=dict(),print_check_t=1):

107 108 """ 109 Initiate an instance of the STAR class. 110 111 @keyword path_gastronoom: path in ~/GASTRoNOoM/ for modeling out/input 112 113 (default: None) 114 @type path_gastronoom: string 115 @keyword path_mcmax: the folder in ~/MCMax/ for modeling out/input 116 117 (default: None) 118 @type path_mcmax: string 119 @keyword example_star: if not None the STAR object is exact duplicate 120 of example_star. Can be a normal dictionary as 121 well. Paths are not copied and need to be given 122 explicitly. 123 124 (default: None) 125 @type example_star: dict or Star() 126 @keyword extra_input: extra input that you wish to add to the dict 127 128 (default: None) 129 @type extra_input: dict or Star() 130 @keyword print_check_t: Print the dust temperature check automatically. 131 Can still be done manually by running s.checkT 132 or asking for any T_DES_%s or R_DES_%s keyword 133 where %s is a dust species. 134 135 (default: 1) 136 @type print_check_t: bool 137 138 @return: STAR object in the shape of a dictionary which includes all 139 stellar data available, if None are passed for both options it 140 is an empty dictionary; if an example star is passed it has a 141 dict that is an exact duplicate of the example star's dict 142 @rtype: Star() 143 144 """ 145 146 super(Star, self).__init__(example_star) 147 if not extra_input is None: self.update(extra_input) 148 self.Rsun = cst.R_sun.cgs.value #in cm Harmanec & Prsa 2011 149 self.Msun = 1.98547e33 #in g Harmanec & Prsa 2011 150 self.Mearth = cst.M_earth.cgs.value # in g 151 self.Tsun = 5779.5747 #in K Harmanec & Psra 2011 152 self.Lsun = cst.L_sun.cgs.value #in erg/s 153 self.au = 149598.0e8 #in cm 154 self.c = cst.c.cgs.value #in cm/s 155 self.h = cst.h.cgs.value #in erg*s, Planck constant 156 self.k = cst.k_B.cgs.value #in erg/K, Boltzmann constant 157 self.sigma = cst.sigma_sb.cgs.value #in erg/cm^2/s/K^4 = g / (K^4 s^3) Stefan_boltzmann constant 158 self.mh = cst.m_p.cgs.value #in g, mass hydrogen atom 159 self.G = cst.G.cgs.value # in cm^3 g^-1 s^-2 160 161 self.path_gastronoom = path_gastronoom 162 self.path_mcmax = path_mcmax 163 self.print_check_t = print_check_t 164 165 #-- Convenience paths 166 cc.path.mout = os.path.join(cc.path.mcmax,self.path_mcmax) 167 cc.path.gout = os.path.join(cc.path.gastronoom,self.path_gastronoom) 168 169 self.dust_list = None

170 171 172

173 - def __getitem__(self,key):

174 175 """ 176 Overriding the standard dictionary __getitem__ method. 177 178 @param key: Star()[key] where key is a string for which a corresponding 179 dictionary value is searched. If the key is not present in 180 the dictionary, an attempt is made to calculate it from 181 already present data; if it fails a KeyError is still 182 raised. 183 @type key: string 184 185 @return: The value from the Star() dict for key 186 @rtype: any 187 188 """ 189 190 if not self.has_key(key): 191 self.missingInput(key) 192 return super(Star,self).__getitem__(key) 193 elif super(Star,self).__getitem__(key) == '%': 194 del self[key] 195 self.missingInput(key) 196 value = super(Star,self).__getitem__(key) 197 self[key] = '%' 198 return value 199 else: 200 return super(Star,self).__getitem__(key)

201 202 203

204 - def __cmp__(self,star):

205 206 """ 207 Overriding the standard dictionary __cmp__ method. 208 209 A parameter set (dictionary of any type) is compared with this instance 210 of Star(). 211 212 An attempt is made to create keys with values in each dict, if the 213 other has keys that are not present in the first. If this fails, False 214 is returned. 215 216 @param star: A different parameter set. 217 @type star: dict or Star() 218 219 @return: The comparison between this object and star 220 @rtype: bool 221 222 """ 223 224 try: 225 all_keys = set(self.keys() + star.keys()) 226 for k in all_keys: 227 if not self.has_key(): 228 self[k] 229 if not star.has_key(): 230 star[k] 231 #if len(self) > len(star): 232 #changed_keys = [star[k] 233 #for k in self.keys() 234 #if not star.has_key(k)] 235 #print "Initialized keys in STAR2 from STAR1 == STAR2 comparison : \n" + str(changed_keys) 236 #if len(self) < len(star): 237 #changed_keys = [self[k] for k in star.keys() if not self.has_key(k)] 238 #print "Initialized keys in STAR1 from STAR1 == STAR2 comparison : \n" + str(changed_keys) 239 except KeyError: 240 print 'Comparison error: Either STAR1 or STAR2 contains a key ' + \ 241 'that cannot be initialized for the other.' 242 print 'Both STAR instances are considered to be unequal.' 243 finally: 244 if isinstance(star, super(Star)): 245 return cmp(super(Star,self), super(Star,star)) 246 else: 247 return cmp(super(Star,self), star)

248 249 250

251 - def readDustProperties(self):

252 253 ''' 254 When requesting the list of dust species, all dust species with nonzero 255 abundance are returned. The list is created here, and remembered. 256 257 When this is done, the dust properties are also read for those species. 258 259 The info is saved in self.dust. 260 261 The dust list itself can be accessed through getDustList(). This method 262 also calls the readDustProperties method to ensure all info is always 263 available. self.dust_list has a fixed order of appearance, according to 264 Dust.dat. Important for, e.g., MCMax.py. 265 266 ''' 267 268 #-- If already read, don't read again 269 if not self.dust_list is None: return 270 271 #-- Read the info 272 dust_info = dict() 273 for k,s in zip(['SPECIES_SHORT','SPEC_DENS','MOLAR_WEIGHT','T_DES',\ 274 'T_DESA','T_DESB','PART_FILE'],\ 275 ['species','sd','molar','tdes','tdesa','tdesb','fn']): 276 dust_info[s] = DataIO.getInputData(keyword=k,filename='Dust.dat') 277 278 #-- Check which dust species are requested and exist in Dust.dat 279 dust_list = [species 280 for species in dust_info['species'] 281 if self.has_key('A_' + species)] 282 dust_list = [species 283 for species in dust_list 284 if float(self['A_' + species]) != 0] 285 286 #-- order in which the species appear is fixed according to the 287 # order of the species in the Dust.dat input file. Note that this 288 # order does not matter for the database. 289 # rgrain species are put first, following the dust_list order. 290 # Then the rest. This is the same as what is used in MCMax.py. 291 dust_list = sorted(dust_list,\ 292 key=lambda x: (not self.has_key('RGRAIN_%s'%x),\ 293 dust_list.index(x))) 294 self.dust_list = tuple(dust_list) 295 296 print '==========' 297 print 'Dust species taken into account during modeling: %s'\ 298 %(', '.join(self.dust_list)) 299 300 #-- Save all relevant info based on which species are requested here 301 self.dust = dict() 302 for species in self.dust_list: 303 index = dust_info['species'].index(species) 304 self.dust[species] = dict() 305 for key in ['sd','molar','tdes','tdesa','tdesb','fn']: 306 self.dust[species][key] = dust_info[key][index]

307 308 309

310 - def addCoolingPars(self):

311 312 ''' 313 Add Star parameters from the cooling database. 314 315 Any existing parameters are overwritten! 316 317 ''' 318 319 #-- Check if the Star instance actually has a model id available. 320 if not self['LAST_GASTRONOOM_MODEL']: 321 print('WARNING! Star() does not have a cooling model id. ' + \ 322 'Cannot add parameters from the database.') 323 return 324 325 #-- Set the path, and extract the parameters dictionary from the db 326 cooling_path = os.path.join(cc.path.gout,\ 327 'GASTRoNOoM_cooling_models.db') 328 cool_db = Database.Database(cooling_path) 329 db_dict = cool_db[self['LAST_GASTRONOOM_MODEL']] 330 331 #-- Only set a small subset of keys for now. These will always be in the 332 # cooling dictionary. 333 cool_keys = ['T_STAR','R_STAR','TEMDUST_FILENAME','MDOT_GAS'] 334 cool_dict = {k: DataIO.convertFloat(db_dict[k]) for k in cool_keys} 335 336 #-- Set the correct units for the stellar radius. 337 cool_dict['R_STAR'] = cool_dict['R_STAR']/self.Rsun 338 self.update(cool_dict)

339 340 341

342 - def writeDensity(self):

343 344 ''' 345 Write dust mass density and n(h2) profile (in Rstar). 346 347 Only if MCMax or GASTRoNOoM model_id is available! 348 349 ''' 350 351 mcmid = self['LAST_MCMAX_MODEL'] 352 gasid = self['LAST_GASTRONOOM_MODEL'] 353 if mcmid: 354 rad = self.getDustRad(unit='rstar') 355 dens = self.getDustDensity() 356 fn = os.path.join(cc.path.mout,'models',mcmid,\ 357 'density_profile_%s.dat'%mcmid) 358 DataIO.writeCols(fn,[rad,dens]) 359 if gasid: 360 rad = self.getGasRad(unit='rstar') 361 nh2 = self.getGasNumberDensity() 362 fn = os.path.join(cc.path.gout,'models',gasid,\ 363 'nh2_density_profile_%s.dat'%gasid) 364 DataIO.writeCols(fn,[rad,nh2])

365 366

367 - def readKappas(self):

368 369 ''' 370 Read the kappas.dat file of an MCMax model. 371 372 ''' 373 374 opas = DataIO.readCols(os.path.join(cc.path.mout,'models',\ 375 self['LAST_MCMAX_MODEL'],\ 376 'kappas.dat')) 377 return opas.pop(0),opas

378 379 380

381 - def removeMutableMCMax(self,mutable,var_pars):

382 383 """ 384 Remove mutable parameters after an MCMax run. 385 386 This method is accessed by the ComboCode package whenever necessary. The 387 user should not have to do this. 388 389 @param mutable: mutable keywords. Listed in 390 aux/Mutable_Parameters_MCMax.dat. 391 392 @type mutable: list[string] 393 @param var_pars: parameters that are varied during gridding, these will 394 not be removed and are kept constant throughout the 395 iteration 396 @type var_pars: list[string] 397 398 """ 399 400 #- remove keys which should be changed by output of new mcmax model, 401 #- but only the mutable input!!! 402 for key in self.keys(): 403 if key in mutable \ 404 + ['R_MAX_' + species for species in self.getDustList()]\ 405 + ['T_DES_' + species for species in self.getDustList()]\ 406 + ['R_DES_' + species for species in self.getDustList()]\ 407 and key not in var_pars: 408 del self[key] 409 410 #- Check the effective destruction temperature of every species, and 411 #- see if max and min T's are as requested. 412 if self.print_check_t: 413 self.checkT()

414 415 416

417 - def removeMutableGastronoom(self,mutable,var_pars):

418 419 """ 420 Remove mutable parameters after a GASTRoNOoM run. 421 422 This method is accessed by the ComboCode package whenever necessary. The 423 user should not have to do this. 424 425 @param mutable: mutable keywords. Listed in 426 aux/Mutable_Parameters_GASTRoNOoM.dat. 427 428 @param mutable: mutable parameters 429 @type mutable: list[string] 430 @param var_pars: parameters that are varied during gridding, these will 431 not be removed and are kept constant throughout the 432 iteration 433 @type var_pars: list[string] 434 435 """ 436 437 for key in self.keys(): 438 if key in mutable and key not in var_pars: 439 del self[key]

440 441 442

443 - def updateMolecules(self,parlist):

444 445 ''' 446 Update variable information in the molecule instances of this star. 447 448 @param parlist: parameters that have to be updated. 449 @type parlist: list[string] 450 451 ''' 452 453 for molec in self['GAS_LIST']: 454 molec.updateParameters(pardict=dict([(k,self[k]) 455 for k in parlist]))

456 457 458

459 - def normalizeDustAbundances(self):

460 461 """ 462 Normalize the dust abundances such that they add up to a total of 1. 463 464 If the MRN_DUST keyword for MCMax is nonzero, all nonzero abundances 465 are set to 1. The abundance given in the inputfile does not matter in 466 this case. 467 468 """ 469 470 abun_ori = [self['A_%s'%sp] for sp in self.getDustList()] 471 self['A_DUST_ORIGINAL'] = abun_ori 472 if int(self['MRN_DUST']): 473 self['A_NO_NORM'] = 1 474 print 'WARNING! Take care when extracting output in MCMax using '+\ 475 'these scripts, if MRN_DUST == 1! Especially if some ' + \ 476 'abundances are set manually and some according to MRN: ' + \ 477 'these are not normalized to 1, since this depends on the '+\ 478 'MRN distributed dust species.' 479 for sp in self.getDustList(): 480 mrn_count = 0 481 if self['MRN_NGRAINS'] != len(self.getDustList()) \ 482 and self.has_key('RGRAIN_%s'%sp): 483 self.__setitem__('A_%s'%sp,2) 484 mrn_count += 1 485 elif self['MRN_NGRAINS'] == len(self.getDustList()): 486 self.__setitem__('A_%s'%sp,2) 487 mrn_count += 1 488 if mrn_count != self['MRN_NGRAINS']: 489 raise IOError('MRN_NGRAINS not equal to amount of RGRAIN_sp keywords.') 490 total = sum(abun_ori) 491 if not int(self['A_NO_NORM']) and '%.3f'%total != '1.000': 492 print 'Normalizing dust abundances to 1, from a total of %f.'%total 493 abun_new = [a/total for a in abun_ori] 494 print ', '.join(['%.2f'%a for a in abun_ori]), ' is changed to ', \ 495 ', '.join(['%.2f'%a for a in abun_new]), ' for ', \ 496 ', '.join(self.getDustList()), '.' 497 [self.__setitem__('A_%s'%sp,a) for a,sp in zip(abun_new,\ 498 self.getDustList())]

499 500 501

502 - def __addLineList(self):

503 504 """ 505 Take molecular transitions from a line list and wavelength range. 506 507 Based on the GASTRoNOoM radiat and indices data files. See Molecule.py 508 for more info. 509 510 """ 511 512 gas_list = [] 513 if type(self['LS_TELESCOPE']) is types.StringType: 514 self['LS_TELESCOPE'] = [self['LS_TELESCOPE']] 515 if not self['LS_NO_VIB']: 516 self['LS_NO_VIB'] = [] 517 elif type(self['LS_NO_VIB']) is types.StringType: 518 self['LS_NO_VIB'] = [self['LS_NO_VIB']] 519 ctrl = [tr.getInputString(include_nquad=0) for tr in self['GAS_LINES']] 520 for molec in self['GAS_LIST']: 521 for telescope in self['LS_TELESCOPE']: 522 if telescope == 'PACS': 523 ls_min = 50 524 ls_max = 200 525 elif telescope == 'SPIRE': 526 ls_min = 180 527 ls_max = 700 528 extra_pars = {'fraction_tau_step':self['FRACTION_TAU_STEP'],\ 529 'min_tau_step':self['MIN_TAU_STEP'],\ 530 'write_intensities':self['WRITE_INTENSITIES'],\ 531 'tau_max':self['TAU_MAX'],\ 532 'n_quad':self['N_QUAD'],\ 533 'offset':self['LS_OFFSET'],\ 534 'tau_min':self['TAU_MIN'],\ 535 'check_tau_step':self['CHECK_TAU_STEP']} 536 nl = Transition.makeTransitionsFromRadiat(molec=molec,\ 537 telescope=telescope,ls_min=ls_min,ls_max=ls_max,\ 538 path_gastronoom=self.path_gastronoom,\ 539 no_vib=molec.molecule in self['LS_NO_VIB'],\ 540 ls_unit='micron',**extra_pars) 541 #-- exclude transitions if they are already included. 542 # This is to avoid adding a double with a different n_quad, in 543 # case it was included manually. Can use GAS_LINES key, as it 544 # contains both manual and other lines, but the latter are 545 # also assigned self['N_QUAD'] anyway 546 nl = [tr for tr in nl 547 if tr.getInputString(include_nquad=0) not in ctrl] 548 gas_list.extend(nl) 549 self['GAS_LINES'].extend(gas_list)

550 551 552

553 - def checkT(self):

554 555 """ 556 Search input list for minimum temperature. 557 558 Method prints the actual minimum T for which the model was calculated. 559 560 Note that the density drops to zero gradually and that the criterium 561 has to be sudden change of slope. Check criterium if the printed T is 562 not good enough as determination of max radius IS correct. 563 564 """ 565 566 coldens = ColumnDensity.ColumnDensity(self) 567 self.calcT_INNER_DUST() 568 for index,species in enumerate(self.getDustList()): 569 self['T_DES_%s'%species] = coldens.t_des[species] 570 self['R_DES_%s'%species] = coldens.r_des[species]\ 571 /self.Rsun/self['R_STAR'] 572 print 'The EFFECTIVE maximum temperature for species %s '%species+\ 573 'is %.2f K, at radius %.2f R_STAR.'\ 574 %(self['T_DES_%s'%species],self['R_DES_%s'%species]) 575 576 species_list_min = [species 577 for species in self.getDustList() 578 if self.has_key('T_MIN_%s'%species) \ 579 or self.has_key('R_MAX_%s'%species)] 580 for species in species_list_min: 581 print 'The EFFECTIVE minimum temperature for species'+\ 582 ' %s is %.2f K at maximum radius %.2f R_STAR.'\ 583 %(species,coldens.t_min[species],\ 584 coldens.r_max[species]/self.Rsun/self['R_STAR']) 585 if self.has_key('T_MIN_%s'%species): 586 print 'The REQUESTED minimum temperature for species '+\ 587 '%s is %.2f K.'%(species,self['T_MIN_%s'%species]) 588 if self.has_key('R_MAX_%s'%species): 589 print 'The REQUESTED maximum radius for species'+\ 590 '%s is %.2f R_STAR.'%(species,self['R_MAX_%s'%species]) 591 print 'The EFFECTIVE outer radius of the shell is %.2f R_STAR.'\ 592 %(coldens.r_outer/self.Rsun/self['R_STAR']) 593 print 'Note that if R_MAX is ~ the effective outer radius, the ' +\ 594 'requested minimum temperature may not be reached.' 595 596 #- No point in keeping zeroes around for T_DES or T_MIN 597 for species in self.getDustList(): 598 for par in ('T_DES_' + species, 'T_MIN_' + species): 599 if self.has_key(par): 600 if not float(self[par]): 601 del self[par]

602 603 604

605 - def getFullNameMolecule(self,short_name):

606 607 ''' 608 Get the full name of a molecule, based on it's short name, 609 if it is present in the GAS_LIST. 610 611 @return: Name the name. None if not available. 612 @rtype: string 613 614 ''' 615 616 molecule = [molec.molecule_full 617 for molec in self['GAS_LIST'] 618 if molec.molecule == short_name] 619 if not molecule: 620 return None 621 #- Return first only, if multiple there's multiple requested molecules 622 #- of the same type (fi different abundance) 623 else: 624 return molecule[0]

625 626 627

628 - def getShortNameMolecule(self,full_name):

629 630 ''' 631 Get the short name of a molecule, based on it's full name, 632 if it is present in the GAS_LIST. 633 634 @param full_name: The full name of a molecule from Molecule.dat 635 @type full_name: string 636 637 @return: None if not available, otherwise the short hand name. 638 @rtype: string 639 640 ''' 641 642 molecule = [molec.molecule 643 for molec in self['GAS_LIST'] 644 if molec.molecule_full == full_name] 645 if not molecule: 646 return None 647 #- Return first only, if multiple there's multiple requested molecules 648 #- of the same type (fi different abundance) 649 else: 650 return molecule[0]

651 652 653

654 - def getMolecule(self,molec_name):

655 656 ''' 657 Get a Molecule() object based on the molecule name. 658 659 A Star() object always has only one version of one molecule. 660 661 @param molec_name: short name of the molecule 662 @type molec_name: string 663 664 @return: The molecule 665 @rtype: Molecule() 666 667 ''' 668 669 try: 670 return [molec 671 for molec in self['GAS_LIST'] 672 if molec.molecule == molec_name][0] 673 except IndexError: 674 return None

675 676 677

678 - def getTransition(self,sample):

679 680 ''' 681 Return a Transition() object that has the same parameters as sample. 682 683 The comparison is done based on the str representation of the trans. 684 This excludes the dictionary entries of the transition! 685 686 The actual model ids or data are not included in this comparison! 687 688 None is returned if no match is found. 689 690 @param sample: A sample transition to be cross referenced with the 691 transitions in this Star() object. If a match is found, 692 it is returned. 693 @type sample: Transition() 694 695 @return: If a match is found, this transition is returned. 696 @rtype: Transition() 697 698 ''' 699 700 try: 701 i = self['GAS_LINES'].index(sample) 702 return self['GAS_LINES'][i] 703 except ValueError: 704 return None

705 706 707

708 - def getTransList(self,**kwargs):

709 710 ''' 711 Return a list of (transmodelid, molecmodelid, dictionary) for every 712 transition in the Star model. 713 714 ''' 715 716 trl = Transition.extractTransFromStars([self],**kwargs) 717 trl_info = [(trans.getModelId(),\ 718 trans.molecule.getModelId(),\ 719 trans.makeDict()) 720 for trans in trl] 721 return trl_info

722 723 724

725 - def getTransitions(self,molec):

726 727 ''' 728 Return a list of all transitions associated with a single molecule. 729 730 @param molec: the shorthand notation of the molecule 731 @type molec: string 732 733 @return: All transitions for this molecule 734 @rtype: list[Transition()] 735 736 ''' 737 738 return [trans 739 for trans in self['GAS_LINES'] 740 if trans.molecule.molecule==molec]

741 742 743

744 - def getDustFn(self,species=''):

745 746 ''' 747 Return the dust output filename for density and temperature. 748 749 You can choose the species for which the output file is retrieved. 750 Unless thermal contact for this model is on. 751 752 @keyword ftype: The species for which dust info is read. Default if the 753 average profiles are requested. Can be any species in 754 Dust.dat as long as the model calculation included it. 755 756 (default: '') 757 @type ftype: str 758 759 @return: The filename of the requested dust output file. 760 @rtype: str 761 762 ''' 763 764 fp = os.path.join(cc.path.mout,'models',self['LAST_MCMAX_MODEL']) 765 766 #- if T_CONTACT: no specific species denstemp files, 767 #- so denstemp.dat is taken 768 if species and not self['T_CONTACT']: 769 ispecies = self.getDustList().index(species) 770 fn = os.path.join(fp,'denstempP%.2i.dat'%(ispecies+1)) 771 else: 772 fn = os.path.join(fp,'denstemp.dat') 773 774 return fn

775 776 777

778 - def getDustRad(self,species='',unit='cm'):

779 780 ''' 781 Return the dust radial grid in cm, au or Rstar. 782 783 @keyword species: The dust species for which to return the grid. If 784 default or if T_CONTACT is on, the file is denstemp 785 otherwise it is the species specific file. 786 787 (default: '') 788 @type species: str 789 @keyword unit: The unit of the returned grid. Can be 'cm','rstar','au', 790 'm' 791 792 (default: 'cm') 793 @type unit: str 794 795 @return: array giving the radial grid. 796 @rtype: array 797 798 ''' 799 800 if not self['LAST_MCMAX_MODEL']: return empty(0) 801 802 fn = self.getDustFn(species) 803 rad = array(DataIO.getKeyData(incr=int(self['NRAD']),\ 804 keyword='RADIUS',filename=fn)) 805 806 unit = str(unit).lower() 807 if unit == 'au': 808 rad = rad/self.au 809 elif unit == 'rstar': 810 rad = rad/self.Rsun/self['R_STAR'] 811 elif unit == 'm': 812 rad = rad*10**-2 813 814 return rad

815 816

817 - def getDustTheta(self,species=''):

818 819 ''' 820 Return the dust theta grid in radians. 821 822 This is the angular grid goin from pole at zero radians to equator at 823 pi/2 radians. 824 825 @keyword species: The dust species for which to return the grid. If 826 default or if T_CONTACT is on, the file is denstemp 827 otherwise it is the species specific file. 828 829 (default: '') 830 @type species: str 831 832 @return: array giving the theta grid. 833 @rtype: array 834 835 ''' 836 837 if not self['LAST_MCMAX_MODEL']: return empty(0) 838 839 fn = self.getDustFn(species) 840 theta = array(DataIO.getKeyData(incr=int(self['NTHETA']),\ 841 keyword='THETA',filename=fn)) 842 return theta

843 844 845

846 - def getDustDensity(self,species='',avg_theta=1):

847 848 ''' 849 Return the dust density profile from the file made for MCMax. 850 851 self.getDustRad returns the radial grid associated with this profile. 852 self.getDustTheta returns the angular grid associated with this profile 853 if angular averaging is off. 854 855 An empty array is returned if the model does not exist. 856 857 Note that MCMax is a 2D code. By default, the theta coordinate (angle 858 pole-equator) is averaged out. This can be turned off, in which case 859 the full density grid is returned giving rho at (r_0,t_0),...,(r_0,t_n), 860 (r_1,t_0),...,(r_1,t_n),...,(r_n,t_0),...,(r_n,t_n). 861 862 @keyword species: The dust species for which to return the grid. If 863 default or if T_CONTACT is on, the file is denstemp 864 otherwise it is the species specific file. 865 866 (default: '') 867 @type species: str 868 @keyword avg_theta: Average out the angular grid in the density 869 distribution. At a given radial point, the densities 870 at all thetas are averaged. On by default. 871 872 (default: 1) 873 @type avg_theta: bool 874 875 @return: the density profile (g/cm3) 876 @rtype: array 877 878 ''' 879 880 if not self['LAST_MCMAX_MODEL']: return empty(0) 881 882 #-- Read the dust density profile and reduce the array by averaging 883 # over the theta coordinate, if requested. 884 fn = self.getDustFn(species) 885 incr = int(self['NRAD'])*int(self['NTHETA']) 886 dens_ori = DataIO.getKeyData(filename=fn,incr=incr,\ 887 keyword='DENSITY') 888 if avg_theta: dens = Data.reduceArray(dens_ori,self['NTHETA']) 889 else: dens = dens_ori 890 891 return dens

892 893 894

895 - def getDustTemperature(self,add_key=0,species='',avg_theta=1):

896 897 ''' 898 Return the dust temperature profile from MCMax model. 899 900 self.getDustRad returns the radial grid associated with this profile. 901 self.getDustTheta returns the angular grid associated with this profile 902 if angular averaging is off. 903 904 An empty array is returned if the model does not exist. 905 906 The total dust temperature without separate components for the 907 different dust species is returned if no species is requested or if 908 thermal contact is on. 909 910 Note that MCMax is a 2D code. By default, the theta coordinate (angle 911 pole-equator) is averaged out. This can be turned off, in which case 912 the full temperature grid is returned giving T at (r_0,t_0),..., 913 (r_0,t_n),(r_1,t_0),...,(r_1,t_n),...,(r_n,t_0),...,(r_n,t_n). 914 915 @keyword species: The dust species for which to return the grid. If 916 default or if T_CONTACT is on, the file is denstemp 917 otherwise it is the species specific file. 918 919 (default: '') 920 @type species: str 921 @keyword add_key: Add a key for a legend to the ouput as third tuple 922 element. 923 924 (default: 0) 925 @type add_key: bool 926 @keyword avg_theta: Average out the angular grid in the temperature 927 distribution. At a given radial point, the temps 928 at all thetas are averaged. On by default. 929 930 (default: 1) 931 @type avg_theta: bool 932 933 @return: The temperature profile (K) as well as a key, if requested. 934 @rtype: array or (array,string) 935 936 ''' 937 938 if not self['LAST_MCMAX_MODEL']: 939 return add_key and (empty(0),'') or empty(0) 940 941 #-- Read the dust temperature profile and reduce the array by averaging 942 # over the theta coordinate, if requested. 943 fn = self.getDustFn(species) 944 incr = int(self['NRAD'])*int(self['NTHETA']) 945 temp_ori = DataIO.getKeyData(incr=incr,keyword='TEMPERATURE',\ 946 filename=fn) 947 if avg_theta: temp = Data.reduceArray(temp_ori,self['NTHETA']) 948 else: temp = temp_ori 949 950 if add_key: 951 key = '$T_{\mathrm{d, avg}}$ for %s'\ 952 %self['LAST_MCMAX_MODEL'].replace('_','\_') 953 return temp,key 954 else: 955 return temp

956 957 958

959 - def getGasNumberDensity(self,molecule=False,**kwargs):

960 961 ''' 962 Give the n(h2) number density profile of the gas read from a GASTRoNOoM 963 model. 964 965 Additional input keywords for self.getCoolFn() can be passed along. 966 967 An empty is returned in case no model is available. 968 969 @keyword molecule: Return the molecular number density instead of the 970 H_2 number density. 971 972 (default: False) 973 @type molecule: bool 974 975 @return: The number density (in cm-3) profile of n(H_2) 976 @rtype: array 977 978 ''' 979 980 if not self['LAST_GASTRONOOM_MODEL']: return empty(0) 981 fgr_file = self.getCoolFn(**kwargs) 982 983 kws = dict() 984 if molecule: 985 kws['keyword'] = 'N(MOLEC)' 986 kws['key_index'] = 8 987 else: 988 kws['keyword'] = 'N(H2)' 989 990 nmol = DataIO.getGastronoomOutput(filename=fgr_file,return_array=1,\ 991 **kws) 992 993 return nmol

994 995 996

997 - def getGasVelocity(self,**kwargs):

998 999 ''' 1000 Give the velocity profile of the gas read from a GASTRoNOoM model. 1001 1002 Additional input keywords for self.getCoolFn() can be passed along. 1003 1004 An empty is returned in case no model is available. 1005 1006 @return: The velocity (in cm/s) profile 1007 @rtype: array 1008 1009 ''' 1010 1011 if not self['LAST_GASTRONOOM_MODEL']: return empty(0) 1012 fgr_file = self.getCoolFn(**kwargs) 1013 vel = DataIO.getGastronoomOutput(filename=fgr_file,keyword='VEL',\ 1014 return_array=1) 1015 return vel

1016 1017 1018

1019 - def getGasTemperature(self,**kwargs):

1020 1021 ''' 1022 Give the temperature profile of the gas read from a GASTRoNOoM model. 1023 1024 Additional input keywords for self.getCoolFn() can be passed along. 1025 1026 An empty is returned in case no model is available. 1027 1028 @return: The temperature (in K) profile 1029 @rtype: array 1030 1031 ''' 1032 1033 if not self['LAST_GASTRONOOM_MODEL']: return empty(0) 1034 fgr_file = self.getCoolFn(**kwargs) 1035 temp = DataIO.getGastronoomOutput(filename=fgr_file,keyword='TEMP',\ 1036 return_array=1) 1037 return temp

1038 1039 1040

1041 - def getGasRad(self,unit='cm',ftype='fgr_all',**kwargs):

1042 1043 ''' 1044 Return the gas radial grid in cm, AU or Rstar. 1045 1046 Additional input keywords for self.getCoolFn() can be passed along. 1047 1048 An empty is returned in case no model is available. 1049 1050 @keyword unit: The unit of the returned grid. Can be 'cm','rstar','au', 1051 'm' 1052 1053 (default: 'cm') 1054 @type unit: str 1055 @keyword ftype: The cooling output file type. Either '1', '2', '3', 1056 'fgr', 'fgr_all', or 'rate'. 1057 1058 (default: 'fgr_all') 1059 @type ftype: str 1060 1061 @return: the radial grid 1062 @rtype: array 1063 1064 ''' 1065 1066 if not self['LAST_GASTRONOOM_MODEL']: return empty(0) 1067 1068 #-- coolrate***.dat does not give radii. 1069 if ftype == 'rate': return empty(0) 1070 1071 unit = str(unit).lower() 1072 fgr_file = self.getCoolFn(ftype=ftype,**kwargs) 1073 rad = DataIO.getGastronoomOutput(filename=fgr_file,keyword='RADIUS',\ 1074 return_array=1) 1075 #-- fgr_all gives radius in cm. Others in rstar. Convert others to cm 1076 if ftype != 'fgr_all': 1077 rad = rad*self['R_STAR']*self.Rsun 1078 1079 if unit == 'au': 1080 rad = rad/self.au 1081 elif unit == 'rstar': 1082 rad = rad/self.Rsun/self['R_STAR'] 1083 elif unit == 'm': 1084 rad = rad*10**-2 1085 return rad

1086 1087 1088

1089 - def getCoolFn(self,ftype='fgr_all',mstr='',modelid=''):

1090 1091 ''' 1092 Return the cooling output filename. 1093 1094 You can define the type of cooling file you want, as well as an 1095 additional identification string for the molecule/sampling. 1096 1097 @keyword ftype: The cooling output file type. Either '1', '2', '3', 1098 'fgr', 'fgr_all', or 'rate'. 1099 1100 (default: 'fgr_all') 1101 @type ftype: str 1102 @keyword mstr: The additional identication string. Not applicable to 1103 'fgr' or 'fgr_all'. Can be any molecule, or 'sampling'. 1104 File must exist to be used further! 1105 1106 (default: '') 1107 @type mstr: str 1108 @keyword modelid: The model id to be used. If default, the cooling id 1109 is used. If defined, it could refer to an id 1110 different from the cooling id such as an mline id. 1111 1112 (default: '') 1113 @type modelid: str 1114 1115 @return: The filename of the requested cooling output file. 1116 @rtype: str 1117 1118 ''' 1119 1120 if not modelid: 1121 modelid = self['LAST_GASTRONOOM_MODEL'] 1122 ftype = str(ftype) 1123 if ftype == 'fgr' or ftype == 'fgr_all': 1124 mstr = '' 1125 if mstr: 1126 mstr = '_' + mstr 1127 fn = os.path.join(cc.path.gout,'models',modelid,\ 1128 'cool%s%s%s.dat'%(ftype,modelid,mstr)) 1129 return fn

1130 1131 1132

1133 - def getWindDensity(self,denstype='gas'):

1134 1135 ''' 1136 Read the n_h2 number density profile and calculate the total gas or dust 1137 density profile from that. This requires a GASTRoNOoM cooling model to 1138 have been calculated. 1139 1140 The gas density is simply calculated by multiplying with the mass of 1141 molecular hydrogen. 1142 1143 The dust density is calculated by multiplying the H_2 density with the 1144 dust-to-gas ratio. The radial dependence of the velocity profiles is 1145 taken into account. Not the radial dependence of the mass-loss rate, but 1146 CC assumes the d2g ratio remains constant in case of variable mass-loss 1147 rate. 1148 1149 @keyword denstype: The type of density profile, either 'dust' or 'gas'. 1150 1151 (default: 'gas') 1152 @type denstype: str 1153 1154 @return: The density profile (g/cm3) 1155 @rtype: array 1156 1157 ''' 1158 1159 denstype = denstype.lower() 1160 if denstype not in ['gas','dust']: 1161 print "Density type not known. Must be gas or dust." 1162 return 1163 #-- GASTRoNOoM calculates smoother density profiles than 1164 # this formula ever can accomplish 1165 #dens = float(self['MDOT_DUST'])*self.Msun\ 1166 # /((vg+drift)*4.*pi*rad**2.*self.year) 1167 1168 #-- Use the H2 density profile to take into account any 1169 # type of variable mass loss (including exponents in 1170 # r_points_mass_loss. 1171 nh2 = self.getGasNumberDensity() 1172 gas_dens = nh2*self.mh*2. 1173 1174 if denstype == 'gas': return gas_dens 1175 1176 #-- The gas velocity profile 1177 vg = self.getGasVelocity() 1178 #-- The drift profile, corrected for the average grain size 1179 drift = self.getAverageDrift() 1180 1181 #-- Calc dust density based on md/mg instead of d2g to take 1182 # into account velocity profiles instead of terminal vels 1183 dust_dens = gas_dens*self['MDOT_DUST']/self['MDOT_GAS']*vg/(vg+drift) 1184 1185 return dust_dens

1186 1187 1188

1189 - def getAverageDrift(self):

1190 1191 ''' 1192 Return an array with the average drift velocity as a function of 1193 radius, from coolfgr_all, in cm/s. 1194 1195 ''' 1196 1197 inputfile = self.getCoolFn(ftype='fgr_all') 1198 drift = DataIO.getGastronoomOutput(inputfile,keyword='VDRIFT') 1199 opa_gs_max = 2.5e-1 1200 opa_gs_min = 5.0e-3 1201 return array(drift)/sqrt(0.25)*1.25\ 1202 *(opa_gs_max**(-2.)-opa_gs_min**(-2.))\ 1203 /(opa_gs_max**(-2.5)-opa_gs_min**(-2.5))

1204 1205 1206

1207 - def getOpticalDepth(self,wavelength=0):

1208 1209 ''' 1210 Calculate the optical depth. 1211 1212 If wavelength keyword is given, tau at wavelength is returned. 1213 1214 Otherwise, the full wavelength array is returned. 1215 1216 @keyword wavelength: the wavelength in micron. If 0, the whole 1217 wavelength array is returned. 1218 1219 (default: 0) 1220 @type wavelength: float 1221 1222 @return: The optical depth at requested wavelength or the full 1223 wavelength and optical depth arrays 1224 @rtype: float or (array,array) 1225 1226 ''' 1227 1228 wavelength = float(wavelength) 1229 rad = self.getDustRad() 1230 dens = self.getDustDensity() 1231 wave_list,kappas = self.getWeightedKappas() 1232 if wavelength: 1233 wave_index = argmin(abs(wave_list-wavelength)) 1234 return integrate.trapz(y=dens*kappas[wave_index],x=rad) 1235 else: 1236 return (wave_list,array([integrate.trapz(y=dens*kappas[i],x=rad) 1237 for i in xrange(len(wave_list))]))

1238 1239 1240

1241 - def getWeightedKappas(self):

1242 1243 ''' 1244 Return the wavelength and kappas weighted with their respective dust 1245 mass fractions. 1246 1247 Typically you only want the absorption coefficients because GASTRoNOoM 1248 does not take into account scattering. You could try approximating 1249 the effect of scattering on the acceleration, but at this point this is 1250 not taken into account accurately. 1251 1252 @return: The wavelength and weighted kappas grid 1253 @rtype: (array,array) 1254 1255 ''' 1256 1257 wave_list,kappas = self.readKappas() 1258 if self['INCLUDE_SCAT_GAS']: 1259 #-- First the absorption coefficients of all dust species are given 1260 # Then the scattering coefficients. So iterate twice over the 1261 # dust list. 1262 wkappas = [sum([float(self['A_%s'%(species)])*float(kappas[i][j]) 1263 for i,species in enumerate(self.getDustList()*2)]) 1264 for j in xrange(len(kappas[0]))] 1265 else: 1266 #-- Only iterate once over the dust list to just take the 1267 # absorption coefficients. 1268 wkappas = [sum([float(self['A_%s'%(species)])*float(kappas[i][j]) 1269 for i,species in enumerate(self.getDustList())]) 1270 for j in xrange(len(kappas[0]))] 1271 return array(wave_list),array(wkappas)

1272 1273 1274

1275 - def getDustList(self):

1276 1277 ''' 1278 Return (and initialize) the list of nonzero abundance dust species in 1279 the model. 1280 1281 @return: The dust species 1282 @rtype: tuple(str) 1283 1284 ''' 1285 1286 if self.dust_list is None: self.readDustProperties() 1287 return self.dust_list

1288 1289 1290

1291 - def calcA_NO_NORM(self):

1292 1293 """ 1294 Set the default value of A_NO_NORM to 0. 1295 1296 """ 1297 1298 if not self.has_key('A_NO_NORM'): 1299 self['A_NO_NORM'] = 0 1300 else: 1301 pass

1302 1303 1304

1305 - def calcLS_NO_VIB(self):

1306 1307 """ 1308 Set the default value of LS_NO_VIB (remove vibrational states from the 1309 calculation and plots) to zero. 1310 1311 """ 1312 1313 if not self.has_key('LS_NO_VIB'): 1314 self['LS_NO_VIB'] = [] 1315 else: 1316 pass

1317 1318 1319

1320 - def calcN_QUAD(self):

1321 1322 """ 1323 Set the default value of N_QUAD to 100. 1324 1325 Only used when auto selecting transition based on a wavelength range. 1326 1327 """ 1328 1329 if not self.has_key('N_QUAD'): 1330 self['N_QUAD'] = 100 1331 else: 1332 pass

1333 1334 1335

1336 - def calcLS_OFFSET(self):

1337 1338 """ 1339 Set the default value of LS_OFFSET to 0.0. 1340 1341 Only used when auto selecting transitions based on a wavelength range. 1342 1343 """ 1344 1345 if not self.has_key('LS_OFFSET'): 1346 self['LS_OFFSET'] = 0.0 1347 else: 1348 pass

1349 1350 1351

1352 - def calcTLR(self):

1353 1354 """ 1355 Stefan-Boltzmann's law. 1356 1357 Star() object needs to have at least 2 out of 3 parameters (T,L,R), 1358 with in L and R in solar values and T in K. 1359 1360 The one missing parameter is calculated. 1361 1362 This method does nothing if all three are present. 1363 1364 """ 1365 1366 if not self.has_key('T_STAR'): 1367 self['T_STAR']=(float(self['L_STAR'])/float(self['R_STAR'])**2.)\ 1368 **(1/4.)*self.Tsun 1369 elif not self.has_key('L_STAR'): 1370 self['L_STAR']=(float(self['R_STAR']))**2.*\ 1371 (float(self['T_STAR'])/self.Tsun)**4. 1372 elif not self.has_key('R_STAR'): 1373 self['R_STAR']=(float(self['L_STAR'])*\ 1374 (self.Tsun/float(self['T_STAR']))**4)**(1/2.) 1375 else: 1376 pass

1377 1378 1379

1380 - def calcLL_GAS_LIST(self):

1381 1382 ''' 1383 Define Molecule() objects for the molecules requested in the LineList 1384 mode. 1385 1386 ''' 1387 1388 if not self.has_key('LL_GAS_LIST'): 1389 if type(self['LL_MOLECULES']) is types.StringType: 1390 self['LL_GAS_LIST'] = [Molecule.Molecule(linelist=1,\ 1391 molecule=self['LL_MOLECULES'])] 1392 else: 1393 self['LL_GAS_LIST'] = [Molecule.Molecule(molecule=molec, 1394 linelist=1) 1395 for molec in self['LL_MOLECULES']] 1396 keys = ['LL_CAT','LL_MAX_EXC','LL_MIN_STRENGTH'] 1397 for k in keys: 1398 if not isinstance(self[k],tuple): 1399 self[k] = [self[k]]*len(self['LL_GAS_LIST']) 1400 else: 1401 pass

1402 1403 1404

1405 - def calcUSE_MASER_IN_SPHINX(self):

1406 1407 ''' 1408 Set the default value of USE_MASER_IN_SPHINX parameter. 1409 1410 ''' 1411 1412 if not self.has_key('USE_MASER_IN_SPHINX'): 1413 self['USE_MASER_IN_SPHINX'] = 1 1414 else: 1415 pass

1416 1417 1418

1419 - def calcUSE_NO_MASER_OPTION(self):

1420 1421 ''' 1422 Set the default value of USE_NO_MASER_OPTION parameter. 1423 1424 ''' 1425 1426 if not self.has_key('USE_NO_MASER_OPTION'): 1427 self['USE_NO_MASER_OPTION'] = 0 1428 else: 1429 pass

1430 1431 1432

1433 - def calcXDEX(self):

1434 1435 ''' 1436 Set the default value of XDEX parameter. 1437 1438 ''' 1439 1440 if not self.has_key('XDEX'): 1441 self['XDEX'] = 2. 1442 else: 1443 pass

1444 1445 1446

1447 - def calcSTOCHASTIC_VEL_MODE(self):

1448 1449 ''' 1450 Set the default value of STOCHASTIC_VEL_MODE parameter. 1451 1452 ''' 1453 1454 if not self.has_key('STOCHASTIC_VEL_MODE'): 1455 self['STOCHASTIC_VEL_MODE'] = 'CONSTANT' 1456 else: 1457 pass

1458 1459 1460

1461 - def calcSTOCHASTIC_VEL(self):

1462 1463 ''' 1464 Set the default value of STOCHASTIC_VEL parameter. 1465 1466 ''' 1467 1468 if not self.has_key('STOCHASTIC_VEL'): 1469 self['STOCHASTIC_VEL'] = 0. 1470 else: 1471 pass

1472 1473 1474

1475 - def calcSTOCHASTIC_VEL_POWER(self):

1476 1477 ''' 1478 Set the default value of STOCHASTIC_VEL_POWER parameter. 1479 1480 ''' 1481 1482 if not self.has_key('STOCHASTIC_VEL_POWER'): 1483 self['STOCHASTIC_VEL_POWER'] = 0. 1484 else: 1485 pass

1486 1487 1488

1489 - def calcSTOCHASTIC_VEL_INNER(self):

1490 1491 ''' 1492 Set the default value of STOCHASTIC_VEL_INNER parameter. 1493 1494 ''' 1495 1496 if not self.has_key('STOCHASTIC_VEL_INNER'): 1497 self['STOCHASTIC_VEL_INNER'] = 1.0 1498 else: 1499 pass

1500 1501 1502

1503 - def calcSTOCHASTIC_VEL_FILENAME(self):

1504 1505 ''' 1506 Set the default value of STOCHASTIC_VEL_FILENAME parameter. 1507 1508 ''' 1509 1510 if not self.has_key('STOCHASTIC_VEL_FILENAME'): 1511 self['STOCHASTIC_VEL_FILENAME'] = '' 1512 else: 1513 pass

1514 1515 1516

1517 - def calcFEHLER(self):

1518 1519 ''' 1520 Set the default value of FEHLER parameter. 1521 1522 ''' 1523 1524 if not self.has_key('FEHLER'): 1525 self['FEHLER'] = 1e-4 1526 else: 1527 pass

1528 1529 1530

1531 - def calcN_FREQ(self):

1532 1533 ''' 1534 Set the default value of N_FREQ parameter. 1535 1536 ''' 1537 1538 if not self.has_key('N_FREQ'): 1539 self['N_FREQ'] = 30 1540 else: 1541 pass

1542 1543 1544

1545 - def calcSTART_APPROX(self):

1546 1547 ''' 1548 Set the default value of START_APPROX parameter. 1549 1550 ''' 1551 1552 if not self.has_key('START_APPROX'): 1553 self['START_APPROX'] = 0 1554 else: 1555 pass

1556 1557 1558

1559 - def calcUSE_FRACTION_LEVEL_CORR(self):

1560 1561 ''' 1562 Set the default value of USE_FRACTION_LEVEL_CORR parameter. 1563 1564 ''' 1565 1566 if not self.has_key('USE_FRACTION_LEVEL_CORR'): 1567 self['USE_FRACTION_LEVEL_CORR'] = 1 1568 else: 1569 pass

1570 1571 1572

1573 - def calcFRACTION_LEVEL_CORR(self):

1574 1575 ''' 1576 Set the default value of FRACTION_LEVEL_CORR parameter. 1577 1578 ''' 1579 1580 if not self.has_key('FRACTION_LEVEL_CORR'): 1581 self['FRACTION_LEVEL_CORR'] = 0.8 1582 else: 1583 pass

1584 1585 1586

1587 - def calcNUMBER_LEVEL_MAX_CORR(self):

1588 1589 ''' 1590 Set the default value of NUMBER_LEVEL_MAX_CORR parameter. 1591 1592 ''' 1593 1594 if not self.has_key('NUMBER_LEVEL_MAX_CORR'): 1595 self['NUMBER_LEVEL_MAX_CORR'] = 1e-12 1596 else: 1597 pass

1598 1599 1600

1601 - def calcFRACTION_TAU_STEP(self):

1602 1603 ''' 1604 Set the default value of FRACTION_TAU_STEP parameter. 1605 1606 ''' 1607 1608 if not self.has_key('FRACTION_TAU_STEP'): 1609 self['FRACTION_TAU_STEP'] = 1e-2 1610 else: 1611 pass

1612 1613 1614

1615 - def calcMIN_TAU_STEP(self):

1616 1617 ''' 1618 Set the default value of MIN_TAU_STEP parameter. 1619 1620 ''' 1621 1622 if not self.has_key('MIN_TAU_STEP'): 1623 self['MIN_TAU_STEP'] = 1e-4 1624 else: 1625 pass

1626 1627 1628

1629 - def calcWRITE_INTENSITIES(self):

1630 1631 ''' 1632 Set the default value of WRITE_INTENSITIES parameter. 1633 1634 ''' 1635 1636 if not self.has_key('WRITE_INTENSITIES'): 1637 self['WRITE_INTENSITIES'] = 0 1638 else: 1639 pass

1640 1641 1642

1643 - def calcTAU_MAX(self):

1644 1645 ''' 1646 Set the default value of TAU_MAX parameter. 1647 1648 ''' 1649 1650 if not self.has_key('TAU_MAX'): 1651 self['TAU_MAX'] = 12. 1652 else: 1653 pass

1654 1655 1656

1657 - def calcTAU_MIN(self):

1658 1659 ''' 1660 Set the default value of TAU_MIN parameter. 1661 1662 ''' 1663 1664 if not self.has_key('TAU_MIN'): 1665 self['TAU_MIN'] = -6. 1666 else: 1667 pass

1668 1669 1670

1671 - def calcCHECK_TAU_STEP(self):

1672 1673 ''' 1674 Set the default value of CHECK_TAU_STEP parameter. 1675 1676 ''' 1677 1678 if not self.has_key('CHECK_TAU_STEP'): 1679 self['CHECK_TAU_STEP'] = 1e-2 1680 else: 1681 pass

1682 1683 1684

1685 - def calcLOGG(self):

1686 1687 """ 1688 Set the default value of LOGG to 0. 1689 1690 """ 1691 1692 if not self.has_key('LOGG'): 1693 self['LOGG']=0 1694 else: 1695 pass

1696 1697

1698 - def calcF_CONT_TYPE(self):

1699 1700 """ 1701 Set the default value of F_CONT_TYPE to MCMax. This is the type of 1702 derivation of measured continuum fluxes. Can be: ISO, MSX, PHOT, MCMax 1703 1704 """ 1705 1706 if not self.has_key('F_CONT_TYPE'): 1707 self['F_CONT_TYPE'] = 'MCMax' 1708 else: 1709 pass

1710 1711

1712 - def calcAH2O_RATE(self):

1713 1714 ''' 1715 Calculate the outflow rate of H2O, by multiplying the H2O abundance 1716 with the mass-loss rate. 1717 1718 Value is set in units of Msun/yr 1719 1720 ''' 1721 1722 if not self.has_key('AH2O_RATE'): 1723 self['AH2O_RATE'] = self['F_H2O'] * self['MDOT_GAS'] 1724 else: 1725 pass

1726 1727 1728

1729 - def calcT_INNER_DUST(self):

1730 1731 """ 1732 Find the dust temperature at the inner radius in Kelvin. 1733 1734 Taken from last mcmax model, and defined by the dust species able to 1735 exist at the highest temperature; if no mcmax model is present, the 1736 temperature is taken to be zero, indicating no inner radius T is 1737 available. 1738 1739 """ 1740 1741 if not self.has_key('T_INNER_DUST'): 1742 rad = self.getDustRad(unit='rstar') 1743 temp = self.getDustTemperature() 1744 if not rad.size: 1745 self['T_INNER_DUST'] = 0 1746 return 1747 rin = self['R_INNER_DUST'] 1748 self['T_INNER_DUST'] = temp[argmin(abs(rad-rin))] 1749 else: 1750 pass

1751 1752 1753

1754 - def calcTEMPERATURE_EPSILON_GAS(self):

1755 1756 """ 1757 If not present in input, TEMPERATURE_EPSILON_GAS is equal to 0.5. 1758 1759 """ 1760 1761 if not self.has_key('TEMPERATURE_EPSILON_GAS'): 1762 self['TEMPERATURE_EPSILON_GAS'] = 0.5 1763 else: 1764 pass

1765 1766 1767

1768 - def calcTEMPERATURE_EPSILON2_GAS(self):

1769 1770 """ 1771 If not present in input, TEMPERATURE_EPSILON2_GAS is equal to 0, 1772 in which case it will be ignored when making input file. 1773 1774 """ 1775 1776 if not self.has_key('TEMPERATURE_EPSILON2_GAS'): 1777 self['TEMPERATURE_EPSILON2_GAS'] = 0 1778 else: 1779 pass

1780 1781 1782

1783 - def calcRADIUS_EPSILON2_GAS(self):

1784 1785 """ 1786 If not present in input, RADIUS_EPSILON2_GAS is equal to 0, \ 1787 in which case it will be ignored when making input file. 1788 1789 """ 1790 1791 if not self.has_key('RADIUS_EPSILON2_GAS'): 1792 self['RADIUS_EPSILON2_GAS'] = 0 1793 else: 1794 pass

1795 1796 1797

1798 - def calcTEMPERATURE_EPSILON3_GAS(self):

1799 1800 """ 1801 If not present in input, TEMPERATURE_EPSILON3_GAS is equal to 0, 1802 in which case it will be ignored when making input file. 1803 1804 """ 1805 1806 if not self.has_key('TEMPERATURE_EPSILON3_GAS'): 1807 self['TEMPERATURE_EPSILON3_GAS'] = 0 1808 else: 1809 pass

1810 1811 1812

1813 - def calcRADIUS_EPSILON3_GAS(self):

1814 1815 """ 1816 If not present in input, RADIUS_EPSILON3_GAS is equal to 0, 1817 in which case it will be ignored when making input file. 1818 1819 """ 1820 1821 if not self.has_key('RADIUS_EPSILON3_GAS'): 1822 self['RADIUS_EPSILON3_GAS'] = 0 1823 else: 1824 pass

1825 1826 1827

1828 - def calcDUST_TO_GAS_CHANGE_ML_SP(self):

1829 1830 """ 1831 Set default value of sphinx/mline specific d2g ratio to the 1832 semi-empirical d2g ratio, ie based on MDOT_DUST and MDOT_GAS. 1833 1834 In order to turn this off, set this parameter to 0 in the input file, 1835 in which case the iterated acceleration d2g ratio is used. 1836 1837 Both MDOT_GAS and MDOT_DUST have to be defined explicitly if this 1838 parameter is not. 1839 1840 This parameter has to be defined explicitly if one of MDOT_GAS and 1841 MDOT_DUST is not defined explicitly. 1842 1843 Note that the DUST_TO_GAS keyword is the internal representation of the 1844 dust_to_gas ratio and should never be explicitly defined. For all 1845 practical purposes, use DUST_TO_GAS_CHANGE_ML_SP. 1846 1847 """ 1848 1849 if not self.has_key('DUST_TO_GAS_CHANGE_ML_SP'): 1850 if not self.has_key('MDOT_DUST'): 1851 raise IOError('Both MDOT_DUST and DUST_TO_GAS_CHANGE_ML_SP '+\ 1852 'are undefined.') 1853 if not self.has_key('MDOT_GAS'): 1854 raise IOError('Both MDOT_GAS and DUST_TO_GAS_CHANGE_ML_SP '+\ 1855 'are undefined.') 1856 self['DUST_TO_GAS_CHANGE_ML_SP'] = self['DUST_TO_GAS'] 1857 else: 1858 pass

1859 1860 1861

1862 - def calcR_INNER_GAS(self):

1863 1864 """ 1865 If not present in input, R_INNER_GAS is equal to R_INNER_DUST 1866 1867 """ 1868 1869 if not self.has_key('R_INNER_GAS'): 1870 self['R_INNER_GAS'] = self['R_INNER_DUST'] 1871 else: 1872 pass

1873 1874 1875

1876 - def calcUSE_DENSITY_NON_CONSISTENT(self):

1877 1878 """ 1879 Set USE_DENSITY_NON_CONSISTENT off by default. 1880 1881 """ 1882 1883 if not self.has_key('USE_DENSITY_NON_CONSISTENT'): 1884 self['USE_DENSITY_NON_CONSISTENT'] = 0 1885 else: 1886 pass

1887 1888 1889

1890 - def calcR_OUTER_DUST(self):

1891 1892 """ 1893 If not present in input, R_OUTER_DUST is calculated from 1894 R_OUTER_DUST_AU. 1895 1896 """ 1897 1898 if not self.has_key('R_OUTER_DUST'): 1899 if self.has_key('R_OUTER_DUST_AU'): 1900 self['R_OUTER_DUST'] = self['R_OUTER_DUST_AU']*self.au\ 1901 /self['R_STAR']/self.Rsun 1902 elif self.has_key('R_OUTER_MULTIPLY'): 1903 self['R_OUTER_DUST'] = self['R_INNER_DUST']\ 1904 *self['R_OUTER_MULTIPLY'] 1905 else: 1906 pass

1907 1908 1909

1910 - def calcR_INNER_DUST(self):

1911 1912 """ 1913 Calculate the inner radius from MCMax output in stellar radii. 1914 1915 If no MCMax model is calculated yet, R_{i,d} is the stellar radius. 1916 1917 Else, the inner dust radius is taken where the density reaches a 1918 threshold, defined by R_INNER_DUST_MODE: 1919 1920 - MAX: Density reaches a maximum value, depends on the different 1921 condensation temperatures of the dust species taken into 1922 account 1923 - ABSOLUTE: Density becomes larger than 10**(-30) 1924 - RELATIVE: Density becomes larger than 1% of maximum density 1925 1926 Unless defined in the CC input file, the dust radius is updated every 1927 time a new iteration starts. 1928 1929 If no MCMax model is known, and destruction temperature iteration is 1930 off, the inner radius is 2 stellar radii for calculation time reasons. 1931 1932 """ 1933 1934 if not self.has_key('R_INNER_DUST'): 1935 if self.has_key('R_INNER_DUST_AU'): 1936 self['R_INNER_DUST'] = self['R_INNER_DUST_AU']*self.au\ 1937 /self['R_STAR']/self.Rsun 1938 else: 1939 rad = self.getDustRad() 1940 if rad.size == 0: 1941 self['R_INNER_DUST'] = 1.0 1942 return 1943 dens = self.getDustDensity() 1944 if self['R_INNER_DUST_MODE'] == 'MAX': 1945 ri_cm = rad[argmax(dens)] 1946 elif self['R_INNER_DUST_MODE'] == 'ABSOLUTE': 1947 ri_cm = rad[dens>10**(-30)][0] 1948 else: 1949 ri_cm = rad[dens>0.01*max(dens)][0] 1950 self['R_INNER_DUST'] = ri_cm/self.Rsun/self['R_STAR'] 1951 else: 1952 pass

1953 1954 1955

1956 - def calcR_INNER_DUST_MODE(self):

1957 1958 """ 1959 The mode of calculating the inner radius from MCMax output, if present. 1960 1961 Can be ABSOLUTE (dens>10**-50) or RELATIVE (dens[i]>dens[i+1]*0.01). 1962 1963 CLASSIC reproduces the old method. 1964 1965 Set here to the default value of ABSOLUTE. 1966 1967 """ 1968 1969 if not self.has_key('R_INNER_DUST_MODE'): 1970 self['R_INNER_DUST_MODE'] = 'ABSOLUTE' 1971 else: 1972 pass

1973 1974 1975

1976 - def calcRID_TEST(self):

1977 1978 """ 1979 The mode of determining the dust temp profile. 1980 1981 Only for testing purposes. 1982 1983 - Default is 'R_STAR', ie the temperature is taken from the stellar 1984 radius onward, regardless of what the inner radius is. 1985 1986 - 'R_INNER_GAS' is used for taking the dust temperature from the 1987 inner gas radius onward. 1988 1989 - 'BUGGED_CASE' is the old version where r [R*] > R_STAR [Rsun]. 1990 1991 """ 1992 1993 if not self.has_key('RID_TEST'): 1994 self['RID_TEST'] = 'R_STAR' 1995 else: 1996 pass

1997 1998 1999

2000 - def calcSPEC_DENS_DUST(self):

2001 2002 """ 2003 Calculating average specific density of dust shell. 2004 2005 This is based on the input dust species abundances and their specific 2006 densities. 2007 2008 """ 2009 2010 if not self.has_key('SPEC_DENS_DUST'): 2011 if int(self['MRN_DUST']): 2012 these_sd = [self.dust[sp]['sd'] 2013 for sp in self.getDustList() 2014 if self.has_key('RGRAIN_%s'%sp)] 2015 sd_mrn = sum(these_sd)/len(these_sd) 2016 a_sd = sum([self['A_' + sp]*self.dust[sp]['sd'] 2017 for sp in self.getDustList() 2018 if self['A_%s'%sp] != 2.]) 2019 self['SPEC_DENS_DUST'] = (sd_mrn + a_sd)/2. 2020 else: 2021 these_sd = [self['A_' + sp]*self.dust[sp]['sd'] 2022 for sp in self.getDustList()] 2023 self['SPEC_DENS_DUST'] = sum(these_sd) 2024 else: 2025 pass

2026 2027 2028

2029 - def calcLAST_MCMAX_MODEL(self):

2030 2031 """ 2032 Creates empty string if not present yet. 2033 2034 """ 2035 2036 if not self.has_key('LAST_MCMAX_MODEL'): 2037 self['LAST_MCMAX_MODEL'] = '' 2038 else: 2039 pass

2040 2041 2042

2043 - def calcLAST_PACS_MODEL(self):

2044 2045 """ 2046 Sets to None if not present yet. 2047 2048 """ 2049 2050 if not self.has_key('LAST_PACS_MODEL'): 2051 self['LAST_PACS_MODEL'] = None 2052 else: 2053 pass

2054 2055 2056

2057 - def calcLAST_SPIRE_MODEL(self):

2058 2059 """ 2060 Sets to None if not present yet. 2061 2062 Note that this is an index if it IS present, and can be zero. Always 2063 check with None instead of boolean. 2064 2065 """ 2066 2067 if not self.has_key('LAST_SPIRE_MODEL'): 2068 self['LAST_SPIRE_MODEL'] = None 2069 else: 2070 pass

2071 2072 2073

2074 - def calcLAST_GASTRONOOM_MODEL(self):

2075 2076 """ 2077 Creates empty string if not present yet. 2078 2079 """ 2080 2081 if not self.has_key('LAST_GASTRONOOM_MODEL'): 2082 self['LAST_GASTRONOOM_MODEL'] = '' 2083 else: 2084 pass

2085 2086

2087 - def calcDRIFT_TYPE(self):

2088 2089 """ 2090 Set the type of drift between dust and gas taken into account. 2091 2092 Is either consistent (from momentum transfer calculation or zero). 2093 2094 """ 2095 2096 if not self.has_key('DRIFT_TYPE'): 2097 if self['V_EXP_DUST'] == self['VEL_INFINITY_GAS']: 2098 self['DRIFT_TYPE'] = 'ZERO' 2099 else: 2100 self['DRIFT_TYPE'] = 'CONSISTENT' 2101 else: 2102 pass

2103 2104 2105

2106 - def calcDRIFT(self):

2107 2108 """ 2109 Find terminal drift velocity from last calculated GASTRoNOoM model. 2110 2111 Units are km/s and is given for grain size 0.25 micron. 2112 2113 If no GASTRoNOoM model exists, the drift is taken to be 0. 2114 2115 """ 2116 2117 if not self.has_key('DRIFT'): 2118 try: 2119 #- The last 10 entries should be roughly constant anyway, 2120 #- ~terminal values 2121 self['DRIFT'] = self.getAverageDrift()[-5]/10.**5 2122 except IOError: 2123 self['DRIFT'] = 0 2124 print 'No GASTRoNOoM model has been calculated yet, drift ' + \ 2125 'is unknown and set to the default of 0.' 2126 else: 2127 pass

2128 2129 2130

2131 - def calcM_DUST(self):

2132 2133 """ 2134 Calculate the total dust mass based on the requested denstype 2135 prescription. 2136 2137 The dust mass is set in solar masses. 2138 2139 If Denstype == 'POW': 2140 Find total dust mass, based on sigma_0 in the case of a power law. 2141 2142 In all other cases, the dust mass is calculated from the density profile 2143 of the GASTRoNOoM model after correcting for d2g ratio, for now. 2144 2145 Note that in the latter case the dust mass will be an upper limit 2146 when compared to the actual dust mass determined by MCMax, if tdesiter 2147 is on: Some of the dust in the density profile will be destroyed based 2148 on the condensation temperature. 2149 2150 Only in the case of denstypes POW and MEIXNER this dust mass is actually 2151 used as input for MCMax. This never takes into account the MCMax output. 2152 To retrieve the dust mass calculated with MCMax, check the MCMax log 2153 file. 2154 2155 Currently GASTRoNOoM is never calculated first, so giving this as input 2156 does not yet work for any denstype other than POW. The value returned is 2157 correct, however. When implemented, 2 iterations will be required. 2158 2159 """ 2160 2161 if not self.has_key('M_DUST'): 2162 #-- M_DUST is an input parameter calculated from an equation 2163 if self['DENSTYPE'] == 'POW': 2164 if self['DENSPOW'] == 2: 2165 self['M_DUST'] \ 2166 = 2*pi*self['DENSSIGMA_0']\ 2167 *(self['R_INNER_DUST']*self['R_STAR']*self.Rsun)**2\ 2168 *log(self['R_OUTER_DUST']/float(self['R_INNER_DUST']))\ 2169 /self.Msun 2170 else: 2171 self['M_DUST'] \ 2172 = 2*pi*self['DENSSIGMA_0']\ 2173 *(self['R_INNER_DUST']*self['R_STAR']*self.Rsun)**2\ 2174 /(2.-self['DENSPOW'])/self.Msun\ 2175 *((self['R_OUTER_DUST']/float(self['R_INNER_DUST']))\ 2176 **(2.-self['DENSPOW'])-1.) 2177 #-- M_DUST is not an input parameter and can be calculated from the 2178 # MCMax density profile. Note that for DENSTYPE==MEIXNER M_DUST 2179 # has to be given in the CC inputfile, so this is never calculated 2180 #-- If MCMax was not yet ran, calculate the dust mass from the gas 2181 # density profile times the dust-to-gas ratio. 2182 else: 2183 if self['LAST_MCMAX_MODEL']: 2184 rad = self.getDustRad(unit='cm') 2185 dens = self.getDustDensity() 2186 m_dust = integrate.trapz(x=rad,y=dens*pi*4*rad**2) 2187 self['M_DUST'] = m_dust/self.Msun 2188 elif self['LAST_GASTRONOOM_MODEL']: 2189 rad = self.getGasRad(unit='rstar') 2190 dens = self.getWindDensity(denstype='dust') 2191 selection = (rad>self["R_INNER_DUST"])\ 2192 *(rad<self['R_OUTER_DUST']) 2193 dens = dens[selection] 2194 rad = rad[selection]*self['R_STAR']*self.Rsun 2195 m_dust = integrate.trapz(x=rad,y=dens*pi*4*rad**2) 2196 self['M_DUST'] = m_dust/self.Msun 2197 else: 2198 pass 2199 else: 2200 pass

2201 2202 2203

2204 - def calcMDOT_DUST(self):

2205 2206 ''' 2207 Calculate the value of MDOT_DUST from the DUST_TO_GAS_RATIO_ML_SP. 2208 2209 Requires MDOT_GAS and VEL_INFINITY_GAS to be defined. 2210 2211 This parameter is recalculated after every iteration and updates 2212 V_EXP_DUST in the equation. 2213 2214 MDOT_DUST can be given explicitly in the inputfile in which case it 2215 remains unchanged. 2216 2217 MDOT_DUST is used to calculate the real DUST_TO_GAS ratio parameter. So 2218 through explicit definition of 2 parameters out of MDOT_GAS, MDOT_DUST 2219 and DUST_TO_GAS_CHANGE_ML_SP you can control what the internal 2220 dust-to-gas ratio should be. 2221 2222 If DUST_TO_GAS_CHANGE_ML_SP is not given, MDOT_DUST and MDOT_GAS have 2223 to be defined explicitly. 2224 2225 ''' 2226 2227 if not self.has_key('MDOT_DUST'): 2228 if not self.has_key('DUST_TO_GAS_CHANGE_ML_SP'): 2229 raise IOError('Both MDOT_DUST and DUST_TO_GAS_CHANGE_ML_SP '+\ 2230 'are undefined.') 2231 if not self.has_key('MDOT_GAS'): 2232 raise IOError('Both MDOT_DUST and MDOT_GAS are undefined.') 2233 self['MDOT_DUST'] = float(self['DUST_TO_GAS_CHANGE_ML_SP'])\ 2234 /float(self['VEL_INFINITY_GAS'])\ 2235 *float(self['V_EXP_DUST'])\ 2236 *float(self['MDOT_GAS']) 2237 else: 2238 pass 2239 2240

2241 - def calcMDOT_GAS(self):

2242 2243 ''' 2244 Calculate the value of MDOT_GAS from the DUST_TO_GAS_RATIO_ML_SP. 2245 2246 Requires MDOT_DUST and VEL_INFINITY_GAS to be defined. 2247 2248 This parameter is recalculated after every iteration and updates 2249 V_EXP_DUST in the equation. 2250 2251 MDOT_GAS can be given explicitly in the inputfile in which case it 2252 remains unchanged. 2253 2254 MDOT_GAS is used to calculate the real DUST_TO_GAS ratio parameter. So 2255 through explicit definition of 2 parameters out of MDOT_GAS, MDOT_DUST 2256 and DUST_TO_GAS_CHANGE_ML_SP you can control what the internal 2257 dust-to-gas ratio should be. 2258 2259 If DUST_TO_GAS_CHANGE_ML_SP is not given, MDOT_GAS has to be defined 2260 explicitly. 2261 2262 ''' 2263 2264 if not self.has_key('MDOT_GAS'): 2265 if not self.has_key('DUST_TO_GAS_CHANGE_ML_SP'): 2266 raise IOError('Both MDOT_GAS and DUST_TO_GAS_CHANGE_ML_SP '+\ 2267 'are undefined.') 2268 if not self.has_key('MDOT_DUST'): 2269 raise IOError('Both MDOT_DUST and MDOT_GAS are undefined.') 2270 self['MDOT_GAS'] = float(self['VEL_INFINITY_GAS'])\ 2271 /float(self['V_EXP_DUST'])\ 2272 *float(self['MDOT_DUST'])\ 2273 /float(self['DUST_TO_GAS_CHANGE_ML_SP']) 2274 else: 2275 pass 2276 2277

2278 - def calcMDOT_MODE(self):

2279 2280 ''' 2281 Set the default value of MDOT_MODE to constant. 2282 2283 ''' 2284 2285 if not self.has_key('MDOT_MODE'): 2286 self['MDOT_MODE'] = 'CONSTANT' 2287 else: 2288 pass

2289 2290 2291

2292 - def calcMDOT_GAS_START(self):

2293 2294 ''' 2295 Set the default value of MDOT_GAS_START equal to MDOT_GAS. 2296 2297 ''' 2298 2299 if not self.has_key('MDOT_GAS_START'): 2300 self['MDOT_GAS_START'] = self['MDOT_GAS'] 2301 else: 2302 pass

2303 2304

2305 - def calcSHELLMASS_CLASS(self):

2306 2307 ''' 2308 Set the order of magnitude of SHELLMASS = Mdot/v_inf. 2309 0: Mdot/v_inf < 5e-8 2310 1: 5e-8 <= Mdot/v_inf < 2e-7 2311 2: 2e-7 <= Mdot/v_inf < 5e-7 2312 3: 5e-7 <= Mdot/v_inf 2313 2314 ''' 2315 2316 if not self.has_key('SHELLMASS_CLASS'): 2317 if self['SHELLMASS'] < 5e-8: 2318 #self['SHELLMASS_CLASS'] = (0,r'$\dot{M}_\mathrm{g}/v_{\infty\mathrm{,g}} < 5 \times 10^{-8}\ \mathrm{M}_\odot\ \mathrm{yr}^{-1}\ \mathrm{km}^{-1}\ \mathrm{s}$') 2319 self['SHELLMASS_CLASS'] = (0,r'$\dot{M}_\mathrm{g}/v_{\infty\mathrm{,g}} < 5 \times 10^{-8}$') 2320 elif self['SHELLMASS'] >= 3e-7: 2321 #self['SHELLMASS_CLASS'] = (3,r'$\dot{M}_\mathrm{g}/v_{\infty\mathrm{,g}} \geq 3 \times 10^{-7}\ \mathrm{M}_\odot\ \mathrm{yr}^{-1}\ \mathrm{km}^{-1}\ \mathrm{s}$') 2322 self['SHELLMASS_CLASS'] = (3,r'$\dot{M}_\mathrm{g}/v_{\infty\mathrm{,g}} \geq 3 \times 10^{-7}$') 2323 elif self['SHELLMASS'] >= 5e-8 and self['SHELLMASS'] < 1.0e-7: 2324 self['SHELLMASS_CLASS'] = (1,r'$5 \times 10^{-8}$ $\leq \dot{M}_\mathrm{g}/v_{\infty\mathrm{,g}} < 1 \times 10^{-7}$') 2325 #self['SHELLMASS_CLASS'] = (1,r'$5 \times 10^{-8}\ \mathrm{M}_\odot\ \mathrm{yr}^{-1}\ \mathrm{km}^{-1}\ \mathrm{s}$ $\leq \dot{M}_\mathrm{g}/v_{\infty\mathrm{,g}} < 1 \times 10^{-7}\ \mathrm{M}_\odot\ \mathrm{yr}^{-1}\ \mathrm{km}^{-1}\ \mathrm{s}$') 2326 else: 2327 #self['SHELLMASS_CLASS'] = (2,r'$1 \times 10^{-7}\ \mathrm{M}_\odot\ \mathrm{yr}^{-1}\ \mathrm{km}^{-1}\ \mathrm{s}$ $\leq \dot{M}_\mathrm{g}/v_{\infty\mathrm{,g}} < 3 \times 10^{-7}\ \mathrm{M}_\odot\ \mathrm{yr}^{-1}\ \mathrm{km}^{-1}\ \mathrm{s}$') 2328 self['SHELLMASS_CLASS'] = (2,r'$1 \times 10^{-7}$ $\leq \dot{M}_\mathrm{g}/v_{\infty\mathrm{,g}} < 3 \times 10^{-7}$') 2329 else: 2330 pass

2331 2332 2333

2334 - def calcMDOT_CLASS(self):

2335 2336 ''' 2337 Set the order of magnitude of MDOT. 2338 0: Mdot < 3e-7 2339 1: 3e-7 <= Mdot < 3e-6 2340 2: 3e-6 <= Mdot < 1e-5 2341 3: 1e-5 <= Mdot 2342 2343 ''' 2344 2345 if not self.has_key('MDOT_CLASS'): 2346 if self['MDOT_GAS'] < 3e-7: 2347 self['MDOT_CLASS'] = (0,r'$\dot{M}_\mathrm{g} < 3 \times 10^{-7}\ \mathrm{M}_\odot\ \mathrm{yr}^{-1}$') 2348 elif self['MDOT_GAS'] >= 1e-5: 2349 self['MDOT_CLASS'] = (3,r'$\dot{M}_\mathrm{g} \geq 1 \times 10^{-5}\ \mathrm{M}_\odot\ \mathrm{yr}^{-1}$') 2350 elif self['MDOT_GAS'] >= 3e-7 and self['MDOT_GAS'] < 3e-6: 2351 self['MDOT_CLASS'] = (1,r'$3 \times 10^{-7}\ \mathrm{M}_\odot\ \mathrm{yr}^{-1}$ $\leq \dot{M}_\mathrm{g} < 3 \times 10^{-6}\ \mathrm{M}_\odot\ \mathrm{yr}^{-1}$') 2352 else: 2353 self['MDOT_CLASS'] = (2,r'$3 \times 10^{-6}\ \mathrm{M}_\odot\ \mathrm{yr}^{-1}$ $\leq \dot{M}_\mathrm{g} < 1 \times 10^{-5}\ \mathrm{M}_\odot\ \mathrm{yr}^{-1}$') 2354 else: 2355 pass

2356 2357

2358 - def calcQ_STAR(self):

2359 2360 ''' 2361 Set the stellar pulsation constant (che1992). 2362 2363 ''' 2364 2365 if not self.has_key('Q_STAR'): 2366 self['Q_STAR'] = self['P_STAR']*self['M_STAR']**0.5\ 2367 *self['R_STAR']**(-3/2.) 2368 else: 2369 pass

2370 2371

2372 - def calcSCD_CLASS(self):

2373 2374 ''' 2375 Set the order of magnitude of SHELLCOLDENS. 2376 0: scd < 0.06 2377 1: 0.06 <= scd < 0.15 2378 2: 0.15 <= scd < 0.4 2379 3: 0.4 <= scd 2380 2381 ''' 2382 2383 if not self.has_key('SCD_CLASS'): 2384 if self['SHELLCOLDENS'] < 0.07: 2385 self['SCD_CLASS'] = (0,r'$\bar{m} < 0.07\ \mathrm{g\;cm}^{-2}$') 2386 elif self['SHELLCOLDENS'] >= 0.07 and self['SHELLCOLDENS'] < 0.15: 2387 self['SCD_CLASS'] = (1,r'$0.07\ \mathrm{g\;cm}^{-2}$ $\leq \bar{m} < 0.15\ \mathrm{g\;cm}^{-2}$') 2388 elif self['SHELLCOLDENS'] >=0.4: 2389 self['SCD_CLASS'] = (3,r'$\bar{m} \geq 0.4\ \mathrm{g\;cm}^{-2}$') 2390 else: 2391 self['SCD_CLASS'] = (2,r'$0.15\ \mathrm{g\;cm}^{-2}$ $\leq \bar{m} < 0.4\ \mathrm{g\;cm}^{-2}$') 2392 else: 2393 pass

2394 2395

2396 - def calcL_CLASS(self):

2397 2398 ''' 2399 Set the order of magnitude of L_STAR. 2400 2401 0: lstar < 6000 2402 1: 6000 <= lstar < 8000 2403 2: 8000 <= lstar < 10000 2404 3: 10000 <= lstar 2405 2406 ''' 2407 2408 if not self.has_key('L_CLASS'): 2409 if self['L_STAR'] < 6000: 2410 self['L_CLASS'] = (0,r'$L_\star < 6000$ $\mathrm{L}_\odot$') 2411 elif self['L_STAR'] >= 10000: 2412 self['L_CLASS'] = (3,r'$L_\star \geq 10000$ $\mathrm{L}_\odot$') 2413 elif self['L_STAR'] >= 8000 and self['L_STAR'] < 10000: 2414 self['L_CLASS'] = (2,r'$8000$ $\mathrm{L}_\odot$ $\leq L_\star < 10000$ $\mathrm{L}_\odot$') 2415 else: 2416 self['L_CLASS'] = (1,r'$6000$ $\mathrm{L}_\odot$ $\leq L_\star < 8000$ $\mathrm{L}_\odot$') 2417 else: 2418 pass

2419 2420 2421

2422 - def calcSCATTYPE(self):

2423 2424 ''' 2425 Set the default scattering type for MCMax to 'ISOTROPIC'. 2426 2427 Can also be 'NONE', or 'FULL'. 2428 2429 'FULL' requires dust opacity .particle files with full scattering 2430 matrices. 2431 2432 ''' 2433 2434 if not self.has_key('SCATTYPE'): 2435 self['SCATTYPE'] = 'ISOTROPIC' 2436 else: 2437 pass

2438 2439 2440

2441 - def calcT_CONTACT(self):

2442 2443 ''' 2444 Set thermal contact off. 2445 2446 ''' 2447 2448 if not self.has_key('T_CONTACT'): 2449 self['T_CONTACT'] = 0 2450 else: 2451 pass

2452 2453 2454

2455 - def calcT_CLASS(self):

2456 2457 ''' 2458 Set the order of magnitude of T_STAR. 2459 2460 0: tstar < 2000 2461 1: 2000 <= tstar < 2200 2462 2: 2200 <= tstar < 2500 2463 3: 2500 <= tstar 2464 2465 ''' 2466 2467 if not self.has_key('T_CLASS'): 2468 if self['T_STAR'] < 2000: 2469 self['T_CLASS'] = (0,r'$T_\star < 2000\ \mathrm{K}$') 2470 elif self['T_STAR'] >= 2500: 2471 self['T_CLASS'] = (3,r'$T_\star \geq 2500\ \mathrm{K}$') 2472 elif self['T_STAR'] >= 2250 and self['T_STAR'] < 2500: 2473 self['T_CLASS'] = (2,r'$2250\ \mathrm{K}$ $\leq T_\star < 2500\ \mathrm{K}$') 2474 else: 2475 self['T_CLASS'] = (1,r'$2000\ \mathrm{K}$ $\leq T_\star < 2250\ \mathrm{K}$') 2476 else: 2477 pass

2478 2479

2480 - def calcVG_CLASS(self):

2481 2482 ''' 2483 Set the order of magnitude of VEL_INFINITY_GAS 2484 2485 0: vg < 10 2486 1: 10 <= vg < 15 2487 2: 15 <= vg < 20 2488 3: 20 <= vg 2489 2490 ''' 2491 2492 if not self.has_key('VG_CLASS'): 2493 if self['VEL_INFINITY_GAS'] < 10.: 2494 self['VG_CLASS'] = (0,r'$v_{\infty\mathrm{,g}} < 10\ \mathrm{km\;s}^{-1}$') 2495 elif self['VEL_INFINITY_GAS'] >= 20.: 2496 self['VG_CLASS'] = (3,r'$v_{\infty\mathrm{,g}} \geq 20\ \mathrm{km\;s}^{-1}$') 2497 elif self['VEL_INFINITY_GAS'] >= 15. and self['VEL_INFINITY_GAS'] < 20.: 2498 self['VG_CLASS'] = (2,r'$15\ \mathrm{km\;s}^{-1}$ $\leq v_{\infty\mathrm{,g}} < 20\ \mathrm{km\;s}^{-1}$') 2499 else: 2500 self['VG_CLASS'] = (1,r'$10\ \mathrm{km\;s}^{-1}$ $\leq v_{\infty\mathrm{,g}} < 15\ \mathrm{km\;s}^{-1}$') 2501 else: 2502 pass

2503 2504 2505

2506 - def calcV_EXP_DUST(self):

2507 2508 """ 2509 Calculate dust terminal velocity from gas terminal velocity and drift. 2510 2511 Given in km/s. 2512 2513 """ 2514 2515 if not self.has_key('V_EXP_DUST'): 2516 self['V_EXP_DUST']= float(self['VEL_INFINITY_GAS']) \ 2517 + float(self['DRIFT']) 2518 else: 2519 pass

2520 2521 2522

2523 - def calcREDDENING(self):

2524 2525 ''' 2526 A boolean flag for applying interstellar reddening or not. This is 2527 model (read: distance) dependent, hence belongs in Star() objects. 2528 2529 Having this available here makes it possible to compare using reddening 2530 or not. 2531 2532 Default value is set to 0. 2533 2534 ''' 2535 2536 if not self.has_key('REDDENING'): 2537 self['REDDENING']= 0 2538 else: 2539 pass

2540 2541 2542

2543 - def calcREDDENING_MAP(self):

2544 2545 ''' 2546 The interstellar extinction map used for determining the interstellar 2547 extinction in K-band at a given distance, in the direction of given 2548 longitude and latitude (set in Star.dat). 2549 2550 Default is Marshall et al. 2006 (marshall), but is replaced by Drimmel 2551 et al. 2003 (drimmel) in case ll and bb are outside the range of 2552 availability in Marshall. 2553 2554 Alternatives are Arenou et al. 1992 (arenou) and Schlegel et al. 1998 2555 (schlegel). 2556 2557 ''' 2558 2559 if not self.has_key('REDDENING_MAP'): 2560 self['REDDENING_MAP']= 'marshall' 2561 else: 2562 pass

2563 2564 2565

2566 - def calcREDDENING_LAW(self):

2567 2568 ''' 2569 The extinction law used to redden model spectra. 2570 2571 Default is the combination of the laws by Fitzpatrick et al. 2004 2572 (Optical) and Chiar & Tielens 2006 (IR), see IvS repo for more details 2573 at cc.ivs.sed.reddening. 2574 2575 Alternatives include cardelli1989, donnell1994, fitzpatrick1999, 2576 fitzpatrick2004, chiar2006. 2577 2578 ''' 2579 2580 if not self.has_key('REDDENING_LAW'): 2581 self['REDDENING_LAW']= 'fitz2004chiar2006' 2582 else: 2583 pass

2584 2585

2586 - def calcRT_INCLINATION(self):

2587 2588 ''' 2589 Set the default value of MCMax ray-tracing inclination to 45 degrees. 2590 2591 ''' 2592 2593 if not self.has_key('RT_INCLINATION'): 2594 self['RT_INCLINATION']= 45.0 2595 else: 2596 pass

2597 2598

2599 - def calcRT_NOSOURCE(self):

2600 2601 ''' 2602 Set the default value of excluding the central source from the ray 2603 tracing of the MCMax model to False. 2604 2605 ''' 2606 2607 if not self.has_key('RT_NOSOURCE'): 2608 self['RT_NOSOURCE']= 0 2609 else: 2610 pass

2611 2612

2613 - def calcRT_SPEC(self):

2614 2615 ''' 2616 Set the default value of MCMax ray-tracing of the spectrum to False. 2617 2618 ''' 2619 2620 if not self.has_key('RT_SPEC'): 2621 self['RT_SPEC']= 0 2622 else: 2623 pass

2624 2625 2626

2627 - def calcRT_IMAGE(self):

2628 2629 ''' 2630 Set the default value of MCMax image to False. 2631 2632 ''' 2633 2634 if not self.has_key('RT_IMAGE'): 2635 self['RT_IMAGE']= 0 2636 else: 2637 pass

2638 2639

2640 - def calcRT_VIS(self):

2641 2642 ''' 2643 Set the default value of MCMax visibilities to False. 2644 2645 ''' 2646 2647 if not self.has_key('RT_VIS'): 2648 self['RT_VIS']= 0 2649 else: 2650 pass

2651 2652

2653 - def calcRT_BASEVIS(self):

2654 2655 ''' 2656 Set the default value of MCMax visibilities versus baseline to False. 2657 2658 ''' 2659 2660 if not self.has_key('RT_BASEVIS'): 2661 self['RT_BASEVIS']= 0 2662 else: 2663 pass

2664 2665

2666 - def calcRT_REDO(self):

2667 2668 ''' 2669 Set the default value of redoing the MCMax model observation to False. 2670 2671 ''' 2672 2673 if not self.has_key('RT_REDO'): 2674 self['RT_REDO']= 0 2675 else: 2676 pass

2677 2678

2679 - def calcRT_OUTPUTFOLDER(self):

2680 2681 ''' 2682 Set the default value of MCMax ray-tracing outputfolder to empty. The 2683 output model observations are then saved in the model output folder. 2684 2685 ''' 2686 2687 if not self.has_key('RT_OUTPUTFOLDER'): 2688 self['RT_OUTPUTFOLDER']= '' 2689 else: 2690 pass

2691 2692

2693 - def calcR_MAX(self,missing_key):

2694 2695 """ 2696 Calculate the maximum existence radii for dust species. 2697 2698 Based on T_MIN_SPECIES for the species, and derived from mcmax output. 2699 2700 If not MCMax model available, a power law is assumed. If T_MIN is not 2701 given, no boundaries are assumed. 2702 2703 Is given in solar radii. 2704 2705 @param missing_key: the missing max radius for a species that is needed 2706 Of the format R_MAX_SPECIES. 2707 @type missing_key: string 2708 2709 """ 2710 2711 if not self.has_key(missing_key): 2712 #- R_MAX is for T_MIN 2713 try: 2714 tmin = float(self[missing_key.replace('R_MAX','T_MIN',1)]) 2715 if self['LAST_MCMAX_MODEL']: 2716 species = missing_key[6:] 2717 rad = self.getDustRad(species=species,unit='rstar') 2718 temp = self.getDustTemperature(species=species) 2719 temp_interp = interp1d(rad,temp) 2720 try: 2721 #-- interp1d raises a ValueError when outside of bounds 2722 rmax = temp_interp(tmin) 2723 except ValueError: 2724 #- if T_MIN (for R_MAX) > temp[0] then no dust can 2725 #- be present of the species 2726 #- RMAX should be made very small, ie R* 2727 if tmin > temp[0]: 2728 rmax = self['R_STAR'] 2729 2730 #- on the other hand, if TMIN < temp[-1] 2731 #- all dust is allowed, so raise KeyError, such that 2732 #- R_MAX is set to None 2733 elif tmin < temp[-1]: 2734 raise KeyError 2735 2736 #- In the alternative cases, something is seriously 2737 # wrong! 2738 else: 2739 raise IOError('Something went wrong when '+\ 2740 'searching for R_MAX corresponding'+\ 2741 ' to T_MIN... Debug!') 2742 else: 2743 #-- Grab the power typical for optically thin temperature 2744 # profiles in the Rayleigh limit, i.e. s=1 2745 power = -2./(4+1) 2746 #-- Take the reciprocal relation of the dust temperature 2747 # profile (as in profilers.Temperature.Tdust()) 2748 # rmax is in Rstar! 2749 rmax = (tmin/self['T_STAR'])**(1/power)/2. 2750 self[missing_key] = rmax 2751 except KeyError: 2752 self[missing_key] = None 2753 else: 2754 pass

2755 2756 2757

2758 - def calcT_DES(self,sp):

2759 2760 """ 2761 Find the max temperature at which a dust species can exist. 2762 2763 First, the CC inputfile is searched for T_MAX_SPECIES, in which case 2764 the sublimation temperature is constant. T_MAX is never made by Star()! 2765 2766 If not present, Dust.dat info is taken, being either a sublimation 2767 temperature, or the coefficients to calculate a pressure dependent 2768 sublimation temperature. These are set using T_DESA_ and T_DESB_SPECIES 2769 2770 Note that tdesa and tdesb from Dust.dat are the coefficients given in 2771 Kama et al 2009. MCMax uses a slightly different definition, and the 2772 notation has crossed that of the paper. In any case, MCMax's definition 2773 of tdesa and tdesb is defined as shown here. 2774 2775 This assumes TDESITER to be on. 2776 2777 @param sp: The dust species 2778 @type sp: string 2779 2780 """ 2781 2782 if not self.has_key('T_DESA_' + sp) \ 2783 or not self.has_key('T_DESB_' + sp): 2784 try: 2785 #-- Check if a max temperature was given for the dust species 2786 if not self['T_MAX_' + sp]: 2787 del self['T_MAX_' + sp] 2788 #-- If yes, then set the coefficients so that a constant 2789 # condensation temperature is used 2790 self['T_DESA_'+sp] = 10.0**(4)/self['T_MAX_' + sp] 2791 self['T_DESB_'+sp] = 10.0**(-4) 2792 except KeyError: 2793 if self.dust[sp]['tdesa']: 2794 #-- Set the MCMax coefficients based on the Kama 2009 vals 2795 self['T_DESA_'+sp] = 1e4*self.dust[sp]['tdesb']\ 2796 /self.dust[sp]['tdesa'] 2797 self['T_DESB_'+sp] = 1e4/self.dust[sp]['tdesa'] 2798 else: 2799 #-- If no coefficients are available, use a more or less 2800 # constant condensation temperature based on TDES in 2801 # Dust.dat. 2802 self['T_DESA_'+sp] = 1e4/self.dust[sp]['tdes'] 2803 self['T_DESB_'+sp] = 1e-4 2804 else: 2805 pass

2806 2807 2808 2809

2810 - def calcR_SHELL_UNIT(self):

2811 2812 ''' 2813 Set default value of R_SHELL_UNIT to R_STAR. 2814 2815 ''' 2816 2817 if not self.has_key('R_SHELL_UNIT'): 2818 self['R_SHELL_UNIT'] = 'R_STAR' 2819 else: 2820 pass

2821 2822 2823

2824 - def calcDENSTYPE(self):

2825 2826 """ 2827 Define the type of density distribution. 2828 2829 Default is 'MASSLOSS' for first iteration, otherwise SHELLFILE. 2830 2831 If second iteration, a DENSFILE is created taking into account the 2832 acceleration zone. This file is only created if not already present. 2833 2834 The dust density profile is calculated from the h2 number density, 2835 after scaling to the dust mass-loss rate and correcting for the dust 2836 velocity profile. 2837 2838 """ 2839 2840 if not self.has_key('DENSTYPE') or not self.has_key('DENSFILE'): 2841 if self['MDOT_MODE'] != 'CONSTANT': 2842 exstr = '_var' 2843 else: 2844 exstr = '' 2845 filename = os.path.join(cc.path.gout,'data_for_mcmax',\ 2846 '_'.join(['dens',\ 2847 self['LAST_GASTRONOOM_MODEL'],\ 2848 'mdotd%s%.2e.dat'%(exstr,\ 2849 self['MDOT_DUST'])])) 2850 if os.path.isfile(filename): 2851 self['DENSFILE'] = filename 2852 self['DENSTYPE'] = "SHELLFILE" 2853 else: 2854 if self['LAST_GASTRONOOM_MODEL']: 2855 if self.has_key('DENSTYPE'): 2856 if self['DENSTYPE'] == "MASSLOSS": 2857 raise IOError 2858 self['DENSTYPE'] = "SHELLFILE" 2859 self['DENSFILE'] = filename 2860 rad = self.getGasRad() 2861 dens = self.getWindDensity(denstype='dust') 2862 DataIO.writeCols(filename,[rad/self.au,dens]) 2863 print '** Made MCMax density input file at %s.'%filename 2864 else: 2865 print '** Writing and/or reading DENSFILE output and/or '+\ 2866 'input failed. Assuming standard mass-loss density'+\ 2867 ' distribution.' 2868 self['DENSTYPE'] = "MASSLOSS" 2869 self['DENSFILE'] = '' 2870 else: 2871 pass

2872 2873 2874

2875 - def calcSHELLMASS(self):

2876 2877 """ 2878 Calculate the average mass per shell of the circumstellar envelope. 2879 2880 Calculated by Mdot_gas/vexp. 2881 2882 """ 2883 2884 if not self.has_key('SHELLMASS'): 2885 self['SHELLMASS'] = self['MDOT_GAS']/self['VEL_INFINITY_GAS'] 2886 else: 2887 pass

2888 2889

2890 - def calcSHELLDENS(self):

2891 2892 """ 2893 Calculate the average density of the circumstellar envelope. 2894 2895 """ 2896 2897 if not self.has_key('SHELLDENS'): 2898 2899 mdot_cgs = (float(self['MDOT_GAS'])*u.M_sun/u.yr).to(u.g/u.s).value 2900 self['SHELLDENS'] = mdot_cgs/((self['VEL_INFINITY_GAS']*10**5)\ 2901 *(self['R_STAR']*self.Rsun)**2*4.*pi) 2902 else: 2903 pass

2904 2905 2906

2907 - def calcSHELLCOLDENS(self):

2908 2909 """ 2910 Calculate a proxy for the average column density of the circumstellar 2911 shell. 2912 2913 This is (intuitively) rho * R_STAR, which is important for radiative 2914 excitation (density tracing the source of the radiation, R_STAR setting 2915 the scale of the envelope). Two types of radiative excitation can be 2916 related to this: direct stellar light, and thermal dust emission. 2917 2918 Especially important for water, but in a balance with other excitation 2919 mechanisms. 2920 2921 Note that this quantity is also related to the optical depth through 2922 tau = kappa*coldens. 2923 2924 """ 2925 2926 if not self.has_key('SHELLCOLDENS'): 2927 self['SHELLCOLDENS'] = self['SHELLDENS']\ 2928 *self['R_STAR']*self.Rsun 2929 else: 2930 pass

2931 2932

2933 - def calcSHELLDENS2(self):

2934 2935 """ 2936 Calculate a proxy for the average degree of collisional excitation in 2937 the circumstellar shell. 2938 2939 This is (intuitively) sqrt(rho * rho * R_STAR): two density factors 2940 tracing both collisional partners, and R_STAR setting the scale of the 2941 envelope. 2942 2943 Sqrt is taken for easy comparison between this and the mass-loss rate 2944 to which it is directly proportional. 2945 2946 Especially important for CO, but also to some degree for water where it 2947 is in balance with other excitation mechanisms. 2948 2949 Calculated by taking SHELLDENS**2*R_STAR ~ R_STAR^3/2. 2950 2951 """ 2952 2953 if not self.has_key('SHELLDENS2'): 2954 self['SHELLDENS2'] = sqrt(self['SHELLDENS']**2\ 2955 *self['R_STAR']*self.Rsun) 2956 else: 2957 pass

2958 2959

2960 - def calcDENSFILE(self):

2961 2962 """ 2963 Pointer to the calcDENSTYPE method in case DENSFILE is missing. 2964 2965 """ 2966 2967 self.calcDENSTYPE()

2968 2969 2970

2971 - def calcMRN_DUST(self):

2972 2973 ''' 2974 Set the default value for MRN_DUST to 0. 2975 2976 ''' 2977 2978 if not self.has_key('MRN_DUST'): 2979 self['MRN_DUST'] = 0 2980 else: 2981 pass

2982 2983 2984

2985 - def calcMRN_INDEX(self):

2986 2987 ''' 2988 Set the default value for MRN_INDEX to 3.5 (standard power law in ISM 2989 dust grain size distribution). 2990 2991 ''' 2992 2993 if not self.has_key('MRN_INDEX'): 2994 self['MRN_INDEX'] = 3.5 2995 else: 2996 pass

2997 2998 2999

3000 - def calcMRN_NGRAINS(self):

3001 3002 ''' 3003 Set the default balue for MRN_NGRAINS to the max number of dust species 3004 involved. 3005 3006 This means that all dust species are treated in the mrn treatment of 3007 MCMax. 3008 3009 If the max is set to less species, then the extra species are treated 3010 as normal, with manually set abundances. 3011 3012 ''' 3013 3014 if not self.has_key('MRN_NGRAINS'): 3015 self['MRN_NGRAINS'] = len(self.getDustList()) 3016 else: 3017 pass

3018 3019 3020

3021 - def calcMRN_RMAX(self):

3022 3023 ''' 3024 Set the default value for the maximum grain size in micron. 3025 3026 Abundances of bigger grains will be set to 0. 3027 3028 ''' 3029 3030 if not self.has_key('MRN_RMAX'): 3031 self['MRN_RMAX'] = 1000. 3032 else: 3033 pass

3034 3035 3036

3037 - def calcMRN_RMIN(self):

3038 3039 ''' 3040 Set the default value for the minimum grain size in micron. 3041 3042 Abundances of smaller grains will be set to 0. 3043 3044 ''' 3045 3046 if not self.has_key('MRN_RMIN'): 3047 self['MRN_RMIN'] = 0.01 3048 else: 3049 pass

3050 3051

3052 - def calcSCSET(self):

3053 3054 ''' 3055 Set default of self-consistent settling to False. 3056 3057 ''' 3058 3059 if not self.has_key('SCSET'): 3060 self['SCSET'] = 0 3061 else: 3062 pass

3063 3064

3065 - def calcSCSETEQ(self):

3066 3067 ''' 3068 Set default of self-consistent settling to True. 3069 3070 Only relevant if SCSET == 1. 3071 3072 ''' 3073 3074 if not self.has_key('SCSETEQ'): 3075 self['SCSETEQ'] = 1 3076 else: 3077 pass

3078 3079 3080

3081 - def calcALPHATURB(self):

3082 3083 ''' 3084 Set default of the turbulent mixing strenght to 1e-4. 3085 3086 ''' 3087 3088 if not self.has_key('ALPHATURB'): 3089 self['ALPHATURB'] = 1e-4 3090 else: 3091 pass

3092 3093 3094

3095 - def calcMOLECULE(self):

3096 3097 ''' 3098 Set the MOLECULE keyword to empty list if not given in the input. 3099 3100 ''' 3101 3102 if not self.has_key('MOLECULE'): 3103 self['MOLECULE'] = [] 3104 else: 3105 pass

3106 3107 3108

3109 - def calcGAS_LIST(self):

3110 3111 """ 3112 Set the GAS_LIST keyword based on the MOLECULE keyword. 3113 3114 The input MOLECULE format from the CC input is converted into 3115 Molecule() objects. 3116 3117 """ 3118 3119 if not self.has_key('GAS_LIST') and self['MOLECULE']: 3120 if len(self['MOLECULE'][0]) == 2: 3121 #- First convert the long GASTRoNOoM input molecule names to 3122 #- the short names, since if len() is 2, it comes from 3123 #- PlottingSession.setPacsFromDb 3124 molec_indices \ 3125 = [DataIO.getInputData(keyword='MOLEC_TYPE',make_float=0,\ 3126 filename='Molecule.dat')\ 3127 .index(molec[0]) 3128 for molec in self['MOLECULE']] 3129 molecules_long = [molec[0] for molec in self['MOLECULE']] 3130 self['MOLECULE'] \ 3131 = [[DataIO.getInputData(keyword='TYPE_SHORT',\ 3132 filename='Molecule.dat')[index]] \ 3133 + [molec[1]] 3134 for molec,index in zip(self['MOLECULE'],molec_indices)] 3135 self['TRANSITION'] \ 3136 = [[DataIO.getInputData(keyword='TYPE_SHORT',\ 3137 filename='Molecule.dat')\ 3138 [molec_indices[molecules_long.index(trans[0])]]] \ 3139 + trans[1:] 3140 for trans in self['TRANSITION']] 3141 #- Pull the info from the db 3142 self['GAS_LIST'] = [] 3143 for molec,model_id in self['MOLECULE']: 3144 self['GAS_LIST'].append(Molecule.makeMoleculeFromDb(\ 3145 model_id=model_id,molecule=molec,\ 3146 path_gastronoom=self.path_gastronoom)) 3147 else: 3148 for key,index in zip(['R_OUTER','CHANGE_FRACTION_FILENAME',\ 3149 'SET_KEYWORD_CHANGE_ABUNDANCE',\ 3150 'NEW_TEMPERATURE_FILENAME',\ 3151 'SET_KEYWORD_CHANGE_TEMPERATURE',\ 3152 'ENHANCE_ABUNDANCE_FACTOR',\ 3153 'ABUNDANCE_FILENAME'],\ 3154 [13,16,17,18,19,15,14]): 3155 if self['%s_H2O'%key]: 3156 self['MOLECULE'] = \ 3157 [[(i==index and molec[0] in ['1H1H16O','p1H1H16O',\ 3158 '1H1H17O','p1H1H17O',\ 3159 '1H1H18O','p1H1H18O']) 3160 and self['%s_H2O'%key] 3161 or str(entry) 3162 for i,entry in enumerate(molec)] 3163 for molec in self['MOLECULE']] 3164 #-- Check if startype is not BB, because if starfile is given 3165 # while BB is requested, the starfile cannot be given to the 3166 # Molecule class. 3167 starfile = self['STARTYPE'] != 'BB' and self['STARFILE'] or '' 3168 self['GAS_LIST'] = \ 3169 [Molecule.Molecule(\ 3170 molecule=molec[0],ny_low=int(molec[1]),\ 3171 ny_up=int(molec[2]),nline=int(molec[3]),\ 3172 n_impact=int(molec[4]),n_impact_extra=int(molec[5]),\ 3173 abun_molec=float(molec[6]),\ 3174 abun_molec_rinner=float(molec[7]),\ 3175 abun_molec_re=float(molec[8]),\ 3176 rmax_molec=float(molec[9]),itera=int(molec[10]),\ 3177 lte_request=int(molec[11]),\ 3178 use_collis_radiat_switch=int(molec[12]),\ 3179 abundance_filename=molec[14],\ 3180 enhance_abundance_factor=float(molec[15]),\ 3181 opr=self['OPR'],\ 3182 ratio_12c_to_13c=self['RATIO_12C_TO_13C'],\ 3183 ratio_16o_to_18o=self['RATIO_16O_TO_18O'],\ 3184 ratio_16o_to_17o=self['RATIO_16O_TO_17O'],\ 3185 r_outer=float(molec[13]) \ 3186 and float(molec[13]) \ 3187 or self['R_OUTER_GAS'],\ 3188 outer_r_mode=float(molec[13]) \ 3189 and 'FIXED' \ 3190 or self['OUTER_R_MODE'],\ 3191 dust_to_gas_change_ml_sp=self\ 3192 ['DUST_TO_GAS_CHANGE_ML_SP'],\ 3193 set_keyword_change_abundance=int(molec[17]),\ 3194 change_fraction_filename=molec[16],\ 3195 set_keyword_change_temperature=int(molec[19]),\ 3196 new_temperature_filename=molec[18],\ 3197 use_maser_in_sphinx=self['USE_MASER_IN_SPHINX'],\ 3198 use_no_maser_option=self['USE_NO_MASER_OPTION'],\ 3199 fehler=self['FEHLER'],n_freq=self['N_FREQ'],\ 3200 start_approx=self['START_APPROX'],xdex=self['XDEX'],\ 3201 use_fraction_level_corr=self['USE_FRACTION_LEVEL_CORR'],\ 3202 fraction_level_corr=self['FRACTION_LEVEL_CORR'],\ 3203 number_level_max_corr=self['NUMBER_LEVEL_MAX_CORR'],\ 3204 starfile=starfile,path_gastronoom=self.path_gastronoom) 3205 3206 for molec in self['MOLECULE']] 3207 3208 #- safety check 3209 requested_molecules = set([molec.molecule 3210 for molec in self['GAS_LIST']]) 3211 if not len(self['GAS_LIST']) == len(requested_molecules): 3212 raise IOError('Multiple parameter sets for a single molecule'+\ 3213 ' passed. This is impossible! Contact Robin...') 3214 print 'Gas molecules that are taken into account are ' + \ 3215 ', '.join(sorted([molec[0] for molec in self['MOLECULE']]))+\ 3216 '.' 3217 elif not self.has_key('GAS_LIST') and not self['MOLECULE']: 3218 self['GAS_LIST'] = [] 3219 else: 3220 pass

3221 3222 3223

3224 - def calcR_OUTER_H2O(self):

3225 3226 ''' 3227 Set default value of R_OUTER_H2O to 0. 3228 3229 ''' 3230 3231 if not self.has_key('R_OUTER_H2O'): 3232 self['R_OUTER_H2O'] = 0 3233 else: 3234 pass

3235 3236 3237

3238 - def calcNEW_TEMPERATURE_FILENAME_H2O(self):

3239 3240 ''' 3241 Set default value of NEW_TEMPERATURE_FILENAME_H2O to ''. 3242 3243 ''' 3244 3245 if not self.has_key('NEW_TEMPERATURE_FILENAME_H2O'): 3246 self['NEW_TEMPERATURE_FILENAME_H2O'] = '' 3247 else: 3248 pass

3249 3250 3251

3252 - def calcCHANGE_FRACTION_FILENAME_H2O(self):

3253 3254 ''' 3255 Set default value of CHANGE_FRACTION_FILENAME_H2O to ''. 3256 3257 ''' 3258 3259 if not self.has_key('CHANGE_FRACTION_FILENAME_H2O'): 3260 self['CHANGE_FRACTION_FILENAME_H2O'] = '' 3261 else: 3262 pass

3263 3264 3265

3266 - def calcSET_KEYWORD_CHANGE_TEMPERATURE_H2O(self):

3267 3268 ''' 3269 Set default value of SET_KEYWORD_CHANGE_TEMPERATURE_H2O to ''. 3270 3271 ''' 3272 3273 if not self.has_key('SET_KEYWORD_CHANGE_TEMPERATURE_H2O'): 3274 self['SET_KEYWORD_CHANGE_TEMPERATURE_H2O'] = 0 3275 else: 3276 pass

3277 3278 3279

3280 - def calcSET_KEYWORD_CHANGE_ABUNDANCE_H2O(self):

3281 3282 ''' 3283 Set default value of SET_KEYWORD_CHANGE_ABUNDANCE_H2O to ''. 3284 3285 ''' 3286 3287 if not self.has_key('SET_KEYWORD_CHANGE_ABUNDANCE_H2O'): 3288 self['SET_KEYWORD_CHANGE_ABUNDANCE_H2O'] = 0 3289 else: 3290 pass

3291 3292 3293

3294 - def calcENHANCE_ABUNDANCE_FACTOR_H2O(self):

3295 3296 ''' 3297 Set default value of ENHANCE_ABUNDANCE_FACTOR_H2O to ''. 3298 3299 ''' 3300 3301 if not self.has_key('ENHANCE_ABUNDANCE_FACTOR_H2O'): 3302 self['ENHANCE_ABUNDANCE_FACTOR_H2O'] = 0 3303 else: 3304 pass

3305 3306 3307

3308 - def calcABUNDANCE_FILENAME_H2O(self):

3309 3310 ''' 3311 Set default value of ABUNDANCE_FILENAME_H2O to ''. 3312 3313 ''' 3314 3315 if not self.has_key('ABUNDANCE_FILENAME_H2O'): 3316 self['ABUNDANCE_FILENAME_H2O'] = '' 3317 else: 3318 pass

3319 3320 3321

3322 - def calcR_OUTER_EFFECTIVE(self):

3323 3324 ''' 3325 Get the effective outer radius (either from Mamon, or a fixed value). 3326 3327 ''' 3328 3329 filename = os.path.join(cc.path.gout,'models',\ 3330 self['LAST_GASTRONOOM_MODEL'],\ 3331 'input%s.dat'%self['LAST_GASTRONOOM_MODEL']) 3332 3333 if not self.has_key('R_OUTER_EFFECTIVE'): 3334 self['R_OUTER_EFFECTIVE'] \ 3335 = float(DataIO.readFile(filename=filename,delimiter=' ')[0][4]) 3336 else: 3337 pass

3338 3339 3340

3341 - def calcKEYWORD_DUST_TEMPERATURE_TABLE(self):

3342 3343 ''' 3344 Set KEYWORD_DUST_TEMPERATURE_TABLE to False for now. 3345 3346 If it was not yet defined, there is not ftemperature file anyway. 3347 3348 ''' 3349 3350 if not self.has_key('KEYWORD_DUST_TEMPERATURE_TABLE'): 3351 if self['DUST_TEMPERATURE_FILENAME']: 3352 self['KEYWORD_DUST_TEMPERATURE_TABLE'] = 1 3353 else: 3354 self['KEYWORD_DUST_TEMPERATURE_TABLE'] = 0 3355 else: 3356 pass

3357 3358 3359

3360 - def calcNUMBER_INPUT_DUST_TEMP_VALUES(self):

3361 3362 ''' 3363 Set NUMBER_INPUT_DUST_TEMP_VALUES to length of input file for dust temp 3364 3365 If it does not exist set to 0. 3366 3367 ''' 3368 3369 if not self.has_key('NUMBER_INPUT_DUST_TEMP_VALUES'): 3370 if self['DUST_TEMPERATURE_FILENAME']: 3371 self['NUMBER_INPUT_DUST_TEMP_VALUES'] \ 3372 = len([1 3373 for line in DataIO.readFile(\ 3374 self['DUST_TEMPERATURE_FILENAME']) 3375 if line]) 3376 else: 3377 self['NUMBER_INPUT_DUST_TEMP_VALUES'] = 0 3378 else: 3379 pass

3380 3381 3382

3383 - def calcDUST_TEMPERATURE_FILENAME(self):

3384 3385 """ 3386 Look for the temperature stratification of the star. 3387 3388 If a last mcmax model is available, the filename is given, (for now 2d). 3389 3390 Else an empty string is given, and a power law is used in GASTRoNOoM. 3391 3392 """ 3393 3394 if not self.has_key('DUST_TEMPERATURE_FILENAME'): 3395 filename = self['RID_TEST'] != 'R_STAR' \ 3396 and os.path.join(cc.path.mout,\ 3397 'data_for_gastronoom',\ 3398 '_'.join(['Td',\ 3399 self['LAST_MCMAX_MODEL'],\ 3400 self['RID_TEST']\ 3401 +'.dat']))\ 3402 or os.path.join(cc.path.mout,\ 3403 'data_for_gastronoom',\ 3404 '_'.join(['Td',\ 3405 self['LAST_MCMAX_MODEL']\ 3406 + '.dat'])) 3407 if os.path.isfile(filename): 3408 self['DUST_TEMPERATURE_FILENAME'] = filename 3409 else: 3410 if self['LAST_MCMAX_MODEL']: 3411 rad = self.getDustRad(unit='rstar') 3412 temp = self.getDustTemperature() 3413 if self['RID_TEST'] == 'R_STAR': 3414 temp = temp[rad > 1] 3415 rad = rad[rad > 1] 3416 elif self['RID_TEST'] == 'R_INNER_GAS': 3417 temp = temp[rad > self['R_INNER_GAS']] 3418 rad = rad[rad > self['R_INNER_GAS']] 3419 elif self['RID_TEST'] == 'BUGGED_CASE': 3420 temp = temp[rad > self['R_STAR']] 3421 rad = rad[rad > self['R_STAR']] 3422 self['DUST_TEMPERATURE_FILENAME'] = filename 3423 DataIO.writeCols(filename,[rad,temp]) 3424 self['KEYWORD_DUST_TEMPERATURE_TABLE'] = 1 3425 self['NUMBER_INPUT_DUST_TEMP_VALUES'] = len(rad) 3426 print '** Made dust temperature stratifaction file at %s.'\ 3427 %filename 3428 if self['NUMBER_INPUT_DUST_TEMP_VALUES'] > 999: 3429 ss = '** WARNING! The dust temperature file contains'+\ 3430 ' more than 999 grid points. GASTRoNOoM will '+\ 3431 'fail because it requires less than 1000 points.' 3432 print ss 3433 else: 3434 self['DUST_TEMPERATURE_FILENAME'] = '' 3435 else: 3436 pass

3437 3438 3439

3440 - def calcGAS_LINES(self):

3441 3442 """ 3443 Making transition line input for gas data (auto search) 3444 and additional no-data lines. 3445 3446 The Transition() objects are created then for these lines and added 3447 to the GAS_LINES list. 3448 3449 """ 3450 3451 if not self.has_key('GAS_LINES'): 3452 self['GAS_LINES'] = list() 3453 #-- To make sure the GAS_LIST is done, and the conversion of 3454 # TRANSITION to the right molecule names is done 3455 # (in case of PlottingSession.setPacsFromDb is used) 3456 self.calcGAS_LIST() 3457 3458 #-- Check if specific transition were requested in addition to data 3459 # Note that these include autosearch transitions if requested 3460 # (See ComboCode.py) 3461 if self.has_key('TRANSITION'): 3462 #-- Keep if molecule is available in this model 3463 molecules = [m.molecule for m in self['GAS_LIST']] 3464 self['TRANSITION'] = [trans 3465 for trans in self['TRANSITION'] 3466 if trans[0] in molecules] 3467 3468 #-- Check if identical transitions both with and without n_quad 3469 # are requested: Only keep the one with. (ie in the case of 3470 # radio_autosearch returns a manually requested transition) 3471 trl_ras = [tr for tr in self['TRANSITION'] if len(tr) == 11] 3472 trl_man = [tr for tr in self['TRANSITION'] if len(tr) == 12] 3473 trl_filt = [tr for tr in trl_ras 3474 if tr not in [tt[:-1] for tt in trl_man]] 3475 self['TRANSITION'] = trl_man + trl_filt 3476 3477 #-- transitions taken from radio db have only 11 keys, and thus 3478 # have no n_quad. This is added from the Star() object. 3479 nl = [Transition.makeTransition(star=self,trans=trans) 3480 for trans in self['TRANSITION']] 3481 nl = [trans for trans in nl if trans] 3482 self['GAS_LINES'].extend(nl) 3483 3484 #- Check if molecular line catalogues have to be browsed to create 3485 #- line lists in addition to the data 3486 if self['LINE_SELECT']: 3487 if self['LINE_SELECT'] == 1: 3488 self.__addLineList() 3489 elif self['LINE_SELECT'] == 2: 3490 trl = DataIO.readDict(self['LS_FILE'],\ 3491 multi_keys=['TRANSITION']) 3492 trl_sorted = DataIO.checkEntryInfo(trl['TRANSITION'],11,\ 3493 'TRANSITION') 3494 nl = [Transition.makeTransition(trans=trans,star=self) 3495 for trans in trl_sorted] 3496 nl = [trans for trans in nl if trans] 3497 3498 #-- Kick out those transitions that were requested manually 3499 # Can use GAS_LINES key, as it contains both manual and 3500 # other lines, but the latter are also assigned 3501 # self['N_QUAD'] anyway 3502 ctrl = [tr.getInputString(include_nquad=0) 3503 for tr in self['GAS_LINES']] 3504 nl = [tr for tr in nl 3505 if tr.getInputString(include_nquad=0) not in ctrl] 3506 self['GAS_LINES'].extend(nl) 3507 3508 #-- Sort the transitions. 3509 self['GAS_LINES'] = sorted(list(self['GAS_LINES']),\ 3510 key=lambda x: str(x)) 3511 #-- Check uniqueness. Same N_QUAD double transitions can still occur 3512 self['GAS_LINES'] = Transition.checkUniqueness(self['GAS_LINES']) 3513 3514 #-- Is this still needed? 3515 requested_transitions = set([str(trans) 3516 for trans in self['GAS_LINES']]) 3517 if not len(self['GAS_LINES']) == len(requested_transitions): 3518 print 'Length of the requested transition list: %i'\ 3519 %len(self['GAS_LINES']) 3520 print 'Length of the requested transition list with only ' + \ 3521 'the "transition string" parameters: %i'\ 3522 %len(requested_transitions) 3523 print 'Guilty transitions:' 3524 trans_strings = [str(trans) for trans in self['GAS_LINES']] 3525 print '\n'.join([str(trans) 3526 for trans in self['GAS_LINES'] 3527 if trans_strings.count(str(trans))>1]) 3528 raise IOError('Multiple parameter sets for a single ' + \ 3529 'transition requested. This is impossible! '+ \ 3530 'Check if N_QUAD and TRANSITION definitions ' + \ 3531 'have the same value in your inputfile.' + \ 3532 'If yes, check code/contact Robin.') 3533 else: 3534 pass

3535 3536 3537

3538 - def calcSTARTYPE(self):

3539 3540 """ 3541 Set the default value for STARTYPE, which is the blackbody 3542 assumption (BB). 3543 3544 """ 3545 3546 if not self.has_key('STARTYPE'): 3547 self['STARTYPE'] = 'BB' 3548 else: 3549 pass

3550 3551 3552

3553 - def calcSTARFILE(self):

3554 3555 """ 3556 Set the default value for STARFILE, which is an empty string 3557 (ie STARTYPE is BB, no inputfile). 3558 3559 """ 3560 3561 if not self.has_key('STARFILE'): 3562 if self['STARTYPE'] == 'BB': 3563 self['STARFILE'] = '' 3564 elif self['STARTYPE'] == 'ATMOSPHERE': 3565 modeltypes = ['comarcs','marcs','kurucz'] 3566 modeltype = None 3567 for mt in modeltypes: 3568 if mt in self['ATM_FILENAME']: 3569 modeltype = mt 3570 continue 3571 if modeltype is None: 3572 raise IOError('Atmosphere model type is unknown.') 3573 DataIO.testFolderExistence(cc.path.starf) 3574 atmfile = self['ATM_FILENAME'] 3575 atmos = Atmosphere.Atmosphere(modeltype,filename=atmfile) 3576 atmosmodel = atmos.getModel(teff=self['T_STAR'],\ 3577 logg=self['LOGG']) 3578 starfile = os.path.join(cc.path.starf,'%s_teff%s_logg%s.dat'\ 3579 %(os.path.splitext(atmos.filename)[0],\ 3580 str(atmos.teff_actual),\ 3581 str(atmos.logg_actual))) 3582 if not os.path.isfile(starfile): 3583 savetxt(starfile,atmosmodel,fmt=('%.8e')) 3584 print 'Using input model atmosphere at ' 3585 print starfile 3586 self['STARFILE'] = starfile 3587 elif self['STARTYPE'] == 'TABLE': 3588 self['STARTABLE'] = self['STARTABLE'].strip('"').strip("'") 3589 if not (os.path.isfile(self['STARTABLE']) \ 3590 and os.path.split(self['STARTABLE'])[0]): 3591 self['STARFILE'] = os.path.join(cc.path.starf,\ 3592 self['STARTABLE']) 3593 else: 3594 self['STARFILE'] = self['STARTABLE'] 3595 print 'Using input star spectrum at ' 3596 print self['STARFILE'] 3597 else: 3598 pass

3599 3600 3601

3602 - def calcLINE_SELECT(self):

3603 3604 ''' 3605 If the LINE_SELECT keyword is not present, set to False. 3606 3607 ''' 3608 3609 if not self.has_key('LINE_SELECT'): 3610 self['LINE_SELECT'] = 0 3611 else: 3612 pass

3613 3614 3615

3616 - def calcDUST_TO_GAS(self):

3617 3618 ''' 3619 Calculate the empirical value oft he dust to gas ratio. 3620 3621 ''' 3622 3623 if not self.has_key('DUST_TO_GAS'): 3624 self['DUST_TO_GAS'] = float(self['MDOT_DUST'])\ 3625 *float(self['VEL_INFINITY_GAS'])\ 3626 /float(self['MDOT_GAS'])\ 3627 /float(self['V_EXP_DUST']) 3628 else: 3629 pass

3630 3631 3632

3633 - def calcDUST_TO_GAS_INITIAL(self):

3634 3635 ''' 3636 Set a default value for the initial dust-to-gas ratio at 0.002. 3637 3638 ''' 3639 3640 if not self.has_key('DUST_TO_GAS_INITIAL'): 3641 self['DUST_TO_GAS_INITIAL'] = 0.002 3642 else: 3643 pass

3644 3645 3646

3647 - def calcDUST_TO_GAS_ITERATED(self):

3648 3649 ''' 3650 Fetch the iterated result of the dust-to-gas ratio from cooling. 3651 3652 ''' 3653 3654 if not self.has_key('DUST_TO_GAS_ITERATED'): 3655 try: 3656 filename = os.path.join(cc.path.gout,'models',\ 3657 self['LAST_GASTRONOOM_MODEL'],\ 3658 'input%s.dat'\ 3659 %self['LAST_GASTRONOOM_MODEL']) 3660 self['DUST_TO_GAS_ITERATED'] = float(DataIO.readFile(\ 3661 filename=filename,\ 3662 delimiter=' ')[0][6]) 3663 except IOError: 3664 self['DUST_TO_GAS_ITERATED'] = None 3665 else: 3666 pass

3667 3668 3669

3670 - def calcINCLUDE_SCAT_GAS(self):

3671 3672 ''' 3673 Set the keyword INCLUDE_SCAT_GAS to 0. 3674 3675 The keyword decides whether to take into account the scattering 3676 coefficients in GASTRoNOoM as if they contributed to the absorption 3677 coefficients. 3678 3679 ''' 3680 3681 if not self.has_key('INCLUDE_SCAT_GAS'): 3682 self['INCLUDE_SCAT_GAS'] = 0 3683 else: 3684 pass

3685 3686 3687

3688 - def calcRATIO_12C_TO_13C(self):

3689 3690 ''' 3691 Set default value for ratio_12c_to_13c to 0. 3692 3693 ''' 3694 3695 if not self.has_key('RATIO_12C_TO_13C'): 3696 self['RATIO_12C_TO_13C'] = 0 3697 else: 3698 pass

3699 3700 3701

3702 - def calcRATIO_16O_TO_17O(self):

3703 3704 ''' 3705 Set default value for ratio_16o_to_17o to 0. 3706 3707 ''' 3708 3709 if not self.has_key('RATIO_16O_TO_17O'): 3710 self['RATIO_16O_TO_17O'] = 0 3711 else: 3712 pass

3713 3714 3715

3716 - def calcRATIO_16O_TO_18O(self):

3717 3718 ''' 3719 Set default value for ratio_16o_to_18o to 0. 3720 3721 ''' 3722 3723 if not self.has_key('RATIO_16O_TO_18O'): 3724 self['RATIO_16O_TO_18O'] = 0 3725 else: 3726 pass

3727 3728 3729

3730 - def calcOPR(self):

3731 3732 ''' 3733 Set default value for opr to 0. 3734 3735 ''' 3736 3737 if not self.has_key('OPR'): 3738 self['OPR'] = 0 3739 else: 3740 pass

3741 3742 3743

3744 - def calcUSE_NEW_DUST_KAPPA_FILES(self):

3745 3746 ''' 3747 Set the default value of USE_NEW_DUST_KAPPA_FILES to 1. 3748 3749 ''' 3750 3751 if not self.has_key('USE_NEW_DUST_KAPPA_FILES'): 3752 self['USE_NEW_DUST_KAPPA_FILES'] = 1 3753 else: 3754 pass

3755 3756 3757

3758 - def calcTEMDUST_FILENAME(self):

3759 3760 """ 3761 Making extinction efficiency input files for GASTRoNOoM from MCMax 3762 output mass extinction coefficients. 3763 3764 If no MCMax output available, this file is temdust.kappa, the standard. 3765 3766 In units of cm^-1, Q_ext/a. 3767 3768 """ 3769 3770 if not self.has_key('TEMDUST_FILENAME'): 3771 if self['NLAM'] > 2000: 3772 #- For now not supported, GASTRoNOoM cannot take more than 2000 3773 #- wavelength opacity points 3774 raise IOError('NLAM > 2000 not supported due to GASTRoNOoM '+\ 3775 'opacities!') 3776 filename = os.path.join(cc.path.gdata,'temdust_%s.dat'\ 3777 %self['LAST_MCMAX_MODEL']) 3778 if not int(self['USE_NEW_DUST_KAPPA_FILES']) \ 3779 or not self['LAST_MCMAX_MODEL']: 3780 self['TEMDUST_FILENAME'] = 'temdust.kappa' 3781 elif os.path.isfile(filename): 3782 self['TEMDUST_FILENAME'] = os.path.split(filename)[1] 3783 else: 3784 try: 3785 wavelength,q_ext = self.getWeightedKappas() 3786 q_ext *= self['SPEC_DENS_DUST']*(4.0/3) 3787 wavelength = list(wavelength) 3788 wavelength.reverse() 3789 q_ext = list(q_ext) 3790 q_ext.reverse() 3791 self['TEMDUST_FILENAME'] = os.path.split(filename)[1] 3792 DataIO.writeCols(filename,[wavelength,q_ext]) 3793 print '** Made opacity file at ' + filename + '.' 3794 except IOError: 3795 self['TEMDUST_FILENAME'] = 'temdust.kappa' 3796 else: 3797 pass

3798 3799 3800

3801 - def calcR_OH1612_AS(self):

3802 3803 ''' 3804 Set the R_OH1612_AS to the default value of 0 as. 3805 3806 ''' 3807 3808 if not self.has_key('R_OH1612_AS'): 3809 self['R_OH1612_AS'] = 0 3810 else: 3811 pass

3812 3813 3814

3815 - def calcR_OH1612(self):

3816 3817 ''' 3818 Calculate the R_OH1612 in R_STAR. 3819 3820 ''' 3821 3822 if not self.has_key('R_OH1612'): 3823 #-- Get the angular size equivalency 3824 equiv = eq.angularSize(self['DISTANCE']) 3825 #-- Convert as to cm 3826 rad = (self['R_OH1612_AS']*u.arcsec).to(u.cm,equivalencies=equiv) 3827 self['R_OH1612'] = rad.value/self['R_STAR']/self.Rsun 3828 else: 3829 pass

3830 3831 3832

3833 - def calcR_OH1612_NETZER(self):

3834 3835 ''' 3836 Calculate the radial OH maser peak distance in cm. 3837 3838 Taken from Netzer & Knapp 1987, eq. 29. 3839 3840 The interstellar radiation factor is taken as A = 5.4 3841 (avg Habing Field) 3842 3843 ''' 3844 3845 if not self.has_key('R_OH1612_NETZER'): 3846 mg = self['MDOT_GAS']/1e-5 3847 vg = self['VEL_INFINITY_GAS'] 3848 self['R_OH1612_NETZER'] = ((5.4*mg**0.7/vg**0.4)**-4.8\ 3849 + (74.*mg/vg)**-4.8)**(-1/4.8)*1e16\ 3850 /self['R_STAR']/self.Rsun 3851 else: 3852 pass

3853 3854 3855

3856 - def missingInput(self,missing_key):

3857 3858 """ 3859 Try to resolve a missing key. 3860 3861 @param missing_key: the missing key for which an attempt will be made 3862 to calculate its value based on already present 3863 parameters 3864 @type missing_key: string 3865 3866 """ 3867 3868 if missing_key in ('T_STAR','L_STAR','R_STAR'): 3869 self.calcTLR() 3870 elif missing_key in ['R_MAX_' + species 3871 for species in self.getDustList()]: 3872 self.calcR_MAX(missing_key) 3873 elif missing_key in ['R_DES_' + species 3874 for species in self.getDustList()]: 3875 self.checkT() 3876 elif missing_key in ['T_DESA_' + species 3877 for species in self.getDustList()] + \ 3878 ['T_DESB_' + species 3879 for species in self.getDustList()]: 3880 self.calcT_DES(missing_key[7:]) 3881 elif missing_key in ['T_DES_' + species 3882 for species in self.getDustList()]: 3883 self.checkT() 3884 elif hasattr(self,'calc' + missing_key): 3885 getattr(self,'calc' + missing_key)() 3886 else: 3887 pass

3888

Source Code for Module ComboCode.cc.modeling.objects.Star