ComboCode.cc.plotting.objects.PlotDust

31 32 """ 33 Plotting environment for SEDs and all dust parameters. 34 35 """ 36

37 - def __init__(self,star_name='model',sed=None,path_mcmax='',\ 38 inputfilename=None,fn_add_star=0):

39 40 ''' 41 Initializing PlotDust session. 42 43 @keyword star_name: name of the star from Star.dat, use default only 44 when never using any star model specific things 45 46 (default: "model") 47 @type star_name: string 48 @keyword path_mcmax: Output modeling folder in MCMax home folder 49 50 (default: '') 51 @type path_mcmax: string 52 @keyword inputfilename: name of inputfile that is also copied to the 53 output folder of the plots, 54 if None nothing is copied 55 56 (default: None) 57 @type inputfilename: string 58 @keyword sed: an Sed object needed for plotting the SED. None if not 59 plotting an Sed, but just dust parameters such as \ 60 temperature 61 62 (default: None) 63 @type sed: Sed() 64 @keyword fn_add_star: Add the star name to the requested plot filename. 65 Only relevant if fn_plt is given in a sub method. 66 67 (default: 0) 68 @type fn_add_star: bool 69 70 ''' 71 72 super(PlotDust, self).__init__(star_name=star_name,\ 73 path=path_mcmax,\ 74 code='MCMax',\ 75 inputfilename=inputfilename,\ 76 fn_add_star=fn_add_star) 77 #-- Convenience path 78 cc.path.mout = os.path.join(cc.path.mcmax,self.path) 79 self.sed = sed

80 81 82

83 - def plotSed(self,star_grid=[],fn_plt='',cfg='',iterative=0,no_models=0,\ 84 show_phot_filter=0,**kwargs):

85 86 """ 87 Creating an SED with 0, 1 or more models and data. 88 89 Includes data preparation on the spot. 90 91 Additional plotCols keywords can be passed through kwargs. They 92 overwrite cfg input. 93 94 @keyword star_grid: list of Star() models to plot. If star_grid is [], 95 only data are plotted. 96 97 (default: []) 98 @type star_grid: list[Star()] 99 @keyword fn_plt: A base plot filename. Includes folder. If not, a 100 default is added 101 102 (default: '') 103 @type fn_plt: string 104 @keyword cfg: path to the Plotting2.plotCols config file. If default, 105 the hard-coded default plotting options are used. 106 107 (default: '') 108 @type cfg: string 109 @keyword iterative: add an extra suffix to the filename for each 110 iteratively calculated model, with this number 111 giving the model muber (index in star_grid), 112 0 if not used. 113 114 (default: 0) 115 @type iterative: int 116 @keyword no_models: Only show data. 117 118 (default: 0) 119 @type no_models: bool 120 @keyword show_phot_filter: Show the wavelength band of the photometric 121 filters as an x error bar on the model phot 122 123 (default: 0) 124 @type show_phot_filter: bool 125 126 """ 127 128 if self.sed is None: 129 print 'No dsed given in Path.dat. Cannot plot SED. Aborting...' 130 return 131 print '***********************************' 132 print '** Creating SED plot.' 133 134 cfg_dict = Plotting2.readCfg(cfg) 135 if cfg_dict.has_key('no_models'): 136 no_models = cfg_dict['no_models'] 137 if cfg_dict.has_key('show_phot_filter'): 138 show_phot_filter = bool(cfg_dict['show_phot_filter']) 139 if cfg_dict.has_key('filename'): 140 fn_plt = cfg_dict.pop('filename') 141 142 data_labels = dict([(dt,(n,ls)) 143 for n,dt,ls in zip(DataIO.getInputData(path=cc.path.usr,\ 144 keyword='PLOT_NAMES',\ 145 filename='Sed.dat',\ 146 remove_underscore=1),\ 147 DataIO.getInputData(path=cc.path.usr,\ 148 keyword='DATA_TYPES',\ 149 filename='Sed.dat'),\ 150 DataIO.getInputData(path=cc.path.usr,\ 151 keyword='LINE_TYPES',\ 152 filename='Sed.dat'))]) 153 154 plot_title='SED %s'%self.star_name_plots 155 156 #- prepare and collect data, keytags and line types 157 keytags = [] 158 data_x = [] 159 data_y = [] 160 data_xerr = [] 161 data_yerr = [] 162 line_types = [] 163 for (dt,fn),tdata in sorted([dset 164 for dset in self.sed.data.items() 165 if 'PHOT' not in dset[0][0].upper()]): 166 keytags.append(data_labels[dt][0]) 167 data_x.append(tdata[0]) 168 data_y.append(tdata[1]) 169 #data_err.append(tdata[2]) 170 #-- For now, no error bars for spectra. 171 data_xerr.append(None) 172 data_yerr.append(None) 173 line_types.append(data_labels[dt][1]) 174 175 for (dt,fn),(w,f,err) in sorted([dset 176 for dset in self.sed.data.items() 177 if 'PHOT' in dset[0][0].upper()]): 178 keytags.append(data_labels[dt][0]) 179 data_x.append(w) 180 data_xerr.append(None) 181 data_y.append(f) 182 data_yerr.append(err) 183 line_types.append(data_labels[dt][1]) 184 185 #-- Set line types manually, so the photometry can be plotted with the 186 # right colors. 187 elp = Plotting2.getLineTypes() 188 elp = [lp for lp in elp if lp not in line_types] 189 190 #- Collect model data as well as keytags and set line types 191 model_ids_mcm = [s['LAST_MCMAX_MODEL'] 192 for s in star_grid 193 if s['LAST_MCMAX_MODEL']] 194 195 #- Only if the model_ids list is not empty, MCMax models are available 196 #- Otherwise the ray tracing keyword is unnecessary. 197 if no_models: 198 model_ids_mcm = [] 199 200 wmodels = [] 201 fmodels = [] 202 for model_id,s in zip(model_ids_mcm,star_grid): 203 dpath = os.path.join(cc.path.mout,'models',model_id) 204 fn_spec = 'spectrum{:04.1f}'.format(s['RT_INCLINATION']) 205 if s['RT_NOSOURCE']: fn_spec += 'NOSTAR' 206 fn_spec += '.dat' 207 w,f = MCMax.readModelSpectrum(dpath,s['RT_SPEC'],fn_spec) 208 if s['REDDENING']: 209 print 'Reddening models to correct for interstellar extinction.' 210 ak = self.sed.getAk(s['DISTANCE'],s['REDDENING_MAP'],\ 211 s['REDDENING_LAW']) 212 f = Reddening.redden(w,f,ak,law=s['REDDENING_LAW']) 213 wmodels.append(w) 214 fmodels.append(f) 215 216 for i,(w,f,model_id) in enumerate(zip(wmodels,fmodels,model_ids_mcm)): 217 data_x.append(w) 218 data_y.append(f) 219 data_xerr.append(None) 220 data_yerr.append(None) 221 keytags.append(model_id.replace('_','\_')) 222 line_types.append(elp[i]) 223 224 if self.sed.photbands.size: 225 filts = self.sed.filter_info 226 for i,(w,f) in enumerate(zip(wmodels,fmodels)): 227 mphot = Sed.calcPhotometry(w,f,self.sed.photbands) 228 data_x.append(self.sed.photwave) 229 data_y.append(mphot) 230 if show_phot_filter: 231 data_xerr.append([self.sed.photwave-filts.wlower,\ 232 filts.wupper-self.sed.photwave]) 233 else: 234 data_xerr.append(None) 235 data_yerr.append(None) 236 line_types.append('o%s'%Plotting2.splitLineType(elp[i])[1]) 237 238 keytags = [tag.replace('#','') for tag in keytags] 239 extra_pars = dict() 240 extra_pars['line_types'] = line_types 241 extra_pars['keytags'] = keytags 242 extra_pars['ymax'] = 1.3*max([max(dy[np.isfinite(dy)]) 243 for dy in data_y]) 244 extra_pars['ymin'] = 0.5*min([min(dy[np.isfinite(dy)]) 245 for dy in data_y]) 246 extra_pars['xmin'] = 2 247 extra_pars['xmax'] = 200 248 extra_pars['fontsize_key'] = 16 249 extra_pars['xlogscale'] = 1 250 extra_pars['ylogscale'] = 0 251 extra_pars['xerr_capsize'] = 0 252 extra_pars['xerr_linewidth'] = 4 253 extra_pars.update(cfg_dict) 254 extra_pars.update(kwargs) 255 256 #-- Set filename plot 257 pfn = fn_plt if fn_plt else 'SED' 258 suff = 'iterative_{}'.format(iterative) if iterative else '' 259 pfn = self.setFnPlt(pfn,fn_suffix=suff) 260 261 pfn = Plotting2.plotCols(x=data_x,y=data_y,xerr=data_xerr,\ 262 yerr=data_yerr,filename=pfn,**extra_pars) 263 print '** Your SED plots can be found at:' 264 print pfn 265 print '***********************************'

