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 cc.ivs.sigproc import funclib
214 >>> from cc.ivs.sed import model
215 >>> from cc.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 from astropy.io import fits as pyfits
322 import logging
323 import numpy as np
324
325 from cc.ivs.aux.decorators import memoized
326 from cc.ivs.aux import decorators
327 from cc.ivs.aux import loggers
328 from cc.ivs.io import ascii
329
330 basedir = os.path.dirname(__file__)
331
332 logger = logging.getLogger("SED.FILT")
333 logger.addHandler(loggers.NullHandler())
334
335 custom_filters = {'_prefer_file':True}
340 """
341 Retrieve the response curve of a photometric system 'SYSTEM.FILTER'
342
343 OPEN.BOL represents a bolometric open filter.
344
345 Example usage:
346
347 >>> p = pl.figure()
348 >>> for band in ['J','H','KS']:
349 ... p = pl.plot(*get_response('2MASS.%s'%(band)))
350
351 If you defined a custom filter with the same name as an existing one and
352 you want to use that one in the future, set C{prefer_file=False} in the
353 C{custom_filters} module dictionary.
354
355 @param photband: photometric passband
356 @type photband: str ('SYSTEM.FILTER')
357 @return: (wavelength [A], response)
358 @rtype: (array, array)
359 """
360 photband = photband.upper()
361 prefer_file = custom_filters['_prefer_file']
362 if photband=='OPEN.BOL':
363 return np.array([1,1e10]),np.array([1/(1e10-1),1/(1e10-1)])
364
365 photfile = os.path.join(basedir,'filters',photband)
366 photfile_is_file = os.path.isfile(photfile)
367
368 if photfile_is_file and prefer_file:
369 wave, response = ascii.read2array(photfile).T[:2]
370
371 elif photband in custom_filters:
372 wave, response = custom_filters[photband]['response']
373
374 elif photfile_is_file:
375 wave, response = ascii.read2array(photfile).T[:2]
376 else:
377 raise IOError,('{0} does not exist {1}'.format(photband,custom_filters.keys()))
378 sa = np.argsort(wave)
379 return wave[sa],response[sa]
380
382 """
383 Create a custom filter as a sum of Gaussians.
384
385 @param wave: wavelength to evaluate the profile on
386 @type wave: ndarray
387 @param peaks: heights of the peaks
388 @type peaks: ndarray of length N, with N peaks
389 @param range: wavelength range of the peaks
390 @type range: tuple
391 @param sigma: width of the peaks in units of (range/N)
392 @type sigma: float
393 @return: filter profile
394 @rtype: ndarray
395 """
396 wpnts = np.linspace(range[0],range[1],len(peaks)+2)[1:-1]
397 sigma = (range[1]-range[0])/(sigma*len(peaks))
398 gauss = lambda x,p: p[0] * np.exp( -(x-p[1])**2 / (2.0*p[2]**2))
399 els = [gauss(wave,[pk,mu,sigma]) for pk,mu in zip(peaks,wpnts)]
400 profile = np.array(els).sum(axis=0)
401 return profile
402
405 """
406 Add a custom filter to the set of predefined filters.
407
408 Extra keywords are:
409 'eff_wave', 'type',
410 'vegamag', 'vegamag_lit',
411 'ABmag', 'ABmag_lit',
412 'STmag', 'STmag_lit',
413 'Flam0', 'Flam0_units', 'Flam0_lit',
414 'Fnu0', 'Fnu0_units', 'Fnu0_lit',
415 'source'
416
417 default C{type} is 'CCD'.
418 default C{photband} is 'CUSTOM.FILTER'
419
420 @param wave: wavelength (angstrom)
421 @type wave: ndarray
422 @param response: response
423 @type response: ndarray
424 @param photband: photometric passband
425 @type photband: str ('SYSTEM.FILTER')
426 """
427 kwargs.setdefault('photband','CUSTOM.FILTER')
428 kwargs.setdefault('copy_from','JOHNSON.V')
429 kwargs.setdefault('force',False)
430 photband = kwargs['photband']
431
432 photfile = os.path.join(basedir,'filters',photband)
433 if os.path.isfile(photfile) and not kwargs['force']:
434 raise ValueError,'bandpass {0} already exists'.format(photfile)
435 elif photband in custom_filters:
436 logger.debug('Overwriting previous definition of {0}'.format(photband))
437 custom_filters[photband] = dict(response=(wave,response))
438
439 kwargs.setdefault('type','CCD')
440 kwargs.setdefault('eff_wave',eff_wave(photband,det_type=kwargs['type']))
441
442
443
444 myrow = get_info([kwargs['copy_from']])
445 for name in myrow.dtype.names:
446 if 'lit' in name:
447 myrow[name] = 0
448 myrow[name] = kwargs.pop(name,myrow[name])
449 del decorators.memory[__name__]
450
451 custom_filters[photband]['zp'] = myrow
452 logger.debug('Added photband {0} to the predefined set'.format(photband))
453
455 """
456 Set whether to prefer files or custom filters when both exist.
457
458 @param prefer_file: boolean
459 @type prefer_file: bool
460 """
461 custom_filters['_prefer_file'] = prefer_file
462 logger.info("Prefering {}".format(prefer_file and 'files' or 'custom filters'))
463
466
467 Delta = np.log10(1.+1./R)
468 x = np.arange(np.log10(lambda0),np.log10(lambdan)+Delta,Delta)
469 x = 10**x
470 photbands = []
471 for i in range(len(x)-1):
472 wave = np.linspace(x[i],x[i+1],100)
473 resp = np.ones_like(wave)
474 dw = wave[1]-wave[0]
475 wave = np.hstack([wave[0]-dw,wave,wave[-1]+dw])
476 resp = np.hstack([0,resp,0])
477 photband = 'BOXCAR_R{0:d}.{1:d}'.format(int(R),int(x[i]))
478 try:
479 add_custom_filter(wave,resp,photband=photband)
480 except ValueError:
481 logger.info('{0} already exists, skipping'.format(photband))
482 photbands.append(photband)
483 logger.info('Added spectrophotometric filters')
484 return photbands
485
488 """
489 List available response curves.
490
491 Specify a glob string C{name} and/or a wavelength range to make a selection
492 of all available curves. If nothing is supplied, all curves will be returned.
493
494 @param name: list all curves containing this string
495 @type name: str
496 @param wave_range: list all curves within this wavelength range (A)
497 @type wave_range: (float, float)
498 @return: list of curve files
499 @rtype: list of str
500 """
501
502 if not '*' in name:
503 name_ = '*' + name + '*'
504 else:
505 name_ = name
506 curve_files = sorted(glob.glob(os.path.join(basedir,'filters',name_.upper())))
507 curve_files = sorted(curve_files+[key for key in custom_filters.keys() if ((name in key) and not (key=='_prefer_file'))])
508 curve_files = [cf for cf in curve_files if not ('HUMAN' in cf or 'EYE' in cf) ]
509
510 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])]
511
512 for curve_file in curve_files: logger.info(curve_file)
513 return curve_files
514
518 """
519 Return true if the photometric passband is actually a color.
520
521 @param photband: name of the photometric passband
522 @type photband: string
523 @return: True or False
524 @rtype: bool
525 """
526 if '-' in photband.split('.')[1]:
527 return True
528 elif photband.split('.')[1].upper() in ['M1','C1']:
529 return True
530 else:
531 return False
532
535 """
536 Retrieve the photometric bands from color
537
538 @param photband: name of the photometric passband
539 @type photband: string
540 @return: tuple of strings
541 @rtype: tuple
542 """
543 system,band = photband.split('.')
544 band = band.strip()
545 if '-' in band:
546 bands = tuple(['%s.%s'%(system,iband) for iband in band.split('-')])
547 elif band.upper()=='M1':
548 bands = tuple(['%s.%s'%(system,iband) for iband in ['V','B','Y']])
549 elif band.upper()=='C1':
550 bands = tuple(['%s.%s'%(system,iband) for iband in ['V','B','U']])
551 else:
552 raise ValueError('cannot recognize color {}'.format(photband))
553 return bands
554
556 """
557 Make a color from a color name and fluxes.
558
559 You get two things: a list of photbands that are need to construct the color,
560 and a function which you need to pass fluxes to compute the color.
561
562 >>> bands, func = make_color('JOHNSON.B-V')
563 >>> print(bands)
564 ('JOHNSON.B', 'JOHNSON.V')
565 >>> print(func(2,3.))
566 0.666666666667
567
568 @return: photbands, function to construct color
569 @rtype: tuple,callable
570 """
571 system,band = photband.split('.')
572 band = band.strip()
573 photbands = get_color_photband(photband)
574 if len(band.split('-'))==2:
575 func = lambda f0,f1: f0/f1
576 elif band=='M1':
577 func = lambda fv,fb,fy: fv*fy/fb**2
578 elif band=='C1':
579 func = lambda fv,fb,fu: fu*fb/fv**2
580 else:
581 raise ValueError('cannot recognize color {}'.format(photband))
582 return photbands,func
583
584
585 -def eff_wave(photband,model=None,det_type=None):
586 """
587 Return the effective wavelength of a photometric passband.
588
589 The effective wavelength is defined as the average wavelength weighed with
590 the response curve.
591
592 >>> eff_wave('2MASS.J')
593 12393.093155655277
594
595 If you give model fluxes as an extra argument, the wavelengths will take
596 these into account to calculate the `true' effective wavelength (e.g.,
597 Van Der Bliek, 1996), eq 2.
598
599 @param photband: photometric passband
600 @type photband: str ('SYSTEM.FILTER') or array/list of str
601 @param model: model wavelength and fluxes
602 @type model: tuple of 1D arrays (wave,flux)
603 @return: effective wavelength [A]
604 @rtype: float or numpy array
605 """
606
607
608
609 if isinstance(photband,unicode):
610 photband = str(photband)
611 if isinstance(photband,str):
612 single_band = True
613 photband = [photband]
614
615 else:
616 single_band = False
617
618 my_eff_wave = []
619 for iphotband in photband:
620 try:
621 wave,response = get_response(iphotband)
622
623 if det_type is None and len(get_info([iphotband])):
624 det_type = get_info([iphotband])['type'][0]
625 elif det_type is None:
626 det_type = 'CCD'
627 if model is None:
628
629 if det_type=='BOL':
630 this_eff_wave = np.sqrt(np.trapz(response,x=wave)/np.trapz(response/wave**2,x=wave))
631 else:
632 this_eff_wave = np.sqrt(np.trapz(wave*response,x=wave)/np.trapz(response/wave,x=wave))
633 else:
634
635
636 is_response = response>1e-10
637 start_response,end_response = wave[is_response].min(),wave[is_response].max()
638 fluxm = np.sqrt(10**np.interp(np.log10(wave),np.log10(model[0]),np.log10(model[1])))
639
640 if det_type=='CCD':
641 this_eff_wave = np.sqrt(np.trapz(wave*fluxm*response,x=wave) / np.trapz(fluxm*response/wave,x=wave))
642 elif det_type=='BOL':
643 this_eff_wave = np.sqrt(np.trapz(fluxm*response,x=wave) / np.trapz(fluxm*response/wave**2,x=wave))
644
645 except IOError:
646 this_eff_wave = np.nan
647 my_eff_wave.append(this_eff_wave)
648
649 if single_band:
650 my_eff_wave = my_eff_wave[0]
651 else:
652 my_eff_wave = np.array(my_eff_wave,float)
653
654 return my_eff_wave
655
656 @memoized
657 -def get_info(photbands=None):
658 """
659 Return a record array containing all filter information.
660
661 The record arrays contains following columns:
662 - photband
663 - eff_wave
664 - type
665 - vegamag, vegamag_lit
666 - ABmag, ABmag_lit
667 - STmag, STmag_lit
668 - Flam0, Flam0_units, Flam0_lit
669 - Fnu0, Fnu0_units, Fnu0_lit,
670 - source
671
672 @param photbands: list of photbands to get the information from. The input
673 order is equal to the output order. If C{None}, all filters are returned.
674 @type photbands: iterable container (list, tuple, 1Darray)
675 @return: record array containing all information on the requested photbands.
676 @rtype: record array
677 """
678 zp_file = os.path.join(os.path.dirname(os.path.abspath(__file__)),'zeropoints.dat')
679 zp = ascii.read2recarray(zp_file)
680 for iph in custom_filters:
681 if iph=='_prefer_file': continue
682 if 'zp' in custom_filters[iph]:
683 zp = np.hstack([zp,custom_filters[iph]['zp']])
684 zp = zp[np.argsort(zp['photband'])]
685
686
687
688 if photbands is not None:
689 order = np.searchsorted(zp['photband'],photbands)
690 zp = zp[order]
691 keep = (zp['photband']==photbands)
692 zp = zp[keep]
693
694 return zp
695
701 """
702 Update information in zeropoint file, e.g. after calibration.
703
704 Call first L{cc.ivs.sed.model.calibrate} without arguments, and pass the output
705 to this function.
706
707 @param zp: updated contents from C{zeropoints.dat}
708 @type zp: recarray
709 """
710 zp_file = os.path.join(os.path.dirname(os.path.abspath(__file__)),'zeropoints.dat')
711 zp_,comms = ascii.read2recarray(zp_file,return_comments=True)
712 existing = [str(i.strip()) for i in zp_['photband']]
713 resp_files = sorted(glob.glob(os.path.join(os.path.dirname(os.path.abspath(__file__)),'filters/*')))
714 resp_files = [os.path.basename(ff) for ff in resp_files if not os.path.basename(ff) in existing]
715 resp_files.remove('HUMAN.EYE')
716 resp_files.remove('HUMAN.CONES')
717 resp_files.remove('CONES.EYE')
718 if zp is None:
719 zp = zp_
720 logger.info('No new calibrations; previous information on existing response curves is copied')
721 else:
722 logger.info('Received new calibrations contents of zeropoints.dat will be updated')
723
724
725 new_zp = np.zeros(len(resp_files),dtype=zp.dtype)
726 logger.info('Found {} new response curves, adding them with default information'.format(len(resp_files)))
727 for i,respfile in enumerate(resp_files):
728 new_zp[i]['photband'] = respfile
729 new_zp[i]['eff_wave'] = float(eff_wave(respfile))
730 new_zp[i]['type'] = 'CCD'
731 new_zp[i]['vegamag'] = np.nan
732 new_zp[i]['ABmag'] = np.nan
733 new_zp[i]['STmag'] = np.nan
734 new_zp[i]['Flam0_units'] = 'erg/s/cm2/AA'
735 new_zp[i]['Fnu0_units'] = 'erg/s/cm2/AA'
736 new_zp[i]['source'] = 'nan'
737 zp = np.hstack([zp,new_zp])
738 sa = np.argsort(zp['photband'])
739 ascii.write_array(zp[sa],'zeropoints.dat',header=True,auto_width=True,comments=['#'+line for line in comms[:-2]],use_float='%g')
740
741
742
743 if __name__=="__main__":
744 import sys
745 import pylab as pl
746 if not sys.argv[1:]:
747 import doctest
748 doctest.testmod()
749 pl.show()
750
751 else:
752 import itertools
753 responses = list_response()
754 systems = [response.split('.')[0] for response in responses]
755 set_responses = sorted(set([response.split('.')[0] for response in systems]))
756 this_filter = 0
757 for i,resp in enumerate(responses):
758
759 this_system = resp.split('.')[0]
760 nr_filters = systems.count(this_system)
761
762
763
764 p = pl.figure(set_responses.index(this_system),figsize=(10,4.5))
765 if not hasattr(pl.gca(),'color_cycle'):
766 color_cycle = itertools.cycle([pl.cm.spectral(j) for j in np.linspace(0, 1.0, nr_filters)])
767 p = pl.gca().color_cycle = color_cycle
768 color = pl.gca().color_cycle.next()
769 p = pl.title(resp.split('.')[0])
770
771 wave,trans = get_response(resp)
772 p = pl.plot(wave/1e4,trans,label=resp,color=color)
773
774 p = pl.xlabel('Wavelength [micron]')
775 p = pl.ylabel('Transmission')
776
777
778 this_filter+=1
779 if this_filter==nr_filters:
780 this_filter = 0
781 p = pl.legend(prop=dict(size='small'))
782 p = pl.savefig('/home/pieterd/python/ivs/doc/images/ivs_sed_filters_%s'%(this_system));p = pl.close()
783