ivs.sed.limbdark

1 # -*- coding: utf-8 -*- 2 """ 3 Interface to the limb-darkening library. 4 5 Section 1. Basic interface 6 ========================== 7 8 Retrieve limb darkening coeffcients and goodness of fit parameters: 9 10 >>> coefs,fac,ssr,idiff = get_coeffs(teff=10000,logg=4.0,'JOHNSON.V',law='claret') 11 12 When you need this for a large grid of parameters, consider using L{get_itable}. 13 The available laws are C{claret}, C{linear}, C{logarithmic}, C{quadratic} and 14 C{power}. 15 16 Retrieve a normalised passband-integrated limb darkening profile via: 17 18 >>> mu,intensities = get_limbdarkening(teff=10000,logg=4.0,photbands=['JOHNSON.V'],normalised=True) 19 20 and fit a particular law via: 21 22 >>> coeffs,ssr,idiff = fit_law(mu,intensities[:,0],law='claret',fitmethod='equidist_mu_leastsq') 23 >>> print(coeffs) 24 [ 2.00943266 -3.71261279 4.25469623 -1.70466934] 25 26 You can evaluate this law on a particular angle-sampling grid: 27 28 >>> mus = np.linspace(0,1,1000) 29 >>> mu_law = ld_claret(mus,coeffs) 30 31 Or on a disk-sampling grid: 32 33 >>> rs = np.linspace(0,1,1000) 34 >>> r_law = ld_claret(_mu(rs),coeffs) 35 36 Make a simple plot: 37 38 >>> p = pl.figure() 39 >>> p,q = pl.subplot(121),pl.title('Angle-sampled') 40 >>> p = pl.plot(mu,intensities,'ko-') 41 >>> p = pl.plot(mus,mu_law,'r-',lw=2) 42 >>> p,q = pl.subplot(122),pl.title('Disk-sampled') 43 >>> p = pl.plot(_r(mu),intensities,'ko-') 44 >>> p = pl.plot(rs,r_law,'r-',lw=2) 45 46 ]include figure]]ivs_limbdark_basic.png] 47 48 Section 2. Fit comparison 49 ========================= 50 51 You can fit a limb darkening law in various ways: using different laws, different 52 optimizers and different x-coordinates (limb angle mu or radius). In the figure 53 below, you can see the influence of some of these choices: 54 55 >>> p = pl.figure() 56 57 >>> laws = ['claret','linear','logarithmic','quadratic','power'] 58 >>> fitmethods = ['equidist_mu_leastsq','equidist_r_leastsq','leastsq'] 59 60 >>> r,rs = _r(mu),_r(mus) 61 >>> integ_mu = np.trapz(intensities[:,0],x=mu) 62 >>> integ_r = np.trapz(intensities[:,0],x=r) 63 >>> for i,fitmethod in enumerate(fitmethods): 64 ... p = pl.subplot(3,2,2*i+1) 65 ... p = pl.title(fitmethod) 66 ... p = pl.plot(mu,intensities[:,0],'ko-') 67 ... p = pl.xlabel('$\mu$ = cos(limb angle)') 68 ... p = pl.ylabel('Normalised intensity') 69 ... p = pl.subplot(3,2,2*i+2) 70 ... p = pl.title(fitmethod) 71 ... p = pl.plot(r,intensities[:,0],'ko-') 72 ... p = pl.xlabel('r') 73 ... p = pl.ylabel('Normalised intensity') 74 ... for law in laws: 75 ... cl,ssr,idiff = fit_law(mu,intensities[:,0],law=law,fitmethod=fitmethod) 76 ... print("{:12s} ({:19s}): {}".format(law,fitmethod,cl)) 77 ... Imus = globals()['ld_{}'.format(law)](mus,cl) 78 ... p = pl.subplot(3,2,2*i+1) 79 ... idiff = (np.trapz(Imus,x=mus)-integ_mu)/integ_mu 80 ... p = pl.plot(mus,Imus,'-',lw=2,label="{} ({:.3f},{:.3f}%)".format(law,ssr,idiff*100)) 81 ... p = pl.subplot(3,2,2*i+2) 82 ... idiff = (np.trapz(Imus,x=rs)-integ_r)/integ_r 83 ... p = pl.plot(rs,Imus,'-',lw=2,label="{} ({:.3f},{:.3}%)".format(law,ssr,idiff*100)) 84 ... p = pl.subplot(3,2,2*i+1) 85 ... leg = pl.legend(loc='best',fancybox=True,prop=dict(size='small')) 86 ... leg.get_frame().set_alpha(0.5) 87 ... p = pl.subplot(3,2,2*i+2) 88 ... leg = pl.legend(loc='best',fancybox=True,prop=dict(size='small')) 89 ... leg.get_frame().set_alpha(0.5) 90 claret (equidist_mu_leastsq): [ 2.00943266 -3.71261279 4.25469623 -1.70466934] 91 linear (equidist_mu_leastsq): [ 0.48655636] 92 divide by zero encountered in log 93 invalid value encountered in multiply 94 logarithmic (equidist_mu_leastsq): [ 0.6 0. ] 95 quadratic (equidist_mu_leastsq): [ 0.21209044 0.36591797] 96 power (equidist_mu_leastsq): [ 0.29025864] 97 claret (equidist_r_leastsq ): [ 1.99805089 -3.04504952 3.02780048 -1.09105603] 98 linear (equidist_r_leastsq ): [ 0.42641628] 99 logarithmic (equidist_r_leastsq ): [ 0.6 0. ] 100 quadratic (equidist_r_leastsq ): [ 0.28650391 0.24476527] 101 power (equidist_r_leastsq ): [ 0.31195553] 102 claret (leastsq ): [ 3.49169665 -8.96767096 11.16873314 -4.72869922] 103 linear (leastsq ): [ 0.57180245] 104 logarithmic (leastsq ): [ 0.6 0. ] 105 quadratic (leastsq ): [ 0.00717934 0.65416204] 106 power (leastsq ): [ 0.27032565] 107 108 ]include figure]]ivs_limbdark_fitcomp.png] 109 110 Author: Pieter Degroote, with thanks to Steven Bloemen. 111 """ 112 import logging 113 import os 114 import itertools 115 import astropy.io.fits as pf 116 import numpy as np 117 from scipy.optimize import leastsq,fmin 118 from scipy.interpolate import splrep, splev 119 from scipy.interpolate import LinearNDInterpolator 120 from Scientific.Functions.Interpolation import InterpolatingFunction 121 from ivs.aux import loggers 122 from ivs.sed import reddening 123 from ivs.sed import model 124 from ivs.spectra import tools 125 from ivs.units import constants 126 from ivs.aux.decorators import memoized,clear_memoization 127 from ivs import config 128 129 logger = logging.getLogger("SED.LIMBDARK") 130 logger.addHandler(loggers.NullHandler()) 131 132 #-- default values for grids 133 defaults = dict(grid='kurucz',odfnew=False,z=+0.0,vturb=2, 134 alpha=False,nover=False) 135 #-- relative location of the grids 136 basedir = 'ldtables'

