1
2 """
3 Functions relevant for photometric calibration
4
5 Table of contents:
6
7 1. Available response functions
8 2. Adding filters on the fly
9 - Defining a new filter
10 - Temporarily modifying an existing filter
11 3. Adding filters permanently
12
13 Section 1. Available response functions
14 =======================================
15
16 Short list of available systems:
17
18 >>> responses = list_response()
19 >>> systems = [response.split('.')[0] for response in responses]
20 >>> set_responses = sorted(set([response.split('.')[0] for response in systems]))
21 >>> for i,resp in enumerate(set_responses):
22 ... print '%10s (%3d filters)'%(resp,systems.count(resp))
23 2MASS ( 3 filters)
24 ACSHRC ( 17 filters)
25 ACSSBC ( 6 filters)
26 ACSWFC ( 12 filters)
27 AKARI ( 13 filters)
28 ANS ( 6 filters)
29 APEX ( 1 filters)
30 ARGUE ( 3 filters)
31 BESSEL ( 6 filters)
32 BESSELL ( 6 filters)
33 COROT ( 2 filters)
34 COUSINS ( 3 filters)
35 DDO ( 7 filters)
36 DENIS ( 3 filters)
37 DIRBE ( 10 filters)
38 EEV4280 ( 1 filters)
39 ESOIR ( 10 filters)
40 GAIA ( 4 filters)
41 GALEX ( 2 filters)
42 GENEVA ( 7 filters)
43 HIPPARCOS ( 1 filters)
44 IPHAS ( 3 filters)
45 IRAC ( 4 filters)
46 IRAS ( 4 filters)
47 ISOCAM ( 21 filters)
48 JOHNSON ( 25 filters)
49 KEPLER ( 43 filters)
50 KRON ( 2 filters)
51 LANDOLT ( 6 filters)
52 MIPS ( 3 filters)
53 MOST ( 1 filters)
54 MSX ( 6 filters)
55 NARROW ( 1 filters)
56 NICMOS ( 6 filters)
57 OAO2 ( 12 filters)
58 PACS ( 3 filters)
59 SAAO ( 13 filters)
60 SCUBA ( 6 filters)
61 SDSS ( 10 filters)
62 SLOAN ( 2 filters)
63 SPIRE ( 3 filters)
64 STEBBINS ( 6 filters)
65 STISCCD ( 2 filters)
66 STISFUV ( 4 filters)
67 STISNUV ( 7 filters)
68 STROMGREN ( 6 filters)
69 TD1 ( 4 filters)
70 TYCHO ( 2 filters)
71 TYCHO2 ( 2 filters)
72 ULTRACAM ( 5 filters)
73 USNOB1 ( 2 filters)
74 UVEX ( 5 filters)
75 VILNIUS ( 7 filters)
76 VISIR ( 13 filters)
77 WALRAVEN ( 5 filters)
78 WFCAM ( 5 filters)
79 WFPC2 ( 21 filters)
80 WISE ( 4 filters)
81 WOOD ( 12 filters)
82
83 Plots of all passbands of all systems:
84
85 ]include figure]]ivs_sed_filters_2MASS.png]
86
87 ]include figure]]ivs_sed_filters_ACSHRC.png]
88
89 ]include figure]]ivs_sed_filters_ACSSBC.png]
90
91 ]include figure]]ivs_sed_filters_ACSWFC.png]
92
93 ]include figure]]ivs_sed_filters_AKARI.png]
94
95 ]include figure]]ivs_sed_filters_ANS.png]
96
97 ]include figure]]ivs_sed_filters_APEX.png]
98
99 ]include figure]]ivs_sed_filters_ARGUE.png]
100
101 ]include figure]]ivs_sed_filters_BESSEL.png]
102
103 ]include figure]]ivs_sed_filters_BESSELL.png]
104
105 ]include figure]]ivs_sed_filters_COROT.png]
106
107 ]include figure]]ivs_sed_filters_COUSINS.png]
108
109 ]include figure]]ivs_sed_filters_DDO.png]
110
111 ]include figure]]ivs_sed_filters_DENIS.png]
112
113 ]include figure]]ivs_sed_filters_DIRBE.png]
114
115 ]include figure]]ivs_sed_filters_ESOIR.png]
116
117 ]include figure]]ivs_sed_filters_EEV4280.png]
118
119 ]include figure]]ivs_sed_filters_GAIA.png]
120
121 ]include figure]]ivs_sed_filters_GALEX.png]
122
123 ]include figure]]ivs_sed_filters_GENEVA.png]
124
125 ]include figure]]ivs_sed_filters_HIPPARCOS.png]
126
127 ]include figure]]ivs_sed_filters_IPHAS.png]
128
129 ]include figure]]ivs_sed_filters_IRAC.png]
130
131 ]include figure]]ivs_sed_filters_IRAS.png]
132
133 ]include figure]]ivs_sed_filters_ISOCAM.png]
134
135 ]include figure]]ivs_sed_filters_JOHNSON.png]
136
137 ]include figure]]ivs_sed_filters_KEPLER.png]
138
139 ]include figure]]ivs_sed_filters_KRON.png]
140
141 ]include figure]]ivs_sed_filters_LANDOLT.png]
142
143 ]include figure]]ivs_sed_filters_MIPS.png]
144
145 ]include figure]]ivs_sed_filters_MOST.png]
146
147 ]include figure]]ivs_sed_filters_MSX.png]
148
149 ]include figure]]ivs_sed_filters_NARROW.png]
150
151 ]include figure]]ivs_sed_filters_NICMOS.png]
152
153 ]include figure]]ivs_sed_filters_OAO2.png]
154
155 ]include figure]]ivs_sed_filters_PACS.png]
156
157 ]include figure]]ivs_sed_filters_SAAO.png]
158
159 ]include figure]]ivs_sed_filters_SCUBA.png]
160
161 ]include figure]]ivs_sed_filters_SDSS.png]
162
163 ]include figure]]ivs_sed_filters_SLOAN.png]
164
165 ]include figure]]ivs_sed_filters_SPIRE.png]
166
167 ]include figure]]ivs_sed_filters_STEBBINS.png]
168
169 ]include figure]]ivs_sed_filters_STISCCD.png]
170
171 ]include figure]]ivs_sed_filters_STISFUV.png]
172
173 ]include figure]]ivs_sed_filters_STISNUV.png]
174
175 ]include figure]]ivs_sed_filters_STROMGREN.png]
176
177 ]include figure]]ivs_sed_filters_TD1.png]
178
179 ]include figure]]ivs_sed_filters_TYCHO.png]
180
181 ]include figure]]ivs_sed_filters_TYCHO2.png]
182
183 ]include figure]]ivs_sed_filters_ULTRACAM.png]
184
185 ]include figure]]ivs_sed_filters_USNOB1.png]
186
187 ]include figure]]ivs_sed_filters_UVEX.png]
188
189 ]include figure]]ivs_sed_filters_VILNIUS.png]
190
191 ]include figure]]ivs_sed_filters_VISIR.png]
192
193 ]include figure]]ivs_sed_filters_WALRAVEN.png]
194
195 ]include figure]]ivs_sed_filters_WFPC2.png]
196
197 ]include figure]]ivs_sed_filters_WISE.png]
198
199 ]include figure]]ivs_sed_filters_WOOD.png]
200
201 Section 2: Adding filters on the fly
202 ====================================
203
204 Section 2.1: Defining a new filter
205 ----------------------------------
206
207 You can add custom filters on the fly using L{add_custom_filter}. In this
208 example we add a weird-looking filter and check the definition of Flambda and
209 Fnu and its relation to the effective wavelength of a passband:
210
211 Prerequisites: some modules that come in handy:
212
213 >>> from ivs.sigproc import funclib
214 >>> from ivs.sed import model
215 >>> from ivs.units import conversions
216
217 First, we'll define a double peakd Gaussian profile on the wavelength grid of
218 the WISE.W3 response curve:
219
220 >>> wave = get_response('WISE.W3')[0]
221 >>> trans = funclib.evaluate('gauss',wave,[1.5,76000.,10000.,0.])
222 >>> trans+= funclib.evaluate('gauss',wave,[1.0,160000.,25000.,0.])
223
224 This is what it looks like:
225
226 >>> p = pl.figure()
227 >>> p = pl.plot(wave/1e4,trans,'k-')
228 >>> p = pl.xlabel("Wavelength [micron]")
229 >>> p = pl.ylabel("Transmission [arbitrary units]")
230
231 ]include figure]]ivs_sed_filters_weird_trans.png]
232
233 We can add this filter to the list of predefined filters in the following way
234 (for the doctests to work, we have to do a little work around and call
235 filters via that module, this is not needed in a normal workflow):
236
237 >>> model.filters.add_custom_filter(wave,trans,photband='LAMBDA.CCD',type='CCD')
238 >>> model.filters.add_custom_filter(wave,trans,photband='LAMBDA.BOL',type='BOL')
239
240 Note that we add the filter twice, once assuming that it is mounted on a
241 bolometer, and once on a CCD device. We'll call the filter C{LAMBDA.CCD} and
242 C{LAMBDA.BOL}. From now on, they are available within functions as L{get_info}
243 and L{get_response}. For example, what is the effective (actually pivot)
244 wavelength?
245
246 >>> effwave_ccd = model.filters.eff_wave('LAMBDA.CCD')
247 >>> effwave_bol = model.filters.eff_wave('LAMBDA.BOL')
248 >>> print(effwave_ccd,effwave_bol)
249 (119263.54911400242, 102544.27931275869)
250
251 Let's do some synthetic photometry now. Suppose we have a black body atmosphere:
252
253 >>> bb = model.blackbody(wave,5777.)
254
255 We now calculate the synthetic flux, assuming the CCD and BOL device. We
256 compute the synthetic flux both in Flambda and Fnu:
257
258 >>> flam_ccd,flam_bol = model.synthetic_flux(wave,bb,['LAMBDA.CCD','LAMBDA.BOL'])
259 >>> fnu_ccd,fnu_bol = model.synthetic_flux(wave,bb,['LAMBDA.CCD','LAMBDA.BOL'],units=['FNU','FNU'])
260
261 You can see that the fluxes can be quite different when you assume photon or
262 energy counting devices!
263
264 >>> flam_ccd,flam_bol
265 (897.68536911320564, 1495.248213834755)
266 >>> fnu_ccd,fnu_bol
267 (4.2591095543803019e-06, 5.2446332430111098e-06)
268
269 Can we now readily convert Flambda to Fnu with assuming the pivot wavelength?
270
271 >>> fnu_fromflam_ccd = conversions.convert('erg/s/cm2/AA','erg/s/cm2/Hz',flam_ccd,wave=(effwave_ccd,'A'))
272 >>> fnu_fromflam_bol = conversions.convert('erg/s/cm2/AA','erg/s/cm2/Hz',flam_bol,wave=(effwave_bol,'A'))
273
274 Which is equivalent with:
275
276 >>> fnu_fromflam_ccd = conversions.convert('erg/s/cm2/AA','erg/s/cm2/Hz',flam_ccd,photband='LAMBDA.CCD')
277 >>> fnu_fromflam_bol = conversions.convert('erg/s/cm2/AA','erg/s/cm2/Hz',flam_bol,photband='LAMBDA.BOL')
278
279 Apparently, with the definition of pivot wavelength, you can safely convert from
280 Fnu to Flambda:
281
282 >>> print(fnu_ccd,fnu_fromflam_ccd)
283 (4.2591095543803019e-06, 4.259110447428463e-06)
284 >>> print(fnu_bol,fnu_fromflam_bol)
285 (5.2446332430111098e-06, 5.2446373530017525e-06)
286
287 The slight difference you see is numerical.
288
289 Section 2.2: Temporarily modifying an existing filter
290 -----------------------------------------------------
291
292 Under usual conditions, you are prohibited from overwriting an existing predefined
293 response curve. That is, if you try to L{add_custom_filter} with a C{photband}
294 that already exists as a file, a C{ValueError} will be raised (this is not the
295 case for a custom defined filter, which you can overwrite without problems!).
296 If, for testing purposes, you want to use another definition of a predefined
297 response curve, you need to set C{force=True} in L{add_custom_filter}, and then
298 call
299
300 >>> set_prefer_file(False)
301
302 To reset and use the original definitions again, do
303
304 >>> set_prefer_file(True)
305
306 Section 3.: Adding filters permanently
307 --------------------------------------
308
309 Add a new response curve file to the ivs/sed/filters directory. The file should
310 contain two columns, the first column is the wavelength in angstrom, the second
311 column is the transmission curve. The units of the later are not important.
312
313 Then, call L{update_info}. The contents of C{zeropoints.dat} will automatically
314 be updated. Make sure to add any additional information on the new filters
315 manually in that file (e.g. is t a CCD or bolometer, what are the zeropoint
316 magnitudes etc).
317
318 """
319 import os
320 import glob
321
322 import logging
323 import numpy as np
324
325 from ivs import config
326 from ivs.aux.decorators import memoized
327 from ivs.aux import decorators
328 from ivs.aux import loggers
329 from ivs.inout import ascii
330
331 basedir = os.path.dirname(__file__)
332
333 logger = logging.getLogger("SED.FILT")
334 logger.addHandler(loggers.NullHandler())
335
336 custom_filters = {'_prefer_file':True}
341 """
342 Retrieve the response curve of a photometric system 'SYSTEM.FILTER'
343
344 OPEN.BOL represents a bolometric open filter.
345
346 Example usage:
347
348 >>> p = pl.figure()
349 >>> for band in ['J','H','KS']:
350 ... p = pl.plot(*get_response('2MASS.%s'%(band)))
351
352 If you defined a custom filter with the same name as an existing one and
353 you want to use that one in the future, set C{prefer_file=False} in the
354 C{custom_filters} module dictionary.
355
356 @param photband: photometric passband
357 @type photband: str ('SYSTEM.FILTER')
358 @return: (wavelength [A], response)
359 @rtype: (array, array)
360 """
361 photband = photband.upper()
362 prefer_file = custom_filters['_prefer_file']
363 if photband=='OPEN.BOL':
364 return np.array([1,1e10]),np.array([1/(1e10-1),1/(1e10-1)])
365
366 photfile = os.path.join(basedir,'filters',photband)
367 photfile_is_file = os.path.isfile(photfile)
368
369 if photfile_is_file and prefer_file:
370 wave, response = ascii.read2array(photfile).T[:2]
371
372 elif photband in custom_filters:
373 wave, response = custom_filters[photband]['response']
374
375 elif photfile_is_file:
376 wave, response = ascii.read2array(photfile).T[:2]
377 else:
378 raise IOError,('{0} does not exist {1}'.format(photband,custom_filters.keys()))
379 sa = np.argsort(wave)
380 return wave[sa],response[sa]
381
383 """
384 Create a custom filter as a sum of Gaussians.
385
386 @param wave: wavelength to evaluate the profile on
387 @type wave: ndarray
388 @param peaks: heights of the peaks
389 @type peaks: ndarray of length N, with N peaks
390 @param range: wavelength range of the peaks
391 @type range: tuple
392 @param sigma: width of the peaks in units of (range/N)
393 @type sigma: float
394 @return: filter profile
395 @rtype: ndarray
396 """
397 wpnts = np.linspace(range[0],range[1],len(peaks)+2)[1:-1]
398 sigma = (range[1]-range[0])/(sigma*len(peaks))
399 gauss = lambda x,p: p[0] * np.exp( -(x-p[1])**2 / (2.0*p[2]**2))
400 els = [gauss(wave,[pk,mu,sigma]) for pk,mu in zip(peaks,wpnts)]
401 profile = np.array(els).sum(axis=0)
402 return profile
403
406 """
407 Add a custom filter to the set of predefined filters.
408
409 Extra keywords are:
410 'eff_wave', 'type',
411 'vegamag', 'vegamag_lit',
412 'ABmag', 'ABmag_lit',
413 'STmag', 'STmag_lit',
414 'Flam0', 'Flam0_units', 'Flam0_lit',
415 'Fnu0', 'Fnu0_units', 'Fnu0_lit',
416 'source'
417
418 default C{type} is 'CCD'.
419 default C{photband} is 'CUSTOM.FILTER'
420
421 @param wave: wavelength (angstrom)
422 @type wave: ndarray
423 @param response: response
424 @type response: ndarray
425 @param photband: photometric passband
426 @type photband: str ('SYSTEM.FILTER')
427 """
428 kwargs.setdefault('photband','CUSTOM.FILTER')
429 kwargs.setdefault('copy_from','JOHNSON.V')
430 kwargs.setdefault('force',False)
431 photband = kwargs['photband']
432
433 photfile = os.path.join(basedir,'filters',photband)
434 if os.path.isfile(photfile) and not kwargs['force']:
435 raise ValueError,'bandpass {0} already exists'.format(photfile)
436 elif photband in custom_filters:
437 logger.debug('Overwriting previous definition of {0}'.format(photband))
438 custom_filters[photband] = dict(response=(wave,response))
439
440 kwargs.setdefault('type','CCD')
441 kwargs.setdefault('eff_wave',eff_wave(photband,det_type=kwargs['type']))
442
443
444
445 myrow = get_info([kwargs['copy_from']])
446 for name in myrow.dtype.names:
447 if 'lit' in name:
448 myrow[name] = 0
449 myrow[name] = kwargs.pop(name,myrow[name])
450 del decorators.memory[__name__]
451
452 custom_filters[photband]['zp'] = myrow
453 logger.debug('Added photband {0} to the predefined set'.format(photband))
454
456 """
457 Set whether to prefer files or custom filters when both exist.
458
459 @param prefer_file: boolean
460 @type prefer_file: bool
461 """
462 custom_filters['_prefer_file'] = prefer_file
463 logger.info("Prefering {}".format(prefer_file and 'files' or 'custom filters'))
464
467
468 Delta = np.log10(1.+1./R)
469 x = np.arange(np.log10(lambda0),np.log10(lambdan)+Delta,Delta)
470 x = 10**x
471 photbands = []
472 for i in range(len(x)-1):
473 wave = np.linspace(x[i],x[i+1],100)
474 resp = np.ones_like(wave)
475 dw = wave[1]-wave[0]
476 wave = np.hstack([wave[0]-dw,wave,wave[-1]+dw])
477 resp = np.hstack([0,resp,0])
478 photband = 'BOXCAR_R{0:d}.{1:d}'.format(int(R),int(x[i]))
479 try:
480 add_custom_filter(wave,resp,photband=photband)
481 except ValueError:
482 logger.info('{0} already exists, skipping'.format(photband))
483 photbands.append(photband)
484 logger.info('Added spectrophotometric filters')
485 return photbands
486
489 """
490 List available response curves.
491
492 Specify a glob string C{name} and/or a wavelength range to make a selection
493 of all available curves. If nothing is supplied, all curves will be returned.
494
495 @param name: list all curves containing this string
496 @type name: str
497 @param wave_range: list all curves within this wavelength range (A)
498 @type wave_range: (float, float)
499 @return: list of curve files
500 @rtype: list of str
501 """
502
503 if not '*' in name:
504 name_ = '*' + name + '*'
505 else:
506 name_ = name
507 curve_files = sorted(glob.glob(os.path.join(basedir,'filters',name_.upper())))
508 curve_files = sorted(curve_files+[key for key in custom_filters.keys() if ((name in key) and not (key=='_prefer_file'))])
509 curve_files = [cf for cf in curve_files if not ('HUMAN' in cf or 'EYE' in cf) ]
510
511 curve_files = [os.path.basename(curve_file) for curve_file in curve_files if (wave_range[0]<=eff_wave(os.path.basename(curve_file))<=wave_range[1])]
512
513 for curve_file in curve_files: logger.info(curve_file)
514 return curve_files
515
519 """
520 Return true if the photometric passband is actually a color.
521
522 @param photband: name of the photometric passband
523 @type photband: string
524 @return: True or False
525 @rtype: bool
526 """
527 if '-' in photband.split('.')[1]:
528 return True
529 elif photband.split('.')[1].upper() in ['M1','C1']:
530 return True
531 else:
532 return False
533
536 """
537 Retrieve the photometric bands from color
538
539 @param photband: name of the photometric passband
540 @type photband: string
541 @return: tuple of strings
542 @rtype: tuple
543 """
544 system,band = photband.split('.')
545 band = band.strip()
546 if '-' in band:
547 bands = tuple(['%s.%s'%(system,iband) for iband in band.split('-')])
548 elif band.upper()=='M1':
549 bands = tuple(['%s.%s'%(system,iband) for iband in ['V','B','Y']])
550 elif band.upper()=='C1':
551 bands = tuple(['%s.%s'%(system,iband) for iband in ['V','B','U']])
552 else:
553 raise ValueError('cannot recognize color {}'.format(photband))
554 return bands
555
557 """
558 Make a color from a color name and fluxes.
559
560 You get two things: a list of photbands that are need to construct the color,
561 and a function which you need to pass fluxes to compute the color.
562
563 >>> bands, func = make_color('JOHNSON.B-V')
564 >>> print(bands)
565 ('JOHNSON.B', 'JOHNSON.V')
566 >>> print(func(2,3.))
567 0.666666666667
568
569 @return: photbands, function to construct color
570 @rtype: tuple,callable
571 """
572 system,band = photband.split('.')
573 band = band.strip()
574 photbands = get_color_photband(photband)
575 if len(band.split('-'))==2:
576 func = lambda f0,f1: f0/f1
577 elif band=='M1':
578 func = lambda fv,fb,fy: fv*fy/fb**2
579 elif band=='C1':
580 func = lambda fv,fb,fu: fu*fb/fv**2
581 else:
582 raise ValueError('cannot recognize color {}'.format(photband))
583 return photbands,func
584
585
586 -def eff_wave(photband,model=None,det_type=None):
587 """
588 Return the effective wavelength of a photometric passband.
589
590 The effective wavelength is defined as the average wavelength weighed with
591 the response curve.
592
593 >>> eff_wave('2MASS.J')
594 12393.093155655277
595
596 If you give model fluxes as an extra argument, the wavelengths will take
597 these into account to calculate the `true' effective wavelength (e.g.,
598 Van Der Bliek, 1996), eq 2.
599
600 @param photband: photometric passband
601 @type photband: str ('SYSTEM.FILTER') or array/list of str
602 @param model: model wavelength and fluxes
603 @type model: tuple of 1D arrays (wave,flux)
604 @return: effective wavelength [A]
605 @rtype: float or numpy array
606 """
607
608
609
610 if isinstance(photband,unicode):
611 photband = str(photband)
612 if isinstance(photband,str):
613 single_band = True
614 photband = [photband]
615
616 else:
617 single_band = False
618
619 my_eff_wave = []
620 for iphotband in photband:
621 try:
622 wave,response = get_response(iphotband)
623
624 if det_type is None and len(get_info([iphotband])):
625 det_type = get_info([iphotband])['type'][0]
626 elif det_type is None:
627 det_type = 'CCD'
628 if model is None:
629
630 if det_type=='BOL':
631 this_eff_wave = np.sqrt(np.trapz(response,x=wave)/np.trapz(response/wave**2,x=wave))
632 else:
633 this_eff_wave = np.sqrt(np.trapz(wave*response,x=wave)/np.trapz(response/wave,x=wave))
634 else:
635
636
637 is_response = response>1e-10
638 start_response,end_response = wave[is_response].min(),wave[is_response].max()
639 fluxm = np.sqrt(10**np.interp(np.log10(wave),np.log10(model[0]),np.log10(model[1])))
640
641 if det_type=='CCD':
642 this_eff_wave = np.sqrt(np.trapz(wave*fluxm*response,x=wave) / np.trapz(fluxm*response/wave,x=wave))
643 elif det_type=='BOL':
644 this_eff_wave = np.sqrt(np.trapz(fluxm*response,x=wave) / np.trapz(fluxm*response/wave**2,x=wave))
645
646 except IOError:
647 this_eff_wave = np.nan
648 my_eff_wave.append(this_eff_wave)
649
650 if single_band:
651 my_eff_wave = my_eff_wave[0]
652 else:
653 my_eff_wave = np.array(my_eff_wave,float)
654
655 return my_eff_wave
656
657 @memoized
658 -def get_info(photbands=None):
659 """
660 Return a record array containing all filter information.
661
662 The record arrays contains following columns:
663 - photband
664 - eff_wave
665 - type
666 - vegamag, vegamag_lit
667 - ABmag, ABmag_lit
668 - STmag, STmag_lit
669 - Flam0, Flam0_units, Flam0_lit
670 - Fnu0, Fnu0_units, Fnu0_lit,
671 - source
672
673 @param photbands: list of photbands to get the information from. The input
674 order is equal to the output order. If C{None}, all filters are returned.
675 @type photbands: iterable container (list, tuple, 1Darray)
676 @return: record array containing all information on the requested photbands.
677 @rtype: record array
678 """
679 zp_file = os.path.join(os.path.dirname(os.path.abspath(__file__)),'zeropoints.dat')
680 zp = ascii.read2recarray(zp_file)
681 for iph in custom_filters:
682 if iph=='_prefer_file': continue
683 if 'zp' in custom_filters[iph]:
684 zp = np.hstack([zp,custom_filters[iph]['zp']])
685 zp = zp[np.argsort(zp['photband'])]
686
687
688
689 if photbands is not None:
690 order = np.searchsorted(zp['photband'],photbands)
691 zp = zp[order]
692 keep = (zp['photband']==photbands)
693 zp = zp[keep]
694
695 return zp
696
702 """
703 Update information in zeropoint file, e.g. after calibration.
704
705 Call first L{ivs.sed.model.calibrate} without arguments, and pass the output
706 to this function.
707
708 @param zp: updated contents from C{zeropoints.dat}
709 @type zp: recarray
710 """
711 zp_file = os.path.join(os.path.dirname(os.path.abspath(__file__)),'zeropoints.dat')
712 zp_,comms = ascii.read2recarray(zp_file,return_comments=True)
713 existing = [str(i.strip()) for i in zp_['photband']]
714 resp_files = sorted(glob.glob(os.path.join(os.path.dirname(os.path.abspath(__file__)),'filters/*')))
715 resp_files = [os.path.basename(ff) for ff in resp_files if not os.path.basename(ff) in existing]
716 resp_files.remove('HUMAN.EYE')
717 resp_files.remove('HUMAN.CONES')
718 resp_files.remove('CONES.EYE')
719 if zp is None:
720 zp = zp_
721 logger.info('No new calibrations; previous information on existing response curves is copied')
722 else:
723 logger.info('Received new calibrations contents of zeropoints.dat will be updated')
724
725
726 new_zp = np.zeros(len(resp_files),dtype=zp.dtype)
727 logger.info('Found {} new response curves, adding them with default information'.format(len(resp_files)))
728 for i,respfile in enumerate(resp_files):
729 new_zp[i]['photband'] = respfile
730 new_zp[i]['eff_wave'] = float(eff_wave(respfile))
731 new_zp[i]['type'] = 'CCD'
732 new_zp[i]['vegamag'] = np.nan
733 new_zp[i]['ABmag'] = np.nan
734 new_zp[i]['STmag'] = np.nan
735 new_zp[i]['Flam0_units'] = 'erg/s/cm2/AA'
736 new_zp[i]['Fnu0_units'] = 'erg/s/cm2/AA'
737 new_zp[i]['source'] = 'nan'
738 zp = np.hstack([zp,new_zp])
739 sa = np.argsort(zp['photband'])
740 ascii.write_array(zp[sa],'zeropoints.dat',header=True,auto_width=True,comments=['#'+line for line in comms[:-2]],use_float='%g')
741
742
743
744 if __name__=="__main__":
745 import sys
746 import pylab as pl
747 if not sys.argv[1:]:
748 import doctest
749 doctest.testmod()
750 pl.show()
751
752 else:
753 import itertools
754 responses = list_response()
755 systems = [response.split('.')[0] for response in responses]
756 set_responses = sorted(set([response.split('.')[0] for response in systems]))
757 this_filter = 0
758 for i,resp in enumerate(responses):
759
760 this_system = resp.split('.')[0]
761 nr_filters = systems.count(this_system)
762
763
764
765 p = pl.figure(set_responses.index(this_system),figsize=(10,4.5))
766 if not hasattr(pl.gca(),'color_cycle'):
767 color_cycle = itertools.cycle([pl.cm.spectral(j) for j in np.linspace(0, 1.0, nr_filters)])
768 p = pl.gca().color_cycle = color_cycle
769 color = pl.gca().color_cycle.next()
770 p = pl.title(resp.split('.')[0])
771
772 wave,trans = get_response(resp)
773 p = pl.plot(wave/1e4,trans,label=resp,color=color)
774
775 p = pl.xlabel('Wavelength [micron]')
776 p = pl.ylabel('Transmission')
777
778
779 this_filter+=1
780 if this_filter==nr_filters:
781 this_filter = 0
782 p = pl.legend(prop=dict(size='small'))
783 p = pl.savefig('/home/pieterd/python/ivs/doc/images/ivs_sed_filters_%s'%(this_system));p = pl.close()
784