266 267 268

269 - def plotCorrflux(self,star_grid=[],fn_plt='',cfg='',no_models=0):

270 271 """ 272 Plot correlated fluxes with 0, 1 or more models and data. 273 274 Includes data preparation on the spot. 275 276 @keyword star_grid: list of Star() models to plot. If star_grid is [], 277 only data are plotted. 278 279 (default: []) 280 @type star_grid: list[Star()] 281 @keyword fn_plt: A base plot filename. Includes folder. If not, a 282 default is added 283 284 (default: '') 285 @type fn_plt: string 286 @keyword cfg: path to the Plotting2.plotCols config file. If default, 287 the hard-coded default plotting options are used. 288 289 (default: '') 290 @type cfg: string 291 @keyword no_models: Only show data. 292 293 (default: 0) 294 @type no_models: bool 295 296 """ 297 298 print '***********************************' 299 print '** Creating Correlated Fluxes plot.' 300 if not cc.path.dcflux: 301 print 'No dcflux given in Path.dat. No data are plotted.' 302 303 cfg_dict = Plotting2.readCfg(cfg) 304 if cfg_dict.has_key('no_models'): 305 no_models = cfg_dict['no_models'] 306 if cfg_dict.has_key('filename'): 307 fn_plt = cfg_dict.pop('filename') 308 309 #-- Select MIDI data. Assumes baseline at the end of the filename. 310 ssd = os.path.join(cc.path.dcflux,self.star_name,\ 311 '_'.join([self.star_name,'MIDI','*.fits'])) 312 files = [os.path.splitext(gi)[0] for gi in glob.glob(ssd)] 313 ggd = dict([(float(gi.split('_')[-1].strip('m')),gi+'.fits') 314 for gi in files]) 315 316 #-- Read the models 317 models = [] 318 if not no_models: 319 for s in star_grid: 320 model_id = s['LAST_MCMAX_MODEL'] 321 dpath = os.path.join(cc.path.mout,'models',model_id) 322 fn_vis = 'visibility{:04.1f}'.format(s['RT_INCLINATION']) 323 if s['RT_NOSOURCE']: fn_vis += 'NOSTAR' 324 fn_vis += '.dat' 325 model = MCMax.readVisibilities(dpath=dpath,fn_vis=fn_vis) 326 models.append(model) 327 real_models = [model for model in models if model] 328 if not real_models: 329 no_models = 1 330 baselines = sorted(ggd.keys()) 331 else: 332 baselines = sorted(real_models[0]['baseline'].keys()) 333 334 #-- prepare and collect data, keytags and line types 335 data = [] 336 for bl in baselines: 337 ddict = dict() 338 data.append(ddict) 339 340 #-- Extract data from the fits file 341 if ggd.has_key(bl): 342 dfits = pyfits.open(ggd[bl]) 343 x = 1e6*dfits['OI_WAVELENGTH'].data['EFF_WAVE'][::-1] 344 ddict['x'] = [x] 345 ddict['y'] = [dfits['OI_VIS'].data['VISAMP'][0]][::-1] 346 ddict['yerr'] = [dfits['OI_VIS'].data['VISAMPERR'][0]][::-1] 347 dfits.close() 348 else: 349 ddict['x'] = [] 350 ddict['y'] = [] 351 ddict['yerr'] = [] 352 353 if no_models: 354 continue 355 356 #-- Extract models from the model folders 357 for model in models: 358 ddict['yerr'].append(None) 359 if not model: 360 ddict['x'].append(np.empty(0)) 361 ddict['y'].append(np.empty(0)) 362 continue 363 ddict['x'].append(model['wavelength']) 364 ddict['y'].append(model['flux']*model['baseline'][bl]) 365 366 #-- Set some plot limits 367 ddict['xmin'] = 8 368 ddict['xmax'] = 13 369 ddict['ymin'] = -0.1 370 ddict['labels'] = [('MIDI %.1f m'%bl,0.05,0.9)] 371 #-- Wavelength limits between 8 and 13 micron, limits of the N band 372 # atmospheric transmission. Outside these ranges, the flux is not 373 # relevant 374 ddict['ymax'] = 1.1*max([max(iy[(ix<=13.)*(ix>=8.)]) 375 for ix,iy in zip(ddict['x'],ddict['y']) 376 if iy.size]) 377 kwargs = dict() 378 kwargs['keytags'] = ['MIDI'] if ggd else [] 379 if not no_models: 380 kwargs['keytags'].extend([s['LAST_MCMAX_MODEL'].replace('_','\_') 381 for s in star_grid]) 382 kwargs['xaxis'] = '$\lambda$ ($\mu$m)' 383 kwargs['yaxis'] = 'Corr.~FLux (Jy)' 384 kwargs['dimensions'] = (1,len(data)+1) 385 kwargs['figsize'] = (10,15) 386 kwargs['fontsize_axis'] = 20 387 kwargs['fontsize_ticklabels'] = 20 388 kwargs['fontsize_key'] = 18 389 kwargs['fontsize_label'] = 14 390 kwargs['linewidth'] = 3 391 kwargs['cfg'] = cfg_dict 392 kwargs['extension'] = '.pdf' 393 kwargs['hspace'] = 0.3 394 kwargs['ws_bot'] = 0.01 395 kwargs['ws_top'] = 0.99 396 kwargs['ws_left'] = 0.10 397 kwargs['ws_right'] = 0.98 398 399 #-- Set filename plot 400 pfn = fn_plt if fn_plt else 'CorrFlux' 401 pfn = self.setFnPlt(pfn) 402 403 pfn = Plotting2.plotTiles(data=data,filename=pfn,**kwargs) 404 print '** Your Correlated Flux plots can be found at:' 405 print pfn 406 print '***********************************'

