Shift a spectrum with towards the red or blue side with some radial velocity.

You can give units with the extra keywords vrad_units (units of wavelengths are not important). The shifted wavelengths will be in the same units as the input wave array.

If units are not supplied, the radial velocity is assumed to be in km/s.

If you want to apply a barycentric (or orbital) correction, you'd probabl want to reverse the sign of the radial velocity!

When the keyword flux is set, the spectrum will be interpolated onto the original wavelength grid (so the original wavelength grid will not change). When the keyword flux is not set, the wavelength array will be changed (but the fluxes not, obviously).

When flux is set, fluxes will be returned.

When flux is not set, wavelengths will be returned.

Example usage: shift a spectrum to the blue ('left') with 20 km/s.

>>> wave = np.linspace(3000,8000,1000)
>>> wave_shift1 = doppler_shift(wave,20.)
>>> wave_shift2 = doppler_shift(wave,20000.,vrad_units='m/s')
>>> print(wave_shift1[0],wave_shift1[-1])
(3000.200138457119, 8000.5337025523177)
>>> print(wave_shift2[0],wave_shift2[-1])
(3000.200138457119, 8000.5337025523177)

Parameters:

• wave (ndarray) - wavelength array
• vrad (float (units: km/s) or tuple (float,'units')) - radial velocity (negative shifts towards blue, positive towards red)
• vrad_units (str (interpretable for units.conversions.convert)) - units of radial velocity (default: km/s)

Returns: ndarray

shifted wavelength array/shifted flux

Deterimine vsini of an observed spectrum via the Fourier transform method.

According to Simon-Diaz (2006) and Carroll (1933):

vsini = 0.660 * c/ (lambda * f1)

But more general (see Reiners 2001, Dravins 1990)

vsini = q1 * c/ (lambda*f1)

where f1 is the first minimum of the Fourier transform.

The error is estimated as the Rayleigh limit of the Fourier Transform

Example usage and tests: Generate some data. We need a central wavelength (A), the speed of light in angstrom/s, limb darkening coefficients and test vsinis:

>>> clam = 4480.
>>> c = conversions.convert('m/s','A/s',constants.cc)
>>> epsilons = np.linspace(0.,1.0,10)
>>> vsinis = np.linspace(50,300,10)

We analytically compute the shape of the Fourier transform in the following domain (and need the scipy.special.j1 for this)

>>> x = np.linspace(0,30,1000)[1:]
>>> from scipy.special import j1

Keep track of the calculated and predicted q1 values:

>>> q1s = np.zeros((len(epsilons),len(vsinis)))
>>> q1s_pred = np.zeros((len(epsilons),len(vsinis)))

Start a figure and set the color cycle

>>> p= pl.figure()
>>> p=pl.subplot(131)
>>> color_cycle = [pl.cm.spectral(j) for j in np.linspace(0, 1.0, len(epsilons))]
>>> p = pl.gca().set_color_cycle(color_cycle)
>>> p=pl.subplot(133);p=pl.title('Broadening kernel')
>>> p = pl.gca().set_color_cycle(color_cycle)

Now run over all epsilons and vsinis and determine the q1 constant:

>>> for j,epsilon in enumerate(epsilons):
...    for i,vsini in enumerate(vsinis):
...       vsini = conversions.convert('km/s','A/s',vsini)
...       delta = clam*vsini/c
...       lambdas = np.linspace(-5.,+5.,20000)
...       #-- analytical rotational broadening kernel
...       y = 1-(lambdas/delta)**2
...       G = (2*(1-epsilon)*np.sqrt(y)+pi*epsilon/2.*y)/(pi*delta*(1-epsilon/3))
...       lambdas = lambdas[-np.isnan(G)]
...       G = G[-np.isnan(G)]
...       G /= max(G)
...       #-- analytical FT of rotational broadening kernel
...       g = 2. / (x*(1-epsilon/3.)) * ( (1-epsilon)*j1(x) +  epsilon* (sin(x)/x**2 - cos(x)/x))
...       #-- numerical FT of rotational broadening kernel
...       sigma,g_ = pergrams.deeming(lambdas-lambdas[0],G,fn=2e0,df=1e-3,norm='power')
...       myx = 2*pi*delta*sigma
...       #-- get the minima
...       rise = np.diff(g_[1:])>=0
...       fall = np.diff(g_[:-1])<=0
...       minima = sigma[1:-1][rise & fall]
...       minvals = g_[1:-1][rise & fall]
...       q1s[j,i] =  vsini / (c/clam/minima[0])
...       q1s_pred[j,i] = 0.610 + 0.062*epsilon + 0.027*epsilon**2 + 0.012*epsilon**3 + 0.004*epsilon**4
...    p= pl.subplot(131)
...    p= pl.plot(vsinis,q1s[j],'o',label='$\epsilon$=%.2f'%(epsilon));pl.gca()._get_lines.count -= 1
...    p= pl.plot(vsinis,q1s_pred[j],'-')
...    p= pl.subplot(133)
...    p= pl.plot(lambdas,G,'-')

And plot the results:

>>> p= pl.subplot(131)
>>> p= pl.xlabel('v sini [km/s]');p = pl.ylabel('q1')
>>> p= pl.legend(prop=dict(size='small'))

>>> p= pl.subplot(132);p=pl.title('Fourier transform')
>>> p= pl.plot(x,g**2,'k-',label='Analytical FT')
>>> p= pl.plot(myx,g_/max(g_),'r-',label='Computed FT')
>>> p= pl.plot(minima*2*pi*delta,minvals/max(g_),'bo',label='Minima')
>>> p= pl.legend(prop=dict(size='small'))
>>> p= pl.gca().set_yscale('log')

]include figure]]ivs_spectra_tools_vsini_kernel.png]

