ivs.spectra.model

1 # -*- coding: utf-8 -*- 2 """ 3 Interface to stellar spectra library and functions to manipulate them. 4 5 Make a plot of the domains of all spectral grids. First collect all the grid 6 names 7 8 >>> grids = get_gridnames() 9 10 Prepare the plot 11 12 >>> p = pl.figure() 13 >>> color_cycle = [pl.cm.spectral(j) for j in np.linspace(0, 1.0, len(grids))] 14 >>> p = pl.gca().set_color_cycle(color_cycle) 15 16 And plot all the grid points. We have to set some custom default values for 17 some grids. 18 19 >>> for grid in grids: 20 ... vturb = 'ostar' in grid and 10 or 2 21 ... t = 0. 22 ... if 'marcs' in grid: t = 1. 23 ... teffs,loggs = get_grid_dimensions(grid=grid,vturb=vturb,t=t) 24 ... p = pl.plot(np.log10(teffs),loggs,'o',ms=7,label=grid) 25 26 Now take care of the plot details 27 28 >>> p = pl.xlim(pl.xlim()[::-1]) 29 >>> p = pl.ylim(pl.ylim()[::-1]) 30 >>> p = pl.xlabel('Effective temperature [K]') 31 >>> p = pl.ylabel('log( Surface gravity [cm s$^{-1}$]) [dex]') 32 >>> xticks = [3000,5000,7000,10000,15000,25000,35000,50000,65000] 33 >>> p = pl.xticks([np.log10(i) for i in xticks],['%d'%(i) for i in xticks]) 34 >>> p = pl.legend(loc='upper left',prop=dict(size='small')) 35 >>> p = pl.grid() 36 37 ]include figure]]ivs_spectra_model_grid.png] 38 39 """ 40 import os 41 import logging 42 import copy 43 import astropy.io.fits as pf 44 import inspect 45 from ivs import config 46 from ivs.units import constants 47 from ivs.units import conversions 48 from ivs.aux.decorators import memoized 49 from ivs.spectra import tools 50 from ivs.aux import loggers 51 52 import numpy as np 53 # from Scientific.Functions.Interpolation import InterpolatingFunction 54 from scipy.interpolate import RegularGridInterpolator as InterpolatingFunction 55 56 logger = logging.getLogger("SPEC.MODEL") 57 logger.addHandler(loggers.NullHandler) 58 59 defaults = dict(grid='atlas',z=+0.0,vturb=2,band='vis', 60 t=0.0,a=0.0,c=0.5,atm='p') # MARCS and COMARCS 61 basedir = 'spectables' 62 63 #{ Interface to library 64

65 -def set_defaults(**kwargs):

66 """ 67 Set defaults of this module. 68 69 If you give no keyword arguments, the default values will be reset. 70 """ 71 if not kwargs: 72 kwargs = dict(grid='atlas',z=+0.0,vturb=2,band='vis', 73 t=0.0,a=0.0,c=0.5,atm='p') # MARCS and COMARCS 74 75 for key in kwargs: 76 if key in defaults: 77 defaults[key] = kwargs[key]

78 79

80 -def get_gridnames():

81 """ 82 Return a list of available grid names. 83 84 @return: list of grid names 85 @rtype: list of str 86 """ 87 return ['cmfgen','ostar2002','bstar2006','atlas','marcs', 'heberb', 'hebersdb','tmapsdb']

88 89 90

91 -def get_file(**kwargs):