407 408 409

410 - def plotVisibilities(self,star_grid=[],fn_plt='',cfg='',no_models=0):

411 412 """ 413 Plot visibilities as a function of baseline. 414 415 Wavelengths plotted are what is requested in the ray tracing 416 417 Includes data preparation on the spot. 418 419 Data location is that of correlated flux (for the visibilities), but 420 also requires an sed object to retrieve the MIDI spectrum. If one of 421 them is not available, models are be plotted without data. 422 423 @keyword star_grid: list of Star() models to plot. If star_grid is [], 424 only data are plotted. 425 426 (default: []) 427 @type star_grid: list[Star()] 428 @keyword fn_plt: A base plot filename. Includes folder. If not, a 429 default is added 430 431 (default: '') 432 @type fn_plt: string 433 @keyword cfg: path to the Plotting2.plotCols config file. If default, 434 the hard-coded default plotting options are used. 435 436 (default: '') 437 @type cfg: string 438 @keyword no_models: Only show data. 439 440 (default: 0) 441 @type no_models: bool 442 @keyword fn_add_star: Add the star name to the requested plot filename. 443 444 (default: 1) 445 @type fn_add_star: bool 446 447 """ 448 449 if not cc.path.dcflux: 450 print 'No dcflux given in Path.dat. Aborting...' 451 return 452 453 print '***********************************' 454 print '** Creating Visibilities plot.' 455 if not self.sed or 'MIDI' not in self.sed.data_types: 456 print 'No dsed given in Path.dat or no MIDI spectral data found.' 457 458 cfg_dict = Plotting2.readCfg(cfg) 459 if cfg_dict.has_key('no_models'): 460 no_models = cfg_dict['no_models'] 461 if cfg_dict.has_key('filename'): 462 fn_plt = cfg_dict.pop('filename') 463 464 #-- Read the models. Wavelengths are taken from the ray-tracing output 465 models = [] 466 if not no_models: 467 for s in star_grid: 468 model_id = s['LAST_MCMAX_MODEL'] 469 dpath = os.path.join(cc.path.mout,'models',model_id) 470 fn_vis = 'basevis{:04.1f}.dat'.format(s['RT_INCLINATION']) 471 model = MCMax.readVisibilities(dpath=dpath,fn_vis=fn_vis) 472 models.append(model) 473 real_models = [model for model in models if model] 474 if not real_models: 475 no_models = 1 476 else: 477 wavelengths = sorted(real_models[0]['wavelength'].keys()) 478 479 if no_models: wavelengths = (8.,10.,13.) 480 481 #-- Grab the MIDI spectrum 482 if self.sed and 'MIDI' in self.sed.data_types: 483 fn = self.sed.data_filenames[self.sed.data_types.index('MIDI')] 484 midi_flux = self.sed.data[('MIDI',fn)][1] 485 midi_err = self.sed.data[('MIDI',fn)][2] 486 midi_relerr = (midi_err/midi_flux)**2 487 488 #-- Select MIDI data. Assumes baseline at the end of the filename. 489 ssd = os.path.join(cc.path.dcflux,self.star_name,\ 490 '_'.join([self.star_name,'MIDI','*.fits'])) 491 files = [os.path.splitext(gi)[0] for gi in glob.glob(ssd)] 492 ggd = dict([(float(gi.split('_')[-1].strip('m')),gi+'.fits') 493 for gi in files]) 494 else: 495 ggd = dict() 496 497 #-- Collect MIDI data from the fits file and calculate visibilities 498 ddf = dict() 499 for k,v in sorted(ggd.items()): 500 ddf[k] = dict() 501 dfits = pyfits.open(v) 502 503 #-- Read the wavelength 504 cwave = 1e6*dfits['OI_WAVELENGTH'].data['EFF_WAVE'][::-1] 505 506 #-- Read flux + err and select the right range 507 cflux = dfits['OI_VIS'].data['VISAMP'][0][::-1] 508 cflux = cflux[(cwave<=13.)*(cwave>=8.)] 509 cflux_err = dfits['OI_VIS'].data['VISAMPERR'][0][::-1] 510 cflux_err = cflux_err[(cwave<=13.)*(cwave>=8.)] 511 512 #-- The visibilities are correlated flux divided by real flux 513 ddf[k]['y'] = cflux/midi_flux 514 515 #-- Error propagation 516 cflux_relerr = (cflux_err/cflux)**2 517 yerr = np.sqrt(midi_relerr + cflux_relerr)*cflux/midi_flux 518 ddf[k]['yerr'] = yerr 519 520 #-- Wavelength grid 521 ddf[k]['x'] = cwave[(cwave<=13.)*(cwave>=8.)] 522 dfits.close() 523 524 #-- prepare and collect plot data, keytags and line types 525 data = [] 526 for w in wavelengths: 527 ddict = dict() 528 data.append(ddict) 529 530 #-- Set the plot x and y if data are available 531 if ddf: 532 bls = [k for k in sorted(ddf.keys())] 533 ddict['x'] = [[bl for bl in bls]] 534 ddict['y'] = [[ddf[bl]['y'][np.argmin(abs(ddf[bl]['x']-w))] 535 for bl in bls]] 536 ddict['yerr'] = [[ddf[bl]['yerr'][np.argmin(abs(ddf[bl]['x']-w))] 537 for bl in bls]] 538 else: 539 ddict['x'] = [] 540 ddict['y'] = [] 541 ddict['yerr'] = [] 542 543 #-- Set labels 544 ddict['labels'] = [('MIDI %s $\\mu$m'%w,0.85,0.9)] 545 546 if no_models: 547 continue 548 549 #-- Extract models from the model folders 550 for model in models: 551 ddict['yerr'].append(None) 552 if not model: 553 ddict['x'].append(np.empty(0)) 554 ddict['y'].append(np.empty(0)) 555 continue 556 ddict['x'].append(model['baseline']) 557 ddict['y'].append(model['wavelength'][w]) 558 559 kwargs = dict() 560 kwargs['keytags'] = ['MIDI'] if ddf else [] 561 if not no_models: 562 kwargs['keytags'].extend([s['LAST_MCMAX_MODEL'].replace('_','\_') 563 for s in star_grid]) 564 kwargs['xaxis'] = 'Baseline (m)' 565 kwargs['yaxis'] = 'Visibility' 566 kwargs['dimensions'] = (1,len(data)+1) 567 kwargs['figsize'] = (10,15) 568 kwargs['fontsize_axis'] = 20 569 kwargs['fontsize_ticklabels'] = 20 570 kwargs['fontsize_key'] = 18 571 kwargs['fontsize_label'] = 14 572 kwargs['linewidth'] = 3 573 kwargs['cfg'] = cfg_dict 574 kwargs['extension'] = '.pdf' 575 kwargs['hspace'] = 0.3 576 kwargs['ws_bot'] = 0.01 577 kwargs['ws_top'] = 0.99 578 kwargs['ws_left'] = 0.10 579 kwargs['ws_right'] = 0.98 580 581 #-- Set filename plot 582 pfn = fn_plt if fn_plt else 'Visibilities' 583 pfn = self.setFnPlt(pfn) 584 585 pfn = Plotting2.plotTiles(data=data,filename=pfn,**kwargs) 586 print '** Your Correlated Flux plots can be found at:' 587 print pfn 588 print '***********************************'