Extra keyword arguments are passed to pergrams.deeming

Parameters:

• wave (ndarray) - wavelength array in Angstrom
• flux (ndarray) - normalised flux of the profile

Returns: (array,array),(array,array),array,float

periodogram (freqs in s/km), extrema (weird units), vsini values (km/s), error (km/s)

Apply rotational broadening to a spectrum assuming a linear limb darkening law.

This function is based on the ROTIN program. See Fortran file for explanations of parameters.

Limb darkening law is linear, default value is epsilon=0.6

Possibility to normalize as well by giving continuum in 'cont' parameter.

Warning: method='python' is still experimental, but should work.

Section 1. Parameters for rotational convolution

VROT: v sin i (in km/s):

• if VROT=0 - rotational convolution is a) either not calculated, b) or, if simultaneously FWHM is rather large (vrot/c*lambda < FWHM/20.), vrot is set to FWHM/20*c/lambda;
• if VROT >0 but the previous condition b) applies, the value of VROT is changed as in the previous case
• if VROT<0 - the value of abs(VROT) is used regardless of how small compared to FWHM it is

CHARD: characteristic scale of the variations of unconvolved stellar spectrum (basically, characteristic distance between two neighbouring wavelength points) - in A:

• if =0 - program sets up default (0.01 A)

STEPR: wavelength step for evaluation rotational convolution;

• if =0, the program sets up default (the wavelength interval corresponding to the rotational velocity devided by 3.)
• if <0, convolved spectrum calculated on the original (detailed) SYNSPEC wavelength mesh

Section 2. parameters for instrumental convolution

FWHM: WARNING: this is not the full width at half maximum for Gaussian instrumental profile, but the sigma (FWHM = 2.3548 sigma).

STEPI: wavelength step for evaluating instrumental convolution

• if =0, the program sets up default (FWHM/10.)
• if <0, convolved spectrum calculated with the previous wavelength mesh: either the original (SYNSPEC) one if vrot=0, or the one used in rotational convolution (vrot > 0)

Section 3. wavelength interval and normalization of spectra

ALAM0: initial wavelength ALAM1: final wavelength IREL: for =1 relative spectrum, =0 absolute spectrum

Returns: array, array

wavelength,flux

Combine and weight-average spectra on a common wavelength grid.

list_of_spectra should be a list of lists/arrays. Each element in the main list should be (wavelength,flux,error).

If you have FUSE fits files, use ivs.fits.read_fuse. If you have IUE FITS files, use ivs.fits.read_iue.

After Peter Woitke.

Parameters:

• R (float) - resolution
• lambda0 (tuple (float,str)) - start wavelength, unit
• lambdan (tuple (float,str)) - end wavelength, unit

Returns: array, array, array

binned spectrum (wavelengths,flux, error)

Method to combine a set of spectra while removing cosmic rays by comparing the spectra with each other and removing the outliers.

Algorithm:

For each spectrum a very rough continuum is determined by using a median filter. Then there are multiple passes through the spectra. In one pass, outliers are identified by compairing the flux point with the median of all normalized spectra plus sigma times the standard deviation of all points. The standard deviation is the median of the std for all flux points in 1 spectrum (spec_std) plus the std of all fluxes of a given wavelength point over all spectra.

spec_std = np.median( np.std(fluxes) ) fn > np.median(fn) + sigma * (np.std(fn) + spec_std)

Where fn are the roughly normalized spectra, NOT the original spectra. In the next passes those outliers are not used to calculat the median and std of the spectrum. After the last pass, the flux points that were rejected are replaced by the rough continuum value, and they are summed to get the final flux.

The value for sigma strongly depends on the number of spectra and the quality of the spectra. General, the more spectra, the higher sigma should be. Using a too low value for sigma will throw away to much information. The effect of the window size and the number of runs is limited.

Returns the wavelengths and fluxes of the merged spectra, and if full_output is True, also a list of accepted and rejected points, produced by np.where()

Parameters:

• waves - list of wavelengths
• fluxes - list of fluxes
• vrads - list of radial velocities (optional)
• vrad_units - units of the radial velocities
• sigma - value used for sigma clipping
• window - window size used in median filter
• runs - number of iterations through the spectra
• full_output - True is need to return accepted and rejected

Returns: array, array (, tuple, tuple)

wavelenght and flux of merged spectrum (, accepted and rejected points)

Cross correlate a spectrum with a template, working in velocity space. The velocity range is controlled by using step, nsteps and start_dev as:

velocity = np.arange(start_dev - nsteps * step , start_dev + nsteps * step , step)

If two_step is set to True, then it will run twice, and in the second run focus on the velocity where the correlation is at its maximum.

Returns the velocity and the normalized correlation function

Returns the response curve of the given instrument. Up till now only a HERMES response cruve is available. This response curve is based on 25 spectra of the single sdB star Feige 66, and has a wavelenght range of 3800 to 8000 A.

Parameters:

• instrument (string) - the instrument of which you want the response curve

Divides the provided spectrum by the response curve of the given instrument. Up till now only a HERMES response curve is available.

Parameters:

• wave (numpy array) - wavelenght array
• flux (numpy array) - flux array
• instrument (string) - the instrument of which you want the response curve

Returns: (array, array)

the new spectrum (wave, flux)