92 """ 93 Retrieve the filename containing the specified SED grid. 94 95 Available grids and example keywords: 96 - grid='fastwind': no options 97 - grid='cmfgen': no options 98 - grid='marcs': options c, atm (p or s) 99 - grid='ostar2002': options z,v 100 - grid='bstar2006a': options z,v 101 - grid='bstar2006b': options z,v 102 - grid='bstar2006': options z,v,a 103 - grid='atlas': options z 104 - grid='heberb': no options 105 - grid='hebersdb': no options 106 - grid='tmapsdb': no options 107 108 Details for grid 'bstar2006': 109 - metallicity in Z/Z0 with Z0 solar. z=0,0.001,0.01,0.033,0.1,0.2,0.5,1.,2 110 - microturbulent velocity: v=2 or v=10 km/s 111 - abundance: a='' or a='CN'. In the latter, the Helium abundance is 112 increased to He/H=0.2, the nitrogen abundance is increased 113 by a factor of 5, and the carbon abundance is halved (CNO cycle processed 114 material brought to the stellar surface) 115 116 Details for grid 'heberb': LTE Grid computed for B-type stars by Uli Heber, 117 reff: Heber et al. 2000 118 119 Details for grid 'hebersdb': LTE Grid computed for sdB stars by Uli Heber, 120 reff: Heber et al. 2000 121 122 Details for grid 'tmapsdb': NLTE Grid computed for sdB stars using the TMAP 123 (TUEBINGEN NLTE MODEL ATMOSPHERE PACKAGE) code. reff: Werner K., et al. 2003 124 and Rauch T., Deetjen J.L. 2003 125 126 @param grid: gridname, or path to grid. 127 @type grid: string 128 """ 129 #-- possibly you give a filename 130 grid = kwargs.get('grid',defaults['grid'])#.lower() 131 if os.path.isfile(grid): 132 return grid 133 134 #-- general 135 z = kwargs.get('z',defaults['z']) 136 #-- only for Kurucz 137 vturb = int(kwargs.get('vturb',defaults['vturb'])) 138 #-- only for Marcs 139 t = kwargs.get('t',defaults['t']) 140 a = kwargs.get('a',defaults['a']) 141 c = kwargs.get('c',defaults['c']) 142 atm = kwargs.get('atm',defaults['atm']) 143 #-- only for TLUSTY 144 band = kwargs.get('band',defaults['band']) 145 146 #-- figure out what grid to use 147 if grid=='cmfgen': 148 basename = 'cmfgen_spectra.fits' 149 elif grid=='marcs': 150 basename = 'marcsp_%sz%.1ft%.1f_a%.2f_c%.2f_spectra.fits'%(atm,z,t,a,c) 151 elif grid=='ostar2002': 152 basename = 'OSTAR2002_z%.3fv%d_%s_spectra.fits'%(z,vturb,band) 153 elif grid=='bstar2006': 154 basename = 'BSTAR2006_z%.3fv%d_%s_spectra.fits'%(z,vturb,band) 155 elif grid=='atlas': 156 basename = 'ATLASp_z%.1ft%.1f_a%.2f_spectra.fits'%(z,t,a) 157 elif grid=='heberb': 158 basename = 'Heber2000_B_h909.fits' 159 elif grid=='hebersdb': 160 basename = 'Heber2000_sdB_h909.fits' 161 elif grid=='tmapsdb': 162 basename = 'TMAP2011_sdB.fits' 163 else: 164 raise ValueError, "grid %s does not exist"%(grid) 165 166 #-- retrieve the absolute path of the file and check if it exists: 167 grid = config.get_datafile(basedir,basename) 168 169 return grid

170 171

172 -def get_table(teff=None,logg=None,vrad=0,vrot=0,**kwargs):