589 590 591

592 - def plotDensMeixner(self,star_grid=[],fn_plt='',cfg='',**kwargs):

593 594 """ 595 A 2D color-scale plot to visualize the density structure of the Meixner 596 model. Requires the Star()[DENSTYPE] to be Meixner. 597 598 The method used for plotting is available in cc.plotting.PlotMeixner.py 599 and can be used standalone. It is called by this method as well, 600 prepping the star_grid information to be used with PlotMeixner. 601 602 Alternatively, if star_grid is not defined, you can also pass the 603 plotting input for PlotMeixner.plot() through this method. args is only 604 passed to the plotDens method in this case. 605 606 For now the density is plotted in cgs with respect to cm. The density is 607 normalized to the maximum density in the grid. (rho=rho/rho_max) 608 609 @keyword star_grid: parameter sets, if [], the parameter 610 sets are determined from the model ids 611 612 (default: []) 613 @type star_grid: list[Star()] 614 @keyword fn_plt: A base plot filename. Includes folder. If not, a 615 default is added 616 617 (default: '') 618 @type fn_plt: string 619 @keyword cfg: Path to config file for this method. PlotMeixner.plotDens 620 settings can be changed in cfg. Can be used to update the 621 kwargs in case the density grid has to be adapted compared 622 to what the star_grid indicates or what the kwargs 623 contain. Give fn_plt as filename in cfg. 624 625 (default: '') 626 @type cfg: string 627 628 """ 629 630 print '***********************************' 631 print '** Your Meixner density plots can be found at:' 632 cfg_dict = Plotting2.readCfg(cfg) 633 if cfg_dict.has_key('filename'): 634 fn_plt = cfg_dict.pop('filename') 635 kwargs.update(cfg_dict) 636 637 pfns = [] 638 if not star_grid: 639 ppars = dict() 640 ppars['A'] = 100 641 ppars['B'] = 2 642 ppars['C'] = 1 643 ppars['D'] = 0.5 644 ppars['E'] = 0 645 ppars['F'] = 2 646 ppars['rSw'] = 5e13 647 ppars['nTheta'] = 50 648 ppars['nRad'] = 50 649 ppars['rPlot'] = 1e16 650 ppars['rMin'] = 1e13 651 ppars.update(kwargs) 652 653 #-- Define a filename suffix based on combination of 654 # Meixner parameters from ppars. Only update ppars with kwargs 655 # after that. 656 suff = '_'.join(['-'.join([k,str(v)]) 657 for k,v in sorted(ppars.items())]) 658 659 #-- Set the filename 660 fn_plt = fn_plt if fn_plt else 'dens_meixner' 661 ppars['filename'] = self.setFnPlt(fn_plt,fn_suffix=suff) 662 663 #-- Make the plot 664 pfn = PlotMeixner.plotDens(**ppars) 665 pfns.append(pfn) 666 667 alldds = [] 668 star_grid = [s for s in star_grid if s['DENSTYPE'] == 'MEIXNER'] 669 for star in star_grid: 670 #-- Set the Meixner density parameters from the Star() object. 671 ppars = dict() 672 ppars['A'] = star['MEIXA'] 673 ppars['B'] = star['MEIXB'] 674 ppars['C'] = star['MEIXC'] 675 ppars['D'] = star['MEIXD'] 676 ppars['E'] = star['MEIXE'] 677 ppars['F'] = star['MEIXF'] 678 ppars['rSw'] = star['MEIXRSW']*star.au 679 ppars['nTheta'] = star['NTHETA'] 680 ppars['nRad'] = star['NRAD'] 681 ppars['rPlot'] = star['R_OUTER_DUST']*star['R_STAR']*star.Rsun 682 ppars['rMin'] = star['R_INNER_DUST']*star['R_STAR']*star.Rsun 683 684 #-- Define a filename suffix based on combination of 685 # Meixner parameters from ppars. Only update ppars with kwargs 686 # after that. 687 suff = '_'.join(['-'.join([k,str(v)]) 688 for k,v in sorted(ppars.items())]) 689 ppars.update(kwargs) 690 691 #-- Set the filename 692 fn_plt = fn_plt if fn_plt else 'dens_meixner' 693 ppars['filename'] = self.setFnPlt(fn_plt,fn_suffix=suff) 694 695 alldds.append(ppars) 696 697 #-- Only make plots for unique Meixner density models. 698 pfns_noext = [] 699 for ppars in alldds: 700 if ppars['filename'] not in pfns_noext: 701 #-- Make the plot 702 pfn = PlotMeixner.plotDens(**ppars) 703 pfns.append(pfn) 704 pfns_noext.append(os.path.splitext(pfn)[0]) 705 706 print '\n'.join(pfns) 707 print '***********************************'

