ComboCode.cc.modeling.objects.Molecule

88 89 ''' 90 A class to deal with molecules in GASTRoNOoM. 91 92 ''' 93

94 - def __init__(self,molecule,ny_up=0,ny_low=0,nline=0,n_impact=0,\ 95 n_impact_extra=0,abun_molec=1.0e-10,abun_molec_rinner=1.0e-10,\ 96 abun_molec_re=1.0e-10,rmax_molec=1.,itera=0,lte_request=None,\ 97 use_collis_radiat_switch=0,dust_to_gas_change_ml_sp=0,\ 98 use_no_maser_option=0,use_maser_in_sphinx=1,\ 99 fehler=1e-4,xdex=2.,n_freq=30,start_approx=0,\ 100 use_fraction_level_corr=1,fraction_level_corr=0.8,\ 101 number_level_max_corr=1e-12,\ 102 ratio_12c_to_13c=0,ratio_16o_to_17o=0,ratio_16o_to_18o=0,\ 103 opr=0,r_outer=0,outer_r_mode='MAMON',abundance_filename=None,\ 104 change_fraction_filename=None,set_keyword_change_abundance=0,\ 105 set_keyword_change_temperature=0,enhance_abundance_factor=0,\ 106 new_temperature_filename=None,linelist=0,starfile='',\ 107 path_gastronoom=None):

