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