708 709 710 711 712

713 - def plotDens(self,star_grid=[],models=[],fn_plt='',unit='cm',cfg=''):

714 715 """ 716 Plotting the temperature stratification of the dust. 717 718 All models are shown in one plot. 719 720 @keyword star_grid: parameter sets, if [], the parameter 721 sets are determined from the model ids 722 723 (default: []) 724 @type star_grid: list[Star()] 725 @keyword models: The model_ids, if [], the parameter sets are expected 726 in star_grid 727 728 (default: []) 729 @type models: list[string] 730 @keyword fn_plt: A base plot filename. Includes folder. If not, a 731 default is added 732 733 (default: '') 734 @type fn_plt: string 735 @keyword unit: The unit of the plotted radial grid. Can be 'cm','rstar', 736 'au', 'm' 737 738 (default: 'cm') 739 @type unit: str 740 @keyword cfg: path to the Plotting2.plotCols config file. If default, 741 the hard-coded default plotting options are used. 742 743 (default: '') 744 @type cfg: string 745 746 """ 747 748 print '***********************************' 749 print '** Starting to plot the dust density profile.' 750 if not star_grid and not models: 751 print 'Input is undefined. Aborting.' 752 return 753 elif not star_grid and models: 754 star_grid = self.makeMCMaxStars(models=models) 755 756 cfg_dict = Plotting2.readCfg(cfg) 757 if cfg_dict.has_key('unit'): 758 unit = cfg_dict['unit'] 759 if cfg_dict.has_key('filename'): 760 fn_plt = cfg_dict.pop('filename') 761 762 rads = [s.getDustRad(unit=unit) for s in star_grid] 763 denss = [s.getDustDensity() for s in star_grid] 764 keys = [s['LAST_MCMAX_MODEL'].replace('_','\_') for s in star_grid] 765 766 ppars = dict() 767 ppars['yaxis'] = r'$\rho_\mathrm{d}\ \mathrm{(g cm}^{-3}\mathrm{)}$' 768 if unit == 'rstar': ppars['xaxis'] = '$r\ \mathrm{(R}_\star\mathrm{)}$' 769 elif unit =='au': ppars['xaxis'] = '$r\ \mathrm{(AU)}$' 770 elif unit == 'm': ppars['xaxis'] = '$r\ \mathrm{(m)}$' 771 else: ppars['xaxis'] = '$r\ \mathrm{(cm)}$' 772 ppars['xlogscale'] = 1 773 ppars['ylogscale'] = 1 774 775 #-- Set plot filename 776 pfn = fn_plt if fn_plt else 'dens' 777 pfn = self.setFnPlt(pfn) 778 779 pfn = Plotting2.plotCols(x=rads,y=denss,filename=pfn,\ 780 keytags=keys,cfg=cfg_dict,**ppars) 781 print '** Your plot can be found at:' 782 print pfn 783 print '***********************************'

784 785 786

787 - def plotTemp(self,star_grid=[],models=[],fn_plt='',power=[],unit='cm',\ 788 cfg=''):