108 109 ''' 110 Initiate a Molecule class, setting all values for the allowed 111 transition parameters to zero (by default). 112 113 @param molecule: shorthand name of the molecule 114 @type molecule: string 115 116 @keyword ny_up: number of levels in first vibration state v=1 117 118 (default: 0) 119 @type ny_up: int 120 @keyword ny_low: number of levels in the ground vibration state v=0 121 122 (default: 0) 123 @type ny_low: int 124 @keyword nline: number of allowed transitions in molecule 125 126 (default: 0) 127 @type nline: int 128 @keyword n_impact: number of depth points in radius mesh 129 130 (default: 0) 131 @type n_impact: int 132 @keyword n_impact_extra: number of depth points in radius_mesh 133 (< n_impact) used for variable mass-loss 134 (0 if constant mdot) 135 136 (default: 0) 137 @type n_impact_extra: int 138 @keyword itera: number of iterations in mline for molecule, LTE 139 approximation if zero 140 141 (default: 0) 142 @type itera: string 143 @keyword abun_molec: molecular abundance at the stellar radius. 144 Default is arbitrary, and used if molecule is co 145 or h2o. 146 147 (default: 1.0e-10) 148 @type abun_molec: float 149 @keyword abun_molec_rinner: molecular abundance at inner shell radius. 150 Default is arbitrary, and used if molecule 151 is co or h2o. 152 153 (default: 1.0e-10) 154 @type abun_molec_rinner: float 155 @keyword abun_molec_re: molecular abundance at rmax_molec. Default is 156 arbitrary, and used if molecule is co or h2o. 157 158 (default: 1.0e-10) 159 @type abun_molec_re: float 160 @keyword rmax_molec: The radius from which the Willacy abundance 161 profiles are used. They are rescaled to 162 abun_molec_re. Default is arbitrary, and used if 163 molecule is co or h2o. 164 165 (default: 1.) 166 @type rmax_molec: float 167 @keyword use_collis_radiat_switch: in case of unstable mline, such as 168 for para h2o sometimes 169 170 (default: 0) 171 @type use_collis_radiat_switch: bool 172 @keyword fehler: convergence criterium in mline 173 174 (default: 1e-4) 175 @type fehler: float 176 @keyword xdex: Controls the distribution of the impact parameters in the 177 interval between R_STAR and R_OUTER. 178 179 (default: 2.) 180 @type xdex: float 181 @keyword n_freq: Number of frequency points in line profile 182 183 (default: 30) 184 @type n_freq: int 185 @keyword start_approx: set to 0 when one wants to start with LTE-approx 186 as starting n(NY,N_IMPACT); set to 1 when 187 starting from another model - with same NY, 188 N_IMPACT, ... 189 190 (default: 0) 191 @type start_approx: bool 192 @keyword use_fraction_level_corr: set to 1 if one wants to put a limit 193 on the level-population correction 194 (BES3). 195 196 (default: 1) 197 @type use_fraction_level_corr: bool 198 @keyword fraction_level_corr: user-defined fraction for maximum change 199 in level-population correction; useful in 200 case of H2O 201 202 (default: 0.8) 203 @type fraction_level_corr: float 204 @keyword number_level_max_corr: user-defined level population. Only the 205 level corrections for levels with a 206 higher level population will be used to 207 determine convergence criterion 208 209 (default: 1e-12) 210 @type number_level_max_corr: float 211 @keyword ratio_12c_to_13c: 12c/13c ratio, only relevant for 13co and 212 other molecules with that isotope 213 214 (default: 0) 215 @type ratio_12c_to_13c: int 216 @keyword ratio_16o_to_17o: 16o/17o ratio, only relevant for h2_17o and 217 other molecules with that isotope 218 219 (default: 0) 220 @type ratio_16o_to_17o: int 221 @keyword ratio_16o_to_18o: 16o/18o ratio, only relevant for h2_18o and 222 other molecules with that isotope 223 224 (default: 0) 225 @type ratio_16o_to_18o: int 226 @keyword opr: ortho-to-para water ratio, only relevant for ph2o, 227 ph2_17o,ph2_18o and other molecules with para h2o 228 229 (default: 0) 230 @type opr: int 231 @keyword r_outer: the outer radius of the shell for this molecule, 232 0 if MAMON 233 234 (default: 0) 235 @type r_outer: float 236 @keyword lte_request: using LTE in mline only (with or without 237 collision rates: determined by itera), if default 238 lte_request is 0 if itera != 0 and 1 if itera ==0 239 240 (default: 0) 241 @type lte_request: bool 242 @keyword outer_r_mode: the mode used for calculating r_outer 243 (FIXED or MAMON) 244 245 (default: 'MAMON') 246 @type outer_r_mode: string 247 @keyword dust_to_gas_change_ml_sp: if 0 not used, otherwise this is an 248 alternative value for the 249 dust-to-gas ratio in mline/sphinx 250 for this molecule and its 251 transitions. 252 253 (default: 0) 254 @type dust_to_gas_change_ml_sp: float 255 @keyword abundance_filename: if enhance_abundance_factor is not zero, 256 this includes the filename and/or path 257 to the file that includes the profile. 258 259 (default: None) 260 @type abundance_filename: string 261 @keyword enhance_abundance_factor: if 0 the Willacy abundance profiles 262 are uses, if not zero the 263 abundance_filename is used and 264 scaled with the factor given here. 265 THIS DOES NOT RESCALE ABUNDANCES BY 266 WILLACY! Only used for filename 267 abundances, hence why this parameter 268 also turns this feature on/off 269 270 (default: 0) 271 @type enhance_abundance_factor: float 272 @keyword set_keyword_change_abundance: Change the abundance calculated 273 in cooling by a radius dependent 274 factor 275 276 (default: 0) 277 @type set_keyword_change_abundance: bool 278 @keyword change_fraction_filename: the filename of the enhancement 279 factors if 280 set_keyword_change_abundance != 0 281 282 (default: None) 283 @type change_fraction_filename: string 284 @keyword set_keyword_change_temperature: Use a different temperature 285 structure in mline and sphinx 286 287 (default: 0) 288 @type set_keyword_change_temperature: bool 289 @keyword new_temperature_filename: the filename for the temperature 290 structure if 291 set_keyword_change_temperature != 0 292 293 (default: None) 294 @type new_temperature_filename: string 295 @keyword use_no_maser_option: Do not allow masers (neg opacs) in mline 296 RT by setting negative line opacs to 1e-60 297 If use_maser_in_sphinx is on, mline 298 will do a final run including masers 299 anyway to see what would happen if they 300 were inluded by allowing negative opacs 301 for the line profile calculations in 302 sphinx (but not for the convergence in 303 mline). 304 305 (default: 0) 306 @type use_no_maser_option: bool 307 @keyword use_maser_in_sphinx: When on, does a final mline run including 308 masers, allowing negative opacities. When 309 off, sets the masing line opacities to 310 1e-60 when writing out the ml3 file. 311 312 (default: 1) 313 @type use_maser_in_sphinx: bool 314 @keyword linelist: The molecule is created for the LineList module. No 315 radiative information is read from GASTRoNOoM input 316 files. 317 318 (default: 0) 319 @type linelist: bool 320 @keyword starfile: input filename for a stellar input spectrum (either 321 user defined or from a model atmosphere spectrum) 322 323 (default: '') 324 @type starfile: str 325 @keyword path_gastronoom: model output folder in the GASTRoNOoM home 326 327 (default: None) 328 @type path_gastronoom: string 329 330 ''' 331 332 self.molecule = str(molecule) 333 self.ny_up = int(ny_up) 334 self.ny_low = int(ny_low) 335 self.nline = int(nline) 336 self.n_impact = int(n_impact) 337 self.n_impact_extra = int(n_impact_extra) 338 self.path_gastronoom = path_gastronoom 339 340 self.molecule_index = DataIO.getInputData(keyword='TYPE_SHORT',\ 341 filename='Molecule.dat')\ 342 .index(self.molecule) 343 mdata = ['MOLEC_TYPE','NAME_SHORT','NAME_PLOT',\ 344 'SPEC_INDICES','USE_INDICES_DAT'] 345 attrs = ['molecule_full','molecule_short','molecule_plot',\ 346 'spec_indices','use_indices_dat'] 347 mfloat = [0,0,0,1,1] 348 for k,a,mf in zip(mdata,attrs,mfloat): 349 setattr(self,a,DataIO.getInputData(keyword=k,make_float=mf,\ 350 filename='Molecule.dat',\ 351 rindex=self.molecule_index,)) 352 353 self.itera = int(itera) 354 #- lte_request may be undefined, but then it would be faulty input, 355 #- where we do want the code to crash... 356 if self.itera==0 and lte_request is None: 357 self.lte_request = 1 358 #- Normally u never use lte if taking into account collision rates 359 elif self.itera != 0: 360 self.lte_request = 0 361 elif self.itera==0 and lte_request is not None: 362 self.lte_request = lte_request 363 self.abun_molec = abun_molec 364 self.abun_molec_rinner = abun_molec_rinner 365 self.abun_molec_re = abun_molec_re 366 self.rmax_molec = rmax_molec 367 self.use_collis_radiat_switch = int(use_collis_radiat_switch) 368 self.ratio_12c_to_13c = ratio_12c_to_13c 369 self.ratio_16o_to_17o = ratio_16o_to_17o 370 self.ratio_16o_to_18o = ratio_16o_to_18o 371 self.opr = opr 372 self.dust_to_gas_change_ml_sp = float(dust_to_gas_change_ml_sp) 373 self.use_no_maser_option = int(use_no_maser_option) 374 self.use_maser_in_sphinx = int(use_maser_in_sphinx) 375 self.fehler = fehler 376 self.xdex = xdex 377 self.n_freq = int(n_freq) 378 self.start_approx = int(start_approx) 379 self.use_fraction_level_corr = int(use_fraction_level_corr) 380 self.fraction_level_corr = fraction_level_corr 381 self.number_level_max_corr = number_level_max_corr 382 383 #-- Set the molecule inputfiles for abundance and temperature, applying 384 # the path from Path.dat if the file does not exist or the path to the 385 # file is not given. (eg a subfolder might be given, but that works) 386 for k in ['abundance_filename','change_fraction_filename',\ 387 'new_temperature_filename']: 388 fn = locals()[k] 389 if fn and not (os.path.isfile(fn) and os.path.split(fn)[0]): 390 fn = os.path.join(cc.path.molf,fn) 391 setattr(self,k,fn) 392 else: 393 setattr(self,k,fn) 394 self.enhance_abundance_factor = float(enhance_abundance_factor) 395 self.set_keyword_change_abundance = int(set_keyword_change_abundance) 396 self.set_keyword_change_temperature = \ 397 int(set_keyword_change_temperature) 398 399 #-- Mainly for plotting purposes: The relative, multiplicative abundance 400 # factor with respect to main isotope (and OPR) is calculated 401 # This does not take into account enhance_abundance_factor! 402 self.abun_factor = self.getAbunFactor() 403 self.outer_r_mode = outer_r_mode 404 if self.outer_r_mode == 'MAMON': self.r_outer = 0 405 else: self.r_outer = float(r_outer) 406 self.__model_id = None 407 if not linelist: 408 if self.use_indices_dat: 409 tag = '_'.join([self.molecule,str(self.ny_low),\ 410 str(self.ny_up),str(self.nline)]) 411 i = DataIO.getInputData(start_index=4,keyword='MOLECULE',\ 412 filename='Indices.dat').index(tag) 413 fn = DataIO.getInputData(path=cc.path.usr,keyword='RADIAT',\ 414 filename='Indices.dat',start_index=4,\ 415 rindex=i) 416 fn = os.path.join(cc.path.gdata,'radiat_backup',fn) 417 else: 418 fn = os.path.join(cc.path.gdata,'%s_radiat.dat'%self.molecule) 419 420 self.radiat = RadiatReader.RadiatReader(fn=fn,nline=self.nline, 421 ny=self.ny_up+self.ny_low) 422 if self.spec_indices: 423 if self.use_indices_dat: 424 f = DataIO.getInputData(path=cc.path.usr,start_index=4,\ 425 keyword='INDICES',rindex=i,\ 426 filename='Indices.dat') 427 filename = os.path.join(cc.path.gdata,'indices_backup',f) 428 else: 429 filename = os.path.join(cc.path.gdata,\ 430 '{}_indices.dat'.format(self.molecule)) 431 rf = DataIO.readFile(filename,' ') 432 self.radiat_indices = [[int(i) for i in line] for line in rf] 433 else: 434 self.radiat = None 435 self.radiat_indices = None 436 self.starfile = starfile 437 self.mline = None

