ComboCode.cc.data.LPTools

1 # -*- coding: utf-8 -*- 2 3 """ 4 Some tools for measuring line profiles of emission lines. 5 6 Author: R. Lombaert 7 8 """ 9 10 import os 11 import numpy as np 12 from scipy import mean,sqrt,log, std,median 13 from scipy import argmin,argmax,array 14 from scipy.integrate import trapz 15 from matplotlib import pyplot as plt 16 from astropy import units as u 17 18 import cc.path 19 from cc.tools.readers import FitsReader, TxtReader 20 from cc.tools.io import DataIO 21 from cc.plotting import Plotting2 22 from cc.tools.units import Equivalency as eq 23 24 from cc.ivs.sigproc import fit, funclib 25 26 27

28 -def readTelescopeProperties(telescope):

29 30 """ 31 Read the telescope properties from Telescope.dat. 32 33 This currently includes the telescope size in m, and the default 34 absolute flux calibration uncertainty. 35 36 @param telescope: The telescope requested 37 @type telescope: str 38 39 @return: The telescope size and absolute flux calibration uncertainty 40 @rtype: (float,float) 41 42 """ 43 44 all_telescopes = DataIO.getInputData(keyword='TELESCOPE',start_index=5,\ 45 filename='Telescope.dat') 46 if 'PACS' in telescope: 47 telescope = 'PACS' 48 else: 49 telescope = telescope 50 try: 51 tel_index = all_telescopes.index(telescope) 52 except ValueError: 53 raise ValueError('%s not found in Telescope.dat.'%telescope) 54 55 size = DataIO.getInputData(keyword='SIZE',start_index=5,\ 56 filename='Telescope.dat',\ 57 rindex=tel_index) 58 abs_err = DataIO.getInputData(keyword='ABS_ERR',start_index=5,\ 59 filename='Telescope.dat',\ 60 rindex=tel_index) 61 return (size,abs_err)

62 63 64

65 -def readLineProfile(filename):

66 67 ''' 68 Read a line profile, independent of the type of extension. 69 70 @keyword filename: The filename to the data file of the line profile. If 71 None a line profile object is expected. 72 73 (default: None) 74 @type filename: string 75 76 @return: The line profile 77 @rtype: LPDataReader() 78 79 ''' 80 81 if filename[-5:] == '.fits': 82 lprof = FitsReader.FitsReader(fn=filename) 83 else: 84 lprof = TxtReader.TxtReader(fn=filename) 85 return lprof

86 87 88

89 -def integrateLPData(vexp,filename=None,lprof=None,window=1.2):

90 91 """ 92 Integrate a line profile read from a fits file or a txt file. 93 94 Requires a terminal expansion velocity to determine the window in which the 95 integration is done. 96 97 If a fits file is read, the source velocity is taken from it (if available) 98 If a text file is read, a source velocity needs to be given as input. 99 100 @param vexp: The terminal gas velocity 101 @type vexp: float 102 103 @keyword filename: The filename to the data file of the line profile. If 104 None a line profile object is expected. 105 106 (default: None) 107 @type filename: string 108 @keyword lprof: A line profile object (LPDataReader or inheriting classes) 109 If None, a filename is expected! 110 111 (default: None) 112 @type lprof: LPDataReader() 113 @keyword window: The factor with which vexp is multiplied when selecting 114 the integration window in the velocity grid. 115 116 (default: 1.2) 117 @type window: float 118 119 @return: The integrated intensity of the profile 120 @rtype: float 121 122 """ 123 124 if lprof is None: 125 lprof = readLineProfile(filename) 126 vel = lprof.getVelocity() 127 flux = lprof.getFlux() 128 vlsr = lprof.getVlsr() 129 vexp = float(vexp) 130 window = float(window) 131 132 vel_sel = vel[(vel > vlsr - window*vexp)*(vel < vlsr + window*vexp)] 133 flux_sel = flux[(vel > vlsr - window*vexp)*(vel < vlsr + window*vexp)] 134 return trapz(x=vel_sel,y=flux_sel)

135 136 137

138 -def getPeakLPData(filename=None,lprof=None,vlsr=None,method='mean',npoints=5):