789 790 """ 791 Plotting the temperature stratification of the dust. 792 793 All models are shown in one plot. 794 795 @keyword star_grid: parameter sets, if [], the parameter 796 sets are determined from the model ids 797 798 (default: []) 799 @type star_grid: list[Star()] 800 @keyword models: The model_ids, if [], the parameter sets are expected 801 in star_grid 802 803 (default: []) 804 @type models: list[string] 805 @keyword fn_plt: A base plot filename. Includes folder. If not, a 806 default is added 807 808 (default: '') 809 @type fn_plt: string 810 @keyword power: A list of values for s in below formula. If [] no power 811 law is included. Power law parameters are taken from 812 star_grid[0]. 813 814 See Thesis p32, where power is s in 815 T(r) = T_eff*(2*r/R_STAR)**(-2/(4+s)). This value is 816 typically 1. 817 818 (default: []) 819 @type power: list 820 @keyword unit: The unit of the plotted radial grid. Can be 'cm','rstar', 821 'au', 'm' 822 823 (default: 'cm') 824 @type unit: str 825 @keyword cfg: path to the Plotting2.plotCols config file. If default, 826 the hard-coded default plotting options are used. 827 828 (default: '') 829 @type cfg: string 830 831 """ 832 833 print '***********************************' 834 print '** Starting to plot dust temperature stratification.' 835 if not star_grid and not models: 836 print 'Input is undefined. Aborting.' 837 return 838 elif not star_grid and models: 839 star_grid = self.makeMCMaxStars(models=models) 840 cfg_dict = Plotting2.readCfg(cfg) 841 if cfg_dict.has_key('power'): 842 power = cfg_dict['power'] 843 if cfg_dict.has_key('unit'): 844 unit = cfg_dict['unit'] 845 if cfg_dict.has_key('filename'): 846 fn_plt = cfg_dict.pop('filename') 847 848 rads = [] 849 temps = [] 850 keytags = [] 851 for star in star_grid: 852 rad = star.getDustRad(unit=unit) 853 temp,key = star.getDustTemperature(add_key=1) 854 rads.append(rad) 855 temps.append(temp) 856 keytags.append(key) 857 858 #-- Add power laws if requested 859 for s in power: 860 rad = star_grid[0].getDustRad(unit=unit) 861 tstar = star_grid[0]['T_STAR'] 862 rstar = star_grid[0]['R_STAR'] 863 key = '$T_\mathrm{d}(r) = %i \\left(\\frac{2r}'%int(tstar) + \ 864 '{\mathrm{R}_\star}\\right)^{\\frac{2}{4+%i}}$'%int(s) 865 temp = Temperature.Tdust(r=rad,T0=tstar,r0=rstar,s=s) 866 rads.append(rad) 867 temps.append(temp) 868 keytags.append(key) 869 870 ppars = dict() 871 ppars['yaxis'] = '$T_\mathrm{d}$ (K)' 872 if unit == 'rstar': ppars['xaxis'] = '$r\ \mathrm{(R}_\star\mathrm{)}$' 873 elif unit =='au': ppars['xaxis'] = '$r\ \mathrm{(AU)}$' 874 elif unit == 'm': ppars['xaxis'] = '$r\ \mathrm{(m)}$' 875 else: ppars['xaxis'] = '$r\ \mathrm{(cm)}$' 876 ppars['plot_title'] = 'Average Dust Temperature Stratification for %s'\ 877 %(self.star_name_plots) 878 879 #-- Set plot filename 880 pfn = fn_plt if fn_plt else 'Td_avg' 881 pfn = self.setFnPlt(pfn) 882 883 pfn = Plotting2.plotCols(x=rads,y=temps,filename=pfn,\ 884 key_location=(0.05,0.05),cfg=cfg_dict,\ 885 xlogscale=1,ylogscale=1,\ 886 keytags=keytags,**ppars) 887 print '** Your plots can be found at:' 888 print pfn 889 print '***********************************'

890 891 892

893 - def plotTempSpecies(self,star_grid=[],models=[],fn_plt='',include_total=1,\ 894 unit='cm',power=[],cfg=''):

895 896 """ 897 Plotting the temperature stratification of the dust for the species 898 separately, per model. 899 900 @keyword star_grid: parameter sets, if [], the parameter 901 sets are determined from the model ids 902 903 (default: []) 904 @type star_grid: list[Star()] 905 @keyword models: The model_ids, if [], the parameter sets are expected 906 in star_grid 907 908 (default: []) 909 @type models: list[string] 910 @keyword fn_plt: A base plot filename. Includes folder. If not, a 911 default is added 912 913 (default: '') 914 @type fn_plt: string 915 @keyword include_total: Include the sum of all temperature profiles as 916 well for comparison. 917 918 (default: 0) 919 @type include_total: bool 920 @keyword unit: The unit of the plotted radial grid. Can be 'cm','rstar', 921 'au', 'm' 922 923 (default: 'cm') 924 @type unit: str 925 @keyword power: A list of values for s in below formula. If [] no power 926 law is included. Power law parameters are taken from 927 star_grid[0]. 928 929 See Thesis p32, where power is s in 930 T(r) = T_eff*(2*r/R_STAR)**(-2/(4+s)). This value is 931 typically 1. 932 933 (default: []) 934 @type power: list 935 @keyword cfg: path to the Plotting2.plotCols config file. If default, 936 the hard-coded default plotting options are used. 937 938 (default: '') 939 @type cfg: string 940 941 """ 942 943 print '***********************************' 944 print '** Starting to plot dust temperature for separate species.' 945 if not star_grid and not models: 946 print 'Input is undefined. Aborting.' 947 return 948 elif not star_grid and models: 949 star_grid = self.makeMCMaxStars(models=models,id_type='MCMax') 950 raise IOError('Reading dust species temperatures from a model id'+\ 951 ' list only, not yet implemented.') 952 #- Requires star.dust_list and T_CONTACT to be taken from the log file. 953 #- It's possible, but needs some programming 954 955 cfg_dict = Plotting2.readCfg(cfg) 956 if cfg_dict.has_key('power'): 957 power = cfg_dict['power'] 958 if cfg_dict.has_key('unit'): 959 unit = cfg_dict['unit'] 960 if cfg_dict.has_key('filename'): 961 fn_plt = cfg_dict.pop('filename') 962 963 plot_filenames = [] 964 for star in star_grid: 965 if not int(star['T_CONTACT']): 966 rads = [star.getDustRad(species=species,unit=unit) 967 for species in star.getDustList()] 968 temps = [star.getDustTemperature(species=species) 969 for species in star.getDustList()] 970 rads = [r[t<=star['T_DES_%s'%sp]] 971 for r,t,sp in zip(rads,temps,star.getDustList())] 972 temps = [t[t<=star['T_DES_%s'%sp]] 973 for t,sp in zip(temps,star.getDustList())] 974 keytags = list(star.getDustList()) 975 else: 976 include_total = 1 977 print 'Thermal contact is on. All dust species share the ' + \ 978 'same temperature profile. Vertical lines indicate ' + \ 979 'inner radii of dust species.' 980 rads, temps, keytags = [], [], [] 981 982 if include_total: 983 rad = star.getDustRad(unit=unit) 984 temp, key = star.getDustTemperature(add_key=1) 985 rads.append(rad[rad>star['R_INNER_DUST']\ 986 *star.Rsun*star['R_STAR']]) 987 temps.append(temp[rad>star['R_INNER_DUST']\ 988 *star.Rsun*star['R_STAR']]) 989 keytags.append(key) 990 991 #-- Add power laws if requested 992 for s in power: 993 rad = star_grid[0].getDustRad(unit=unit) 994 tstar = star_grid[0]['T_STAR'] 995 rstar = star_grid[0]['R_STAR'] 996 key = '$T_\mathrm{d}(r) = %i \\left(\\frac{2r}'%int(tstar) + \ 997 '{\mathrm{R}_\star}\\right)^{\\frac{2}{4+%i}}$'%int(s) 998 temp = Temperature.Tdust(r=rad,T0=tstar,r0=rstar,s=s) 999 rads.append(rad) 1000 temps.append(temp) 1001 keytags.append(key) 1002 1003 ppars = dict() 1004 ppars['yaxis'] = '$T_\mathrm{d}$ (K)' 1005 if unit == 'rstar': ppars['xaxis'] = '$r\ \mathrm{(R}_\star\mathrm{)}$' 1006 elif unit =='au': ppars['xaxis'] = '$r\ \mathrm{(AU)}$' 1007 elif unit == 'm': ppars['xaxis'] = '$r\ \mathrm{(m)}$' 1008 else: ppars['xaxis'] = '$r\ \mathrm{(cm)}$' 1009 1010 #-- Set plot filename 1011 pfn = fn_plt if fn_plt else 'Td_species' 1012 suff = star['LAST_MCMAX_MODEL'] 1013 pfn = self.setFnPlt(pfn,fn_suffix=suff) 1014 1015 plot_filenames.append(Plotting2.plotCols(x=rads,y=temps,\ 1016 cfg=cfg_dict,filename=pfn,\ 1017 keytags=keytags,xlogscale=1,ylogscale=1,\ 1018 fontsize_key=16,**ppars)) 1019 1020 if len(plot_filenames) != len(star_grid): 1021 print 'At least one of the models does not yet have a MCMax model.' 1022 if plot_filenames[0][-4:] == '.pdf': 1023 pfn = fn_plt if fn_plt else 'Td_species' 1024 pfn = self.setFnPlt(pfn) + '.pdf' 1025 DataIO.joinPdf(old=plot_filenames,new=pfn) 1026 print '** Your plots can be found at:' 1027 print pfn 1028 print '***********************************' 1029 else: 1030 print '** Plots can be found at:' 1031 print '\n'.join(plot_filenames) 1032 print '***********************************'