438 439 440

441 - def __str__(self):

442 443 ''' 444 Printing a molecule as it should appear in the GASTRoNOoM input file. 445 446 ''' 447 448 ll = [self.molecule_full,self.ny_low,self.ny_up,self.nline,\ 449 self.n_impact,self.n_impact_extra] 450 return 'MOLECULE={} {:d} {:d} {:d} {:d} {:d}'.format(*ll)

451 452 453

454 - def __eq__(self,other):

455 456 ''' 457 Compare two molecules and return true if equal. 458 459 The condition is the string representation of this Molecule(). Note 460 that the other properties included in the self.makeDict() dictionary are 461 not compared. Those properties do not determine whether a molecule is 462 equal or not. 463 464 In this sense equal refers to the spectroscopy included for this 465 molecule, the number of impact parameters, and the molecule itself, of 466 course. 467 468 @return: The comparison 469 @rtype: bool 470 471 ''' 472 473 try: 474 if str(self) == str(other): 475 return True 476 else: 477 return False 478 except AttributeError: 479 return False

480 481 482

483 - def __ne__(self,other):

484 485 ''' 486 Compare two molecules and return true if not equal. 487 488 The condition is the string representation of this Molecule(). Note 489 that the other properties included in the self.makeDict() dictionary are 490 not compared. Those properties do not determine whether a molecule is 491 equal or not. 492 493 In this sense equal refers to the spectroscopy included for this 494 molecule, the number of impact parameters, and the molecule itself, of 495 course. 496 497 @return: The negative comparison 498 @rtype: bool 499 500 ''' 501 502 try: 503 if str(self) != str(other): 504 return True 505 else: 506 return False 507 except AttributeError: 508 return True