139 140 """ 141 Calculate the peak value of a line profile read from a fits file or a txt 142 file. 143 144 If a fits file is read, the source velocity is taken from it (if available) 145 If a text file is read, a source velocity needs to be given as input. 146 147 The peak value is defined as the mean of central 5 flux points around the 148 source velocity. 149 150 @keyword filename: The filename to the data file of the line profile. If 151 None a line profile object is expected. 152 153 (default: None) 154 @type filename: string 155 @keyword lprof: A line profile object (LPDataReader or inheriting classes) 156 If None, a filename is expected! 157 158 (default: None) 159 @type lprof: LPDataReader() 160 @keyword vlsr: If you want to provide your own vlsr for some reason. Leave 161 as default if you don't have a good reason to do that. 162 163 (default: None) 164 @type vlsr: float 165 @keyword method: The method applied: either 'mean' or 'fit'. Mean derives 166 the peak value from the mean of the npoints flux points 167 around the vlsr. Fit takes the central peak flux at from 168 the fit. 169 170 (default: 'mean') 171 @type method: str 172 @keyword npoints: The number of points around vlsr used for deriving the 173 peak value via the mean method. 174 175 (default: 5) 176 @type npoints: int 177 178 @return: The peak intensity of the profile 179 @rtype: float 180 181 """ 182 183 method = method.lower() 184 if method not in ['mean','fit']: 185 method = 'mean' 186 187 #-- Read line profile 188 if lprof is None: 189 lprof = readLineProfile(filename) 190 191 #-- If fit method, fit the line profile, and estimate peak value at vlsr of 192 # the fit 193 if method == 'fit': 194 fit = fitLP(lprof=lprof) 195 return fit['peak'] 196 197 #-- If mean method, take the mean around vlsr given n_points. 198 vel = lprof.getVelocity() 199 flux = lprof.getFlux() 200 if vlsr is None: 201 vlsr = lprof.getVlsr() 202 else: 203 vlsr = float(vlsr) 204 205 #-- Set number of points on each side of vlsr 206 if npoints%2.==0: 207 npoints = int(npoints/2.) 208 else: 209 npoints = int(npoints/2.) 210 211 #-- Determine mid point, return mean 212 imid = argmin(np.abs(vel-vlsr)) 213 return mean(flux[imid-npoints:imid+npoints+1])

214 215 216

217 -def getLPDataFWHM(lprof):

218 219 ''' 220 Calculate the FWHM of the line profile. 221 222 @param lprof: the line profile object. 223 @type lprof: LPDataReader() 224 225 @return: the fwhm of the line 226 @rtype: float 227 228 ''' 229 230 flux = lprof.getFlux() 231 vel = lprof.getVelocity() 232 vlsr = lprof.getVlsr() 233 maxval = max(flux) 234 i_mid = argmax(flux) 235 flux1 = flux[:i_mid] 236 flux2 = flux[i_mid:] 237 v1 = vel[argmin(abs(flux1-maxval/2.))] 238 v2 = vel[argmin(abs(flux2-maxval/2.))+i_mid] 239 return v2-v1

240 241 242

243 -def varyInitialFit(vel,flux,initial,index,values,vary,\ 244 function=funclib.soft_parabola,vary_window=0):