1033 1034 1035

1036 - def plotOpacities(self,star_grid=[],fn_plt='',scaling=0,species=['AMC'],\ 1037 cfg='',index=0,*args,**kwargs):

1038 1039 """ 1040 Plotting wavelength dependent mass extinction coefficients 1041 (ie opacities). 1042 1043 If based on star_grid or modelslist, they are scaled with abundances if 1044 wanted. 1045 1046 If no model info is given, the input opacities are plotted. 1047 1048 Args and kwargs can be given straight to the plot command. 1049 1050 @keyword star_grid: The input Star() models. If default, the MCMax 1051 input opacities are plotted. 1052 1053 (default: []) 1054 @type star_grid: list(Star()) 1055 @keyword fn_plt: A base plot filename. Includes folder. If not, a 1056 default is added 1057 1058 (default: '') 1059 @type fn_plt: string 1060 @keyword scaling: allow species abundance scaling of opacities 1061 1062 (default: 0) 1063 @type scaling: bool 1064 @keyword species: If no star_grid or model list are given, this gives 1065 the species requested to be plotted from Dust.dat 1066 1067 (default: ['AMC']) 1068 @type species: list(string) 1069 @keyword cfg: path to the Plotting2.plotCols config file. If default, 1070 the hard-coded default plotting options are used. 1071 1072 (default: '') 1073 @type cfg: string 1074 @keyword index: The index of the kappas in the .opacity/.particle file. 1075 0: extinction, 1: absorption, 2: scattering. Only 1076 relevant when plotting without Star() models. 1077 1078 (default: 0) 1079 @type index: int 1080 1081 """ 1082 1083 print '***********************************' 1084 print '** Starting to plot dust opacities.' 1085 1086 #-- Set the filename 1087 cfg_dict = Plotting2.readCfg(cfg) 1088 if cfg_dict.has_key('filename'): 1089 fn_plt = cfg_dict.pop('filename') 1090 if not star_grid and not fn_plt: 1091 fn_plt = os.path.join(cc.path.mopac,\ 1092 'dust_opacities_%s'%'_'.join(species)) 1093 1094 #-- Retrieve index keyword if in cfg. 1095 if cfg_dict.has_key('index'): 1096 index = cfg_dict.pop('index') 1097 1098 #-- Set some plot parameters 1099 ppars = dict() 1100 ppars['xaxis'] = '$\lambda$ ($\mu \mathrm{m}$)' 1101 ppars['yaxis'] = '$\kappa_\lambda$ ($\mathrm{cm}^2\mathrm{/g}$)' 1102 ppars['xlogscale'] = 1 1103 ppars['ylogscale'] = 1 1104 ppars['key_location'] = (0.05,0.05) 1105 ppars.update(kwargs) 1106 ppars.update(cfg_dict) 1107 1108 #-- Check if raw opacities or modeling results are requested 1109 if not star_grid: 1110 kr = KappaReader.KappaReader() 1111 wl_list = [kr.getWavelength(sp) for sp in species] 1112 q_list = [kr.getKappas(sp,index) for sp in species] 1113 if not ppars.has_key('keytags'): 1114 ppars['keytags'] = species 1115 1116 #-- Set plot filename 1117 pfn = fn_plt if fn_plt else 'opacities_species' 1118 pfn = self.setFnPlt(pfn) 1119 pfn = Plotting2.plotCols(x=wl_list,y=q_list,filename=pfn,\ 1120 *args,**ppars) 1121 print '** Your plot can be found at:' 1122 print pfn 1123 else: 1124 fns = [] 1125 for star in star_grid: 1126 try: 1127 wave,opacities = star.readKappas() 1128 except IOError: 1129 continue 1130 if star['INCLUDE_SCAT_GAS']: 1131 n_species = len(star.getDustList()) 1132 opacities = [opacities[i]+opacities[i+n_species] 1133 for i in range(len(star.getDustList()))] 1134 else: 1135 opacities = [opacities[i] 1136 for i in range(len(star.getDustList()))] 1137 1138 if scaling: 1139 opacities = [op*star['A_%s'%sp] 1140 for op,sp in zip(opacities,star.getDustList())] 1141 1142 #-- Set plot filename 1143 pfn = fn_plt if fn_plt else 'opacities_species' 1144 suff = star['LAST_MCMAX_MODEL'] 1145 pfn = self.setFnPlt(pfn,fn_suffix=suff) 1146 1147 keys = ['%s with $A$ = %s and $T_{des} = %i$ K'\ 1148 %(sp,str(star['A_%s'%sp]),int(star['T_DES_%s'%sp])) 1149 for sp in star.getDustList()] 1150 fns.append(Plotting2.plotCols(x=wave,y=opacities,keytags=keys,\ 1151 filename=pfn,*args,**ppars)) 1152 if len(fns) != len(star_grid): 1153 print 'At least one of the models requested does not yet ' + \ 1154 'have a MCMax model.' 1155 print '** Your plots can be found at:' 1156 if fns and fns[-1][-4] == '.pdf': 1157 pfn = fn_plt if fn_plt else 'opacities_species' 1158 pfn = self.setFnPlt(pfn) + '.pdf' 1159 DataIO.joinPdf(old=fns,new=pfn) 1160 print pfn 1161 else: 1162 print '\n'.join(fns) 1163 print '***********************************'