509 510 511

512 - def __hash__(self):

513 514 ''' 515 Return a hash number based on the string of the molecule. 516 517 The condition is the string representation of this Molecule(). Note 518 that the other properties included in the self.makeDict() dictionary are 519 not compared. Those properties do not determine whether a molecule is 520 equal or not. 521 522 In this sense equal refers to the spectroscopy included for this 523 molecule, the number of impact parameters, and the molecule itself, of 524 course. 525 526 @return: The hash number: 527 @rtype: int 528 529 ''' 530 531 return hash(str(self))

532 533 534

535 - def updateParameters(self,pardict):

536 537 ''' 538 Update parameters. 539 540 @param pardict: the parameters with respective values for the update 541 @type pardict: dict() 542 543 ''' 544 545 for k,v in pardict.items(): 546 if hasattr(self,k.lower()): setattr(self,k.lower(),v)

547 548 549

550 - def makeDict(self,path=None,in_progress=0):

551 552 ''' 553 Return a dict with molecule string, and other relevant parameters. 554 555 @keyword path: If a different path is needed, it can be passed here, 556 for files. For instance, when making dictionaries for 557 Molecule() objects in the case of supercomputer copies. 558 559 (default: None) 560 @type path: string 561 @keyword in_progress: add an extra dict entry "IN_PROGRESS" if the 562 molecule is still being calculated somewhere. 563 564 (default: 0) 565 @type in_progress: bool 566 567 @return: The molecule dictionary including all relevant, defining 568 information 569 @rtype: dict() 570 571 ''' 572 573 dd = dict([('MOLECULE',str(self).replace('MOLECULE=','')),\ 574 ('ITERA',self.itera),\ 575 ('ABUN_MOLEC',self.abun_molec),\ 576 ('ABUN_MOLEC_RINNER',self.abun_molec_rinner),\ 577 ('ABUN_MOLEC_RE', self.abun_molec_re),\ 578 ('RMAX_MOLEC',self.rmax_molec),\ 579 ('LTE_REQUEST',self.lte_request),\ 580 ('OUTER_R_MODE',self.outer_r_mode),\ 581 ('USE_COLLIS_RADIAT_SWITCH',self.use_collis_radiat_switch),\ 582 ('R_OUTER',self.r_outer),\ 583 ('USE_NO_MASER_OPTION',self.use_no_maser_option),\ 584 ('USE_MASER_IN_SPHINX',self.use_maser_in_sphinx),\ 585 ('FEHLER',self.fehler),\ 586 ('XDEX',self.xdex),\ 587 ('N_FREQ',self.n_freq),\ 588 ('START_APPROX',self.start_approx),\ 589 ('USE_FRACTION_LEVEL_CORR',self.use_fraction_level_corr),\ 590 ('FRACTION_LEVEL_CORR',self.fraction_level_corr),\ 591 ('NUMBER_LEVEL_MAX_CORR',self.number_level_max_corr) 592 ]) 593 594 for par,isot in [('RATIO_16O_TO_18O','18O'),\ 595 ('RATIO_16O_TO_17O','17O'),\ 596 ('RATIO_12C_TO_13C','13C'),\ 597 ('OPR','p1H')]: 598 if isot in self.molecule: 599 dd[par] = getattr(self,par.lower()) 600 if self.dust_to_gas_change_ml_sp: 601 dd['CHANGE_DUST_TO_GAS_FOR_ML_SP'] = 1 602 dd['DUST_TO_GAS_CHANGE_ML_SP'] = self.dust_to_gas_change_ml_sp 603 604 if self.enhance_abundance_factor: 605 dd['ENHANCE_ABUNDANCE_FACTOR'] = self.enhance_abundance_factor 606 niav = len(DataIO.readCols(filename=self.abundance_filename)[0]) 607 dd['NUMBER_INPUT_ABUNDANCE_VALUES'] = niav 608 dd['KEYWORD_TABLE'] = 1 609 moltab = self.molecule_full[:self.molecule_full.index('.')] 610 dd['MOLECULE_TABLE'] = moltab 611 dd['ISOTOPE_TABLE'] = self.molecule_full 612 if not path is None: 613 afn = os.path.join(path,\ 614 os.path.split(self.abundance_filename)[1]) 615 else: 616 afn = self.abundance_filename 617 dd['ABUNDANCE_FILENAME'] = '"{}"'.format(afn) 618 619 if self.set_keyword_change_abundance: 620 dd['SET_KEYWORD_CHANGE_ABUNDANCE'] \ 621 = self.set_keyword_change_abundance 622 if not path is None: 623 cffn = os.path.join(path,\ 624 os.path.split(self.change_fraction_filename)[1]) 625 else: 626 cffn = self.change_fraction_filename 627 dd['CHANGE_FRACTION_FILENAME'] = '"{}"'.format(cffn) 628 629 if self.set_keyword_change_temperature: 630 dd['SET_KEYWORD_CHANGE_TEMPERATURE'] \ 631 = self.set_keyword_change_temperature 632 if not path is None: 633 tfn = os.path.join(path,\ 634 os.path.split(self.new_temperature_filename)[1]) 635 else: 636 tfn = self.new_temperature_filename 637 dd['NEW_TEMPERATURE_FILENAME'] = '"{}"'.format(tfn) 638 639 if self.starfile: 640 dd['USE_STARFILE'] = 1 641 if not path is None: 642 sfn = os.path.join(path,os.path.split(self.starfile)[1]) 643 644 else: 645 sfn = self.starfile 646 dd['STARFILE'] = '"{}"'.format(sfn) 647 648 if int(in_progress): 649 dd['IN_PROGRESS'] = 1 650 651 return dd