173 """ 174 Retrieve synthetic spectrum. 175 176 Wavelengths in angstrom, fluxes in eddington flux: erg/s/cm2/A. 177 178 It is possible to rotationally broaden the spectrum by supplying at least 179 'vrot' and optionally also other arguments for the rotation.rotin function. 180 Supply vrot in km/s 181 182 It is possibility to include a radial velocity shift in the spectrum. Supply 183 'vrad' in km/s. (+vrad means redshift). Uses spectra.tools.doppler_shift 184 """ 185 gridfile = get_file(**kwargs) 186 template_wave = kwargs.pop('wave',None) 187 188 ff = pf.open(gridfile) 189 190 try: 191 #-- extenstion name as in fits files prepared by Steven 192 mod_name = "T%05d_logg%01.02f" %(teff,logg) 193 mod = ff[mod_name] 194 wave = mod.data.field('wavelength') 195 flux = mod.data.field('flux') 196 cont = mod.data.field('cont') 197 logger.debug('Model spectrum taken directly from file (%s)'%(os.path.basename(gridfile))) 198 if template_wave is not None: 199 flux = np.interp(template_wave,wave,flux) 200 cont = np.interp(template_wave,wave,cont) 201 wave = template_wave 202 #-- if the teff/logg is not present, use the interpolation thing 203 except KeyError: 204 #-- it is possible we first have to set the interpolation function. 205 # This function is memoized, so if it will not be calculated 206 # twice. 207 meshkwargs = copy.copy(kwargs) 208 meshkwargs['wave'] = kwargs.get('wave',None) 209 meshkwargs['teffrange'] = kwargs.get('teffrange',None) 210 meshkwargs['loggrange'] = kwargs.get('loggrange',None) 211 wave,teffs,loggs,flux,flux_grid,cont_grid = get_grid_mesh(wave=template_wave,**kwargs) 212 logger.debug('Model spectrum interpolated from grid %s (%s)'%(os.path.basename(gridfile),meshkwargs)) 213 wave = wave + 0. 214 try: 215 flux = flux_grid(np.log10(teff),logg) + 0. 216 cont = cont_grid(np.log10(teff),logg) + 0. 217 except ValueError: 218 teffs,loggs = get_grid_dimensions(**kwargs) 219 index = np.argmin(np.abs( (teffs-teff)**2 + (loggs-logg)**2 )) 220 #logger.error('teff=%f-->%f, logg=%f-->%f'%(teff,teffs[index],logg,loggs[index])) 221 flux = flux_grid(np.log10(teffs[index]),loggs[index]) + 0. 222 cont = cont_grid(np.log10(teffs[index]),loggs[index]) + 0. 223 224 #-- convert to arrays 225 wave = np.array(wave,float) 226 flux = np.array(flux,float) 227 228 if vrot>0: 229 #-- calculate rotational broadening but reinterpolate on original 230 # wavelength grid. First we need to check which arguments we can pass 231 # through 232 argspec = inspect.getargspec(tools.rotational_broadening) 233 mykwargs = dict(zip(argspec.args[-len(argspec.defaults):],argspec.defaults)) 234 for key in kwargs: 235 if key in mykwargs: 236 mykwargs[key] = kwargs[key] 237 wave_,flux_ = tools.rotational_broadening(wave,flux,vrot,**mykwargs) 238 flux = np.interp(wave,wave_,flux_) 239 if vrad!=0: 240 wave_ = tools.doppler_shift(wave,vrad) 241 flux = np.interp(wave,wave_,flux) 242 243 ff.close() 244 return wave,flux,cont

245 246 247 248 249 250 251 252 253

254 -def get_grid_dimensions(**kwargs):

255 """ 256 Retrieve possible effective temperatures and gravities from a grid. 257 258 @rtype: (ndarray,ndarray) 259 @return: effective temperatures, gravities 260 """ 261 gridfile = get_file(**kwargs) 262 ff = pf.open(gridfile) 263 teffs = [] 264 loggs = [] 265 for mod in ff[1:]: 266 teffs.append(float(mod.header['TEFF'])) 267 loggs.append(float(mod.header['LOGG'])) 268 ff.close() 269 return np.array(teffs),np.array(loggs)

270 271 272 273 274 275 276 #@memoized

277 -def get_grid_mesh(wave=None,teffrange=None,loggrange=None,**kwargs):