1164 1165 1166

1167 - def plotExtinction(self,star_grid=[],models=[],fn_plt='',plot_default=1,\ 1168 cfg=''):

1169 1170 """ 1171 Plotting wavelength dependent extinction efficiencies wrt grain size. 1172 1173 This always depends on a star_grid or one created from a list of MCMax 1174 model ids. 1175 1176 Plotted are the total efficiencies, including relative weights between 1177 the included dust species. This is the input for GASTRoNOoM! 1178 1179 @keyword star_grid: List of Star() instances. If default, model ids 1180 have to be given. 1181 1182 (default: []) 1183 @type star_grid: list[Star()] 1184 @keyword models: The model ids, only required if star_grid is [] 1185 1186 (default: []) 1187 @type models: list[string] 1188 @keyword fn_plt: A base plot filename. Includes folder. If not, a 1189 default is added 1190 1191 (default: '') 1192 @type fn_plt: string 1193 @keyword cfg: path to the Plotting2.plotCols config file. If default, 1194 the hard-coded default plotting options are used. 1195 1196 (default: '') 1197 @type cfg: string 1198 1199 """ 1200 1201 print '***********************************' 1202 print '** Plotting Q_ext/a.' 1203 if not star_grid and not models: 1204 print 'Input is undefined. Aborting.' 1205 return 1206 elif not star_grid and models: 1207 star_grid = self.makeMCMaxStars(models=models) 1208 1209 cfg_dict = Plotting2.readCfg(cfg) 1210 if cfg_dict.has_key('filename'): 1211 fn_plt = cfg_dict.pop('filename') 1212 1213 x = [] 1214 y = [] 1215 keys = [] 1216 for star in star_grid: 1217 try: 1218 inputfile = os.path.join(cc.path.gdata,star['TEMDUST_FILENAME']) 1219 opacities = DataIO.readCols(filename=inputfile) 1220 x.append(opacities[0]) 1221 y.append(opacities[1]) 1222 keys.append('$Q_\mathrm{ext}/a$ for MCMax %s'\ 1223 %star['LAST_MCMAX_MODEL'].replace('_','\_')) 1224 except IOError: 1225 pass 1226 1227 #-- Set plot filename 1228 pfn = fn_plt if fn_plt else 'gastronoom_opacities' 1229 pfn = self.setFnPlt(pfn) 1230 1231 title = 'GASTRoNOoM Extinction Efficiencies in %s'\ 1232 %(self.star_name_plots) 1233 pfn = Plotting2.plotCols(x=x,y=y,cfg=cfg_dict,filename=pfn,\ 1234 xaxis='$\lambda$ ($\mu$m)',keytags=keys,\ 1235 yaxis='$Q_{ext}/a$ (cm$^{-1}$)',\ 1236 plot_title=title,key_location=(0.7,0.6),\ 1237 xlogscale=1,ylogscale=1,fontsize_key=20) 1238 print '** The extinction efficiency plot can be found at:' 1239 print pfn 1240 print '***********************************'

1241 1242 1243

1244 - def makeMCMaxStars(self,models):

1245 1246 ''' 1247 Set parameters for star_list taken from the MCMax database. 1248 1249 Based on the model id of MCMax. 1250 1251 @param models: model_ids for the MCMax db 1252 @type models: list(string) 1253 @return: The model instances 1254 @rtype: list(Star()) 1255 1256 ''' 1257 1258 star_grid = Star.makeStars(models=models,code='MCMax',id_type='MCMax',\ 1259 path=self.path) 1260 for star,model in zip(star_grid,models): 1261 filepath = os.path.join(cc.path.mout,'models',\ 1262 star['LAST_MCMAX_MODEL']) 1263 denstemp = os.path.join(filepath,'denstemp.dat') 1264 logfile = os.path.join(filepath,'log.dat') 1265 grid_shape = DataIO.getKeyData(filename=denstemp,incr=1,\ 1266 keyword='NGRAINS',single=0)[0] 1267 star.update({'NTHETA':int(grid_shape[1]),\ 1268 'NRAD':int(grid_shape[0]),\ 1269 'T_STAR':float(DataIO.getKeyData(filename=logfile,\ 1270 incr=0,\ 1271 keyword='STELLAR TEMPERATURE',\ 1272 single=0)[0][2]),\ 1273 'R_STAR':float(DataIO.getKeyData(filename=logfile,\ 1274 incr=0,\ 1275 keyword='STELLAR RADIUS',\ 1276 single=0)[0][2])}) 1277 return star_grid

Source Code for Module ComboCode.cc.plotting.objects.PlotDust