652 653 654

655 - def makeMlineFilename(self,number='*',include_path=0):

656 657 ''' 658 Return an mline filename for this object. 659 660 @keyword number: the number in the filename (ml*, ml1, ml2, ml3. Hence 661 can be *, 1, 2, 3) 662 663 (default: '*') 664 @type number: string 665 @keyword include_path: Include the full filepath. 666 667 (default: 0) 668 @type include_path: bool 669 670 @return: The sphinx filename for this transition 671 @rtype: string 672 673 ''' 674 675 try: 676 number = str(int(number)) 677 except ValueError: 678 number = str(number) 679 680 fn = 'ml{}{}_{}.dat'.format(number,self.getModelId(),self.molecule) 681 if include_path: 682 fn = os.path.join(cc.path.gastronoom,self.path_gastronoom,'models',\ 683 self.getModelId(),fn) 684 685 return fn

686 687 688

689 - def readMline(self):

690 691 ''' 692 Read the mline output if the model id is valid. 693 694 The mline output is available in the MlineReader object mline, as a 695 property of Molecule(). 696 697 ''' 698 699 if self.mline is None and self.getModelId(): 700 fn = self.makeMlineFilename(include_path=1) 701 self.mline = MlineReader.MlineReader(fn)

702 703 704

705 - def setModelId(self,model_id):