245 246 """ 247 Fit a function to a line profile for different initial guesses of a single 248 parameter. 249 250 The best fit function is returned. Best fit is determined by checking for 251 the smallest relative error for the parameter in question. 252 253 @param vel: The velocity grid 254 @type vel: array 255 @param flux: The intensity grid 256 @type flux: array 257 @param initial: initial parameters (eg [int,vlsr,vexp,gamma] for sp) 258 @type initial: list 259 @param index: Index of initial parameter to be varied 260 @type index: int 261 @param values: The values used for the variable initial parameter 262 @type values: tuple 263 @param vary: Allow initial parameter to be changed in fitting process. 264 Must have same length as initial. 265 @type vary: list[bool] 266 267 @keyword function: The function to be fitted 268 269 (default: funclib.soft_parabola) 270 @type function: funclib.function 271 @keyword vary_window: Vary the velocity window based on the third initial 272 parameter. 1.5 for soft parabola, 3 for other functions 273 centered on the second initial parameter. 274 (mu/vlsr and sigma/vexp respectively) 275 276 (default: 0) 277 @type vary_window: bool 278 279 @return: The model after minimization 280 @rtype: funclib.soft_parabola 281 282 """ 283 284 values, initial, vary_window = list(values), list(initial), int(vary_window) 285 all_init = [[p]*len(values) for i,p in enumerate(initial) if i != index] 286 all_init.insert(index,values) 287 if vary_window: 288 #-- The window for soft parabola can be quite narrow. For Gaussians 289 # it should be wider, taking into account broader wings, and that 290 # sigma/2 < vexp = fwhm/2 291 window = function == funclib.soft_parabola and 1.5 or 3. 292 results = [fitFunction(vel[np.abs(vel-initi[1])<=(initi[2]*window)],\ 293 flux[np.abs(vel-initi[1])<=(initi[2]*window)],\ 294 initi,function,vary=vary) 295 for initi in zip(*all_init)] 296 else: 297 results = [fitFunction(vel,flux,initi,function,vary=vary) 298 for initi in zip(*all_init)] 299 rel_errors = [fg.get_parameters()[1][index]/fg.get_parameters()[0][index] 300 for fg in results] 301 sel_results = [res 302 for res,rel_err in zip(results,rel_errors) 303 if rel_err > 10**(-5)] 304 sel_errors = array([rel_err 305 for rel_err in rel_errors 306 if rel_err > 10**(-5)]) 307 if not sel_results: 308 sel_results = results 309 sel_errors = array(rel_errors) 310 bestresult = sel_results[argmin(abs(sel_errors))] 311 return bestresult 312 313 314

315 -def fitFunction(x,y,initial,function,vary):

316 317 """ 318 Fit a function to a set of x and y values. 319 320 @param x: The x grid 321 @type x: array 322 @param y: The y grid 323 @type y: array 324 @param initial: initial parameters 325 @type initial: list 326 @param function: The function to be fitted 327 @type function: funclib.function (e.g. funclib.soft_parabola,funclib.gauss) 328 @param vary: Allow initial parameter to be changed in fitting process. 329 Must have same length as initial. 330 @type vary: list[bool] 331 332 @return: The model after minimization 333 @rtype: funclib.function() (some function) 334 """ 335 336 #-- fit only soft parabola 337 # 1. setup model 338 #mymodel = funclib.soft_parabola() 339 if function == funclib.gauss and False in vary: 340 mymodel = function(use_jacobian=False) 341 else: 342 mymodel = function() 343 # 2. Initial values: e.g. for SP [int,vlsr,vexp,gamma] 344 mymodel.setup_parameters(values=initial,vary=vary) 345 # 3. minimize and evaluate fit 346 result = fit.minimize(x,y,mymodel) 347 return mymodel

348 349 350

351 -def checkLPShape(vel,flux,vlsr,vexp,window=2.,show=0):

