1
2 """
3 Interface to the SED library.
4
5 The most basic usage of this module is:
6
7 >>> wave,flux = get_table(teff=10000,logg=4.0)
8
9 This will retrieve the model SED with the specified B{effective temperature} and
10 B{logg}, from the standard B{grid}, in standard B{units} and with zero
11 B{reddening}. All these things can be specified though (see below).
12
13 Section 1. Available model grids
14 ================================
15
16 Section 1.1 Available grids
17 ---------------------------
18
19 - kurucz: The Kurucz model grids, (default setting) reference: Kurucz 1993, yCat, 6039, 0
20 - metallicity (z): m01 is -0.1 log metal abundance relative to solar (solar abundances from Anders and Grevesse 1989)
21 - metallicity (z): p01 is +0.1 log metal abundance relative to solar (solar abundances from Anders and Grevesse 1989)
22 - alpha enhancement (alpha): True means alpha enhanced (+0.4)
23 - turbulent velocity (vturb): vturb in km/s
24 - nover= True means no overshoot
25 - odfnew=True means no overshoot but with better opacities and abundances
26
27 - tmap: NLTE grids computed for sdB stars with the Tubingen NLTE Model
28 Atmosphere package. No further parameters are available. Reference:
29 Werner et al. 2003,
30
31
32 Section 1.2 Plotting the domains of all spectral grids
33 ------------------------------------------------------
34
35 We make a plot of the domains of all spectral grids. Therefore, we first collect
36 all the grid names
37
38 >>> grids = get_gridnames()
39
40 Preparation of the plot: set the color cycle of the current axes to the spectral
41 color cycle.
42
43 >>> p = pl.figure(figsize=(10,8))
44 >>> color_cycle = [pl.cm.spectral(j) for j in np.linspace(0, 1.0, len(grids))]
45 >>> p = pl.gca().set_color_cycle(color_cycle)
46
47 To plot all the grid points, we run over all grid names (which are strings), and
48 retrieve their dimensions. The dimensions are just two arrays giving the teff-
49 and logg-coordinate of each SED in the grid. They can thus be easily plot:
50
51 >>> for grid in grids:
52 ... teffs,loggs = get_grid_dimensions(grid=grid)
53 ... p = pl.plot(np.log10(teffs),loggs,'o',ms=7,label=grid)
54
55 And we need to set some of the plotting details to make it look nicer.
56
57 >>> p = pl.xlim(pl.xlim()[::-1])
58 >>> p = pl.ylim(pl.ylim()[::-1])
59 >>> p = pl.xlabel('Effective temperature [K]')
60 >>> p = pl.ylabel('log( Surface gravity [cm s$^{-1}$]) [dex]')
61 >>> xticks = [3000,5000,7000,10000,15000,25000,35000,50000,65000]
62 >>> p = pl.xticks([np.log10(i) for i in xticks],['%d'%(i) for i in xticks])
63 >>> p = pl.legend(loc='upper left',prop=dict(size='small'))
64 >>> p = pl.grid()
65
66 ]include figure]]ivs_sed_model_grid.png]
67
68 Section 2. Retrieval of model SEDs
69 ==================================
70
71 Subsection 2.1 Default settings
72 -------------------------------
73
74 To get information on the grid that is currently defined, you can type the
75 following. Note that not all parameters are relevant for all grids, e.g. the
76 convection theory parameter C{ct} has no influence when the Kurucz grid is
77 chosen.
78
79 >>> print defaults
80 {'use_scratch': False, 'Rv': 3.1, 'co': 1.05, 'c': 0.5, 'grid': 'kurucz', 'alpha': False, 'odfnew': True, 'ct': 'mlt', 'a': 0.0, 'vturb': 2, 'law': 'fitzpatrick2004', 'm': 1.0, 't': 1.0, 'z': 0.0, 'nover': False, 'He': 97}
81
82 or
83
84 >>> print os.path.basename(get_file())
85 kurucz93_z0.0_k2odfnew_sed.fits
86
87 You can change the defaults with the function L{set_defaults}:
88
89 >>> set_defaults(z=0.5)
90 >>> print defaults
91 {'use_scratch': False, 'Rv': 3.1, 'co': 1.05, 'c': 0.5, 'grid': 'kurucz', 'alpha': False, 'odfnew': True, 'ct': 'mlt', 'a': 0.0, 'vturb': 2, 'law': 'fitzpatrick2004', 'm': 1.0, 't': 1.0, 'z': 0.5, 'nover': False, 'He': 97}
92
93 And reset the 'default' default values by calling L{set_defaults} without arguments
94
95 >>> set_defaults()
96 >>> print defaults
97 {'use_scratch': False, 'Rv': 3.1, 'co': 1.05, 'c': 0.5, 'grid': 'kurucz', 'alpha': False, 'odfnew': True, 'ct': 'mlt', 'a': 0.0, 'vturb': 2, 'law': 'fitzpatrick2004', 'm': 1.0, 't': 1.0, 'z': 0.0, 'nover': False, 'He': 97}
98
99 Subsection 2.2 Speeding up
100 --------------------------
101
102 When fitting an sed using the builder class, or repeatedly reading model seds,
103 or integrated photometry, the main bottleneck on the speed will be the disk access
104 This can be circumvented by using the scratch disk. To do this, call the function
105 copy2scratch() after setting the default settings as explained above. f.x.:
106
107 >>> set_defaults(grid='kurucz', z=0.5)
108 >>> copy2scratch()
109
110 You have to do this every time you change a grid setting. This function creates a
111 directory named 'your_username' on the scratch disk and works from there. So you
112 won`t disturbed other users.
113
114 After the fitting process use the function
115
116 >>> clean_scratch()
117
118 to remove the models that you used from the scratch disk. Be carefull with this,
119 because it will remove the models without checking if there is another process
120 using them. So if you have multiple scripts running that are using the same models,
121 only clean the scratch disk after the last process is finished.
122
123 The gain in speed can be up to 70% in single sed fitting, and up to 40% in binary
124 and multiple sed fitting.
125
126 For the sake of the examples, we'll set the defaults back to z=0.0:
127
128 >>> set_defaults()
129
130 Subsection 2.3 Model SEDs
131 -------------------------
132
133 Be careful when you supply parameters: e.g., not all grids are calculated for
134 the same range of metallicities. In L{get_table}, only the effective temperature
135 and logg are 'interpolatable' quantities. You have to set the metallicity to a
136 grid value. The reddening value can take any value: it is not interpolated but
137 calculated. You can thus also specify the type of reddening law (see L{reddening}).
138
139 >>> wave,flux = get_table(teff=12345,logg=4.321,ebv=0.12345,z=0.5)
140
141 but
142
143 >>> try:
144 ... wave,flux = get_table(teff=12345,logg=4.321,ebv=0.12345,z=0.6)
145 ... except IOError,msg:
146 ... print msg
147 File sedtables/modelgrids/kurucz93_z0.6_k2odfnew_sed.fits not found in any of the specified data directories /STER/pieterd/IVSDATA/, /STER/kristofs/IVSdata
148
149 Since the Kurucz model atmospheres have not been calculated for the value of
150 C{z=0.6}.
151
152 Instead of changing the defaults of this module with L{set_defaults}, you can
153 also give extra arguments to L{get_table} to specify the grid you want to use.
154 The default settings will not change in this case.
155
156 >>> wave,flux = get_table(teff=16321,logg=4.321,ebv=0.12345,z=0.3,grid='tlusty')
157
158 The default B{units} of the SEDs are angstrom and erg/s/cm2/AA/sr. To change them,
159 do:
160
161 >>> wave,flux = get_table(teff=16321,logg=4.321,wave_units='micron',flux_units='Jy/sr')
162
163 To B{remove the steradian} from the units when you know the angular diameter of
164 your star in milliarcseconds, you can do (we have to convert diameter to surface):
165
166 >>> ang_diam = 3.21 # mas
167 >>> scale = conversions.convert('mas','sr',ang_diam/2.)
168 >>> wave,flux = get_table(teff=9602,logg=4.1,ebv=0.0,z=0.0,grid='kurucz')
169 >>> flux *= scale
170
171 The example above is representative for the case of Vega. So, if we now calculate
172 the B{synthetic flux} in the GENEVA.V band, we should end up with the zeropoint
173 magnitude of this band, which is close to zero:
174
175 >>> flam = synthetic_flux(wave,flux,photbands=['GENEVA.V'])
176 >>> print '%.3f'%(conversions.convert('erg/s/cm2/AA','mag',flam,photband='GENEVA.V')[0])
177 0.063
178
179 Compare this with the calibrated value
180
181 >>> print filters.get_info(['GENEVA.V'])['vegamag'][0]
182 0.061
183
184 Section 3. Retrieval of integrated photometry
185 =============================================
186
187 Instead of retrieving a model SED, you can immediately retrieve pre-calculated
188 integrated photometry. The benefit of this approach is that it is B{much} faster
189 than retrieving the model SED and then calculating the synthetic flux. Also,
190 you can supply arbitrary metallicities within the grid boundaries, as interpolation
191 is done in effective temperature, surface gravity, reddening B{and} metallicity.
192
193 Note that also the B{reddening law is fixed} now, you need to recalculate the
194 tables for different parameters if you need them.
195
196 The B{massive speed-up} is accomplished the following way: it may take a few tens
197 of seconds to retrieve the first pre-integrated SED, because all available
198 files from the specified grid will be loaded into memory, and a `markerarray'
199 will be made allowing a binary search in the grid. This makes it easy to retrieve
200 all models around the speficied point in N-dimensional space. Next, a linear
201 interpolation method is applied to predict the photometric values of the
202 specified point.
203
204 All defaults set for the retrieval of model SEDs are applicable for the integrated
205 photometry tables as well.
206
207 When retrieving integrated photometry, you also get the B{absolute luminosity}
208 (integration of total SED) as a bonus. This is the absolute luminosity assuming
209 the star has a radius of 1Rsol. Multiply by Rstar**2 to get the true luminosity.
210
211 Because photometric filters cannot trivially be assigned a B{wavelength} to (see
212 L{filters.eff_wave}), by default, no wavelength information is retrieved. If you
213 want to retrieve the effective wavelengths of the filters themselves (not taking
214 into account the model atmospheres), you can give an extra keyword argument
215 C{wave_units}. If you want to take into account the model atmosphere, use
216 L{filters.eff_wave}.
217
218 >>> photbands = ['GENEVA.U','2MASS.J']
219 >>> fluxes,Labs = get_itable(teff=16321,logg=4.321,ebv=0.12345,z=0.123,photbands=photbands)
220 >>> waves,fluxes,Labs = get_itable(teff=16321,logg=4.321,ebv=0.12345,z=0.123,photbands=photbands,wave_units='AA')
221
222 Note that the integration only gives you fluxes, and is thus independent from
223 the zeropoints of the filters (but dependent on the transmission curves). To
224 get the synthetic magnitudes, you can do
225
226 >>> mymags = [conversions.convert('erg/s/cm2/AA','mag',fluxes[i],photband=photbands[i]) for i in range(len(photbands))]
227
228 The mags don't mean anything in this case because they have not been corrected
229 for the distance to the star.
230
231 The retrieval of integrated photometry can go much faster if you want to do
232 it for a whole set of parameters. The L{get_itable_pix} function has a much
233 more flexible, reliable and fast interpolation scheme. It is possible to
234 interpolate also over doppler shift and interstellar Rv, as long as the grids
235 have been computed before. See L{get_itable_pix} for more information.
236
237 Subsection 3. Full example
238 ==========================
239
240 We build an SED of Vega and compute synthetic magnitudes in the GENEVA and
241 2MASS bands.
242
243 These are the relevant parameters of Vega and photometric passbands
244
245 >>> ang_diam = 3.21 # mas
246 >>> teff = 9602
247 >>> logg = 4.1
248 >>> ebv = 0.0
249 >>> z = 0.0
250 >>> photbands = ['GENEVA.U','GENEVA.G','2MASS.J','2MASS.H','2MASS.KS']
251
252 We can compute (R/d) to scale the synthetic flux as
253
254 >>> scale = conversions.convert('mas','sr',ang_diam/2.)
255
256 We retrieve the SED
257
258 >>> wave,flux = get_table(teff=teff,logg=logg,ebv=ebv,z=z,grid='kurucz')
259 >>> flux *= scale
260
261 Then compute the synthetic fluxes, and compare them with the synthetic fluxes as
262 retrieved from the pre-calculated tables
263
264 >>> fluxes_calc = synthetic_flux(wave,flux,photbands)
265 >>> wave_int,fluxes_int,Labs = get_itable(teff=teff,logg=logg,ebv=ebv,z=z,photbands=photbands,wave_units='AA')
266 >>> fluxes_int *= scale
267
268 Convert to magnitudes:
269
270 >>> m1 = [conversions.convert('erg/s/cm2/AA','mag',fluxes_calc[i],photband=photbands[i]) for i in range(len(photbands))]
271 >>> m2 = [conversions.convert('erg/s/cm2/AA','mag',fluxes_int[i],photband=photbands[i]) for i in range(len(photbands))]
272
273 And make a nice plot
274
275 >>> p = pl.figure()
276 >>> p = pl.loglog(wave,flux,'k-',label='Kurucz model')
277 >>> p = pl.plot(wave_int,fluxes_calc,'ro',label='Calculated')
278 >>> p = pl.plot(wave_int,fluxes_int,'bx',ms=10,mew=2,label='Pre-calculated')
279 >>> p = [pl.annotate('%s: %.3f'%(b,m),(w,f),color='r') for b,m,w,f in zip(photbands,m1,wave_int,fluxes_calc)]
280 >>> p = [pl.annotate('%s: %.3f'%(b,m),(w-1000,0.8*f),color='b') for b,m,w,f in zip(photbands,m2,wave_int,fluxes_int)]
281 >>> p = pl.xlabel('Wavelength [Angstrom]')
282 >>> p = pl.ylabel('Flux [erg/s/cm2/AA]')
283
284 ]include figure]]ivs_sed_model_example.png]
285
286 """
287 import re
288 import os
289 import sys
290 import glob
291 import logging
292 import copy
293 import astropy.io.fits as pf
294 import time
295 import numpy as np
296 try:
297 from scipy.interpolate import LinearNDInterpolator
298 from scipy.interpolate import griddata
299 new_scipy = True
300 except ImportError:
301 from Scientific.Functions.Interpolation import InterpolatingFunction
302 new_scipy = False
303 from scipy.interpolate import interp1d
304 from multiprocessing import Process,Manager,cpu_count
305
306 from ivs import config
307 from ivs.units import conversions
308 from ivs.units import constants
309 from ivs.aux import loggers
310 from ivs.aux.decorators import memoized,clear_memoization
311 import itertools
312 import functools
313 from ivs.aux import numpy_ext
314 from ivs.sed import filters
315 from ivs.inout import ascii
316 from ivs.inout import fits
317 from ivs.sigproc import interpol
318 from ivs.catalogs import sesame
319 import reddening
320 import getpass
321 import shutil
322
323 logger = logging.getLogger("SED.MODEL")
324 logger.addHandler(loggers.NullHandler)
325
326 caldir = os.sep.join(['sedtables','calibrators'])
327
328
329 __defaults__ = dict(grid='kurucz',odfnew=True,z=+0.0,vturb=2,
330 alpha=False,nover=False,
331 He=97,
332 ct='mlt',
333 t=1.0,a=0.0,c=0.5,m=1.0,co=1.05,
334 Rv=3.1,law='fitzpatrick2004',
335 use_scratch=False)
336
337 defaults = __defaults__.copy()
338 defaults_multiple = [defaults.copy(),defaults.copy()]
339
340 basedir = 'sedtables/modelgrids/'
341 scratchdir = None
346 """
347 Set defaults of this module
348
349 If you give no keyword arguments, the default values will be reset.
350 """
351 clear_memoization(keys=['ivs.sed.model'])
352
353 if not kwargs:
354 kwargs = __defaults__.copy()
355
356 for key in kwargs:
357 if key in defaults:
358 defaults[key] = kwargs[key]
359 logger.info('Set %s to %s'%(key,kwargs[key]))
360
373
375 """
376 Copy the grids to the scratch directory to speed up the fitting process.
377 Files are placed in the directory: /scratch/uname/ where uname is your username.
378
379 This function checks the grids that are set with the functions set_defaults()
380 and set_defaults_multiple(). Every time a grid setting is changed, this
381 function needs to be called again.
382
383 Don`t forget to remove the files from the scratch directory after the fitting
384 process is completed with clean_scratch()
385
386 It is possible to give z='*' and Rv='*' as an option; when you do that, the grids
387 with all z, Rv values are copied. Don't forget to add that option to clean_scratch too!
388 """
389 global scratchdir
390 uname = getpass.getuser()
391 if not os.path.isdir('/scratch/%s/'%(uname)):
392 os.makedirs('/scratch/%s/'%(uname))
393 scratchdir = '/scratch/%s/'%(uname)
394
395
396 defaults_ = []
397 defaults_.append(defaults)
398 defaults_.extend(defaults_multiple)
399
400
401
402 for default in defaults_:
403 default['use_scratch'] = False
404
405
406 originalDefaults = {}
407 for key in kwargs:
408 if key in default:
409 originalDefaults[key] = default[key]
410 default[key] = kwargs[key]
411 logger.debug('Using provided value for {0:s}={1:s} when copying to scratch'.format(key,str(kwargs[key])))
412
413 fname = get_file(integrated=False,**default)
414
415 if isinstance(fname,str):
416 fname = [fname]
417 for ifname in fname:
418 if not os.path.isfile(scratchdir + os.path.basename(ifname)):
419 shutil.copy(ifname,scratchdir)
420 logger.info('Copied grid: %s to scratch'%(ifname))
421 else:
422 logger.info('Using existing grid: %s from scratch'%(os.path.basename(ifname)))
423
424 fname = get_file(integrated=True,**default)
425 if isinstance(fname,str):
426 fname = [fname]
427 for ifname in fname:
428 if not os.path.isfile(scratchdir + os.path.basename(ifname)):
429 shutil.copy(ifname,scratchdir)
430 logger.info('Copied grid: %s to scratch'%(ifname))
431 else:
432 logger.info('Using existing grid: %s from scratch'%(os.path.basename(ifname)))
433 default['use_scratch'] = True
434 for key in kwargs:
435 if key in default:
436 default[key] = originalDefaults[key]
437
439 """
440 Remove the grids that were copied to the scratch directory by using the
441 function copy2scratch(). Be carefull with this function, as it doesn't check
442 if the models are still in use. If you are running multiple scripts that
443 use the same models, only clean the scratch disk after the last script is
444 finnished.
445 """
446 defaults_ = []
447 defaults_.append(defaults)
448 defaults_.extend(defaults_multiple)
449
450 for default in defaults_:
451 if default['use_scratch']:
452 originalDefaults = {}
453 for key in kwargs:
454 if key in default:
455 originalDefaults[key] = default[key]
456 default[key] = kwargs[key]
457 logger.debug('Using provided value for {0:s}={1:s} when deleting from scratch'.format(key,str(kwargs[key])))
458
459
460
461
462 fname = get_file(integrated=False,**default)
463 if isinstance(fname,str):
464 fname = [fname]
465 for ifname in fname:
466 if os.path.isfile(ifname):
467 logger.info('Removed file: %s'%(ifname))
468 os.remove(ifname)
469
470 fname = get_file(integrated=True,**default)
471 if isinstance(fname,str):
472 fname = [fname]
473 for ifname in fname:
474 if os.path.isfile(ifname):
475 logger.info('Removed file: %s'%(ifname))
476 os.remove(ifname)
477 default['use_scratch'] = False
478 for key in kwargs:
479 if key in default:
480 default[key] = originalDefaults[key]
481
485 """
486 Convert the defaults to a string, e.g. for saving files.
487 """
488 return '_'.join([str(i)+str(defaults[i]) for i in sorted(defaults.keys())])
489
495
498 """
499 Return a list of available grid names.
500
501 If you specificy the grid's name, you get two lists: one with all available
502 original, non-integrated grids, and one with the pre-calculated photometry.
503
504 @parameter grid: name of the type of grid (optional)
505 @type grid: string
506 @return: list of grid names
507 @rtype: list of str
508 """
509 if grid is None:
510 return ['kurucz','fastwind','cmfgen','sdb_uli','wd_boris','wd_da','wd_db',
511 'tlusty','uvblue','atlas12','nemo','tkachenko','marcs','marcs2','tmap',
512 ]
513
514 else:
515 files = config.glob(basedir,'*%s*.fits'%(grid))
516 integrated = [os.path.basename(ff) for ff in files if os.path.basename(ff)[0]=='i']
517 original = [os.path.basename(ff) for ff in files if os.path.basename(ff)[0]!='i']
518 return original,integrated
519
520
521
522 -def get_file(integrated=False,**kwargs):
523 """
524 Retrieve the filename containing the specified SED grid.
525
526 The keyword arguments are specific to the kind of grid you're using.
527
528 Basic keywords are 'grid' for the name of the grid, and 'z' for metallicity.
529 For other keywords, see the source code.
530
531 Available grids and example keywords:
532 - grid='kurucz93':
533 - metallicity (z): m01 is -0.1 log metal abundance relative to solar (solar abundances from Anders and Grevesse 1989)
534 - metallicity (z): p01 is +0.1 log metal abundance relative to solar (solar abundances from Anders and Grevesse 1989)
535 - alpha enhancement (alpha): True means alpha enhanced (+0.4)
536 - turbulent velocity (vturb): vturb in km/s
537 - nover= True means no overshoot
538 - odfnew=True means no overshoot but with better opacities and abundances
539 - grid='tlusty':
540 - z: log10(Z/Z0)
541 - grid='sdb_uli': metallicity and helium fraction (z, he=98)
542 - grid='fastwind': no options
543 - grid='wd_boris': no options
544 - grid='stars': precomputed stars (vega, betelgeuse...)
545 - grid='uvblue'
546 - grid='marcs'
547 - grid='marcs2'
548 - grid='atlas12'
549 - grid='tkachenko': metallicity z
550 - grid='nemo': convection theory and metallicity (CM=Canuto and Mazzitelli 1991),
551 (CGM=Canuto,Goldman,Mazzitelli 1996), (MLT=mixinglengththeory a=0.5)
552 - grid='marcsjorissensp': high resolution spectra from 4000 to 25000 A of (online available) MARCS grid computed by A. Jorissen
553 with turbospectrum v12.1.1 in late 2012, then converted to the Kurucz wavelength grid (by S. Bloemen and M. Hillen).
554
555 @param integrated: choose integrated version of the gridcopy2scratch
556 @type integrated: boolean
557 @keyword grid: gridname (default Kurucz)
558 @type grid: str
559 @return: gridfile
560 @rtype: str
561 """
562
563
564 grid = kwargs.get('grid',defaults['grid'])
565 use_scratch = kwargs.get('use_scratch',defaults['use_scratch'])
566 if os.path.isfile(grid):
567 logger.debug('Selected %s'%(grid))
568 if integrated:
569 basename = os.path.basename(grid)
570 return os.path.join(os.path.dirname(grid),basename[0]=='i' and basename or 'i'+basename)
571 logging.debug('Returning grid path: '+grid)
572 return grid
573
574 grid = grid.lower()
575
576
577 z = kwargs.get('z',defaults['z'])
578 Rv = kwargs.get('Rv',defaults['Rv'])
579
580 vturb = int(kwargs.get('vturb',defaults['vturb']))
581 odfnew = kwargs.get('odfnew',defaults['odfnew'])
582 alpha = kwargs.get('alpha',defaults['alpha'])
583 nover = kwargs.get('nover',defaults['nover'])
584
585 He = int(kwargs.get('He',defaults['He']))
586
587 t = kwargs.get('t',defaults['t'])
588 a = kwargs.get('a',defaults['a'])
589 c = kwargs.get('c',defaults['c'])
590 m = kwargs.get('m',defaults['m'])
591 co= kwargs.get('co',defaults['co'])
592
593 ct = kwargs.get('ct','mlt')
594
595
596 if grid=='fastwind':
597 basename = 'fastwind_sed.fits'
598 elif grid in ['kurucz','kurucz2']:
599 if not isinstance(z,str): z = '%.1f'%(z)
600 if not isinstance(vturb,str): vturb = '%d'%(vturb)
601 if grid=='kurucz2' and integrated:
602 postfix = '_lawfitzpatrick2004_Rv'
603 if not isinstance(Rv,str): Rv = '{:.2f}'.format(Rv)
604 postfix+= Rv
605 else:
606 postfix = ''
607 if not alpha and not nover and not odfnew:
608 basename = 'kurucz93_z%s_k%s_sed%sls.fits'%(z,vturb,postfix)
609 elif alpha and odfnew:
610 basename = 'kurucz93_z%s_ak%sodfnew_sed%s.fits'%(z,vturb,postfix)
611 elif odfnew:
612 basename = 'kurucz93_z%s_k%sodfnew_sed%s.fits'%(z,vturb,postfix)
613 elif nover:
614 basename = 'kurucz93_z%s_k%snover_sed%s.fits'%(z,vturb,postfix)
615 elif grid=='cmfgen':
616 basename = 'cmfgen_sed.fits'
617 elif grid=='sdb_uli':
618 if not isinstance(z,str): z = '%.1f'%(z)
619 if not isinstance(He,str): He = '%d'%(He)
620 basename = 'SED_int_h%s_z%s.fits'%(He,z)
621 elif grid=='wd_boris':
622 basename = 'SED_WD_Gaensicke.fits'
623 elif grid=='wd_da':
624 if integrated:
625 postfix = '_lawfitzpatrick2004_Rv'
626 if not isinstance(Rv,str): Rv = '{:.2f}'.format(Rv)
627 postfix+= Rv
628 else:
629 postfix = ''
630 basename = 'SED_WD_Koester_DA%s.fits'%(postfix)
631 elif grid=='wd_db':
632 basename = 'SED_WD_Koester_DB.fits'
633 elif grid=='marcs':
634 if not isinstance(z,str): z = '%.1f'%(z)
635 if not isinstance(t,str): t = '%.1f'%(t)
636 if not isinstance(a,str): a = '%.2f'%(a)
637 if not isinstance(c,str): c = '%.2f'%(c)
638 basename = 'marcsp_z%st%s_a%s_c%s_sed.fits'%(z,t,a,c)
639 elif grid=='marcs2':
640 basename = 'marcsp2_z0.00t2.0_m.1.0c0.00_sed.fits'
641 elif grid=='comarcs':
642 if not isinstance(z,str): z = '%.2f'%(z)
643 if not isinstance(co,str): co = '%.2f'%(co)
644 if not isinstance(m,str): m = '%.1f'%(m)
645 basename = 'comarcsp_z%sco%sm%sxi2.50_sed.fits'%(z,co,m)
646 elif grid=='stars':
647 basename = 'kurucz_stars_sed.fits'
648 elif grid=='tlusty':
649 if not isinstance(z,str): z = '%.2f'%(z)
650 basename = 'tlusty_z%s_sed.fits'%(z)
651 elif grid=='uvblue':
652 if not isinstance(z,str): z = '%.1f'%(z)
653 basename = 'uvblue_z%s_k2_sed.fits'%(z)
654 elif grid=='atlas12':
655 if not isinstance(z,str): z = '%.1f'%(z)
656 basename = 'atlas12_z%s_sed.fits'%(z)
657 elif grid=='tkachenko':
658 if not isinstance(z,str): z = '%.2f'%(z)
659 basename = 'tkachenko_z%s.fits'%(z)
660 elif grid=='nemo':
661 ct = ct.lower()
662 if ct=='mlt': ct = ct+'072'
663 else: ct = ct+'288'
664 basename = 'nemo_%s_z%.2f_v%d.fits'%(ct,z,vturb)
665 elif grid=='tmap':
666 if integrated:
667 postfix = '_lawfitzpatrick2004_Rv'
668 if not isinstance(Rv,str): Rv = '{:.2f}'.format(Rv)
669 postfix+= Rv
670 else:
671 postfix = ''
672 basename = 'TMAP2012_lowres%s.fits'%(postfix)
673 elif grid=='heberb':
674 basename = 'Heber2000_B_h909_extended.fits'
675 elif grid=='hebersdb':
676 basename = 'Heber2000_sdB_h909_extended.fits'
677
678 elif grid=='tmapsdb':
679
680 if integrated:
681 postfix = '_lawfitzpatrick2004_Rv'
682 if not isinstance(Rv,str): Rv = '{:.2f}'.format(Rv)
683 postfix+= Rv
684 else:
685 postfix = ''
686 basename = 'TMAP2012_sdB_extended%s.fits'%(postfix)
687 elif grid=='kuruczsdb':
688
689 if not isinstance(z,str): z = '%.1f'%(z)
690 if integrated:
691 postfix = '_lawfitzpatrick2004_Rv'
692 if not isinstance(Rv,str): Rv = '{:.2f}'.format(Rv)
693 postfix+= Rv
694 else:
695 postfix = ''
696 basename = 'kurucz_z%s_sdB%s.fits'%(z,postfix)
697 elif grid=='kuruczpagb':
698
699 if not isinstance(z,str): z = '%.1f'%(z)
700 if integrated:
701 postfix = '_lawfitzpatrick2004_Rv'
702 if not isinstance(Rv,str): Rv = '{:.2f}'.format(Rv)
703 postfix+= Rv
704 else:
705 postfix = ''
706 basename = 'kurucz_pAGB_z%s_sed%s.fits'%(z,postfix)
707
708 elif grid=='tmaptest':
709 """ Grids exclusively for testing purposes"""
710 if integrated:
711 postfix = '_lawfitzpatrick2004_Rv'
712 if not isinstance(Rv,str): Rv = '{:.2f}'.format(Rv)
713 postfix+= Rv
714 else:
715 postfix = ''
716 basename = 'TMAP2012_SEDtest%s.fits'%(postfix)
717
718 elif grid=='kurucztest':
719 """ Grids exclusively for testing purposes"""
720 if not isinstance(z,str): z = '%.1f'%(z)
721 if integrated:
722 postfix = '_lawfitzpatrick2004_Rv'
723 if not isinstance(Rv,str): Rv = '{:.2f}'.format(Rv)
724 postfix+= Rv
725 else:
726 postfix = ''
727 basename = 'kurucz_SEDtest_z%s%s.fits'%(z,postfix)
728 elif grid=='marcsjorissensp':
729 if not isinstance(z,str): z = '%.2f'%(z)
730 if not isinstance(a,str): a = '%.2f'%(a)
731 if not isinstance(m,str): m = '%.1f'%(m)
732 basename = 'MARCSJorissenSp_m%s_z%s_a%s_sed.fits'%(m,z,a)
733 else:
734 raise ValueError("Grid {} is not recognized: either give valid descriptive arguments, or give an absolute filepath".format(grid))
735
736 if not '*' in basename:
737 if use_scratch:
738 if integrated:
739 grid = scratchdir+'i'+basename
740 else:
741 grid = scratchdir+basename
742 else:
743 if integrated:
744 grid = config.get_datafile(basedir,'i'+basename)
745 else:
746 grid = config.get_datafile(basedir,basename)
747
748 else:
749 grid = config.glob(basedir,'i'+basename)
750 if use_scratch:
751 if integrated:
752 grid = glob.glob(scratchdir+'i'+basename)
753 else:
754 grid = glob.glob(scratchdir+basename)
755 else:
756 if integrated:
757 grid = config.glob(basedir,'i'+basename)
758 else:
759 grid = config.glob(basedir,basename)
760
761
762 logger.debug('Returning grid path(s): %s'%(grid))
763 return grid
764
767 """
768 Prepare input and output for blackbody-like functions.
769
770 If the user gives wavelength units and Flambda units, we only need to convert
771 everything to SI (and back to the desired units in the end).
772
773 If the user gives frequency units and Fnu units, we only need to convert
774 everything to SI ( and back to the desired units in the end).
775
776 If the user gives wavelength units and Fnu units, we need to convert
777 the wavelengths first to frequency.
778 """
779 @functools.wraps(fctn)
780 def dobb(x,T,**kwargs):
781 wave_units = kwargs.get('wave_units','AA')
782 flux_units = kwargs.get('flux_units','erg/s/cm2/AA')
783 to_kwargs = {}
784
785
786 curr_conv = constants._current_convention
787
788 x_unit_type = conversions.get_type(wave_units)
789 x = conversions.convert(wave_units,curr_conv,x)
790
791 if isinstance(T,tuple):
792 T = conversions.convert(T[1],'K',T[0])
793
794 y_unit_type = conversions.change_convention('SI',flux_units)
795
796 if y_unit_type=='kg1 rad-1 s-2' and x_unit_type=='length':
797 x = conversions.convert(conversions._conventions[curr_conv]['length'],'rad/s',x)
798 x_unit_type = 'frequency'
799 elif y_unit_type=='kg1 m-1 s-3' and x_unit_type=='frequency':
800 x = conversions.convert('rad/s',conversions._conventions[curr_conv]['length'],x)
801 x_unit_type = 'length'
802 elif not y_unit_type in ['kg1 rad-1 s-2','kg1 m-1 s-3']:
803 raise NotImplementedError(flux_units,y_unit_type)
804
805 if x_unit_type=='frequency':
806 x /= (2*np.pi)
807 to_kwargs['freq'] = (x,'Hz')
808 else:
809 to_kwargs['wave'] = (x,conversions._conventions[curr_conv]['length'])
810
811 I = fctn((x,x_unit_type),T)
812
813
814 disc_integrated = kwargs.get('disc_integrated',True)
815 ang_diam = kwargs.get('ang_diam',None)
816 if disc_integrated:
817 I *= np.sqrt(2*np.pi)
818 if ang_diam is not None:
819 scale = conversions.convert(ang_diam[1],'sr',ang_diam[0]/2.)
820 I *= scale
821 I = conversions.convert(curr_conv,flux_units,I,**to_kwargs)
822 return I
823
824 return dobb
825
826
827
828
829
830
831 @_blackbody_input
832 -def blackbody(x,T,wave_units='AA',flux_units='erg/s/cm2/AA',disc_integrated=True,ang_diam=None):
833 """
834 Definition of black body curve.
835
836 To get them into the same units as the Kurucz disc-integrated SEDs, they are
837 multiplied by sqrt(2*pi) (set C{disc_integrated=True}).
838
839 You can only give an angular diameter if disc_integrated is True.
840
841 To convert the scale parameter back to mas, simply do:
842
843 ang_diam = 2*conversions.convert('sr','mas',scale)
844
845 See decorator L{blackbody_input} for details on how the input parameters
846 are handled: the user is free to choose wavelength or frequency units, choose
847 *which* wavelength or frequency units, and can even mix them. To be sure that
848 everything is handled correctly, we need to do some preprocessing and unit
849 conversions.
850
851 Be careful when, e.g. during fitting, scale contains an error: be sure to set
852 the option C{unpack=True} in the L{conversions.convert} function!
853
854 >>> x = np.linspace(2.3595,193.872,500)
855 >>> F1 = blackbody(x,280.,wave_units='AA',flux_units='Jy',ang_diam=(1.,'mas'))
856 >>> F2 = rayleigh_jeans(x,280.,wave_units='micron',flux_units='Jy',ang_diam=(1.,'mas'))
857 >>> F3 = wien(x,280.,wave_units='micron',flux_units='Jy',ang_diam=(1.,'mas'))
858
859
860 >>> p = pl.figure()
861 >>> p = pl.subplot(121)
862 >>> p = pl.plot(x,F1)
863 >>> p = pl.plot(x,F2)
864 >>> p = pl.plot(x,F3)
865
866
867 >>> F1 = blackbody(x,280.,wave_units='AA',flux_units='erg/s/cm2/AA',ang_diam=(1.,'mas'))
868 >>> F2 = rayleigh_jeans(x,280.,wave_units='micron',flux_units='erg/s/cm2/AA',ang_diam=(1.,'mas'))
869 >>> F3 = wien(x,280.,wave_units='micron',flux_units='erg/s/cm2/AA',ang_diam=(1.,'mas'))
870
871
872 >>> p = pl.subplot(122)
873 >>> p = pl.plot(x,F1)
874 >>> p = pl.plot(x,F2)
875 >>> p = pl.plot(x,F3)
876
877
878 @param: wavelength
879 @type: ndarray
880 @param T: temperature, unit
881 @type: tuple (float,str)
882 @param wave_units: wavelength units (frequency or length)
883 @type wave_units: str (units)
884 @param flux_units: flux units (could be in Fnu-units or Flambda-units)
885 @type flux_units: str (units)
886 @param disc_integrated: if True, they are in the same units as Kurucz-disc-integrated SEDs
887 @type disc_integrated: bool
888 @param ang_diam: angular diameter (in mas or rad or something similar)
889 @type ang_diam: (value, unit)
890 @return: intensity
891 @rtype: array
892 """
893 x,x_unit_type = x
894
895 if x_unit_type=='frequency':
896 factor = 2.0 * constants.hh / constants.cc**2
897 expont = constants.hh / (constants.kB*T)
898 I = factor * x**3 * 1. / (np.exp(expont*x) - 1.)
899 elif x_unit_type=='length':
900 factor = 2.0 * constants.hh * constants.cc**2
901 expont = constants.hh*constants.cc / (constants.kB*T)
902 I = factor / x**5. * 1. / (np.exp(expont/x) - 1.)
903 else:
904 raise ValueError(x_unit_type)
905 return I
906
907
908 @_blackbody_input
909 -def rayleigh_jeans(x,T,wave_units='AA',flux_units='erg/s/cm2/AA',disc_integrated=True,ang_diam=None):
910 """
911 Rayleigh-Jeans approximation of a black body.
912
913 Valid at long wavelengths.
914
915 For input details, see L{blackbody}.
916
917 @return: intensity
918 @rtype: array
919 """
920 x,x_unit_type = x
921
922 if x_unit_type=='frequency':
923 factor = 2.0 * constants.kB*T / constants.cc**2
924 I = factor * x**2
925 elif x_unit_type=='length':
926 factor = 2.0 * constants.cc * constants.kB*T
927 I = factor / x**4.
928 else:
929 raise ValueError(unit_type)
930 return I
931
932
933 @_blackbody_input
934 -def wien(x,T,wave_units='AA',flux_units='erg/s/cm2/AA',disc_integrated=True,ang_diam=None):
935 """
936 Wien approximation of a black body.
937
938 Valid at short wavelengths.
939
940 For input details, see L{blackbody}.
941
942 @return: intensity
943 @rtype: array
944 """
945 x,x_unit_type = x
946
947 if x_unit_type=='frequency':
948 factor = 2.0 * constants.hh / constants.cc**2
949 expont = constants.hh / (constants.kB*T)
950 I = factor * x**3 * 1. * np.exp(-expont*x)
951 elif x_unit_type=='length':
952 factor = 2.0 * constants.hh * constants.cc**2
953 expont = constants.hh*constants.cc / (constants.kB*T)
954 I = factor / x**5. * np.exp(-expont/x)
955 else:
956 raise ValueError(unit_type)
957 return I
958
959
960 -def get_table_single(teff=None,logg=None,ebv=None,rad=None,star=None,
961 wave_units='AA',flux_units='erg/s/cm2/AA/sr',**kwargs):
962 """
963 Retrieve the spectral energy distribution of a model atmosphere.
964
965 Wavelengths in A (angstrom)
966 Fluxes in Ilambda = ergs/cm2/s/AA/sr, except specified via 'units',
967
968 If you give 'units', and /sr is not included, you are responsible yourself
969 for giving an extra keyword with the angular diameter C{ang_diam}, or other
970 possibilities offered by the C{units.conversions.convert} function.
971
972 Possibility to redden the fluxes according to the reddening parameter EB_V.
973
974 Extra kwargs can specify the grid type.
975 Extra kwargs can specify constraints on the size of the grid to interpolate.
976 Extra kwargs can specify reddening law types.
977 Extra kwargs can specify information for conversions.
978
979 Example usage:
980
981 >>> from pylab import figure,gca,subplot,title,gcf,loglog
982 >>> p = figure(figsize=(10,6))
983 >>> p=gcf().canvas.set_window_title('Test of <get_table>')
984 >>> p = subplot(131)
985 >>> p = loglog(*get_table(grid='FASTWIND',teff=35000,logg=4.0),**dict(label='Fastwind'))
986 >>> p = loglog(*get_table(grid='KURUCZ',teff=35000,logg=4.0),**dict(label='Kurucz'))
987 >>> p = loglog(*get_table(grid='TLUSTY',teff=35000,logg=4.0),**dict(label='Tlusty'))
988 >>> p = loglog(*get_table(grid='MARCS',teff=5000,logg=2.0),**dict(label='Marcs'))
989 >>> p = loglog(*get_table(grid='KURUCZ',teff=5000,logg=2.0),**dict(label='Kurucz'))
990 >>> p = pl.xlabel('Wavelength [angstrom]');p = pl.ylabel('Flux [erg/s/cm2/AA/sr]')
991 >>> p = pl.legend(loc='upper right',prop=dict(size='small'))
992 >>> p = subplot(132)
993 >>> p = loglog(*get_table(grid='FASTWIND',teff=35000,logg=4.0,wave_units='micron',flux_units='Jy/sr'),**dict(label='Fastwind'))
994 >>> p = loglog(*get_table(grid='KURUCZ',teff=35000,logg=4.0,wave_units='micron',flux_units='Jy/sr'),**dict(label='Kurucz'))
995 >>> p = loglog(*get_table(grid='TLUSTY',teff=35000,logg=4.0,wave_units='micron',flux_units='Jy/sr'),**dict(label='Tlusty'))
996 >>> p = loglog(*get_table(grid='MARCS',teff=5000,logg=2.0,wave_units='micron',flux_units='Jy/sr'),**dict(label='Marcs'))
997 >>> p = loglog(*get_table(grid='KURUCZ',teff=5000,logg=2.0,wave_units='micron',flux_units='Jy/sr'),**dict(label='Kurucz'))
998 >>> p = pl.xlabel('Wavelength [micron]');p = pl.ylabel('Flux [Jy/sr]')
999 >>> p = pl.legend(loc='upper right',prop=dict(size='small'))
1000 >>> p = subplot(133);p = title('Kurucz')
1001 >>> p = loglog(*get_table(grid='KURUCZ',teff=10000,logg=4.0,wave_units='micron',flux_units='Jy/sr'),**dict(label='10000'))
1002 >>> p = loglog(*get_table(grid='KURUCZ',teff=10250,logg=4.0,wave_units='micron',flux_units='Jy/sr'),**dict(label='10250'))
1003 >>> p = loglog(*get_table(grid='KURUCZ',teff=10500,logg=4.0,wave_units='micron',flux_units='Jy/sr'),**dict(label='10500'))
1004 >>> p = loglog(*get_table(grid='KURUCZ',teff=10750,logg=4.0,wave_units='micron',flux_units='Jy/sr'),**dict(label='10750'))
1005 >>> p = loglog(*get_table(grid='KURUCZ',teff=11000,logg=4.0,wave_units='micron',flux_units='Jy/sr'),**dict(label='11000'))
1006 >>> p = pl.xlabel('Wavelength [micron]');p = pl.ylabel('Flux [Jy/sr]')
1007 >>> p = pl.legend(loc='upper right',prop=dict(size='small'))
1008
1009 ]]include figure]]ivs_sed_model_comparison.png]
1010
1011 @param teff: effective temperature
1012 @type teff: float
1013 @param logg: logarithmic gravity (cgs)
1014 @type logg: float
1015 @param ebv: reddening coefficient
1016 @type ebv: float
1017 @param wave_units: units to convert the wavelengths to (if not given, A)
1018 @type wave_units: str
1019 @param flux_units: units to convert the fluxes to (if not given, erg/s/cm2/AA/sr)
1020 @type flux_units: str
1021 @return: wavelength,flux
1022 @rtype: (ndarray,ndarray)
1023 """
1024
1025
1026 gridfile = get_file(**kwargs)
1027
1028
1029 ff = pf.open(gridfile)
1030
1031
1032
1033 if star is not None:
1034 wave = ff[star.upper()].data.field('wavelength')
1035 flux = ff[star.upper()].data.field('flux')
1036 else:
1037 teff = float(teff)
1038 logg = float(logg)
1039
1040 try:
1041
1042 mod_name = "T%05d_logg%01.02f" %(teff,logg)
1043 mod = ff[mod_name]
1044 wave = mod.data.field('wavelength')
1045 flux = mod.data.field('flux')
1046 logger.debug('Model SED taken directly from file (%s)'%(os.path.basename(gridfile)))
1047
1048 except KeyError:
1049
1050
1051
1052 wave,teffs,loggs,flux,flux_grid = get_grid_mesh(**kwargs)
1053 logger.debug('Model SED interpolated from grid %s (%s)'%(os.path.basename(gridfile),kwargs))
1054 wave = wave + 0.
1055 flux = flux_grid(np.log10(teff),logg) + 0.
1056
1057
1058 wave = np.array(wave,float)
1059 flux = np.array(flux,float)
1060
1061 if ebv is not None and ebv>0:
1062 if 'wave' in kwargs.keys():
1063 removed = kwargs.pop('wave')
1064 flux = reddening.redden(flux,wave=wave,ebv=ebv,rtype='flux',**kwargs)
1065 if flux_units!='erg/s/cm2/AA/sr':
1066 flux = conversions.convert('erg/s/cm2/AA/sr',flux_units,flux,wave=(wave,'AA'),**kwargs)
1067 if wave_units!='AA':
1068 wave = conversions.convert('AA',wave_units,wave,**kwargs)
1069
1070 ff.close()
1071
1072 if rad != None:
1073 flux = rad**2 * flux
1074
1075 return wave,flux
1076
1077
1078 -def get_itable_single(teff=None,logg=None,ebv=0,z=0,rad=None,photbands=None,
1079 wave_units=None,flux_units='erg/s/cm2/AA/sr',**kwargs):
1080 """
1081 Retrieve the spectral energy distribution of a model atmosphere in
1082 photometric passbands.
1083
1084 Wavelengths in A (angstrom). If you set 'wavelengths' to None, no effective
1085 wavelengths will be calculated. Otherwise, the effective wavelength is
1086 calculated taking the model flux into account.
1087 Fluxes in Ilambda = ergs/cm2/s/AA/sr, except specified via 'units',
1088
1089 If you give 'units', and /sr is not included, you are responsible yourself
1090 for giving an extra keyword with the angular diameter C{ang_diam}, or other
1091 possibilities offered by the C{units.conversions.convert} function.
1092
1093 Possibility to redden the fluxes according to the reddening parameter EB_V.
1094
1095 Extra kwargs can specify the grid type.
1096 Extra kwargs can specify constraints on the size of the grid to interpolate.
1097 Extra kwargs can specify reddening law types.
1098 Extra kwargs can specify information for conversions.
1099
1100 @param teff: effective temperature
1101 @type teff: float
1102 @param logg: logarithmic gravity (cgs)
1103 @type logg: float
1104 @param ebv: reddening coefficient
1105 @type ebv: float
1106 @param photbands: photometric passbands
1107 @type photbands: list of photometric passbands
1108 @param wave_units: units to convert the wavelengths to (if not given, A)
1109 @type wave_units: str
1110 @param flux_units: units to convert the fluxes to (if not given, erg/s/cm2/AA/sr)
1111 @type flux_units: str
1112 @keyword clear_memory: flag to clear memory from previously loaded SED tables.
1113 If you set it to False, you can easily get an overloaded memory!
1114 @type clear_memory: boolean
1115 @return: (wave,) flux, absolute luminosity
1116 @rtype: (ndarray,)ndarray,float
1117 """
1118 if 'vrad' in kwargs:
1119 logger.debug('vrad is NOT taken into account when interpolating in get_itable()')
1120 if 'rv' in kwargs:
1121 logger.debug('Rv is NOT taken into account when interpolating in get_itable()')
1122
1123 if photbands is None:
1124 raise ValueError('no photometric passbands given')
1125 ebvrange = kwargs.pop('ebvrange',(-np.inf,np.inf))
1126 zrange = kwargs.pop('zrange',(-np.inf,np.inf))
1127 clear_memory = kwargs.pop('clear_memory',True)
1128
1129
1130
1131
1132 markers,(g_teff,g_logg,g_ebv,g_z),gpnts,ext = _get_itable_markers(photbands,ebvrange=ebvrange,zrange=zrange,
1133 include_Labs=True,clear_memory=clear_memory,**kwargs)
1134
1135
1136 try:
1137 input_code = float('%3d%05d%03d%03d'%(int(round((z+5)*100)),int(round(teff)),int(round(logg*100)),int(round(ebv*100))))
1138 index = markers.searchsorted(input_code)
1139 output_code = markers[index]
1140
1141
1142 if not input_code==output_code:
1143 raise KeyError
1144
1145 flux = ext[index]
1146
1147
1148
1149
1150 except KeyError:
1151
1152
1153 if not new_scipy:
1154 teff = teff+1e-2
1155 logg = logg-1e-6
1156 ebv = ebv+1e-6
1157
1158
1159 i_teff = max(1,g_teff.searchsorted(teff))
1160 i_logg = max(1,g_logg.searchsorted(logg))
1161 i_ebv = max(1,g_ebv.searchsorted(ebv))
1162 i_z = max(1,g_z.searchsorted(z))
1163 if i_teff==len(g_teff): i_teff -= 1
1164 if i_logg==len(g_logg): i_logg -= 1
1165 if i_ebv==len(g_ebv): i_ebv -= 1
1166 if i_z==len(g_z): i_z -= 1
1167
1168 teffs_subgrid = g_teff[i_teff-1:i_teff+1]
1169 loggs_subgrid = g_logg[i_logg-1:i_logg+1]
1170 ebvs_subgrid = g_ebv[i_ebv-1:i_ebv+1]
1171 zs_subgrid = g_z[i_z-1:i_z+1]
1172
1173
1174
1175
1176
1177 if not (z in g_z):
1178 fluxes = np.zeros((2,2,2,2,len(photbands)+1))
1179 for i,j,k in itertools.product(xrange(2),xrange(2),xrange(2)):
1180 input_code = float('%3d%05d%03d%03d'%(int(round((zs_subgrid[i]+5)*100)),\
1181 int(round(teffs_subgrid[j])),\
1182 int(round(loggs_subgrid[k]*100)),\
1183 int(round(ebvs_subgrid[1]*100))))
1184 index = markers.searchsorted(input_code)
1185 fluxes[i,j,k] = ext[index-1:index+1]
1186 myf = InterpolatingFunction([zs_subgrid,np.log10(teffs_subgrid),
1187 loggs_subgrid,ebvs_subgrid],np.log10(fluxes),default=-100*np.ones_like(fluxes.shape[1]))
1188 flux = 10**myf(z,np.log10(teff),logg,ebv) + 0.
1189
1190
1191 else:
1192 fluxes = np.zeros((2,2,2,len(photbands)+1))
1193 for i,j in itertools.product(xrange(2),xrange(2)):
1194 input_code = float('%3d%05d%03d%03d'%(int(round((z+5)*100)),\
1195 int(round(teffs_subgrid[i])),\
1196 int(round(loggs_subgrid[j]*100)),\
1197 int(round(ebvs_subgrid[1]*100))))
1198 index = markers.searchsorted(input_code)
1199 fluxes[i,j] = ext[index-1:index+1]
1200 myf = InterpolatingFunction([np.log10(teffs_subgrid),
1201 loggs_subgrid,ebvs_subgrid],np.log10(fluxes),default=-100*np.ones_like(fluxes.shape[1]))
1202 flux = 10**myf(np.log10(teff),logg,ebv) + 0.
1203
1204 else:
1205
1206 i_teff = max(1,g_teff.searchsorted(teff))
1207 i_logg = max(1,g_logg.searchsorted(logg))
1208 i_ebv = max(1,g_ebv.searchsorted(ebv))
1209 i_z = max(1,g_z.searchsorted(z))
1210
1211 if i_teff==len(g_teff): i_teff -= 1
1212 if i_logg==len(g_logg): i_logg -= 1
1213 if i_ebv==len(g_ebv): i_ebv -= 1
1214 if i_z==len(g_z): i_z -= 1
1215 if not (z in g_z):
1216
1217 myflux = np.zeros((16,4+len(photbands)+1))
1218 mygrid = itertools.product(g_teff[i_teff-1:i_teff+1],g_logg[i_logg-1:i_logg+1],g_z[i_z-1:i_z+1])
1219 for i,(t,g,zz) in enumerate(mygrid):
1220 myflux[2*i,:4] = t,g,g_ebv[i_ebv-1],zz
1221 myflux[2*i+1,:4] = t,g,g_ebv[i_ebv],zz
1222 input_code = float('%3d%05d%03d%03d'%(int(round((zz+5)*100)),\
1223 int(round(t)),int(round(g*100)),\
1224 int(round(g_ebv[i_ebv]*100))))
1225 index = markers.searchsorted(input_code)
1226 myflux[2*i,4:] = ext[index-1]
1227 myflux[2*i+1,4:] = ext[index]
1228
1229 myflux[:,0] = np.log10(myflux[:,0])
1230 flux = 10**griddata(myflux[:,:4],np.log10(myflux[:,4:]),(np.log10(teff),logg,ebv,z))
1231 else:
1232
1233 myflux = np.zeros((8,3+len(photbands)+1))
1234 mygrid = itertools.product(g_teff[i_teff-1:i_teff+1],g_logg[i_logg-1:i_logg+1])
1235 for i,(t,g) in enumerate(mygrid):
1236 myflux[2*i,:3] = t,g,g_ebv[i_ebv-1]
1237 myflux[2*i+1,:3] = t,g,g_ebv[i_ebv]
1238 input_code = float('%3d%05d%03d%03d'%(int(round((z+5)*100)),\
1239 int(round(t)),int(round(g*100)),\
1240 int(round(g_ebv[i_ebv]*100))))
1241 index = markers.searchsorted(input_code)
1242 myflux[2*i,3:] = ext[index-1]
1243 myflux[2*i+1,3:] = ext[index]
1244
1245 myflux[:,0] = np.log10(myflux[:,0])
1246 flux = 10**griddata(myflux[:,:3],np.log10(myflux[:,3:]),(np.log10(teff),logg,ebv))
1247 except IndexError:
1248
1249 raise ValueError('point outside of grid (teff={teff}, logg={logg}, ebv={ebv}, z={z}'.format(**locals()))
1250 except ValueError:
1251
1252 raise ValueError('point outside of grid (teff={teff}, logg={logg}, ebv={ebv}, z={z}'.format(**locals()))
1253 if np.any(np.isnan(flux)):
1254
1255 raise ValueError('point outside of grid (teff={teff}, logg={logg}, ebv={ebv}, z={z}'.format(**locals()))
1256 if np.any(np.isinf(flux)):
1257 flux = np.zeros(fluxes.shape[-1])
1258
1259
1260
1261
1262 flux,Labs = np.array(flux[:-1],float),flux[-1]
1263
1264
1265 if rad != None:
1266 flux,Labs = flux*rad**2, Labs*rad**2
1267
1268 if flux_units!='erg/s/cm2/AA/sr':
1269 flux = conversions.nconvert('erg/s/cm2/AA/sr',flux_units,flux,photband=photbands,**kwargs)
1270
1271 if wave_units is not None:
1272 model = get_table(teff=teff,logg=logg,ebv=ebv,**kwargs)
1273 wave = filters.eff_wave(photbands,model=model)
1274 if wave_units !='AA':
1275 wave = wave = conversions.convert('AA',wave_units,wave,**kwargs)
1276
1277 return wave,flux,Labs
1278 else:
1279 return flux,Labs
1280
1281 -def get_itable(photbands=None, wave_units=None, flux_units='erg/s/cm2/AA/sr',
1282 grids=None, **kwargs):
1283 """
1284 Retrieve the integrated spectral energy distribution of a combined model
1285 atmosphere.
1286
1287 >>> teff1,teff2 = 20200,5100
1288 >>> logg1,logg2 = 4.35,2.00
1289 >>> ebv = 0.2,0.2
1290 >>> photbands = ['JOHNSON.U','JOHNSON.V','2MASS.J','2MASS.H','2MASS.KS']
1291
1292 >>> wave1,flux1 = get_table(teff=teff1,logg=logg1,ebv=ebv[0])
1293 >>> wave2,flux2 = get_table(teff=teff2,logg=logg2,ebv=ebv[1])
1294 >>> wave3,flux3 = get_table_multiple(teff=(teff1,teff2),logg=(logg1,logg2),ebv=ebv,radius=[1,20])
1295
1296 >>> iwave1,iflux1,iLabs1 = get_itable(teff=teff1,logg=logg1,ebv=ebv[0],photbands=photbands,wave_units='AA')
1297 >>> iflux2,iLabs2 = get_itable(teff=teff2,logg=logg2,ebv=ebv[1],photbands=photbands)
1298 >>> iflux3,iLabs3 = get_itable_multiple(teff=(teff1,teff2),logg=(logg1,logg2),z=(0,0),ebv=ebv,radius=[1,20.],photbands=photbands)
1299
1300 >>> p = pl.figure()
1301 >>> p = pl.gcf().canvas.set_window_title('Test of <get_itable_multiple>')
1302 >>> p = pl.loglog(wave1,flux1,'r-')
1303 >>> p = pl.loglog(iwave1,iflux1,'ro',ms=10)
1304 >>> p = pl.loglog(wave2,flux2*20**2,'b-')
1305 >>> p = pl.loglog(iwave1,iflux2*20**2,'bo',ms=10)
1306 >>> p = pl.loglog(wave3,flux3,'k-',lw=2)
1307 >>> p = pl.loglog(iwave1,iflux3,'kx',ms=10,mew=2)
1308
1309 @param teff: effective temperature
1310 @type teff: tuple floats
1311 @param logg: logarithmic gravity (cgs)
1312 @type logg: tuple floats
1313 @param ebv: reddening coefficient
1314 @type ebv: tuple floats
1315 @param z: metallicity
1316 @type z: tuple floats
1317 @param radius: ratio of R_i/(R_{i-1})
1318 @type radius: tuple of floats
1319 @param photbands: photometric passbands
1320 @type photbands: list
1321 @param flux_units: units to convert the fluxes to (if not given, erg/s/cm2/AA/sr)
1322 @type flux_units: str
1323 @param grids: specifications for grid1
1324 @type grids: list of dict
1325 @param full_output: return all individual SEDs
1326 @type full_output: boolean
1327 @return: wavelength,flux
1328 @rtype: (ndarray,ndarray)
1329 """
1330
1331 values, parameters, components = {}, set(), set()
1332 for key in kwargs.keys():
1333 if re.search("^(teff|logg|ebv|z|rad)\d?$", key):
1334 par, comp = re.findall("^(teff|logg|ebv|z|rad)(\d?)$", key)[0]
1335 values[key] = kwargs.pop(key)
1336 parameters.add(par)
1337 components.add(comp)
1338
1339
1340 if len(components) == 1:
1341 kwargs.update(values)
1342 return get_itable_single(photbands=photbands,wave_units=wave_units,
1343 flux_units=flux_units,**kwargs)
1344
1345
1346
1347 fluxes, Labs = [],[]
1348 for i, (comp, grid) in enumerate(zip(components,defaults_multiple)):
1349 trash = grid.pop('z',0.0), grid.pop('Rv',0.0)
1350 kwargs_ = kwargs
1351 kwargs_.update(grid)
1352 for par in parameters:
1353 kwargs_[par] = values[par+comp] if par+comp in values else values[par]
1354
1355 f,L = get_itable_single(photbands=photbands,wave_units=None,**kwargs_)
1356
1357 fluxes.append(f)
1358 Labs.append(L)
1359
1360 fluxes = np.sum(fluxes,axis=0)
1361 Labs = np.sum(Labs,axis=0)
1362
1363 if flux_units!='erg/s/cm2/AA/sr':
1364 fluxes = np.array([conversions.convert('erg/s/cm2/AA/sr',flux_units,fluxes[i],photband=photbands[i]) for i in range(len(fluxes))])
1365
1366 if wave_units is not None:
1367 model = get_table_multiple(teff=teff,logg=logg,ebv=ebv, grids=grids,**kwargs)
1368 wave = filters.eff_wave(photbands,model=model)
1369 if wave_units !='AA':
1370 wave = wave = conversions.convert('AA',wave_units,wave)
1371 return wave,fluxes,Labs
1372
1373 return fluxes,Labs
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405 -def get_itable_single_pix(teff=None,logg=None,ebv=None,z=0,rv=3.1,vrad=0,photbands=None,
1406 wave_units=None,flux_units='erg/s/cm2/AA/sr',**kwargs):
1407 """
1408 Super fast grid interpolator.
1409
1410 Possible kwargs are teffrange,loggrange etc.... that are past on to
1411 L{_get_pix_grid}. You should probably use these options when you want to
1412 interpolate in many variables; supplying these ranges will make the grid
1413 smaller and thus decrease memory usage.
1414
1415 It is possible to fix C{teff}, C{logg}, C{ebv}, C{z}, C{rv} and/or C{vrad}
1416 to one value, in which case it B{has} to be a point in the grid. If you want
1417 to retrieve a list of fluxes with the same ebv value that is not in the grid,
1418 you need to give an array with all equal values. The reason is that the
1419 script can try to minimize the number of interpolations, by fixing a
1420 variable on a grid point. The fluxes on the other gridpoints will then not
1421 be interpolated over. These parameter also have to be listed with the
1422 additional C{exc_interpolpar} keyword.
1423
1424 >>> teffs = np.linspace(5000,7000,100)
1425 >>> loggs = np.linspace(4.0,4.5,100)
1426 >>> ebvs = np.linspace(0,1,100)
1427 >>> zs = np.linspace(-0.5,0.5,100)
1428 >>> rvs = np.linspace(2.2,5.0,100)
1429
1430 >>> set_defaults(grid='kurucz2')
1431 >>> flux,labs = get_itable_pix(teffs,loggs,ebvs,zs,rvs,photbands=['JOHNSON.V'])
1432
1433 >>> names = ['teffs','loggs','ebvs','zs','rvs']
1434 >>> p = pl.figure()
1435 >>> for i in range(len(names)):
1436 ... p = pl.subplot(2,3,i+1)
1437 ... p = pl.plot(locals()[names[i]],flux[0],'k-')
1438 ... p = pl.xlabel(names[i])
1439
1440
1441 Thanks to Steven Bloemen for the core implementation of the interpolation
1442 algorithm.
1443
1444 The addition of the exc_interpolpar keyword was done by Michel Hillen (Jan 2016).
1445 """
1446
1447
1448 ebv = np.array([0 for i in teff]) if ebv is None else ebv
1449 z = np.array([0.for i in teff]) if z is None else z
1450 rv = np.array([3.1 for i in teff]) if rv is None else rv
1451 vrad = np.array([0 for i in teff]) if vrad is None else vrad
1452
1453
1454
1455
1456
1457
1458 vrad = 0
1459 N = 1
1460 clear_memory = kwargs.pop('clear_memory',False)
1461
1462 for var in ['teff','logg','ebv','z','rv','vrad']:
1463 if not hasattr(locals()[var],'__iter__'):
1464 kwargs.setdefault(var+'range',(locals()[var],locals()[var]))
1465 else:
1466 N = len(locals()[var])
1467
1468
1469 axis_values,gridpnts,pixelgrid,cols = _get_pix_grid(photbands,
1470 include_Labs=True,clear_memory=clear_memory,**kwargs)
1471
1472
1473
1474
1475
1476
1477
1478 for var in kwargs.get('exc_interpolpar',[]):
1479
1480 var_uniquevalue = np.unique(np.array(locals()[var]))
1481
1482 if len(var_uniquevalue) > 1:
1483 logger.warning('{} is requested to be excluded from interpolation, although fluxes for more than one value are requested!?'.format(var))
1484 else:
1485
1486 var_index = np.where(cols == var)[0]
1487
1488 var_uniquevalue_index = np.where(axis_values[var_index] == var_uniquevalue[0])[0]
1489
1490 if len(var_uniquevalue_index) == 0:
1491 logger.warning('{} can only be excluded from interpolation, as requested, if its values are all equal to an actual grid point!'.format(var))
1492 else:
1493
1494 trash = axis_values.pop(var_index)
1495 cols = np.delete(cols,[var_index])
1496
1497
1498 indices = [x for x in range(pixelgrid.ndim)]
1499 indices.remove(var_index)
1500 indices.insert(0,var_index)
1501 pixelgrid = np.transpose(pixelgrid, indices)
1502
1503 pixelgrid = pixelgrid[var_uniquevalue_index[0]]
1504
1505
1506 values = np.zeros((len(cols),N))
1507 for i,col in enumerate(cols):
1508 values[i] = locals()[col]
1509 pars = 10**interpol.interpolate(values,axis_values,pixelgrid)
1510 flux,Labs = pars[:-1],pars[-1]
1511
1512
1513 if 'rad' in kwargs:
1514 flux,Labs = flux*kwargs['rad']**2, Labs*kwargs['rad']**2
1515
1516
1517 if flux_units!='erg/s/cm2/AA/sr':
1518 flux = conversions.nconvert('erg/s/cm2/AA/sr',flux_units,flux,photband=photbands,**kwargs)
1519 if wave_units is not None:
1520 model = get_table(teff=teff,logg=logg,ebv=ebv,**kwargs)
1521 wave = filters.eff_wave(photbands,model=model)
1522 if wave_units !='AA':
1523 wave = conversions.convert('AA',wave_units,wave,**kwargs)
1524 return wave,flux,Labs
1525 else:
1526 return flux,Labs
1527
1528 -def get_itable_pix(photbands=None, wave_units=None, flux_units='erg/s/cm2/AA/sr',
1529 grids=None, **kwargs):
1530 """
1531 Super fast grid interpolator for multiple tables, completely based on get_itable_pix.
1532 """
1533
1534 values, parameters, components = {}, set(), set()
1535 for key in kwargs.keys():
1536 if re.search("^(teff|logg|ebv|z|rv|vrad|rad)\d?$", key):
1537 par, comp = re.findall("^(teff|logg|ebv|z|rv|vrad|rad)(\d?)$", key)[0]
1538 values[key] = kwargs.pop(key)
1539 parameters.add(par)
1540 components.add(comp)
1541
1542 if len(components) == 1:
1543 kwargs.update(values)
1544 return get_itable_single_pix(photbands=photbands,wave_units=wave_units,
1545 flux_units=flux_units,**kwargs)
1546
1547
1548 fluxes, Labs = [],[]
1549 for i, (comp, grid) in enumerate(zip(components,defaults_multiple)):
1550 trash = grid.pop('z',0.0), grid.pop('Rv',0.0)
1551 kwargs_ = kwargs
1552 kwargs_.update(grid)
1553 for par in parameters:
1554 kwargs_[par] = values[par+comp] if par+comp in values else values[par]
1555
1556 f,L = get_itable_single_pix(photbands=photbands,wave_units=None,**kwargs_)
1557
1558 fluxes.append(f)
1559 Labs.append(L)
1560 fluxes = np.sum(fluxes,axis=0)
1561 Labs = np.sum(Labs,axis=0)
1562
1563 if flux_units!='erg/s/cm2/AA/sr':
1564 fluxes = np.array([conversions.convert('erg/s/cm2/AA/sr',flux_units,fluxes[i],photband=photbands[i]) for i in range(len(fluxes))])
1565
1566 if wave_units is not None:
1567 model = get_table_multiple(teff=teff,logg=logg,ebv=ebv, grids=grids,**kwargs)
1568 wave = filters.eff_wave(photbands,model=model)
1569 if wave_units !='AA':
1570 wave = conversions.convert('AA',wave_units,wave)
1571 return wave,fluxes,Labs
1572 return fluxes,Labs
1573
1574
1575
1576
1577 -def get_table(wave_units='AA',flux_units='erg/cm2/s/AA/sr',grids=None,full_output=False,**kwargs):
1578 """
1579 Retrieve the spectral energy distribution of a combined model atmosphere.
1580
1581 Example usage:
1582
1583 >>> teff1,teff2 = 20200,5100
1584 >>> logg1,logg2 = 4.35,2.00
1585 >>> wave1,flux1 = get_table(teff=teff1,logg=logg1,ebv=0.2)
1586 >>> wave2,flux2 = get_table(teff=teff2,logg=logg2,ebv=0.2)
1587 >>> wave3,flux3 = get_table_multiple(teff=(teff1,teff2),logg=(logg1,logg2),ebv=(0.2,0.2),radius=[1,20])
1588
1589 >>> p = pl.figure()
1590 >>> p = pl.gcf().canvas.set_window_title('Test of <get_table_multiple>')
1591 >>> p = pl.loglog(wave1,flux1,'r-')
1592 >>> p = pl.loglog(wave2,flux2,'b-')
1593 >>> p = pl.loglog(wave2,flux2*20**2,'b--')
1594 >>> p = pl.loglog(wave3,flux3,'k-',lw=2)
1595
1596 @param teff: effective temperature
1597 @type teff: tuple floats
1598 @param logg: logarithmic gravity (cgs)
1599 @type logg: tuple floats
1600 @param ebv: tuple reddening coefficients
1601 @type ebv: tuple floats
1602 @param radius: radii of the stars
1603 @type radius: tuple of floats
1604 @param wave_units: units to convert the wavelengths to (if not given, A)
1605 @type wave_units: str
1606 @param flux_units: units to convert the fluxes to (if not given, erg/s/cm2/AA/sr)
1607 @type flux_units: str
1608 @param grids: specifications for grid1
1609 @type grids: list of dict
1610 @param full_output: return all individual SEDs
1611 @type full_output: boolean
1612 @return: wavelength,flux
1613 @rtype: (ndarray,ndarray)
1614 """
1615 values, parameters, components = {}, set(), set()
1616 for key in kwargs.keys():
1617 if re.search("^(teff|logg|ebv|z|rad)\d?$", key):
1618 par, comp = re.findall("^(teff|logg|ebv|z|rad)(\d?)$", key)[0]
1619 values[key] = kwargs.pop(key)
1620 parameters.add(par)
1621 components.add(comp)
1622
1623
1624 if len(components) == 1:
1625 kwargs.update(values)
1626 return get_table_single(wave_units=wave_units, flux_units=flux_units,
1627 full_output=full_output,**kwargs)
1628
1629
1630
1631 waves, fluxes = [],[]
1632 for i, (comp, grid) in enumerate(zip(components,defaults_multiple)):
1633 trash = grid.pop('z',0.0), grid.pop('Rv',0.0)
1634 kwargs_ = kwargs.copy()
1635 kwargs_.update(grid)
1636 for par in parameters:
1637 kwargs_[par] = values[par+comp] if par+comp in values else values[par]
1638
1639 w,f = get_table_single(**kwargs_)
1640
1641 waves.append(w)
1642 fluxes.append(f)
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660 waves_ = np.sort(np.hstack(waves))
1661 waves_ = np.hstack([waves_[0],waves_[1:][np.diff(waves_)>0]])
1662
1663 wstart = max([w[0] for w in waves])
1664 wend = min([w[-1] for w in waves])
1665 waves_ = waves_[( (wstart<=waves_) & (waves_<=wend))]
1666 if full_output:
1667 fluxes_ = []
1668 else:
1669 fluxes_ = 0.
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679 for i,(wave,flux) in enumerate(zip(waves,fluxes)):
1680 intf = interp1d(np.log10(wave),np.log10(flux),kind='linear')
1681 if full_output:
1682 fluxes_.append(10**intf(np.log10(waves_)))
1683 else:
1684 fluxes_ += 10**intf(np.log10(waves_))
1685
1686 if flux_units!='erg/cm2/s/AA/sr':
1687 fluxes_ = conversions.convert('erg/s/cm2/AA/sr',flux_units,fluxes_,wave=(waves_,'AA'),**kwargs)
1688 if wave_units!='AA':
1689 waves_ = conversions.convert('AA',wave_units,waves_,**kwargs)
1690
1691 if full_output:
1692 fluxes_ = np.vstack(fluxes_)
1693 keep = -np.isnan(np.sum(fluxes_,axis=0))
1694 waves_ = waves_[keep]
1695 fluxes_ = fluxes_[:,keep]
1696 else:
1697 keep = -np.isnan(fluxes_)
1698 waves_ = waves_[keep]
1699 fluxes_ = fluxes_[keep]
1700 return waves_,fluxes_
1701
1704 """
1705 Retrieve possible effective temperatures and gravities from a grid.
1706
1707 E.g. kurucz, sdB, fastwind...
1708
1709 @rtype: (ndarray,ndarray)
1710 @return: effective temperatures, gravities
1711 """
1712 gridfile = get_file(**kwargs)
1713 ff = pf.open(gridfile)
1714 teffs = []
1715 loggs = []
1716 for mod in ff[1:]:
1717 teffs.append(float(mod.header['TEFF']))
1718 loggs.append(float(mod.header['LOGG']))
1719 ff.close()
1720
1721
1722 matrix = np.vstack([np.array(teffs),np.array(loggs)]).T
1723 matrix = numpy_ext.sort_order(matrix,order=[0,1])
1724 teffs,loggs = matrix.T
1725
1726 return teffs,loggs
1727
1732 """
1733 Retrieve possible effective temperatures, surface gravities and reddenings
1734 from an integrated grid.
1735
1736 E.g. kurucz, sdB, fastwind...
1737
1738 @rtype: (ndarray,ndarray,ndarray)
1739 @return: effective temperatures, surface, gravities, E(B-V)s
1740 """
1741 gridfile = get_file(integrated=True,**kwargs)
1742 ff = pf.open(gridfile)
1743 teffs = ff[1].data.field('TEFF')
1744 loggs = ff[1].data.field('LOGG')
1745 ebvs = ff[1].data.field('EBV')
1746 ff.close()
1747
1748
1749
1750
1751 return teffs,loggs,ebvs
1752
1753
1754
1755
1756
1757
1758
1759 @memoized
1760 -def get_grid_mesh(wave=None,teffrange=None,loggrange=None,**kwargs):
1761 """
1762 Return InterpolatingFunction spanning the available grid of atmosphere models.
1763
1764 WARNING: the grid must be entirely defined on a mesh grid, but it does not
1765 need to be equidistant.
1766
1767 It is thus the user's responsibility to know whether the grid is evenly
1768 spaced in logg and teff (e.g. this is not so for the CMFGEN models).
1769
1770 You can supply your own wavelength range, since the grid models'
1771 resolution are not necessarily homogeneous. If not, the first wavelength
1772 array found in the grid will be used as a template.
1773
1774 It might take a long a time and cost a lot of memory if you load the entire
1775 grid. Therefor, you can also set range of temperature and gravity.
1776
1777 WARNING: 30000,50000 did not work out for FASTWIND, since we miss a model!
1778
1779 Example usage:
1780
1781 @param wave: wavelength to define the grid on
1782 @type wave: ndarray
1783 @param teffrange: starting and ending of the grid in teff
1784 @type teffrange: tuple of floats
1785 @param loggrange: starting and ending of the grid in logg
1786 @type loggrange: tuple of floats
1787 @return: wavelengths, teffs, loggs and fluxes of grid, and the interpolating
1788 function
1789 @rtype: (3x1Darray,3Darray,interp_func)
1790 """
1791
1792 teffs,loggs = get_grid_dimensions(**kwargs)
1793
1794
1795 if teffrange is not None:
1796 sa = (teffrange[0]<=teffs) & (teffs<=teffrange[1])
1797 teffs = teffs[sa]
1798 if loggrange is not None:
1799 sa = (loggrange[0]<=loggs) & (loggs<=loggrange[1])
1800 loggs = loggs[sa]
1801
1802
1803 if not new_scipy:
1804 logger.warning('SCIENTIFIC PYTHON')
1805
1806 teffs = list(set(list(teffs)))
1807 loggs = list(set(list(loggs)))
1808 teffs = np.sort(teffs)
1809 loggs = np.sort(loggs)
1810 if wave is not None:
1811 flux = np.ones((len(teffs),len(loggs),len(wave)))
1812
1813
1814 gridfile = get_file(**kwargs)
1815 ff = pf.open(gridfile)
1816 for i,teff in enumerate(teffs):
1817 for j,logg in enumerate(loggs):
1818 try:
1819 mod_name = "T%05d_logg%01.02f" %(teff,logg)
1820 mod = ff[mod_name]
1821 wave_ = mod.data.field('wavelength')
1822 flux_ = mod.data.field('flux')
1823
1824
1825
1826 if wave is None:
1827 wave = wave_
1828 flux = np.ones((len(teffs),len(loggs),len(wave)))
1829 except KeyError:
1830 continue
1831
1832
1833 try:
1834 flux[i,j,:] = flux_
1835 except:
1836 flux[i,j,:] = np.interp(wave,wave_,flux_)
1837 ff.close()
1838 flux_grid = InterpolatingFunction([np.log10(teffs),loggs],flux)
1839 logger.info('Constructed SED interpolation grid')
1840
1841
1842 else:
1843 logger.warning('SCIPY')
1844
1845
1846 gridfile = get_file(**kwargs)
1847 ff = pf.open(gridfile)
1848 if wave is not None:
1849 fluxes = np.zeros((len(teffs),len(wave)))
1850 for i,(teff,logg) in enumerate(zip(teffs,loggs)):
1851 mod_name = "T%05d_logg%01.02f" %(teff,logg)
1852 mod = ff[mod_name]
1853 wave_ = mod.data.field('wavelength')
1854 flux_ = mod.data.field('flux')
1855
1856
1857
1858 if wave is None:
1859 wave = wave_
1860 flux = np.ones((len(teffs),len(wave)))
1861 try:
1862 flux[i] = flux_
1863 except:
1864 flux[i] = np.interp(wave,wave_,flux_)
1865 ff.close()
1866 flux_grid = LinearNDInterpolator(np.array([np.log10(teffs),loggs]).T,flux)
1867 return wave,teffs,loggs,flux,flux_grid
1868
1874 """
1875 Print and return the list of calibrators
1876
1877 @parameter library: name of the library (calspec, ngsl, stelib)
1878 @type library: str
1879 @return: list of calibrator names
1880 @rtype: list of str
1881 """
1882 files = config.glob(os.path.join(caldir,library),'*.fits')
1883 targname = dict(calspec='targetid',ngsl='targname',stelib='object')[library]
1884
1885 names = []
1886 for ff in files:
1887 name = pf.getheader(ff)[targname]
1888
1889
1890
1891
1892 names.append(name)
1893 return names
1894
1895
1896
1897 -def get_calibrator(name='alpha_lyr',version=None,wave_units=None,flux_units=None,library='calspec'):
1898 """
1899 Retrieve a calibration SED
1900
1901 If C{version} is None, get the last version.
1902
1903 Example usage:
1904
1905 >>> wave,flux = get_calibrator(name='alpha_lyr')
1906 >>> wave,flux = get_calibrator(name='alpha_lyr',version='003')
1907
1908 @param name: calibrator name
1909 @type name: str
1910 @param version: version of the calibration file
1911 @type version: str
1912 @param wave_units: units of wavelength arrays (default: AA)
1913 @type wave_units: str (interpretable by C{units.conversions.convert})
1914 @param flux_units: units of flux arrays (default: erg/s/cm2/AA)
1915 @type flux_units: str (interpretable by C{units.conversions.convert})
1916 @return: wavelength and flux arrays of calibrator
1917 @rtype: (ndarray,ndarray)
1918 """
1919
1920 files = config.glob(os.path.join(caldir,library),'*.fits')
1921 targname = dict(calspec='targetid',ngsl='targname',stelib='object')[library]
1922
1923 calfile = None
1924 for ff in files:
1925
1926 fits_file = pf.open(ff)
1927 header = fits_file[0].header
1928 if name in ff or name in header[targname]:
1929
1930 if version is not None and version not in ff:
1931 fits_file.close()
1932 continue
1933
1934 calfile = ff
1935 if library in ['calspec','ngsl']:
1936 wave = fits_file[1].data.field('wavelength')
1937 flux = fits_file[1].data.field('flux')
1938 elif library in ['stelib']:
1939 wave,flux = fits.read_spectrum(ff)
1940 else:
1941 raise ValueError("Don't know what to do with files from library {}".format(library))
1942 fits_file.close()
1943
1944 if calfile is None:
1945 raise ValueError, 'Calibrator %s (version=%s) not found'%(name,version)
1946
1947 if flux_units is not None:
1948 flux = conversions.convert('erg/s/cm2/AA',flux_units,flux,wave=(wave,'AA'))
1949 if wave_units is not None:
1950 wave = conversions.convert('AA',wave_units,wave)
1951
1952
1953 logger.info('Calibrator %s selected'%(calfile))
1954
1955 return wave,flux
1956
1960 filename = config.get_datafile('sedtables/calibrators','{}.ident'.format(library))
1961 names = []
1962 fits_files = []
1963 phot_files = []
1964 with open(filename,'r') as ff:
1965 for line in ff.readlines():
1966 line = line.strip().split(',')
1967 try:
1968 fits_file = config.get_datafile('sedtables/calibrators',line[1])
1969 phot_file = config.get_datafile('sedtables/calibrators',line[2])
1970
1971 except IOError:
1972 continue
1973
1974 names.append(line[0])
1975 fits_files.append(fits_file)
1976 phot_files.append(phot_file)
1977
1978 return names,fits_files,phot_files
1979
1982 """
1983 Calibrate photometry.
1984
1985 Not finished!
1986
1987 ABmag = -2.5 Log F_nu - 48.6 with F_nu in erg/s/cm2/Hz
1988 Flux computed as 10**(-(meas-mag0)/2.5)*F0
1989 Magnitude computed as -2.5*log10(Fmeas/F0)
1990 F0 = 3.6307805477010029e-20 erg/s/cm2/Hz
1991
1992 STmag = -2.5 Log F_lam - 21.10 with F_lam in erg/s/cm2/AA
1993 Flux computed as 10**(-(meas-mag0)/2.5)*F0
1994 Magnitude computed as -2.5*log10(Fmeas/F0)
1995 F0 = 3.6307805477010028e-09 erg/s/cm2/AA
1996
1997 Vegamag = -2.5 Log F_lam - C with F_lam in erg/s/cm2/AA
1998 Flux computed as 10**(-meas/2.5)*F0
1999 Magnitude computed as -2.5*log10(Fmeas/F0)
2000 """
2001 F0ST = 3.6307805477010028e-09
2002 F0AB = 3.6307805477010029e-20
2003
2004 wave,flux = get_calibrator(name='alpha_lyr')
2005 zp = filters.get_info()
2006
2007
2008 syn_flux = synthetic_flux(wave,flux,zp['photband'])
2009 syn_flux_fnu = synthetic_flux(wave,flux,zp['photband'],units='Fnu')
2010 Flam0_lit = conversions.nconvert(zp['Flam0_units'],'erg/s/cm2/AA',zp['Flam0'],photband=zp['photband'])
2011 Fnu0_lit = conversions.nconvert(zp['Fnu0_units'],'erg/s/cm2/Hz',zp['Fnu0'],photband=zp['photband'])
2012
2013
2014 keep = (zp['Flam0_lit']==1) & (zp['Fnu0_lit']==0)
2015 Fnu0 = conversions.nconvert(zp['Flam0_units'],'erg/s/cm2/Hz',zp['Flam0'],photband=zp['photband'])
2016 zp['Fnu0'][keep] = Fnu0[keep]
2017 zp['Fnu0_units'][keep] = 'erg/s/cm2/Hz'
2018
2019
2020 keep = (zp['Flam0_lit']==0) & (zp['Fnu0_lit']==1)
2021 Flam0 = conversions.nconvert(zp['Fnu0_units'],'erg/s/cm2/AA',zp['Fnu0'],photband=zp['photband'])
2022
2023
2024 Flam0 = conversions.nconvert(zp['Flam0_units'],'erg/s/cm2/AA',zp['Flam0'])
2025 Fnu0 = conversions.nconvert(zp['Fnu0_units'],'erg/s/cm2/Hz',zp['Fnu0'])
2026
2027
2028
2029 keep = (zp['Flam0_lit']==0) & (zp['Fnu0_lit']==0)
2030 zp['Flam0'][keep] = syn_flux[keep]
2031 zp['Flam0_units'][keep] = 'erg/s/cm2/AA'
2032 zp['Fnu0'][keep] = syn_flux_fnu[keep]
2033 zp['Fnu0_units'][keep] = 'erg/s/cm2/Hz'
2034
2035 keep = np.array(['DENIS' in photb and True or False for photb in zp['photband']])
2036
2037
2038 keep = (zp['vegamag_lit']==1) & (zp['Flam0_lit']==0)
2039 zp['Flam0'][keep] = syn_flux[keep]
2040 zp['Flam0_units'][keep] = 'erg/s/cm2/AA'
2041
2042
2043 keep = (zp['STmag_lit']==1) & (zp['Flam0_lit']==0)
2044 m_vega = 2.5*np.log10(F0ST/syn_flux) + zp['STmag']
2045 zp['vegamag'][keep] = m_vega[keep]
2046
2047
2048 keep = (zp['ABmag_lit']==1) & (zp['Flam0_lit']==0)
2049 F0AB_lam = conversions.convert('erg/s/cm2/Hz','erg/s/cm2/AA',F0AB,photband=zp['photband'])
2050 m_vega = 2.5*np.log10(F0AB_lam/syn_flux) + zp['ABmag']
2051 zp['vegamag'][keep] = m_vega[keep]
2052
2053
2054 set_wave = np.isnan(zp['eff_wave'])
2055 zp['eff_wave'][set_wave] = filters.eff_wave(zp['photband'][set_wave])
2056
2057 return zp
2058
2059
2060
2061
2062
2063
2064 -def synthetic_flux(wave,flux,photbands,units=None):
2065 """
2066 Extract flux measurements from a synthetic SED (Fnu or Flambda).
2067
2068 The fluxes below 4micron are calculated assuming PHOTON-counting detectors
2069 (e.g. CCDs).
2070
2071 Flam = int(P_lam * f_lam * lam, dlam) / int(P_lam * lam, dlam)
2072
2073 When otherwise specified, we assume ENERGY-counting detectors (e.g. bolometers)
2074
2075 Flam = int(P_lam * f_lam, dlam) / int(P_lam, dlam)
2076
2077 Where P_lam is the total system dimensionless sensitivity function, which
2078 is normalised so that the maximum equals 1. Also, f_lam is the SED of the
2079 object, in units of energy per time per unit area per wavelength.
2080
2081 The PHOTON-counting part of this routine has been thoroughly checked with
2082 respect to johnson UBV, geneva and stromgren filters, and only gives offsets
2083 with respect to the Kurucz integrated files (.geneva and stuff on his websites). These could be
2084 due to different normalisation.
2085
2086 You can also readily integrate in Fnu instead of Flambda by suppling a list
2087 of strings to 'units'. This should have equal length of photbands, and
2088 should contain the strings 'flambda' and 'fnu' corresponding to each filter.
2089 In that case, the above formulas reduce to
2090
2091 Fnu = int(P_nu * f_nu / nu, dnu) / int(P_nu / nu, dnu)
2092
2093 and
2094
2095 Fnu = int(P_nu * f_nu, dnu) / int(P_nu, dnu)
2096
2097 Small note of caution: P_nu is not equal to P_lam according to
2098 Maiz-Apellaniz, he states that P_lam = P_nu/lambda. But in the definition
2099 we use above here, it *is* the same!
2100
2101 The model fluxes should B{always} be given in Flambda (erg/s/cm2/AA). The
2102 program will convert them to Fnu where needed.
2103
2104 The output is a list of numbers, equal in length to the 'photband' inputs.
2105 The units of the output are erg/s/cm2/AA where Flambda was given, and
2106 erg/s/cm2/Hz where Fnu was given.
2107
2108 The difference is only marginal for 'blue' bands. For example, integrating
2109 2MASS in Flambda or Fnu is only different below the 1.1% level:
2110
2111 >>> wave,flux = get_table(teff=10000,logg=4.0)
2112 >>> energys = synthetic_flux(wave,flux,['2MASS.J','2MASS.J'],units=['flambda','fnu'])
2113 >>> e0_conv = conversions.convert('erg/s/cm2/AA','erg/s/cm2/Hz',energys[0],photband='2MASS.J')
2114 >>> np.abs(energys[1]-e0_conv)/energys[1]<0.012
2115 True
2116
2117 But this is not the case for IRAS.F12:
2118
2119 >>> energys = synthetic_flux(wave,flux,['IRAS.F12','IRAS.F12'],units=['flambda','fnu'])
2120 >>> e0_conv = conversions.convert('erg/s/cm2/AA','erg/s/cm2/Hz',energys[0],photband='IRAS.F12')
2121 >>> np.abs(energys[1]-e0_conv)/energys[1]>0.1
2122 True
2123
2124 If you have a spectrum in micron vs Jy and want to calculate the synthetic
2125 fluxes in Jy, a little bit more work is needed to get everything in the
2126 right units. In the following example, we first generate a constant flux
2127 spectrum in micron and Jy. Then, we convert flux to erg/s/cm2/AA using the
2128 wavelengths (this is no approximation) and convert wavelength to angstrom.
2129 Next, we compute the synthetic fluxes in the IRAS band in Fnu, and finally
2130 convert the outcome (in erg/s/cm2/Hz) to Jansky.
2131
2132 >>> wave,flux = np.linspace(0.1,200,10000),np.ones(10000)
2133 >>> flam = conversions.convert('Jy','erg/s/cm2/AA',flux,wave=(wave,'micron'))
2134 >>> lam = conversions.convert('micron','AA',wave)
2135 >>> energys = synthetic_flux(lam,flam,['IRAS.F12','IRAS.F25','IRAS.F60','IRAS.F100'],units=['Fnu','Fnu','Fnu','Fnu'])
2136 >>> energys = conversions.convert('erg/s/cm2/Hz','Jy',energys)
2137
2138 You are responsible yourself for having a response curve covering the
2139 model fluxes!
2140
2141 Now. let's put this all in practice in a more elaborate example: we want
2142 to check if the effective wavelength is well defined. To do that we will:
2143
2144 1. construct a model (black body)
2145 2. make our own weird, double-shaped filter (CCD-type and BOL-type detector)
2146 3. compute fluxes in Flambda, and convert to Fnu via the effective wavelength
2147 4. compute fluxes in Fnu, and compare with step 3.
2148
2149 In an ideal world, the outcome of step (3) and (4) must be equal:
2150
2151 Step (1): We construct a black body model.
2152
2153
2154 WARNING: OPEN.BOL only works in Flambda for now.
2155
2156 See e.g. Maiz-Apellaniz, 2006.
2157
2158 @param wave: model wavelengths (angstrom)
2159 @type wave: ndarray
2160 @param flux: model fluxes (erg/s/cm2/AA)
2161 @type flux: ndarray
2162 @param photbands: list of photometric passbands
2163 @type photbands: list of str
2164 @param units: list containing Flambda or Fnu flag (defaults to all Flambda)
2165 @type units: list of strings or str
2166 @return: model fluxes (erg/s/cm2/AA or erg/s/cm2/Hz)
2167 @rtype: ndarray
2168 """
2169 if isinstance(units,str):
2170 units = [units]*len(photbands)
2171 energys = np.zeros(len(photbands))
2172
2173
2174 filter_info = filters.get_info()
2175 keep = np.searchsorted(filter_info['photband'],photbands)
2176 filter_info = filter_info[keep]
2177
2178 for i,photband in enumerate(photbands):
2179
2180 waver,transr = filters.get_response(photband)
2181
2182
2183
2184 region = ((waver[0]-0.4*waver[0])<=wave) & (wave<=(2*waver[-1]))
2185
2186
2187
2188 if filter_info['eff_wave'][i]>=4e4 and sum(region)<1e5 and sum(region)>1:
2189 logger.debug('%10s: Interpolating model to integrate over response curve'%(photband))
2190 wave_ = np.logspace(np.log10(wave[region][0]),np.log10(wave[region][-1]),1e5)
2191 flux_ = 10**np.interp(np.log10(wave_),np.log10(wave[region]),np.log10(flux[region]),)
2192 else:
2193 wave_ = wave[region]
2194 flux_ = flux[region]
2195 if not len(wave_):
2196 energys[i] = np.nan
2197 continue
2198
2199
2200 if (np.searchsorted(wave_,waver[-1])-np.searchsorted(wave_,waver[0]))<5:
2201 wave__ = np.sort(np.hstack([wave_,waver]))
2202 flux_ = np.interp(wave__,wave_,flux_)
2203 wave_ = wave__
2204
2205 transr = np.interp(wave_,waver,transr,left=0,right=0)
2206
2207
2208
2209 if units is None or ((units is not None) and (units[i].upper()=='FLAMBDA')):
2210 if photband=='OPEN.BOL':
2211 energys[i] = np.trapz(flux_,x=wave_)
2212 elif filter_info['type'][i]=='BOL':
2213 energys[i] = np.trapz(flux_*transr,x=wave_)/np.trapz(transr,x=wave_)
2214 elif filter_info['type'][i]=='CCD':
2215 energys[i] = np.trapz(flux_*transr*wave_,x=wave_)/np.trapz(transr*wave_,x=wave_)
2216
2217
2218 elif units[i].upper()=='FNU':
2219
2220 freq_ = conversions.convert('AA','Hz',wave_)
2221 flux_f = conversions.convert('erg/s/cm2/AA','erg/s/cm2/Hz',flux_,wave=(wave_,'AA'))
2222
2223 sa = np.argsort(freq_)
2224 transr = transr[sa]
2225 freq_ = freq_[sa]
2226 flux_f = flux_f[sa]
2227 if filter_info['type'][i]=='BOL':
2228 energys[i] = np.trapz(flux_f*transr,x=freq_)/np.trapz(transr,x=freq_)
2229 elif filter_info['type'][i]=='CCD':
2230 energys[i] = np.trapz(flux_f*transr/freq_,x=wave_)/np.trapz(transr/freq_,x=wave_)
2231 else:
2232 raise ValueError,'units %s not understood'%(units)
2233
2234
2235 return energys
2236
2239 """
2240 Construct colors from a synthetic SED.
2241
2242 @param wave: model wavelengths (angstrom)
2243 @type wave: ndarray
2244 @param flux: model fluxes (erg/s/cm2/AA)
2245 @type flux: ndarray
2246 @param colors: list of photometric passbands
2247 @type colors: list of str
2248 @param units: list containing Flambda or Fnu flag (defaults to all Flambda)
2249 @type units: list of strings or str
2250 @return: flux ratios or colors
2251 @rtype: ndarray
2252 """
2253 if units is None:
2254 units = [None for color in colors]
2255
2256 syn_colors = np.zeros(len(colors))
2257 for i,(color,unit) in enumerate(zip(colors,units)):
2258
2259
2260 photbands,color_func = filters.make_color(color)
2261
2262 fluxes = synthetic_flux(wave,flux,photbands,units=unit)
2263
2264 syn_colors[i] = color_func(*list(fluxes))
2265
2266 return syn_colors
2267
2271 """
2272 Calculate the bolometric luminosity of a model SED.
2273
2274 Flux should be in cgs per unit wavelength (same unit as wave).
2275 The latter is integrated out, so it is of no importance. After integration,
2276 flux, should have units erg/s/cm2.
2277
2278 Returned luminosity is in solar units.
2279
2280 If you give radius=1 and want to correct afterwards, multiply the obtained
2281 Labs with radius**2.
2282
2283 @param wave: model wavelengths
2284 @type wave: ndarray
2285 @param flux: model fluxes (Flam)
2286 @type flux: ndarray
2287 @param radius: stellar radius in solar units
2288 @type radius: float
2289 @return: total bolometric luminosity
2290 @rtype: float
2291 """
2292 Lint = np.trapz(flux,x=wave)
2293 Labs = Lint*4*np.pi/constants.Lsol_cgs*(radius*constants.Rsol_cgs)**2
2294 return Labs
2295
2296
2297
2298 @memoized
2299 -def _get_itable_markers(photbands,
2300 teffrange=(-np.inf,np.inf),loggrange=(-np.inf,np.inf),
2301 ebvrange=(-np.inf,np.inf),zrange=(-np.inf,np.inf),
2302 include_Labs=True,clear_memory=True,**kwargs):
2303 """
2304 Get a list of markers to more easily retrieve integrated fluxes.
2305 """
2306 if clear_memory:
2307 clear_memoization(keys=['ivs.sed.model'])
2308 gridfiles = get_file(z='*',integrated=True,**kwargs)
2309 if isinstance(gridfiles,str):
2310 gridfiles = [gridfiles]
2311
2312 metals_sa = np.argsort([pf.getheader(ff,1)['z'] for ff in gridfiles])
2313 gridfiles = np.array(gridfiles)[metals_sa]
2314 flux = []
2315 gridpnts = []
2316 grid_z = []
2317 markers = []
2318
2319
2320 for gridfile in gridfiles:
2321 ff = pf.open(gridfile)
2322 ext = ff[1]
2323 z = ff[1].header['z']
2324 if z<zrange[0] or zrange[1]<z:
2325 continue
2326
2327 teffs = ext.data.field('teff')
2328 loggs = ext.data.field('logg')
2329 ebvs = ext.data.field('ebv')
2330 keep = (ebvrange[0]<=ebvs) & (ebvs<=ebvrange[1])
2331
2332
2333
2334
2335
2336
2337 teffs,loggs,ebvs = teffs[keep],loggs[keep],ebvs[keep]
2338 grid_teffs = np.sort(list(set(teffs)))
2339 grid_loggs = np.sort(list(set(loggs)))
2340 grid_ebvs = np.sort(list(set(ebvs)))
2341 grid_z.append(z)
2342
2343
2344
2345
2346
2347 markers.append(np.zeros(len(teffs)))
2348 gridpnts.append(np.zeros((len(teffs),4)))
2349
2350 for i,(it,il,ie) in enumerate(zip(teffs,loggs,ebvs)):
2351 markers[-1][i] = float('%3d%05d%03d%03d'%(int(round((z+5)*100)),int(round(it)),int(round(il*100)),int(round(ie*100))))
2352 gridpnts[-1][i]= it,il,ie,z
2353 flux.append(_get_flux_from_table(ext,photbands,include_Labs=include_Labs))
2354 ff.close()
2355
2356 flux = np.vstack(flux)
2357 markers = np.hstack(markers)
2358 gridpnts = np.vstack(gridpnts)
2359 grid_z = np.sort(grid_z)
2360
2361 return np.array(markers),(grid_teffs,grid_loggs,grid_ebvs,grid_z),gridpnts,flux
2362
2363
2364 @memoized
2365 -def _get_pix_grid(photbands,
2366 teffrange=(-np.inf,np.inf),loggrange=(-np.inf,np.inf),
2367 ebvrange=(-np.inf,np.inf),zrange=(-np.inf,np.inf),
2368 rvrange=(-np.inf,np.inf),vradrange=(-np.inf,np.inf),
2369 include_Labs=True,clear_memory=True,
2370 variables=['teff','logg','ebv','z','rv','vrad'],**kwargs):
2371 """
2372 Prepare the pixalted grid.
2373
2374 In principle, it should be possible to return any number of free parameters
2375 here. I'm thinking about:
2376
2377 teff, logg, ebv, z, Rv, vrad.
2378 """
2379 if clear_memory:
2380 clear_memoization(keys=['ivs.sed.model'])
2381
2382
2383 trash = kwargs.pop('Rv', '*')
2384 trash = kwargs.pop('z', '*')
2385 gridfiles = get_file(z='*',Rv='*',integrated=True,**kwargs)
2386 if isinstance(gridfiles,str):
2387 gridfiles = [gridfiles]
2388 flux = []
2389 grid_pars = []
2390 grid_names = np.array(variables)
2391
2392 for gridfile in gridfiles:
2393 with pf.open(gridfile) as ff:
2394
2395 ext = ff[1]
2396
2397 keep = np.ones(len(ext.data),bool)
2398 for name in variables:
2399
2400 low,high = locals()[name+'range']
2401 in_range = (low<=ext.data.field(name)) & (ext.data.field(name)<=high)
2402 on_edge = np.allclose(ext.data.field(name),low) | np.allclose(ext.data.field(name),high)
2403
2404
2405
2406
2407
2408
2409
2410 keep = keep & (in_range | on_edge)
2411 partial_grid = np.vstack([ext.data.field(name)[keep] for name in variables])
2412 if sum(keep):
2413 grid_pars.append(partial_grid)
2414
2415 flux.append(_get_flux_from_table(ext,photbands,include_Labs=include_Labs)[keep])
2416
2417 flux = np.vstack(flux)
2418 grid_pars = np.hstack(grid_pars)
2419
2420
2421
2422 flux = np.log10(flux)
2423
2424
2425 keep = np.ones(len(grid_names),bool)
2426 for i in range(len(grid_names)):
2427 if np.all(grid_pars[i]==grid_pars[i][0]):
2428 keep[i] = False
2429 grid_pars = grid_pars[keep]
2430
2431
2432 grid_names = grid_names[keep]
2433
2434
2435 axis_values, pixelgrid = interpol.create_pixeltypegrid(grid_pars,flux.T)
2436 return axis_values,grid_pars.T,pixelgrid,grid_names
2437
2442 """
2443 Retrieve flux and flux ratios from an integrated SED table.
2444
2445 @param fits_ext: fits extension containing integrated flux
2446 @type fits_ext: FITS extension
2447 @param photbands: list of photometric passbands
2448 @type photbands: list of str
2449 @param index: slice or index of rows to retrieve
2450 @type index: slice or integer
2451 @return: fluxes or flux ratios
2452 #@rtype: list
2453 """
2454 if index is None:
2455 index = slice(None)
2456 fluxes = []
2457 for photband in photbands:
2458 try:
2459 if not filters.is_color(photband):
2460 fluxes.append(fits_ext.data.field(photband)[index])
2461 else:
2462 system,color = photband.split('.')
2463 if '-' in color:
2464 band0,band1 = color.split('-')
2465 fluxes.append(fits_ext.data.field('%s.%s'%(system,band0))[index]/fits_ext.data.field('%s.%s'%(system,band1))[index])
2466 elif color=='M1':
2467 fv = fits_ext.data.field('STROMGREN.V')[index]
2468 fy = fits_ext.data.field('STROMGREN.Y')[index]
2469 fb = fits_ext.data.field('STROMGREN.B')[index]
2470 fluxes.append(fv*fy/fb**2)
2471 elif color=='C1':
2472 fu = fits_ext.data.field('STROMGREN.U')[index]
2473 fv = fits_ext.data.field('STROMGREN.V')[index]
2474 fb = fits_ext.data.field('STROMGREN.B')[index]
2475 fluxes.append(fu*fb/fv**2)
2476 except KeyError:
2477 logger.warning('Passband %s missing from table'%(photband))
2478 fluxes.append(np.nan*np.ones(len(fits_ext.data)))
2479
2480 if include_Labs:
2481 fluxes.append(fits_ext.data.field("Labs")[index])
2482 fluxes = np.array(fluxes).T
2483 if index is not None:
2484 fluxes = fluxes
2485 return fluxes
2486
2487
2488
2489 if __name__=="__main__":
2490 import doctest
2491 import pylab as pl
2492 doctest.testmod()
2493 pl.show()
2494