706 707 ''' 708 Set a model_id for the molecule, which identifies model_id for MLINE! 709 710 @param model_id: The model_d to be associated with this molecule 711 @type model_id: string 712 713 ''' 714 715 self.__model_id = model_id

716 717 718

719 - def makeLabel(self):

720 721 ''' 722 Return a short-hand label for this particular molecule, 723 taken from the molecule.dat file. 724 725 @return: The label associated with this molecule 726 @rtype: string 727 728 ''' 729 730 return '\ {}\ '.format(self.molecule_plot)

731 732 733

734 - def getModelId(self):

735 736 ''' 737 Return the model_id associated with this molecule. None if not yet set. 738 739 @return: the model id is returned. 740 @rtype: string 741 742 ''' 743 744 return self.__model_id

745 746 747

748 - def isMolecule(self):

749 750 ''' 751 Return True to help the codes know that this is a molecule 752 and not a transition. 753 754 @return: True if molecule or False if transition 755 @rtype: bool 756 757 ''' 758 759 return True

760 761 762

763 - def getAbunFactor(self):

764 765 ''' 766 Return the abundance factor of the molecule, with respect to its main 767 isotope/ortho version. 768 769 ''' 770 771 raw_factors = DataIO.getInputData(keyword='ABUN_FACTOR',\ 772 filename='Molecule.dat',\ 773 rindex=self.molecule_index) 774 factors = [factor == '1' and 1 or float(getattr(self,factor.lower())) 775 for factor in raw_factors.split('*')] 776 #-- abun_factor should only take into account isotope/OPR factors 777 #if float(self.enhance_abundance_factor): 778 # factors.append(float(self.enhance_abundance_factor)) 779 total_factor = 1 780 while factors: 781 total_factor *= factors.pop() 782 return total_factor

783 784 785

786 - def isWater(self):

787 788 ''' 789 Is this molecule a water molecule? 790 791 (ortho, para, isotopologue thereof included) 792 793 @return: True or False 794 @rtype: bool 795 796 ''' 797 798 return self.molecule in ['1H1H16O','p1H1H16O','1H1H17O',\ 799 'p1H1H17O','1H1H18O','p1H1H18O']

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