352 353 """ 354 Check the shape of the line profile, to see if any irregular patterns are 355 present. 356 357 Based on the diff() method, and the std on the result of it. If sharp 358 patterns are present, they should be picked up by this method, and an extra 359 component can be included in fitting the line. 360 361 Detects absorption irregularities, not emission! At least, it tries to 362 filter the emission effects out. If an emission effect is stronger than 363 an absorption effect, it should also not detect any irregularities. 364 365 Currently checks whether a strong downward trend in the flux is associated 366 with a strong upward trend, within a window of 6 km/s. Broader "absorption 367 detections" are unlikely, since ISM absorption lines should not have a line 368 width that is broader than, e.g., 3 km/s (i.e. turbulent velocity). 369 370 Still being tested! 371 372 @param vel: The velocity grid 373 @type vel: array 374 @param flux: The flux grid 375 @type flux: array 376 @param vlsr: the central source velocity 377 @type vlsr: float 378 @param vexp: The expected gas terminal velocity 379 @type vexp: float 380 381 @keyword window: The window for line selection. Value is multiplied with 382 0.6. For the usual window of 2., this is 1.2 (SP). For the 383 Gaussian window of 3., this is 1.8. 384 385 (default: 2.) 386 @type window: float 387 @keyword show: Show the results of the diff and std methods. 388 389 (default: 0) 390 @type show: bool 391 392 @return: The details of the absorption irregularity. Assuming a Gaussian, 393 a depth, mid point velocity, a width and a continuum value are 394 returned. If no irregularity is found, None is returned 395 @rtype: tuple(float) 396 397 """ 398 399 #-- Grab the velocity resolution 400 velres = np.median(np.diff(vel)) 401 402 #-- Get diff of flux to have a measure of change in slope throughout lprof 403 fdf = np.diff(flux) 404 velfdf = vel[1:] 405 fluxfdf = flux[1:] 406 #-- Select emission line. Use 0.8 as factor instead of 0.6 (results in 1.4 407 # for window == 2), to take into account the possible underestimation of 408 # vexp 409 #-- Difference on window factor stems from fact that when window==2, an 0.5 410 # factor really cuts the window too close to the line. For all others, 0.5 411 # is fine. 412 if window == 2: window_factor = 0.8 413 else: window_factor = 0.5 414 fdfwhereline = np.abs(velfdf-vlsr)<=window*window_factor*vexp 415 #-- create window for noise calculation: 416 # - close to the emission line 417 # - contains enough points to be statistically significant 418 # Note that i is determined based on the window: you want enough points, 419 # but not too many, and not less than the window's width. 420 i = window 421 fdfwindow = np.abs(velfdf-vlsr)<=i*vexp 422 while len(velfdf[-fdfwhereline*fdfwindow]) < 60 and i != 6: 423 i += 1 424 fdfwindow = np.abs(velfdf-vlsr)<=i*vexp 425 #-- Check if flux outliers (outside the line) are present: 426 #-- ignore for diff noise calculations 427 noiseflux = std(fluxfdf[-fdfwhereline*fdfwindow]) 428 fdfcleanflux = fluxfdf<3.*noiseflux 429 #-- Calc noise of diff in region where 'not emission line', but also not 430 # too far out, and where no flux outliers are present 431 noisedf = std(fdf[-fdfwhereline*fdfwindow*fdfcleanflux]) 432 433 #-- Show the result 434 if show: 435 plt.clf() 436 velfdfplot = velfdf[fdfwindow] 437 fluxfdfplot = fluxfdf[fdfwindow] 438 plt.step(velfdfplot,fluxfdfplot,'-r',lw=3,where='mid',\ 439 label='Observed profile') 440 cleanflux = fluxfdf.copy() 441 cleanflux[-(-fdfwhereline*fdfwindow*fdfcleanflux)] = np.nan 442 plt.step(velfdf,cleanflux,'-k',lw=2,where='mid',label='Clean noise') 443 plt.plot(velfdf[fdfwindow],fdf[fdfwindow],'b-',lw=3,\ 444 label='diff(profile)') 445 plt.plot([velfdfplot[0],velfdfplot[-1]],[noisedf*3,noisedf*3],'--k',\ 446 label='Upper df limit') 447 plt.plot([velfdfplot[0],velfdfplot[-1]],[-noisedf*3,-noisedf*3],'--k',\ 448 label='Lower df limit') 449 ax = plt.gca() 450 ax.axvline(x=vlsr-window*window_factor*vexp,linewidth=2, ls='--',\ 451 color='k') 452 ax.axvline(x=vlsr+window*window_factor*vexp,linewidth=2, ls='--',\ 453 color='k') 454 leg = plt.legend(loc='best',fancybox=True) 455 leg.get_frame().set_alpha(0.5) 456 plt.show() 457 458 #-- See where the change in profile is significantly large, but only do 459 # this where the line is located 460 fdf_large = np.abs(fdf[fdfwhereline]) >= 3.*noisedf 461 #-- If any large values are detected, check the min and max values 462 if fdf_large.any(): 463 imin = argmin(fdf[fdfwhereline]) 464 imax = argmax(fdf[fdfwhereline]) 465 #-- Check if both the min and max are actually larger than 3*noisedf 466 if not (fdf_large[imin] and fdf_large[imax]): 467 return None 468 #-- If imin > imax, it's likely due to a very sharp profile: You do not 469 # want to change these. 470 # Only select changes in profile when they are reversed! 471 #-- Moreover, interstellar absorption lines are usually narrow, so limit 472 # the maximum allowed line width in velocity to 6 km/s (ism turb vel ~ 473 # 3 km/s max) 474 if imin < imax and (imax-imin)*velres < 5.: 475 #- The FWHM of the irregularity is roughly vel_imax - vel_imin 476 # as the max and min in df give the strongest decrease and increase of 477 # the profile 478 velline = velfdf[fdfwhereline] 479 width = velline[imax] - velline[imin] 480 #- This leads to a gaussian sigma of 481 sigma = width/(2.*sqrt(2.*log(2.))) 482 #- velmid is where the irregularity reaches its most extreme point 483 velmid = (velline[imax] + velline[imin])/2. 484 #- The depth of the irregularity is roughly equal to the difference 485 # in max and min flux in an interval between imin-1/2*width and 486 # imax+1/2*width, ie between velmid-width and velmid+width 487 interval = np.abs(velfdf-velmid) <= width 488 fmax = max(fluxfdf[interval]) 489 fmin = min(fluxfdf[interval]) 490 depth = -abs(fmax - fmin) 491 #-- Return half the sigma, to make sure the gauss fit doesnt try 492 # too wide gaussians (eg whya co10 IRAM) 493 return (depth,velmid,sigma/2.,0) 494 else: 495 #- the found change in slope is likely due to the line itself 496 return None 497 else: 498 return None