137 138 #{ Interface to the library 139 140 -def set_defaults(**kwargs):

141 """ 142 Set defaults of this module 143 144 If you give no keyword arguments, the default values will be reset. 145 """ 146 clear_memoization(keys=['ivs.sed.ld']) 147 if not kwargs: 148 kwargs = dict(grid='kurucz',odfnew=True,z=+0.0,vturb=2, 149 alpha=False,nover=False, # KURUCZ 150 He=97, # WD 151 t=1.0,a=0.0,c=0.5,m=1.0,co=1.05) # MARCS and COMARCS 152 153 for key in kwargs: 154 if key in defaults: 155 defaults[key] = kwargs[key] 156 logger.info('Set %s to %s'%(key,kwargs[key]))

157

158 159 160 -def get_gridnames():

161 """ 162 Return a list of available grid names. 163 164 @return: list of grid names 165 @rtype: list of str 166 """ 167 return ['kurucz']

168

169 170 -def get_grid_dimensions(**kwargs):

171 """ 172 Retrieve the possible effective temperatures and gravities from a grid. 173 174 @rtype: (ndarray,ndarray) 175 @return: effective temperatures, gravities 176 """ 177 gridfile = get_file(**kwargs) 178 ff = pf.open(gridfile) 179 teffs = [] 180 loggs = [] 181 for mod in ff[1:]: 182 teffs.append(float(mod.header['TEFF'])) 183 loggs.append(float(mod.header['LOGG'])) 184 ff.close() 185 return np.array(teffs),np.array(loggs)

186

187 188 @memoized 189 -def get_grid_mesh(wave=None,teffrange=None,loggrange=None,**kwargs):

190 """ 191 Return InterpolatingFunction spanning the available grid of atmosphere models. 192 193 WARNING: the grid must be entirely defined on a mesh grid, but it does not 194 need to be equidistant. 195 196 It is thus the user's responsibility to know whether the grid is evenly 197 spaced in logg and teff (e.g. this is not so for the CMFGEN models). 198 199 You can supply your own wavelength range, since the grid models' 200 resolution are not necessarily homogeneous. If not, the first wavelength 201 array found in the grid will be used as a template. 202 203 It might take a long a time and cost a lot of memory if you load the entire 204 grid. Therefor, you can also set range of temperature and gravity. 205 206 @param wave: wavelength to define the grid on 207 @type wave: ndarray 208 @param teffrange: starting and ending of the grid in teff 209 @type teffrange: tuple of floats 210 @param loggrange: starting and ending of the grid in logg 211 @type loggrange: tuple of floats 212 @return: wavelengths, teffs, loggs and fluxes of grid, and the interpolating 213 function 214 @rtype: (3x1Darray,3Darray,interp_func) 215 """ 216 #-- get the dimensions of the grid 217 teffs,loggs = get_grid_dimensions(**kwargs) 218 #-- build flux grid, assuming a perfectly sampled grid (needs not to be 219 # equidistant) 220 if teffrange is not None: 221 sa = (teffrange[0]<=teffs) & (teffs<=teffrange[1]) 222 teffs = teffs[sa] 223 if loggrange is not None: 224 sa = (loggrange[0]<=loggs) & (loggs<=loggrange[1]) 225 loggs = loggs[sa] 226 227 #-- run over teff and logg, and interpolate the models onto the supplied 228 # wavelength range 229 gridfile = get_file(**kwargs) 230 231 if wave is not None: 232 table = np.zeros((len(teffs),len(wave),18)) 233 ff = pf.open(gridfile) 234 for i,(teff,logg) in enumerate(zip(teffs,loggs)): 235 mod_name = "T%05d_logg%01.02f" %(teff,logg) 236 mod = ff[mod_name] 237 imu = np.array(mod.columns.names[1:], float) 238 itable = np.array(mod.data.tolist())[:,1:] 239 iwave = mod.data.field('wavelength') 240 #-- if there is no wavelength range given, we assume that 241 # the whole grid has the same resolution, and the first 242 # wave-array will be used as a template 243 if wave is None: 244 wave = iwave 245 table = np.zeros((len(teffs),len(wave),len(imu))) 246 try: 247 table[i] = itable 248 except: 249 for j,jimu in enumerate(imu): 250 table[i,:,j] = np.interp(wave,iwave,itable[:,j]) 251 ff.close() 252 table_grid = LinearNDInterpolator(np.array([np.log10(teffs),loggs]).T,table) 253 return imu,wave,teffs,loggs,table,table_grid

254

255 256 257 -def get_file(integrated=False,**kwargs):

258 """ 259 Retrieve the filename containing the specified SED grid. 260 261 The keyword arguments are specific to the kind of grid you're using. 262 263 Basic keywords are 'grid' for the name of the grid, and 'z' for metallicity. 264 For other keywords, see the source code. 265 266 Available grids and example keywords: 267 - grid='kurucz93': 268 * metallicity (z): m01 is -0.1 log metal abundance relative to solar (solar abundances from Anders and Grevesse 1989) 269 * metallicity (z): p01 is +0.1 log metal abundance relative to solar (solar abundances from Anders and Grevesse 1989) 270 * alpha enhancement (alpha): True means alpha enhanced (+0.4) 271 * turbulent velocity (vturb): vturb in km/s 272 * nover= True means no overshoot 273 * odfnew=True means no overshoot but with better opacities and abundances 274 275 @param integrated: choose integrated version of the grid 276 @type integrated: boolean 277 @keyword grid: gridname (default Kurucz) 278 @type grid: str 279 @return: gridfile 280 @rtype: str 281 """ 282 #-- possibly you give a filename 283 grid = kwargs.get('grid',defaults['grid']) 284 if os.path.isfile(grid): 285 return grid 286 287 #-- general 288 z = kwargs.get('z',defaults['z']) 289 #-- only for Kurucz 290 vturb = int(kwargs.get('vturb',defaults['vturb'])) 291 odfnew = kwargs.get('odfnew',defaults['odfnew']) 292 alpha = kwargs.get('alpha',defaults['alpha']) 293 nover = kwargs.get('nover',defaults['nover']) 294 295 #-- figure out what grid to use 296 if grid=='kurucz': 297 if not isinstance(z,str): z = '%.1f'%(z) 298 if not isinstance(vturb,str): vturb = '%d'%(vturb) 299 if not alpha and not nover and not odfnew: 300 basename = 'kurucz93_z%s_k%s_ld.fits'%(z,vturb) 301 elif alpha and odfnew: 302 basename = 'kurucz93_z%s_ak%sodfnew_ld.fits'%(z,vturb) 303 elif odfnew: 304 basename = 'kurucz93_z%s_k%sodfnew_ld.fits'%(z,vturb) 305 elif nover: 306 basename = 'kurucz93_z%s_k%snover_ld.fits'%(z,vturb) 307 else: 308 basename = grid 309 310 #-- retrieve the absolute path of the file and check if it exists: 311 if not '*' in basename: 312 if integrated: 313 grid = config.get_datafile(basedir,'i'+basename) 314 else: 315 grid = config.get_datafile(basedir,basename) 316 #-- we could also ask for a list of files, when wildcards are given: 317 else: 318 grid = config.glob(basedir,'i'+basename) 319 if integrated: 320 grid = config.glob(basedir,'i'+basename) 321 else: 322 grid = config.glob(basedir,basename) 323 logger.debug('Selected %s'%(grid)) 324 325 return grid

326

327 328 329 -def get_table(teff=None,logg=None,ebv=None,vrad=None,star=None, 330 wave_units='AA',flux_units='erg/cm2/s/AA/sr',**kwargs):

331 """ 332 Retrieve the specific intensity of a model atmosphere. 333 334 ebv is reddening 335 vrad is radial velocity: positive is redshift, negative is blueshift (km/s!) 336 337 extra kwargs are for reddening 338 339 You get limb angles, wavelength and a table. The shape of the table is 340 (N_wave,N_mu). 341 342 WARNING: wave and flux units cannot be specificed for the moment. 343 344 >>> mu,wave,table = get_table(10000,4.0) 345 346 >>> p = pl.figure() 347 >>> ax1 = pl.subplot(221) 348 >>> p = pl.title('E(B-V)=0, vrad=0') 349 >>> ax = pl.gca() 350 >>> ax.set_color_cycle([pl.cm.spectral(i) for i in np.linspace(0,1,len(mu))]) 351 >>> p = pl.loglog(wave,table) 352 >>> p = pl.xlim(700,15000) 353 >>> p = pl.ylim(1e3,1e8) 354 355 >>> mu,wave,table = get_table(10000,4.0,ebv=0.5) 356 357 >>> p = pl.subplot(222,sharex=ax1,sharey=ax1) 358 >>> p = pl.title('E(B-V)=0.5, vrad=0') 359 >>> ax = pl.gca() 360 >>> ax.set_color_cycle([pl.cm.spectral(i) for i in np.linspace(0,1,len(mu))]) 361 >>> p = pl.loglog(wave,table) 362 >>> p = pl.xlim(700,15000) 363 >>> p = pl.ylim(1e3,1e8) 364 365 >>> mu,wave,table = get_table(10000,4.0,vrad=-1000.) 366 367 >>> p = pl.subplot(223,sharex=ax1,sharey=ax1) 368 >>> p = pl.title('E(B-V)=0., vrad=-10000.') 369 >>> ax = pl.gca() 370 >>> ax.set_color_cycle([pl.cm.spectral(i) for i in np.linspace(0,1,len(mu))]) 371 >>> p = pl.loglog(wave,table) 372 >>> p = pl.xlim(700,15000) 373 >>> p = pl.ylim(1e3,1e8) 374 375 >>> mu,wave,table = get_table(10000,4.0,vrad=-1000.,ebv=0.5) 376 377 >>> p = pl.subplot(224,sharex=ax1,sharey=ax1) 378 >>> p = pl.title('E(B-V)=0.5, vrad=-10000.') 379 >>> ax = pl.gca() 380 >>> ax.set_color_cycle([pl.cm.spectral(i) for i in np.linspace(0,1,len(mu))]) 381 >>> p = pl.loglog(wave,table) 382 >>> p = pl.xlim(700,15000) 383 >>> p = pl.ylim(1e3,1e8) 384 385 >>> mu,wave,table = get_table(10050,4.12) 386 387 >>> p = pl.figure() 388 >>> pl.gca().set_color_cycle([pl.cm.spectral(i) for i in np.linspace(0,1,len(mu))]) 389 >>> p = pl.loglog(wave,table) 390 391 392 @param teff: effective temperature (K) 393 @type teff: float 394 @param logg: log surface gravity (cgs, dex) 395 @type logg: float 396 @param ebv: reddening (mag) 397 @type ebv: float 398 @param vrad: radial velocity (for doppler shifting) (km/s) 399 @type vrad: float 400 @return: mu angles, wavelengths, table (Nwave x Nmu) 401 @rtype: array, array, array 402 """ 403 #-- get the FITS-file containing the tables 404 gridfile = get_file(**kwargs) 405 406 #-- read the file: 407 ff = pf.open(gridfile) 408 409 teff = float(teff) 410 logg = float(logg) 411 412 #-- if we have a grid model, no need for interpolation 413 try: 414 #-- extenstion name as in fits files prepared by Steven 415 mod_name = "T%05d_logg%01.02f" %(teff,logg) 416 mod = ff[mod_name] 417 mu = np.array(mod.columns.names[1:], float) 418 table = np.array(mod.data.tolist())[:,1:] 419 wave = mod.data.field('wavelength') 420 logger.debug('Model LD taken directly from file (%s)'%(os.path.basename(gridfile))) 421 except KeyError: 422 mu,wave,teffs,loggs,flux,flux_grid = get_grid_mesh(**kwargs) 423 logger.debug('Model LD interpolated from grid %s (%s)'%(os.path.basename(gridfile),kwargs)) 424 wave = wave + 0. 425 table = flux_grid(np.log10(teff),logg) + 0. 426 427 ff.close() 428 429 #-- velocity shift if necessary 430 if vrad is not None and vrad!=0: 431 cc = constants.cc/1000. #speed of light in km/s 432 for i in range(len(mu)): 433 flux_shift = tools.doppler_shift(wave,vrad,flux=table[:,i]) 434 table[:,i] = flux_shift - 5.*vrad/cc*table[:,i] 435 436 #-- redden if necessary 437 if ebv is not None and ebv>0: 438 for i in range(len(mu)): 439 table[:,i] = reddening.redden(table[:,i],wave=wave,ebv=ebv,rtype='flux',**kwargs) 440 441 #-- that's it! 442 return mu,wave,table

443

444 445 446 447 448 449 450 451 -def get_itable2(teff=None,logg=None,theta=None,mu=1,photbands=None,absolute=False,**kwargs):

452 """ 453 mu=1 is center of disk 454 """ 455 if theta is not None: 456 mu = np.cos(theta) 457 458 try: 459 out = get_ld_grid(photbands,integrated=True,**kwargs)(teff,logg) 460 except ValueError: 461 print 'Used teff and logg',teff,logg 462 raise 463 a1x_,a2x_,a3x_,a4x_, I_x1 = out.reshape((len(photbands),5)).T 464 Imu = ld_eval(mu,[a1x_,a2x_,a3x_,a4x_]) 465 if absolute: 466 return Imu*I_x1 467 else: 468 return Imu

469

470 @memoized 471 -def _get_itable_markers(photband,gridfile, 472 teffrange=(-np.inf,np.inf),loggrange=(-np.inf,np.inf)):

473 """ 474 Get a list of markers to more easily retrieve integrated fluxes. 475 """ 476 clear_memoization() 477 478 ff = pf.open(gridfile) 479 ext = ff[photband] 480 columns = ext.columns.names 481 names = ['teff','logg'] 482 grid = [ext.data.field('teff'), 483 ext.data.field('logg')] 484 if 'ebv' in columns: 485 names.append('ebv') 486 grid.append(ext.data.field('ebv')) 487 if 'vrad' in columns: 488 names.append('vrad') 489 grid.append(ext.data.field('vrad')) 490 491 grid_axes = [np.sort(list(set(i))) for i in grid] 492 493 #-- we construct an array representing the teff-logg-ebv content, but 494 # in one number: 50000400 means: 495 # T=50000,logg=4.0 496 N = len(grid[0]) 497 markers = np.zeros(N,float) 498 gridpnts = np.zeros((N,len(grid))) 499 pars = np.zeros((N,5)) 500 501 for i,entries in enumerate(zip(*grid)): 502 markers[i] = encode(**{key:entry for key,entry in zip(names,entries)}) 503 gridpnts[i]= entries 504 pars[i] = list(ext.data[i][-5:]) 505 ff.close() 506 sa = np.argsort(markers) 507 print 'read in gridfile',gridfile 508 pars[:,-1] = np.log10(pars[:,-1]) 509 return markers[sa],grid_axes,gridpnts[sa],pars[sa]

510

511 512 513 514 -def get_limbdarkening(teff=None,logg=None,ebv=None,vrad=None,z=None,photbands=None,normalised=False,**kwargs):

515 """ 516 Retrieve a limb darkening law for a specific star and specific bandpass. 517 518 Possibility to include reddening via EB-V parameter. If not given, 519 reddening will not be performed... 520 521 You choose your own reddening function. 522 523 See e.g. Heyrovsky et al., 2007 524 525 If you specify one angle (mu in radians), it will take the closest match 526 from the grid. 527 528 Mu = cos(theta) where theta is the angle between the surface and the line 529 of sight. mu=1 means theta=0 means center of the star. 530 531 Example usage: 532 533 >>> teff,logg = 5000,4.5 534 >>> mu,intensities = get_limbdarkening(teff=teff,logg=logg,photbands=['JOHNSON.V']) 535 536 @keyword teff: effective temperature 537 @type teff: float 538 @keyword logg: logarithmic gravity (cgs) 539 @type logg: float 540 @keyword system: bandpass system 541 @type system: string 542 @keyword photbands: bandpass filters 543 @type photbands: list of strings 544 @keyword ebv: reddening coefficient 545 @type ebv: float 546 @keyword vrad: radial velocity (+ is redshift, - is blueshift) 547 @type vrad: float 548 @keyword mu: specificy specific angle 549 @type mu: float 550 """ 551 if z is None: 552 z = defaults['z'] 553 #-- retrieve model atmosphere for a given teff and logg 554 mus,wave,table = get_table(teff=teff,logg=logg,ebv=ebv,vrad=vrad,z=z,**kwargs) 555 #-- compute intensity over the stellar disk, and normalise 556 intensities = np.zeros((len(mus),len(photbands))) 557 for i in range(len(mus)): 558 intensities[i] = model.synthetic_flux(wave,table[:,i],photbands) 559 if normalised: 560 intensities/= intensities.max(axis=0) 561 #-- or compute the intensity only for one angle: 562 logger.info('Calculated LD') 563 return mus,intensities

564

565 -def get_coeffs(teff,logg,photband,ebv=None,vrad=None,law='claret',fitmethod='equidist_r_leastsq'):

566 """ 567 Retrieve limb darkening coefficients on the fly. 568 569 This is particularly useful if you only need a few and speed is not really 570 a problem (although this function is nearly instantaneous, looping over it 571 will expose that it's actually pretty slow). 572 573 @keyword teff: effective temperature 574 @type teff: float 575 @keyword logg: logarithmic gravity (cgs) 576 @type logg: float 577 @keyword photbands: bandpass filters 578 @type photbands: list of strings 579 @keyword ebv: reddening coefficient 580 @type ebv: float 581 @keyword vrad: radial velocity (+ is redshift, - is blueshift) 582 @type vrad: float 583 """ 584 mu,intensities = get_limbdarkening(teff=teff,logg=logg,ebv=ebv,vrad=vrad,photbands=[photband]) 585 norm_factor = intensities.max(axis=0) 586 coeffs,ssr,idiff = fit_law(mu,intensities[:,0]/norm_factor,law=law,fitmethod=fitmethod) 587 return coeffs,norm_factor,ssr,idiff

588

589 590 591 592 -def get_ld_grid_dimensions(**kwargs):

593 """ 594 Returns the gridpoints of the limbdarkening grid (not unique values). 595 """ 596 #-- get the FITS-file containing the tables 597 gridfile = get_file(**kwargs) 598 #-- the teff and logg range is the same for every extension, so just 599 # take the first one 600 ff = pf.open(gridfile) 601 teff_,logg_ = ff[1].data.field('Teff'),ff[1].data.field('logg') 602 ff.close() 603 return teff_,logg_

604

605 -def intensity_moment(coeffs,ll=0,**kwargs):

606 """ 607 Calculate the intensity moment (see Townsend 2003, eq. 39). 608 609 Test analytical versus numerical implementation: 610 611 #>>> photband = 'JOHNSON.V' 612 #>>> l,logg = 2.,4.0 613 #>>> gridfile = 'tables/ikurucz93_z0.0_k2_ld.fits' 614 #>>> mu = np.linspace(0,1,100000) 615 #>>> P_l = legendre(l)(mu) 616 #>>> check = [] 617 #>>> for teff in np.linspace(9000,12000,19): 618 #... a1x,a2x,a3x,a4x,I_x1 = get_itable(gridfile,teff=teff,logg=logg,mu=1,photband=photband,evaluate=False) 619 #... coeffs = [a1x,a2x,a3x,a4x,I_x1] 620 #... Ix_mu = ld_claret(mu,coeffs) 621 #... Ilx1 = I_x1*np.trapz(P_l*mu*Ix_mu,x=mu) 622 #... a0x = 1 - a1x - a2x - a3x - a4x 623 #... limb_coeffs = [a0x,a1x,a2x,a3x,a4x] 624 #... Ilx2 = 0 625 #... for r,ar in enumerate(limb_coeffs): 626 #... s = 1 + r/2. 627 #... myIls = _I_ls(l,s) 628 #... Ilx2 += ar*myIls 629 #... Ilx2 = I_x1*Ilx2 630 #... Ilx3 = intensity_moment(teff,logg,photband,coeffs) 631 #... check.append(abs(Ilx1-Ilx2)<0.1) 632 #... check.append(abs(Ilx1-Ilx3)<0.1) 633 #>>> np.all(np.array(check)) 634 #True 635 636 637 @parameter teff: effecitve temperature (K) 638 @type teff: float 639 @parameter logg: log of surface gravity (cgs) 640 @type logg: float 641 @parameter ll: degree of the mode 642 @type ll: float 643 @parameter photband: photometric passbands 644 @type photband: list of strings ( or iterable container) 645 @keyword full_output: if True, returns intensity, coefficients and integrals 646 separately 647 @type full_output: boolean 648 @return: intensity moment 649 """ 650 full_output = kwargs.get('full_output',False) 651 #-- notation of Townsend 2002: coeffs in code are hat coeffs in the paper 652 # (for i=1..4 they are the same) 653 #-- get the LD coefficients at the given temperature and logg 654 a1x_,a2x_,a3x_,a4x_, I_x1 = coeffs 655 a0x_ = 1 - a1x_ - a2x_ - a3x_ - a4x_ 656 limb_coeffs = np.array([a0x_,a1x_,a2x_,a3x_,a4x_]) 657 658 #-- compute intensity moment via helper coefficients 659 int_moms = np.array([_I_ls(ll,1 + r/2.) for r in range(0,5,1)]) 660 #int_moms = np.outer(int_moms,np.ones(1)) 661 I_lx = I_x1 * (limb_coeffs * int_moms).sum(axis=0) 662 #-- return value depends on keyword 663 if full_output: 664 return I_x1,limb_coeffs,int_moms 665 else: 666 return I_lx

667

668 669 670 #} 671 #{ Limbdarkening laws 672 673 -def ld_eval(mu,coeffs):

674 """ 675 Evaluate Claret's limb darkening law. 676 """ 677 return ld_claret(mu,coeffs)

678

679 -def ld_claret(mu,coeffs):

680 """ 681 Claret's limb darkening law. 682 """ 683 a1,a2,a3,a4 = coeffs 684 Imu = 1-a1*(1-mu**0.5)-a2*(1-mu)-a3*(1-mu**1.5)-a4*(1-mu**2.) 685 return Imu

686

687 -def ld_linear(mu,coeffs):

688 """ 689 Linear or linear cosine law 690 """ 691 return 1-coeffs[0]*(1-mu)

692

693 -def ld_nonlinear(mu,coeffs):

694 """ 695 Nonlinear or logarithmic law 696 """ 697 return 1-coeffs[0]*(1-mu)-coeffs[1]*mu*np.log(mu)

698

699 -def ld_logarithmic(mu,coeffs):

700 """ 701 Nonlinear or logarithmic law 702 """ 703 return ld_nonlinear(mu,coeffs)

704

705 -def ld_quadratic(mu,coeffs):

706 """ 707 Quadratic law 708 """ 709 return 1-coeffs[0]*(1-mu)-coeffs[1]*(1-mu)**2.0

710

711 -def ld_uniform(mu,coeffs):

712 """ 713 Uniform law. 714 """ 715 return 1.

716

717 -def ld_power(mu,coeffs):

718 """ 719 Power law. 720 """ 721 return mu**coeffs[0]

722

723 #} 724 725 #{ Fitting routines (thanks to Steven Bloemen) 726 727 -def fit_law(mu,Imu,law='claret',fitmethod='equidist_r_leastsq'):

728 """ 729 Fit an LD law to a sampled set of limb angles/intensities. 730 731 In my (Pieter) experience, C{fitmethod='equidist_r_leastsq' seems 732 appropriate for the Kurucz models. 733 734 Make sure the intensities are normalised! 735 736 >>> mu,intensities = get_limbdarkening(teff=10000,logg=4.0,photbands=['JOHNSON.V'],normalised=True) 737 >>> p = pl.figure() 738 >>> p = pl.plot(mu,intensities[:,0],'ko-') 739 >>> cl,ssr,rdiff = fit_law(mu,intensities[:,0],law='claret') 740 741 >>> mus = np.linspace(0,1,1000) 742 >>> Imus = ld_claret(mus,cl) 743 >>> p = pl.plot(mus,Imus,'r-',lw=2) 744 745 746 747 @return: coefficients, sum of squared residuals,relative flux difference between prediction and model integrated intensity 748 @rtype: array, float, float 749 """ 750 #-- prepare array for coefficients and set the initial guess 751 Ncoeffs = dict(claret=4,linear=1,nonlinear=2,logarithmic=2,quadratic=2, 752 power=1) 753 c0 = np.zeros(Ncoeffs[law]) 754 c0[0] = 0.6 755 #-- do the fitting 756 if fitmethod=='leastsq': 757 (csol, ierr) = leastsq(ldres_leastsq, c0, args=(mu,Imu,law)) 758 elif fitmethod=='fmin': 759 csol = fmin(ldres_fmin, c0, maxiter=1000, maxfun=2000,args=(mu,Imu,law),disp=0) 760 elif fitmethod=='equidist_mu_leastsq': 761 mu_order = np.argsort(mu) 762 tck = splrep(mu[mu_order],Imu[mu_order],s=0.0, k=2) 763 mu_spl = np.linspace(mu[mu_order][0],1,5000) 764 Imu_spl = splev(mu_spl,tck,der=0) 765 (csol, ierr) = leastsq(ldres_leastsq, c0, args=(mu_spl,Imu_spl,law)) 766 elif fitmethod=='equidist_r_leastsq': 767 mu_order = np.argsort(mu) 768 tck = splrep(mu[mu_order],Imu[mu_order],s=0., k=2) 769 r_spl = np.linspace(mu[mu_order][0],1,5000) 770 mu_spl = np.sqrt(1-r_spl**2) 771 Imu_spl = splev(mu_spl,tck,der=0) 772 (csol,ierr) = leastsq(ldres_leastsq, c0, args=(mu_spl,Imu_spl,law)) 773 elif fitmethod=='equidist_mu_fmin': 774 mu_order = np.argsort(mu) 775 tck = splrep(mu[mu_order],Imu[mu_order],k=2, s=0.0) 776 mu_spl = np.linspace(mu[mu_order][0],1,5000) 777 Imu_spl = splev(mu_spl,tck,der=0) 778 csol = fmin(ldres_fmin, c0, maxiter=1000, maxfun=2000,args=(mu_spl,Imu_spl,law),disp=0) 779 elif fitmethod=='equidist_r_fmin': 780 mu_order = np.argsort(mu) 781 tck = splrep(mu[mu_order],Imu[mu_order],k=2, s=0.0) 782 r_spl = np.linspace(mu[mu_order][0],1,5000) 783 mu_spl = np.sqrt(1-r_spl**2) 784 Imu_spl = splev(mu_spl,tck,der=0) 785 csol = fmin(ldres_fmin, c0, maxiter=1000, maxfun=2000,args=(mu_spl,Imu_spl,law),disp=0) 786 else: 787 raise ValueError("Fitmethod {} not recognised".format(fitmethod)) 788 myfit = globals()['ld_%s'%(law)](mu,csol) 789 res = np.sum(Imu - myfit)**2 790 int1,int2 = np.trapz(Imu,x=mu),np.trapz(myfit,x=mu) 791 dflux = (int1-int2)/int1 792 return csol,res,dflux

793

794 795 -def fit_law_to_grid(photband,vrads=[0],ebvs=[0],zs=[0], 796 law='claret',fitmethod='equidist_r_leastsq',**kwargs):

797 """ 798 Gets the grid and fits LD law to all the models. 799 800 Does not use mu=0 point in the fitting process. 801 """ 802 teffs, loggs = get_grid_dimensions(**kwargs) 803 teffs=teffs 804 loggs=loggs 805 grid_pars = [] 806 grid_coeffs = [] 807 Imu1s = [] 808 logger.info('Fitting photband {}'.format(photband)) 809 for teff_, logg_ in zip(teffs, loggs): 810 for ebv_,vrad_,z_ in itertools.product(ebvs,vrads,zs): 811 logger.info('teff={}, logg={}, ebv={}, vrad={}, z={}'.format(teff_, logg_,ebv_,vrad_,z_)) 812 mu, Imu = get_limbdarkening(teff=teff_, logg=logg_, ebv=ebv_,vrad=vrad_,z=z_,photbands=[photband],**kwargs) 813 Imu1 = Imu.max() 814 Imu = Imu[:,0]/Imu1 815 coeffs,res,dflux = fit_law(mu[mu>0], Imu[mu>0],law=law,fitmethod=fitmethod) 816 grid_coeffs.append(coeffs) 817 grid_pars.append([teff_,logg_,ebv_,vrad_,z_]) 818 Imu1s.append([Imu1,res,dflux]) 819 #- wrap up results in nice arrays 820 grid_pars = np.array(grid_pars) 821 grid_coeffs = np.array(grid_coeffs) 822 Imu1s = np.array(Imu1s) 823 824 return grid_pars, grid_coeffs, Imu1s

825

826 -def generate_grid(photbands,vrads=[0],ebvs=[0],zs=[0], 827 law='claret',fitmethod='equidist_r_leastsq',outfile='mygrid.fits',**kwargs):

828 829 if os.path.isfile(outfile): 830 hdulist = pf.open(outfile,mode='update') 831 existing_bands = [ext.header['extname'] for ext in hdulist[1:]] 832 else: 833 hdulist = pf.HDUList([]) 834 hdulist.append(pf.PrimaryHDU(np.array([[0,0]]))) 835 existing_bands = [] 836 837 hd = hdulist[0].header 838 hd.update('FIT', fitmethod, 'FIT ROUTINE') 839 hd.update('LAW', law, 'FITTED LD LAW') 840 hd.update('GRID', kwargs.get('grid',defaults['grid']), 'GRID') 841 842 for photband in photbands: 843 if photband in existing_bands: 844 logger.info('BAND {} already exists: skipping'.format(photband)) 845 continue 846 pars,coeffs,Imu1s = fit_law_to_grid(photband,vrads=vrads,ebvs=ebvs,zs=zs,**kwargs) 847 cols = [] 848 849 cols.append(pf.Column(name='Teff', format='E', array=pars[:,0])) 850 cols.append(pf.Column(name="logg", format='E', array=pars[:,1])) 851 cols.append(pf.Column(name="ebv" , format='E', array=pars[:,2])) 852 cols.append(pf.Column(name="vrad", format='E', array=pars[:,3])) 853 cols.append(pf.Column(name="z" , format='E', array=pars[:,4])) 854 for col in range(coeffs.shape[1]): 855 cols.append(pf.Column(name='a{:d}'.format(col+1), format='E', array=coeffs[:,col])) 856 cols.append(pf.Column(name='Imu1', format='E', array=Imu1s[:,0])) 857 cols.append(pf.Column(name='SRS', format='E', array=Imu1s[:,1])) 858 cols.append(pf.Column(name='dint', format='E', array=Imu1s[:,2])) 859 860 newtable = pf.new_table(pf.ColDefs(cols)) 861 newtable.header.update('EXTNAME', photband, "SYSTEM.FILTER") 862 newtable.header.update('SYSTEM', photband.split('.')[0], 'PASSBAND SYSTEM') 863 newtable.header.update('FILTER', photband.split('.')[1], 'PASSBAND FILTER') 864 hdulist.append(newtable) 865 866 if os.path.isfile(outfile): 867 hdulist.close() 868 else: 869 hdulist.writeto(outfile)

870 #}

871 872 873 #{ Internal functions 874 875 -def _r(mu):

876 """ 877 Convert mu to r coordinates 878 """ 879 return np.sqrt(1.-mu**2.)

880

881 -def _mu(r_):

882 """ 883 Convert r to mu coordinates 884 """ 885 return np.sqrt(1.-r_**2.)

886

887 -def ldres_fmin(coeffs, mu, Imu, law):

888 """ 889 Residual function for Nelder-Mead optimizer. 890 """ 891 return sum((Imu - globals()['ld_%s'%(law)](mu,coeffs))**2)

892

893 -def ldres_leastsq(coeffs, mu, Imu, law):

894 """ 895 Residual function for Levenberg-Marquardt optimizer. 896 """ 897 return Imu - globals()['ld_%s'%(law)](mu,coeffs)

898

899 -def _I_ls(ll,ss):

900 """ 901 Limb darkening moments (Dziembowski 1977, Townsend 2002) 902 903 Recursive implementation. 904 905 >>> _I_ls(0,0.5) 906 0.6666666666666666 907 >>> _I_ls(1,0.5) 908 0.4 909 >>> _I_ls(2,0.5) 910 0.09523809523809523 911 """ 912 if ll==0: 913 return 1./(1.+ss) 914 elif ll==1: 915 return 1./(2.+ss) 916 else: 917 return (ss-ll+2)/(ss+ll-2+3.)*_I_ls(ll-2,ss)

918

919 @memoized 920 -def get_ld_grid(photband,**kwargs):

921 """ 922 Retrieve an interpolating grid for the LD coefficients 923 924 Check outcome: 925 926 #>>> bands = ['GENEVA.U', 'GENEVA.B', 'GENEVA.G', 'GENEVA.V'] 927 #>>> f_ld_grid = get_ld_grid(bands) 928 #>>> ff = pf.open(_atmos['file']) 929 #>>> all(ff['GENEVA.U'].data[257][2:]==f_ld_grid(ff['GENEVA.U'].data[257][0],ff['GENEVA.U'].data[257][1])[0:5]) 930 #True 931 #>>> all(ff['GENEVA.G'].data[257][2:]==f_ld_grid(ff['GENEVA.G'].data[257][0],ff['GENEVA.G'].data[257][1])[10:15]) 932 #True 933 #>>> ff.close() 934 935 #Make some plots: 936 937 #>>> photband = ['GENEVA.V'] 938 #>>> f_ld = get_ld_grid(photband) 939 #>>> logg = 4.0 940 #>>> mu = linspace(0,1,100) 941 #>>> p = figure() 942 #>>> p = gcf().canvas.set_window_title('test of function <get_ld_grid>') 943 #>>> for teff in linspace(9000,12000,19): 944 #... out = f_ld(teff,logg) 945 #... a1x,a2x,a3x,a4x, I_x1 = out.reshape((len(photband),5)).T 946 #... p = subplot(221);p = title('Interpolation of absolute intensities') 947 #... p = plot(teff,I_x1,'ko') 948 #... p = subplot(222);p = title('Interpolation of LD coefficients') 949 #... p = scatter(4*[teff],[a1x,a2x,a3x,a4x],c=range(4),vmin=0,vmax=3,cmap=cm.spectral,edgecolors='none') 950 #... p = subplot(223);p = title('Without absolute intensity') 951 #... p = plot(mu,ld_eval(mu,[a1x,a2x,a3x,a4x]),'-') 952 #... p = subplot(224);p = title('With absolute intensity') 953 #... p = plot(mu,I_x1*ld_eval(mu,[a1x,a2x,a3x,a4x]),'-') 954 955 """ 956 #-- retrieve the grid points (unique values) 957 teffs,loggs = get_ld_grid_dimensions(**kwargs) 958 teffs_grid = np.sort(np.unique1d(teffs)) 959 loggs_grid = np.sort(np.unique1d(loggs)) 960 coeff_grid = np.zeros((len(teffs_grid),len(loggs_grid),5*len(photband))) 961 962 #-- get the FITS-file containing the tables 963 gridfile = get_file(**kwargs) 964 #-- fill the grid 965 ff = pf.open(gridfile) 966 for pp,iband in enumerate(photband): 967 teffs = ff[iband].data.field('Teff') 968 loggs = ff[iband].data.field('logg') 969 for ii,(iteff,ilogg) in enumerate(zip(teffs,loggs)): 970 indext = np.searchsorted(teffs_grid,iteff) 971 indexg = np.searchsorted(loggs_grid,ilogg) 972 #-- array and list are added for backwards compatibility with some 973 # pyfits versions 974 coeff_grid[indext,indexg,5*pp:5*(pp+1)] = np.array(list(ff[iband].data[ii]))[2:] 975 ff.close() 976 #-- make an interpolating function 977 f_ld_grid = InterpolatingFunction([teffs_grid,loggs_grid],coeff_grid) 978 return f_ld_grid

979

980 981 @memoized 982 -def _get_itable_markers(photband,gridfile, 983 teffrange=(-np.inf,np.inf),loggrange=(-np.inf,np.inf)):

984 """ 985 Get a list of markers to more easily retrieve integrated fluxes. 986 """ 987 clear_memoization() 988 989 ff = pf.open(gridfile) 990 ext = ff[photband] 991 columns = ext.columns.names 992 names = ['teff','logg'] 993 grid = [ext.data.field('teff'), 994 ext.data.field('logg')] 995 if 'ebv' in columns: 996 names.append('ebv') 997 grid.append(ext.data.field('ebv')) 998 if 'vrad' in columns: 999 names.append('vrad') 1000 grid.append(ext.data.field('vrad')) 1001 1002 grid_axes = [np.sort(list(set(i))) for i in grid] 1003 1004 #-- we construct an array representing the teff-logg-ebv content, but 1005 # in one number: 50000400 means: 1006 # T=50000,logg=4.0 1007 N = len(grid[0]) 1008 markers = np.zeros(N,float) 1009 gridpnts = np.zeros((N,len(grid))) 1010 pars = np.zeros((N,5)) 1011 1012 for i,entries in enumerate(zip(*grid)): 1013 markers[i] = encode(**{key:entry for key,entry in zip(names,entries)}) 1014 gridpnts[i]= entries 1015 pars[i] = list(ext.data[i][-5:]) 1016 ff.close() 1017 sa = np.argsort(markers) 1018 print 'read in gridfile',gridfile 1019 pars[:,-1] = np.log10(pars[:,-1]) 1020 return markers[sa],grid_axes,gridpnts[sa],pars[sa]

1021 1022 #} 1023 1024 1025 if __name__=="__main__": 1026 import sys 1027 if not sys.argv[1:]: 1028 import pylab as pl 1029 import doctest 1030 doctest.testmod() 1031 pl.show() 1032 else: 1033 #from ivs.aux import argkwargparser 1034 #method,args,kwargs = argkwargparser.parse() 1035 #print("Calling {} with args=({}) and kwargs=({})".format(method,args,kwargs) 1036 photbands = ['JOHNSON.U','JOHNSON.B','JOHNSON.V','KEPLER.V','COROT.SIS','COROT.EXO'] 1037 photbands+= ['2MASS.J','2MASS.H','2MASS.KS'] 1038 #generate_grid(photbands,vrads=np.arange(-500,501,50),ebvs=np.arange(0,2.005,0.01),law='claret',outfile='claret.fits') 1039 #generate_grid(photbands,vrads=np.arange(-500,501,50),ebvs=np.arange(0,2.005,0.01),law='linear',outfile='linear.fits') 1040 generate_grid(['MOST.V','IRAC.36','COROT.EXO'],vrads=[0],ebvs=[0.06],zs=[0,0], 1041 law='claret',fitmethod='equidist_r_leastsq',outfile='HD261903.fits') 1042

Source Code for Module ivs.sed.limbdark