278 """ 279 Return InterpolatingFunction spanning the available grid of spectrum models. 280 281 WARNING: the grid must be entirely defined on a mesh grid, but it does not 282 need to be equidistant. 283 284 It is thus the user's responsibility to know whether the grid is evenly 285 spaced in logg and teff 286 287 You can supply your own wavelength range, since the grid models' 288 resolution are not necessarily homogeneous. If not, the first wavelength 289 array found in the grid will be used as a template. 290 291 It might take a long a time and cost a lot of memory if you load the entire 292 grid. Therefore, you can also set range of temperature and gravity. 293 294 @param wave: wavelength to define the grid on 295 @type wave: ndarray 296 @param teffrange: starting and ending of the grid in teff 297 @type teffrange: tuple of floats 298 @param loggrange: starting and ending of the grid in logg 299 @type loggrange: tuple of floats 300 @return: wavelengths, teffs, loggs and fluxes of grid, and the interpolating 301 function 302 @rtype: (1Darray,1Darray,1Darray,3Darray,InterpolatingFunction) 303 """ 304 #-- get the dimensions of the grid 305 teffs,loggs = get_grid_dimensions(**kwargs) 306 #-- build flux grid, assuming a perfectly sampled grid (needs not to be 307 # equidistant) 308 if teffrange is not None: 309 sa = (teffrange[0]<=teffs) & (teffs<=teffrange[1]) 310 teffs = teffs[sa] 311 if loggrange is not None: 312 sa = (loggrange[0]<=loggs) & (loggs<=loggrange[1]) 313 loggs = loggs[sa] 314 #-- clip if necessary 315 teffs = list(set(list(teffs))) 316 loggs = list(set(list(loggs))) 317 teffs = np.sort(teffs) 318 loggs = np.sort(loggs) 319 if wave is not None: 320 flux = np.ones((len(teffs),len(loggs),len(wave))) 321 cont = np.ones((len(teffs),len(loggs),len(wave))) 322 #-- run over teff and logg, and interpolate the models onto the supplied 323 # wavelength range 324 gridfile = get_file(**kwargs) 325 ff = pf.open(gridfile) 326 for i,teff in enumerate(teffs): 327 for j,logg in enumerate(loggs): 328 try: 329 mod_name = "T%05d_logg%01.02f" %(teff,logg) 330 mod = ff[mod_name] 331 wave_ = mod.data.field('wavelength')#array(mod.data.tolist())[:,0] 332 flux_ = mod.data.field('flux')#array(mod.data.tolist())[:,1] 333 cont_ = mod.data.field('cont')#array(mod.data.tolist())[:,1] 334 #-- if there is no wavelength range given, we assume that 335 # the whole grid has the same resolution, and the first 336 # wave-array will be used as a template 337 if wave is None: 338 wave = wave_ 339 flux = np.ones((len(teffs),len(loggs),len(wave))) 340 cont = np.ones((len(teffs),len(loggs),len(wave))) 341 except KeyError: 342 continue 343 #-- it could be that we're lucky and the grid is completely 344 # homogeneous. In that case, there is no need for interpolation 345 try: 346 flux[i,j,:] = flux_ 347 cont[i,j,:] = cont_ 348 except: 349 flux[i,j,:] = np.interp(wave,wave_,flux_) 350 cont[i,j,:] = np.interp(wave,wave_,cont_) 351 # flux_grid = InterpolatingFunction([np.log10(teffs),loggs],flux) 352 # cont_grid = InterpolatingFunction([np.log10(teffs),loggs],cont) 353 flux_grid = InterpolatingFunction((np.log10(teffs),loggs),flux) 354 cont_grid = InterpolatingFunction((np.log10(teffs),loggs),cont) 355 #logger.info('Constructed spectrum interpolation grid') 356 return wave,teffs,loggs,flux,flux_grid,cont_grid

357 358 359 360 #} 361 362 if __name__=="__main__": 363 from ivs.aux import loggers 364 logger = loggers.get_basic_logger() 365 logger.setLevel(10) 366 #import doctest 367 import pylab as pl 368 import numpy as np 369 #doctest.testmod() 370 371 grids = get_gridnames() 372 373 p = pl.figure() 374 color_cycle = [pl.cm.spectral(j) for j in np.linspace(0, 1.0, len(grids))] 375 p = pl.gca().set_color_cycle(color_cycle) 376 377 for grid in grids: 378 vturb = 'ostar' in grid and 10 or 2 379 t = 0. 380 if 'marcs' in grid: t = 1. 381 teffs,loggs = get_grid_dimensions(grid=grid,vturb=vturb,t=t) 382 p = pl.plot(np.log10(teffs),loggs,'o',ms=7,label=grid) 383 384 p = pl.xlim(pl.xlim()[::-1]) 385 p = pl.ylim(pl.ylim()[::-1]) 386 p = pl.xlabel('Effective temperature [K]') 387 p = pl.ylabel('log( Surface gravity [cm s$^{-1}$]) [dex]') 388 xticks = [3000,5000,7000,10000,15000,25000,35000,50000,65000] 389 p = pl.xticks([np.log10(i) for i in xticks],['%d'%(i) for i in xticks]) 390 p = pl.legend(loc='upper left',prop=dict(size='small')) 391 p = pl.grid() 392 393 pl.show() 394

Source Code for Module ivs.spectra.model