499 500 501 502

503 -def fitLP(filename=None,lprof=None,theory=0,show=0,cfg='',convert_ms_kms=0,\ 504 vary_pars=['vexp'],i_vexp=15.0,i_gamma=1.0,do_gauss=0):

505 506 ''' 507 Fit a line profile with a soft parabola, and a Gaussian component if 508 required. 509 510 The method automatically checks if a second component is needed (eg an 511 extra absorption component). An estimate of the expansion velocity (width 512 of the profile) and an improved guess of the vlsr are given. 513 514 A guess for the gas terminal velocity is returned, as well as its error and 515 the fitted profile (sp/gaussian, and if applicable extra gaussian and the 516 full fit). 517 518 @keyword filename: The filename to the data file of the line profile. If 519 None a line profile object is expected. 520 521 (default: None) 522 @type filename: string 523 @keyword lprof: A line profile object (LPDataReader or inheriting classes) 524 If None, a filename is expected! If not None, the results 525 are saved in this object as well as returned upon method 526 call 527 528 (default: None) 529 @type lprof: LPDataReader() 530 @keyword convert_ms_kms: Convert velocity grid from m/s to km/s. 531 532 (default: 0) 533 @type convert_ms_kms: bool 534 @keyword theory: If theoretical profile, and filename is given, set vlsr to 535 0 and read two columns. lprof not relevant if True. 536 537 (default: 0) 538 @type theory: bool 539 @keyword vary_pars: The soft parabola parameters varied (can only be vexp 540 or gamma for now). The initial values for parameters 541 listed here are not used. If 'gamma' is requested, a 542 reasonable guess for i_vexp when calling the method 543 will improve the fitting results. This is done for the 544 first guess only! If a Gaussian absorption is present 545 improving these first guesses won't make much of a 546 difference. However, the first guess value for gamma 547 is still used. Vexp is always varied if absorption is 548 present. 549 550 (default: ['vexp']) 551 @type vary_pars: list[string] 552 @keyword i_vexp: The initial guess for the expansion velocity. Not relevant 553 if vexp is included in vary_pars. 554 555 (default: 15.0) 556 @type i_vexp: float 557 @keyword i_gamma: The initial guess for the gamma parameter of soft parab. 558 Not relevant if gamma is included in vary_pars. 559 560 (default: 1.0) 561 @type i_gamma: float 562 @keyword do_gauss: Force a Gaussian fit regardless of soft parabola fit 563 results. Still does the soft parabola fit first to allow 564 for comparison of parameters. 565 566 (default: 0) 567 @type do_gauss: bool 568 @keyword show: Show the results of the fit 569 570 (default: 0) 571 @type show: bool 572 573 @return: dictionary including [vexp,evexp,gamma,egamma,fitprof,gaussian,\ 574 fullfit,dintint,fgintint] 575 @rtype: dict[float,float,float,float,funclib.Function(),\ 576 funclib.Function(),funclib.Function()] 577 578 ''' 579 580 print '*******************************************' 581 if theory and not filename is None: 582 d = DataIO.readCols(filename=filename) 583 vel = d[0] 584 flux = d[1] 585 vlsr = 0.0 586 else: 587 if filename is None: 588 filename = lprof.filename 589 print '** Fitting line profile in %s.'%filename 590 if lprof is None: 591 lprof = readLineProfile(filename) 592 vel = lprof.getVelocity() 593 flux = lprof.getFlux() 594 vel = vel[-np.isnan(flux)] 595 flux = flux[-np.isnan(flux)] 596 vlsr = lprof.getVlsr() 597 598 if convert_ms_kms: 599 vel = vel/1000. 600 #-- Initial values: [peak tmb,vlsr,vexp,gamma] 601 # For the central peak value, get a first guess from the data 602 # Attempt multiple vexp values and return the best fitting case. 603 # The initial values are given, with an arbitrary value for the vexp key 604 i_mid = argmin(np.abs(vel-vlsr)) 605 peak = mean(flux[i_mid-2:i_mid+3]) 606 607 #-- Vary vexp or gamma if requested. If not requested i_vexp or i_gamma are 608 # used. 609 # Multiple values for gamma are tried and the best fitting model 610 # is then chosen based on the relative error of the fitted gamma. 611 if 'gamma' in vary_pars: 612 igammas = array([-0.5,-0.1,0.1,0.5,1.0,2.0,4.0]) 613 firstguess = varyInitialFit(vel,flux,[peak,vlsr,i_vexp,0.0],index=3,\ 614 values=igammas,vary_window=1,vary=[1,1,1,1],\ 615 function=funclib.soft_parabola) 616 i_gamma = firstguess.get_parameters()[0][3] 617 #-- varyInitialFit adapts the velocity window itself. No more 618 # assumptions needed for the expansion velocity 619 ivexps = array([50.,40.,30.,25.,20.,15.,10.]) 620 if 'vexp' in vary_pars: 621 firstguess = varyInitialFit(vel,flux,[peak,vlsr,0.,i_gamma],index=2,\ 622 values=ivexps,vary_window=1,vary=[1,1,1,1],\ 623 function=funclib.soft_parabola) 624 625 vexp = abs(firstguess.get_parameters()[0][2]) 626 window = 2. 627 print 'First guess fit, using a soft parabola:' 628 print firstguess.param2str(accuracy=5) 629 630 #-- If vexp > 100, replace with 50. This value is unrealistic, and might be 631 # improved with an extra Gaussian. If not, it will be replaced with a 632 # full Gaussian fit anyway 633 if vexp > 100: vexp = 50. 634 #-- Check whether irregularities are present in the profile. 635 # Initial parameters for a gaussian are returned if true. 636 # For this, a broad selection is made of the profile, to allow for a 637 # decent noise determination outside the line profile 638 keep = np.abs(vel-vlsr)<=(2*window*vexp) 639 velsel,fluxsel = vel[keep],flux[keep] 640 include_gauss = checkLPShape(velsel,fluxsel,vlsr,vexp,window,show=show) 641 642 #-- Do the fit of the line again, including an extra gaussian if 643 # irregularities are present. 644 if not include_gauss is None: 645 #-- fit soft para model + gaussian 646 # 1. Set up new soft parabola for several guesses of vexp 647 ivexps = list(ivexps) 648 initial = [peak,vlsr,0.,i_gamma] 649 all_init = [[p]*len(ivexps) for i,p in enumerate(initial) if i != 2] 650 all_init.insert(2,ivexps) 651 functions = [funclib.soft_parabola() for i in ivexps] 652 [ff.setup_parameters(values=initi) 653 for ff,initi in zip(functions,zip(*all_init))] 654 # 2. setup gaussian 655 gaussians = [funclib.gauss() for ff in functions] 656 # initial guesses assuming an interstellar absorption line from the 657 # checkLPShape method 658 [gg.setup_parameters(values=include_gauss,vary=[True,True,True,False]) 659 for gg in gaussians] 660 # 3. combine soft para + gaussian, and minimize fit 661 mymodels = [fit.Model(functions=[ff,gg]) 662 for ff,gg in zip(functions,gaussians)] 663 [fit.minimize(vel[np.abs(vel-vlsr)<=(init[2]*1.5)],\ 664 flux[np.abs(vel-vlsr)<=(init[2]*1.5)],\ 665 mymodel) 666 for mymodel,init in zip(mymodels,zip(*all_init))] 667 # 4. Select the best fitting result based on the error on vexp 668 mymodels = [fg 669 for fg in mymodels 670 if fg.get_parameters()[1][2] != 0.] 671 functions = [ff 672 for fg,ff in zip(mymodels,functions) 673 if fg.get_parameters()[1][2] != 0.] 674 gaussians = [gg 675 for fg,gg in zip(mymodels,gaussians) 676 if fg.get_parameters()[1][2] != 0.] 677 fitvalues = array([fg.get_parameters()[0][2] 678 for fg in mymodels]) 679 fiterrors = array([fg.get_parameters()[1][2] 680 for fg in mymodels]) 681 mymodel = mymodels[argmin(abs(fiterrors/fitvalues))] 682 finalfit = functions[argmin(abs(fiterrors/fitvalues))] 683 gaussian = gaussians[argmin(abs(fiterrors/fitvalues))] 684 print 'Improved fit, including extra Gaussian:' 685 print mymodel.param2str(accuracy=5) 686 else: 687 #-- if gamma is requested to be varied, allow another iteration on 688 # gamma with the best vexp guess we already have. 689 if 'gamma' in vary_pars: 690 finalfit = varyInitialFit(vel,flux,[peak,vlsr,vexp,0.0],\ 691 index=3,values=igammas,vary_window=1,\ 692 function=funclib.soft_parabola,\ 693 vary=[True,True,True,True]) 694 print 'Final fit with soft parabola, second gamma iteration:' 695 print finalfit.param2str(accuracy=5) 696 #-- firstguess is best we can do at the moment 697 else: 698 finalfit = firstguess 699 700 #-- If the relative error on vexp is larger than 30%, usually something 701 # funky is going on in the emission line. Try a Gaussian instead. 702 703 vexp = abs(finalfit.get_parameters()[0][2]) 704 evexp = abs(finalfit.get_parameters()[1][2]) 705 gamma = finalfit.get_parameters()[0][3] 706 egamma = finalfit.get_parameters()[1][3] 707 708 #-- Gamma has to be positive. If it isnt, dont bother with Gaussian 709 # (double peaked line profile will not be fitted well with a Gaussian!) 710 if (evexp/vexp > 0.40 and gamma > 0) or (evexp/vexp > 0.20 and vexp> 30.) \ 711 or do_gauss: 712 #-- Go back to default window to try a Gaussian fit 713 #keep = np.abs(vel-vlsr)<=(80) 714 #velselg,fluxselg = vel[keep],flux[keep] 715 do_gauss = 1 716 include_gauss = None 717 #-- FWHM is twice vexp! 718 sigmas = 2*ivexps/(2.*sqrt(2.*log(2.))) 719 finalfit = varyInitialFit(vel,flux,[peak,vlsr,0.,0.],index=2,\ 720 values=sigmas,function=funclib.gauss,\ 721 vary_window=1,vary=[True,True,True,False]) 722 vexp = abs(finalfit.get_parameters()[0][2])*(2.*sqrt(2.*log(2.)))/2. 723 evexp = abs(finalfit.get_parameters()[1][2])*(2.*sqrt(2.*log(2.)))/2. 724 gamma, egamma = None,None 725 window = 3. 726 print 'Improved fit, using a gaussian instead of soft parabola:' 727 print finalfit.param2str(accuracy=5) 728 729 #-- Extract some shared parameters between the gauss and sp. 730 peak = finalfit.get_parameters()[0][0] 731 epeak = finalfit.get_parameters()[1][0] 732 fvlsr = finalfit.get_parameters()[0][1] 733 fevlsr = finalfit.get_parameters()[1][1] 734 735 #-- Compute numerical integrations. 736 # After fitting, window for integration should be 0.6*window. vexp is 737 # not expected to be too small anymore as in checkLPShape 738 keep = np.abs(vel-vlsr)<=(0.6*window*vexp) 739 velsel = vel[keep] 740 flux_first = firstguess.evaluate(velsel) 741 flux_final = finalfit.evaluate(velsel) 742 dimb = trapz(y=flux[keep],x=velsel) 743 fi_final = trapz(y=flux_final,x=velsel) 744 print('I_mb (emission line data): %f'\ 745 %dimb) 746 print('I_mb (SP -- initial guess): %f'\ 747 %trapz(y=flux_first,x=velsel)) 748 print('I_mb (SP -- final guess): %f'\ 749 %fi_final) 750 if not include_gauss is None: 751 fitted_flux = mymodel.evaluate(velsel) 752 print('I_mb (SP + Gauss fit): %f'\ 753 %trapz(y=fitted_flux,x=velsel)) 754 print('Final v_exp guess: %.4f +/- %.4f km/s'%(vexp,evexp)) 755 if not gamma is None: 756 print('Final gamma guess: %.4f +/- %.4f'%(gamma,egamma)) 757 print('Final vlsr guess: %.4f +/- %.4f'%(fvlsr,fevlsr)) 758 print('Final peak Tmb guess at v_lsr: %.4f +/- %.4f K'%(peak,epeak)) 759 fwhm = getLPDataFWHM(lprof) 760 print('The FWHM is %.2f km/s.'%(fwhm)) 761 762 #-- plot 763 if show or cfg: 764 plt.clf() 765 #-- improve velocity window for plotting 766 keep = np.abs(vel-vlsr)<=(1.5*window*vexp) 767 velsel,fluxsel = vel[keep],flux[keep] 768 vel_highres = np.linspace(velsel[0],velsel[-1],10000) 769 flux_final_highres = finalfit.evaluate(vel_highres) 770 flux_first_highres = firstguess.evaluate(vel_highres) 771 if not include_gauss is None: 772 flux_full_highres = mymodel.evaluate(vel_highres) 773 if show: 774 plt.step(velsel,fluxsel,'-r',where='mid',lw=3,\ 775 label='Observed profile') 776 plt.plot(vel_highres,flux_first_highres,'b-',lw=3,\ 777 label='First guess') 778 plt.plot(vel_highres,flux_final_highres,'g--',lw=3,\ 779 label='Improved guess') 780 if not include_gauss is None: 781 plt.plot(vel_highres,flux_full_highres,'g-',lw=2,\ 782 label='Full fit (including Gaussian)') 783 leg = plt.legend(loc='best',fancybox=True) 784 leg.get_frame().set_alpha(0.5) 785 plt.show() 786 if cfg: 787 pf = '%s_fitted_%s'%(do_gauss and 'gaussFit' or 'SPFit',\ 788 os.path.split(filename)[1]) 789 keytags = ['Observed profile','Improved guess'] 790 line_types = ['-r','-b',] 791 x = [velsel,vel_highres] 792 y = [fluxsel,flux_final_highres] 793 if not include_gauss is None: 794 line_types.append('g--') 795 x.append(vel_highres) 796 y.append(flux_full_highres) 797 keytags.append('Full fit (including Gaussian)') 798 pf = Plotting2.plotCols(x=x,y=y,filename=pf,cfg=cfg,linewidth=5,\ 799 yaxis='$T_\mathrm{mb}\ (\mathrm{K})$',\ 800 xaxis='$v (\mathrm{km}/\mathrm{s})$',\ 801 keytags=keytags,line_types=line_types,\ 802 histoplot=[0]) 803 print 'Your figure can be found at %s .'%pf 804 #-- Collecting all relevant results and returning. 805 results = dict() 806 #-- If a Gaussian was used for the main profile fit 807 results['do_gauss'] = do_gauss 808 #-- Fitted parameters and errors 809 results['vlsr'] = fvlsr 810 results['evlsr'] = fevlsr 811 results['vexp'] = vexp 812 results['evexp'] = evexp 813 results['fwhm'] = fwhm 814 results['peak'] = peak 815 results['epeak'] = epeak 816 817 #-- Gamma is None if no soft parabola was fitted 818 results['gamma'] = gamma 819 results['egamma'] = egamma 820 821 #-- Integrated line strengths: main line fit, data themselves, fit window 822 results['fgintint'] = fi_final 823 results['dintint'] = dimb 824 results['intwindow'] = window*0.6 825 826 #-- Saving parameters for later evaluation. Full fit is accessible by 827 # making the functions separately and setting pars, then using fit.Model 828 results['fitprof'] = (do_gauss and 'gauss' or 'soft_parabola',\ 829 list(finalfit.get_parameters()[0])) 830 if not include_gauss is None: 831 results['fitabs'] = ('gauss',list(gaussian.get_parameters()[0])) 832 else: 833 results['fitabs'] = None 834 835 return results

836

Source Code for Module ComboCode.cc.data.LPTools