1 """
2 Manipulate or extract information from spectra.
3
4 Subsection 1. Example usage
5 ===========================
6
7 In this example, we:
8
9 1. retrieve a synthetic spectrum
10 2. apply a doppler shift of 20 km/s
11 3. rotationally broaden the synthetic spectrum according to a vsini of 66 km/s,
12 a solar-like limb darkening, and an instrumental profile with FWHM=0.25 angstrom
13 4. add noise to the synthetic spectrum
14 5. calculate the vsini with the Fourier Transform Method and check the influence
15 of the limb darkening coefficients
16
17 Retrieve a synthetic spectrum with effective temperature and surface gravity
18 resembling a main sequence star of spectral type A0. For this purpose, we do
19 not need the whole spectrum but just cut out a piece
20
21 >>> from ivs.spectra import model
22 >>> wave = np.linspace(4570,4574,300)
23 >>> wave,flux,cont = model.get_table(teff=10000,logg=4.0,wave=wave)
24 >>> clam = wave[np.argmin(flux)]
25
26 Apply a doppler shift of 20 km/s (the star travels away from us)
27
28 >>> wave_ = doppler_shift(wave,20.)
29 >>> clam_shift = doppler_shift(clam,20.)
30
31 Rotationally broaden the spectrum and convolve with an instrumental profile
32
33 >>> wave_,flux_ = rotational_broadening(wave_,flux,66.,stepr=-1,fwhm=0.25,stepi=-1,epsilon=0.6)
34
35 Add some noise
36
37 >>> fluxn_ = flux_ + np.random.normal(size=len(flux_),scale=0.01)
38
39 Calculate the vsini with the Fourier transform method. To compare the results,
40 first compute the FT of the synthetic broadened spectrum without noise:
41
42 >>> pergram,minima,vsinis,error = vsini(wave_,flux_,epsilon=0.6,clam=clam_shift,window=(4571,4573.5),df=1e-4)
43
44 Make a plot of what we already have:
45
46 >>> p = pl.figure()
47 >>> p = pl.subplot(121)
48 >>> p = pl.plot(wave,flux,'k-',label='Original template')
49 >>> p = pl.plot(wave_,fluxn_,'b-',label='Spectrum + noise')
50 >>> p = pl.plot(wave_,flux_,'r-',lw=2,label='Broadened')
51 >>> p = pl.legend(loc='lower right',prop=dict(size='small'))
52
53 Set the color cycle of the Fourier Transform plot to spectral
54
55 >>> p = pl.subplot(122)
56 >>> color_cycle = [pl.cm.spectral(j) for j in np.linspace(0, 1.0, 10)]
57 >>> p = pl.gca().set_color_cycle(color_cycle)
58 >>> p = pl.plot(pergram[0],pergram[1],'r-',lw=2,label='Correct')
59 >>> p = pl.gca().set_yscale('log')
60
61 Now compute the vsini of the noisy spectrum, assuming different limb darkening
62 parameters
63
64 >>> for epsilon in np.linspace(0.0,1.0,10):
65 ... pergram,minima,vsinis,error = vsini(wave_,fluxn_,epsilon=epsilon,clam=clam_shift,window=(4571,4573.5),df=1e-4)
66 ... p = pl.plot(pergram[0],pergram[1],label='$\epsilon$=%.2f: vsini = %.1f km/s'%(epsilon,vsinis[0]))
67
68 Set the xticks to vsini values for clarity:
69
70 >>> xticks = np.array([25,35,50.,66,80,100,500][::-1])
71 >>> p = pl.xticks(1/xticks,['%.0f'%(i) for i in xticks])
72 >>> p = pl.xlim(0,0.04);p = pl.grid()
73 >>> p = pl.ylim(5e-4,1)
74 >>> p = pl.legend(loc='lower left',prop=dict(size='small'))
75
76 ]include figure]]ivs_spectra_tools_vsini01.png]
77
78 ]include figure]]ivs_spectra_tools_vsini02.png]
79
80 """
81 from ivs.spectra import pyrotin4
82 import numpy as np
83 import logging
84 from numpy import pi,sin,cos,sqrt
85 import scipy.stats
86 from scipy.signal import fftconvolve, medfilt
87 from ivs.timeseries import pergrams
88 from ivs.units import conversions
89 from ivs.units import constants
90
91 from ivs import config
92 from ivs.inout import fits
93 from scipy.interpolate import interp1d
94
95 logger = logging.getLogger("SPEC.TOOLS")
96
98 """
99 Shift a spectrum with towards the red or blue side with some radial velocity.
100
101 You can give units with the extra keywords C{vrad_units} (units of
102 wavelengths are not important). The shifted wavelengths will be in the same
103 units as the input wave array.
104
105 If units are not supplied, the radial velocity is assumed to be in km/s.
106
107 If you want to apply a barycentric (or orbital) correction, you'd probabl
108 want to reverse the sign of the radial velocity!
109
110 When the keyword C{flux} is set, the spectrum will be interpolated onto
111 the original wavelength grid (so the original wavelength grid will not
112 change). When the keyword C{flux} is not set, the wavelength array will be
113 changed (but the fluxes not, obviously).
114
115 When C{flux} is set, fluxes will be returned.
116
117 When C{flux} is not set, wavelengths will be returned.
118
119 Example usage: shift a spectrum to the blue ('left') with 20 km/s.
120
121 >>> wave = np.linspace(3000,8000,1000)
122 >>> wave_shift1 = doppler_shift(wave,20.)
123 >>> wave_shift2 = doppler_shift(wave,20000.,vrad_units='m/s')
124 >>> print(wave_shift1[0],wave_shift1[-1])
125 (3000.200138457119, 8000.5337025523177)
126 >>> print(wave_shift2[0],wave_shift2[-1])
127 (3000.200138457119, 8000.5337025523177)
128
129 @param wave: wavelength array
130 @type wave: ndarray
131 @param vrad: radial velocity (negative shifts towards blue, positive towards red)
132 @type vrad: float (units: km/s) or tuple (float,'units')
133 @param vrad_units: units of radial velocity (default: km/s)
134 @type vrad_units: str (interpretable for C{units.conversions.convert})
135 @return: shifted wavelength array/shifted flux
136 @rtype: ndarray
137 """
138 cc = constants.cc
139 cc = conversions.convert('m/s',vrad_units,cc)
140 wave_out = wave * (1+vrad/cc)
141 if flux is not None:
142 flux = np.interp(wave,wave_out,flux)
143 return flux
144 else:
145 return wave_out
146
147 -def vsini(wave,flux,epsilon=0.6,clam=None,window=None,**kwargs):
148 """
149 Deterimine vsini of an observed spectrum via the Fourier transform method.
150
151 According to Simon-Diaz (2006) and Carroll (1933):
152
153 vsini = 0.660 * c/ (lambda * f1)
154
155 But more general (see Reiners 2001, Dravins 1990)
156
157 vsini = q1 * c/ (lambda*f1)
158
159 where f1 is the first minimum of the Fourier transform.
160
161 The error is estimated as the Rayleigh limit of the Fourier Transform
162
163 Example usage and tests: Generate some data. We need a central wavelength (A),
164 the speed of light in angstrom/s, limb darkening coefficients and test
165 vsinis:
166
167 >>> clam = 4480.
168 >>> c = conversions.convert('m/s','A/s',constants.cc)
169 >>> epsilons = np.linspace(0.,1.0,10)
170 >>> vsinis = np.linspace(50,300,10)
171
172 We analytically compute the shape of the Fourier transform in the following
173 domain (and need the C{scipy.special.j1} for this)
174
175 >>> x = np.linspace(0,30,1000)[1:]
176 >>> from scipy.special import j1
177
178 Keep track of the calculated and predicted q1 values:
179
180 >>> q1s = np.zeros((len(epsilons),len(vsinis)))
181 >>> q1s_pred = np.zeros((len(epsilons),len(vsinis)))
182
183 Start a figure and set the color cycle
184
185 >>> p= pl.figure()
186 >>> p=pl.subplot(131)
187 >>> color_cycle = [pl.cm.spectral(j) for j in np.linspace(0, 1.0, len(epsilons))]
188 >>> p = pl.gca().set_color_cycle(color_cycle)
189 >>> p=pl.subplot(133);p=pl.title('Broadening kernel')
190 >>> p = pl.gca().set_color_cycle(color_cycle)
191
192 Now run over all epsilons and vsinis and determine the q1 constant:
193
194 >>> for j,epsilon in enumerate(epsilons):
195 ... for i,vsini in enumerate(vsinis):
196 ... vsini = conversions.convert('km/s','A/s',vsini)
197 ... delta = clam*vsini/c
198 ... lambdas = np.linspace(-5.,+5.,20000)
199 ... #-- analytical rotational broadening kernel
200 ... y = 1-(lambdas/delta)**2
201 ... G = (2*(1-epsilon)*np.sqrt(y)+pi*epsilon/2.*y)/(pi*delta*(1-epsilon/3))
202 ... lambdas = lambdas[-np.isnan(G)]
203 ... G = G[-np.isnan(G)]
204 ... G /= max(G)
205 ... #-- analytical FT of rotational broadening kernel
206 ... g = 2. / (x*(1-epsilon/3.)) * ( (1-epsilon)*j1(x) + epsilon* (sin(x)/x**2 - cos(x)/x))
207 ... #-- numerical FT of rotational broadening kernel
208 ... sigma,g_ = pergrams.deeming(lambdas-lambdas[0],G,fn=2e0,df=1e-3,norm='power')
209 ... myx = 2*pi*delta*sigma
210 ... #-- get the minima
211 ... rise = np.diff(g_[1:])>=0
212 ... fall = np.diff(g_[:-1])<=0
213 ... minima = sigma[1:-1][rise & fall]
214 ... minvals = g_[1:-1][rise & fall]
215 ... q1s[j,i] = vsini / (c/clam/minima[0])
216 ... q1s_pred[j,i] = 0.610 + 0.062*epsilon + 0.027*epsilon**2 + 0.012*epsilon**3 + 0.004*epsilon**4
217 ... p= pl.subplot(131)
218 ... p= pl.plot(vsinis,q1s[j],'o',label='$\epsilon$=%.2f'%(epsilon));pl.gca()._get_lines.count -= 1
219 ... p= pl.plot(vsinis,q1s_pred[j],'-')
220 ... p= pl.subplot(133)
221 ... p= pl.plot(lambdas,G,'-')
222
223 And plot the results:
224
225 >>> p= pl.subplot(131)
226 >>> p= pl.xlabel('v sini [km/s]');p = pl.ylabel('q1')
227 >>> p= pl.legend(prop=dict(size='small'))
228
229
230 >>> p= pl.subplot(132);p=pl.title('Fourier transform')
231 >>> p= pl.plot(x,g**2,'k-',label='Analytical FT')
232 >>> p= pl.plot(myx,g_/max(g_),'r-',label='Computed FT')
233 >>> p= pl.plot(minima*2*pi*delta,minvals/max(g_),'bo',label='Minima')
234 >>> p= pl.legend(prop=dict(size='small'))
235 >>> p= pl.gca().set_yscale('log')
236
237 ]include figure]]ivs_spectra_tools_vsini_kernel.png]
238
239 Extra keyword arguments are passed to L{pergrams.deeming}
240
241 @param wave: wavelength array in Angstrom
242 @type wave: ndarray
243 @param flux: normalised flux of the profile
244 @type flux: ndarray
245 @rtype: (array,array),(array,array),array,float
246 @return: periodogram (freqs in s/km), extrema (weird units), vsini values (km/s), error (km/s)
247 """
248 cc = conversions.convert('m/s','AA/s',constants.cc)
249
250 if window is not None:
251 keep = (window[0]<=wave) & (wave<=window[1])
252 wave,flux = wave[keep],flux[keep]
253
254
255 if clam is None: clam = ((wave[0]+wave[-1])/2)
256
257 q1 = 0.610 + 0.062*epsilon + 0.027*epsilon**2 + 0.012*epsilon**3 + 0.004*epsilon**4
258
259
260 freqs,ampls = pergrams.deeming(wave,(1-flux),**kwargs)
261 ampls = ampls/max(ampls)
262 error = 1./wave.ptp()
263
264 rise = np.diff(ampls[1:])>=0
265 fall = np.diff(ampls[:-1])<=0
266 minima = freqs[1:-1][rise & fall]
267 minvals = ampls[1:-1][rise & fall]
268
269 freqs = freqs*clam/q1/cc
270 freqs = conversions.convert('s/AA','s/km',freqs,wave=(clam,'AA'))
271 vsini_values = cc/clam*q1/minima
272 vsini_values = conversions.convert('AA/s','km/s',vsini_values)
273
274 error = error*clam/q1/cc
275 error = conversions.convert('s/AA','s/km',error,wave=(clam,'AA'))
276 return (freqs,ampls),(minima,minvals),vsini_values,error
277
278
279
280
281 -def rotational_broadening(wave_spec,flux_spec,vrot,fwhm=0.25,epsilon=0.6,
282 chard=None,stepr=0,stepi=0,alam0=None,alam1=None,
283 irel=0,cont=None,method='fortran'):
284 """
285 Apply rotational broadening to a spectrum assuming a linear limb darkening
286 law.
287
288 This function is based on the ROTIN program. See Fortran file for
289 explanations of parameters.
290
291 Limb darkening law is linear, default value is epsilon=0.6
292
293 Possibility to normalize as well by giving continuum in 'cont' parameter.
294
295 B{Warning}: C{method='python'} is still experimental, but should work.
296
297 Section 1. Parameters for rotational convolution
298 ================================================
299
300 C{VROT}: v sin i (in km/s):
301
302 - if VROT=0 - rotational convolution is
303 a) either not calculated,
304 b) or, if simultaneously FWHM is rather large
305 (vrot/c*lambda < FWHM/20.),
306 vrot is set to FWHM/20*c/lambda;
307 - if VROT >0 but the previous condition b) applies, the
308 value of VROT is changed as in the previous case
309 - if VROT<0 - the value of abs(VROT) is used regardless of
310 how small compared to FWHM it is
311
312 C{CHARD}: characteristic scale of the variations of unconvolved stellar
313 spectrum (basically, characteristic distance between two neighbouring
314 wavelength points) - in A:
315
316 - if =0 - program sets up default (0.01 A)
317
318 C{STEPR}: wavelength step for evaluation rotational convolution;
319
320 - if =0, the program sets up default (the wavelength
321 interval corresponding to the rotational velocity
322 devided by 3.)
323 - if <0, convolved spectrum calculated on the original
324 (detailed) SYNSPEC wavelength mesh
325
326
327 Section 2. parameters for instrumental convolution
328 ==================================================
329
330 C{FWHM}: WARNING: this is not the full width at half maximum for Gaussian
331 instrumental profile, but the sigma (FWHM = 2.3548 sigma).
332
333 C{STEPI}: wavelength step for evaluating instrumental convolution
334 - if =0, the program sets up default (FWHM/10.)
335 - if <0, convolved spectrum calculated with the previous
336 wavelength mesh:
337 either the original (SYNSPEC) one if vrot=0,
338 or the one used in rotational convolution (vrot > 0)
339
340
341 Section 3. wavelength interval and normalization of spectra
342 ===========================================================
343
344 C{ALAM0}: initial wavelength
345 C{ALAM1}: final wavelength
346 C{IREL}: for =1 relative spectrum, =0 absolute spectrum
347
348 @return: wavelength,flux
349 @rtype: array, array
350 """
351 if method=='fortran':
352
353 if alam0 is None: alam0 = wave_spec[0]
354 if alam1 is None: alam1 = wave_spec[-1]
355 if cont is None: cont = (np.ones(1),np.ones(1))
356 contw,contf = cont
357 if chard is None:
358 chard = np.diff(wave_spec).mean()
359
360
361 logger.info('ROTIN rot.broad. with vrot=%.3f (epsilon=%.2f)'%(vrot,epsilon))
362 w3,f3,ind = pyrotin4.pyrotin(wave_spec,flux_spec,contw,contf,
363 vrot,chard,stepr,fwhm,stepi,alam0,alam1,irel,epsilon)
364
365 return w3[:ind],f3[:ind]
366 elif method=='python':
367 logger.info("PYTHON rot.broad with vrot=%.3f (epsilon=%.2f)"%(vrot,epsilon))
368
369 if fwhm>0:
370 fwhm /= 2.3548
371
372 wave_ = np.linspace(wave_spec[0],wave_spec[-1],len(wave_spec))
373 flux_ = np.interp(wave_,wave_spec,flux_spec)
374 dwave = wave_[1]-wave_[0]
375 n = int(2*4*fwhm/dwave)
376 wave_k = np.arange(n)*dwave
377 wave_k-= wave_k[-1]/2.
378 kernel = np.exp(- (wave_k)**2/(2*fwhm**2))
379 kernel /= sum(kernel)
380 flux_conv = fftconvolve(1-flux_,kernel,mode='same')
381 flux_spec = np.interp(wave_spec+dwave/2,wave_,1-flux_conv,left=1,right=1)
382 if vrot>0:
383
384
385 wave_ = np.log(wave_spec)
386 velo_ = np.linspace(wave_[0],wave_[-1],len(wave_))
387 flux_ = np.interp(velo_,wave_,flux_spec)
388 dvelo = velo_[1]-velo_[0]
389 vrot = vrot/(constants.cc*1e-3)
390
391 n = int(2*vrot/dvelo)
392 velo_k = np.arange(n)*dvelo
393 velo_k -= velo_k[-1]/2.
394 y = 1 - (velo_k/vrot)**2
395 G = (2*(1-epsilon)*sqrt(y)+pi*epsilon/2.*y)/(pi*vrot*(1-epsilon/3.0))
396 G /= G.sum()
397
398 flux_conv = fftconvolve(1-flux_,G,mode='same')
399 velo_ = np.arange(len(flux_conv))*dvelo+velo_[0]
400 wave_conv = np.exp(velo_)
401 return wave_conv,1-flux_conv
402 return wave_spec,flux_spec
403 else:
404 raise ValueError("don't understand method {}".format(method))
405
406 -def combine(list_of_spectra,R=200.,lambda0=(950.,'AA'),lambdan=(3350.,'AA')):
407 """
408 Combine and weight-average spectra on a common wavelength grid.
409
410 C{list_of_spectra} should be a list of lists/arrays. Each element in the
411 main list should be (wavelength,flux,error).
412
413 If you have FUSE fits files, use L{ivs.fits.read_fuse}.
414 If you have IUE FITS files, use L{ivs.fits.read_iue}.
415
416 After Peter Woitke.
417
418 @param R: resolution
419 @type R: float
420 @param lambda0: start wavelength, unit
421 @type lambda0: tuple (float,str)
422 @param lambdan: end wavelength, unit
423 @type lambdan: tuple (float,str)
424 @return: binned spectrum (wavelengths,flux, error)
425 @rtype: array, array, array
426 """
427 l0 = conversions.convert(lambda0[1],'AA',lambda0[0])
428 ln = conversions.convert(lambdan[1],'AA',lambdan[0])
429
430 Delta = np.log10(1.+1./R)
431 x = np.arange(np.log10(l0),np.log10(ln)+Delta,Delta)
432 x = 10**x
433 lamc_j = 0.5*(np.roll(x,1)+x)
434
435
436 Ns = len(list_of_spectra)
437 Nw = len(lamc_j)-1
438 binned_fluxes = np.zeros((Ns,Nw))
439 binned_errors = np.inf*np.ones((Ns,Nw))
440
441 for snr,(wave,flux,err) in enumerate(list_of_spectra):
442 wave0 = np.roll(wave,1)
443 wave1 = np.roll(wave,-1)
444 lam_i0_dc = 0.5*(wave0+wave)
445 lam_i1_dc = 0.5*(wave1+wave)
446 dlam_i = lam_i1_dc-lam_i0_dc
447
448 for j in range(Nw):
449 A = np.min(np.vstack([lamc_j[j+1]*np.ones(len(wave)),lam_i1_dc]),axis=0)
450 B = np.max(np.vstack([lamc_j[j]*np.ones(len(wave)),lam_i0_dc]),axis=0)
451 overlaps = scipy.stats.threshold(A-B,threshmin=0)
452 norm = np.sum(overlaps)
453 binned_fluxes[snr,j] = np.sum(flux*overlaps)/norm
454 binned_errors[snr,j] = np.sqrt(np.sum((err*overlaps)**2))/norm
455
456
457
458 binned_fluxes[np.isnan(binned_fluxes)] = 0
459 binned_errors[np.isnan(binned_errors)] = 1e300
460 weights = 1./binned_errors**2
461
462 totalflux = np.sum(weights*binned_fluxes,axis=0)/np.sum(weights,axis=0)
463 totalerr = np.sqrt(np.sum((weights*binned_errors)**2,axis=0))/np.sum(weights,axis=0)
464 totalspec = np.sum(binned_fluxes>0,axis=0)
465
466
467 return x[:-1],totalflux,totalerr,totalspec
468
469 -def merge_cosmic_clipping(waves, fluxes, vrads=None, vrad_units='km/s', sigma=3.0,
470 base='average', offset='std', window=51, runs=2,
471 full_output=False, **kwargs):
472 """
473 Method to combine a set of spectra while removing cosmic rays by comparing the
474 spectra with each other and removing the outliers.
475
476 Algorithm:
477
478 For each spectrum a very rough continuum is determined by using a median filter.
479 Then there are multiple passes through the spectra. In one pass, outliers are
480 identified by compairing the flux point with the median of all normalized spectra
481 plus sigma times the standard deviation of all points. The standard deviation is
482 the median of the std for all flux points in 1 spectrum (spec_std) plus the std
483 of all fluxes of a given wavelength point over all spectra.
484
485 C{spec_std = np.median( np.std(fluxes) )}
486 C{fn > np.median(fn) + sigma * (np.std(fn) + spec_std)}
487
488 Where fn are the roughly normalized spectra, NOT the original spectra.
489 In the next passes those outliers are not used to calculat the median and std
490 of the spectrum. After the last pass, the flux points that were rejected are
491 replaced by the rough continuum value, and they are summed to get the final
492 flux.
493
494 The value for sigma strongly depends on the number of spectra and the quality
495 of the spectra. General, the more spectra, the higher sigma should be. Using
496 a too low value for sigma will throw away to much information. The effect of
497 the window size and the number of runs is limited.
498
499 Returns the wavelengths and fluxes of the merged spectra, and if full_output
500 is True, also a list of accepted and rejected points, produced by np.where()
501
502 @param waves: list of wavelengths
503 @param fluxes: list of fluxes
504 @param vrads: list of radial velocities (optional)
505 @param vrad_units: units of the radial velocities
506 @param sigma: value used for sigma clipping
507 @param window: window size used in median filter
508 @param runs: number of iterations through the spectra
509 @param full_output: True is need to return accepted and rejected
510
511 @return: wavelenght and flux of merged spectrum (, accepted and rejected points)
512 @rtype: array, array (, tuple, tuple)
513 """
514
515
516 base = getattr(np.ma, base)
517
518
519 if vrads != None:
520 for i, (wave, rv) in enumerate(zip(waves, vrads)):
521 waves[i] = doppler_shift(wave, -rv, vrad_units=vrad_units)
522
523
524 wave = waves[0]
525 flux = np.zeros_like(waves[0])
526
527
528 fc = np.array([np.interp(wave, w_, medfilt(f_, window)) for w_, f_ in zip(waves, fluxes)])
529 fc = np.ma.masked_array( fc, mask = fc == 0. )
530
531
532 fo = np.array([np.interp(wave, w_, f_) for w_, f_ in zip(waves, fluxes)])
533 fo = np.ma.masked_array( fo, mask=np.isfinite(fo) == False )
534
535
536 fn = np.array([f_/c_ for f_, c_ in zip(fo, fc)])
537 fn = np.ma.masked_array( fn, mask=np.isfinite(fn) == False )
538
539 for i in range(runs):
540
541 a = base(fn, axis=0)
542 if offset == 'std': a += np.median(np.ma.std(fn, axis=1))
543 s = np.ma.std(fn, axis=0)
544
545
546 fn.mask = np.ma.mask_or( fn.mask, np.ma.make_mask(fn > a+sigma*s) )
547
548
549 flux = np.sum( np.where(fn.mask, fc, fo), axis=0)
550 logger.debug('Merged %i spectra with sigma = %f and base = %s'%(len(waves), sigma, base))
551
552 if full_output:
553 rejected = np.where(fn.mask)
554 accepted = np.where(fn.mask == False)
555 return wave, flux, accepted, rejected
556 else:
557 return wave, flux
558
559 -def cross_correlate(obj_wave, obj_flux, temp_wave, temp_flux, step=0.3, nsteps=500,
560 start_dev=0.0, two_step=False, verbose=False, **kwargs):
561 """
562 Cross correlate a spectrum with a template, working in velocity space. The velocity
563 range is controlled by using step, nsteps and start_dev as:
564
565 velocity = np.arange(start_dev - nsteps * step , start_dev + nsteps * step , step)
566
567 If two_step is set to True, then it will run twice, and in the second run focus on
568 the velocity where the correlation is at its maximum.
569
570 Returns the velocity and the normalized correlation function
571 """
572
573 def correlate(dvel):
574 rebin_flux = interp1d(temp_wave * ( 1 + 1000. * dvel / constants.cc ), temp_flux)(obj_wave)
575 s2 = np.sqrt(np.sum(rebin_flux**2)/len(rebin_flux))
576 return 1. / ( len(obj_flux) * s1 * s2 ) * np.sum( obj_flux * rebin_flux )
577
578
579 s1 = np.sqrt(sum(obj_flux**2)/len(obj_flux))
580 velocity = np.arange(start_dev - nsteps * step , start_dev + nsteps * step , step)
581 correlation = np.array([correlate(dvel) for dvel in velocity])
582
583
584 if two_step:
585 start_dev = velocity[correlation == np.max(correlation)]
586 velocity = np.arange(start_dev - nsteps * step , start_dev + nsteps * step , step)
587 correlation = np.array([correlate(dvel) for dvel in velocity])
588
589
590 correlation = correlation / correlation[0]
591
592 return velocity, correlation
593
595 """
596 Returns the response curve of the given instrument. Up till now only a HERMES
597 response cruve is available. This response curve is based on 25 spectra of the
598 single sdB star Feige 66, and has a wavelenght range of 3800 to 8000 A.
599
600 @param instrument: the instrument of which you want the response curve
601 @type instrument: string
602 """
603
604 basedir = 'spectables/responses/'
605
606 if instrument == 'hermes':
607 basename = 'response_Feige66.fits'
608
609 response = config.get_datafile(basedir,basename)
610
611 wave, flux = fits.read_spectrum(response)
612
613 return wave, flux
614
616 """
617 Divides the provided spectrum by the response curve of the given instrument. Up till now only
618 a HERMES response curve is available.
619
620 @param wave: wavelenght array
621 @type wave: numpy array
622 @param flux: flux array
623 @type flux: numpy array
624 @param instrument: the instrument of which you want the response curve
625 @type instrument: string
626
627 @return: the new spectrum (wave, flux)
628 @rtype: (array, array)
629 """
630
631
632 wr, fr = get_response(instrument=instrument)
633 fr_ = interp1d(wr, fr, kind='linear')
634
635
636 flux = flux[(wave >= wr[0]) & (wave <= wr[-1])]
637 wave = wave[(wave >= wr[0]) & (wave <= wr[-1])]
638
639 flux = flux / fr_(wave)
640
641 return wave, flux
642
643
644 if __name__=="__main__":
645
646 import pylab as pl
647 from ivs.inout import fits
648 import time
649
650 temp = '/STER/mercator/hermes/%s/reduced/%s_HRF_OBJ_ext_CosmicsRemoved_log_merged_c.fits'
651 objlist = [('20090619','237033'), ('20090701','240226'),
652 ('20090712','241334'), ('20090712','241335'),
653 ('20100107','00268012'), ('20100120','00272619'),
654 ('20100120','00272620'), ('20100203','00273577'),
655 ('20100203','00273578'), ('20100303','00275671'),
656 ('20100303','00275672'), ('20100410','00281505'),
657 ('20100519','00284636'), ('20110222','00334558'),
658 ('20110319','00336547'), ('20110324','00339848'),
659 ('20110401','00342273'), ('20110402','00342363'),
660 ('20110406','00342699'), ('20110408','00342896'),
661 ('20120107','00391289'), ('20120110','00391633'),
662 ('20120116','00392217'), ('20120127','00393151'),
663 ('20120209','00394175'), ('20120330','00399697'),
664 ('20120420','00404769'), ('20120506','00406531'),
665 ('20130106','00445346'), ('20130215','00452556'),
666 ('20130406','00457718'), ('20130530','00474128')]
667 mergeList = [temp%o for o in objlist]
668
669 t1 = time.time()
670 waves, fluxes = [], []
671 wtotal, ftotal = fits.read_spectrum(mergeList[0])
672 wtotal = wtotal[(wtotal>5000) & (wtotal<7000)]
673 ftotal = np.zeros_like(wtotal)
674 for ifile in mergeList:
675 w, f = fits.read_spectrum(ifile)
676 f = f[(w>5000) & (w<7000)]
677 w = w[(w>5000) & (w<7000)]
678 waves.append(w)
679 fluxes.append(f)
680 ftotal += np.interp(wtotal,w,f)
681 t1 = time.time() - t1
682
683 t2 = time.time()
684 wave, flux, accepted, rejected = merge_cosmic_clipping(waves, fluxes, sigma=5.0,
685 window=21, offset='std', base='average', full_output=True)
686 t2 = time.time() - t2
687
688 print 'reading ', t1
689 print 'processing ', t2
690
691 wrej = wtotal[rejected[1]]
692 frej = ftotal[rejected[1]]
693 print len(wrej)
694
695
696 ftotal = ftotal[(wtotal>5000) & (wtotal<7000)]
697 wtotal = wtotal[(wtotal>5000) & (wtotal<7000)]
698 flux = flux[(wave>5000) & (wave<7000)]
699 wave = wave[(wave>5000) & (wave<7000)]
700
701 pl.plot(wtotal, ftotal)
702 pl.plot(wave, flux)
703 pl.plot(wrej, frej, '.r')
704 pl.show()
705
706
707
708
709
710
711
712