Background information spectral translations
===============================================
.. highlight:: none
:linenothreshold: 1000
Introduction
++++++++++++
This background information has been written for two reasons. First we wanted to get
some understanding of the conversions between spectral quantities and second,
we wanted to have some knowledge about legacy FITS headers (of which there must be
a lot) where applying the conversions of WCSLIB in the context of module :mod:`wcs`
without modifications will give wrong results.
.. warning::
One needs to be aware of the fact that WCSLIB converts between frequencies
and velocities in the same reference system while in legacy FITS headers it is
common to give a topocentric reference frequency and a reference velocity in a
different reference system.
.. _spectralbackground-alt:
Alternate headers for a spectral line example
+++++++++++++++++++++++++++++++++++++++++++++
In "Representations of spectral coordinates in FITS" ([Ref3]_ ), section 10.1
deals with an example of a VLA spectral line cube which is regularly sampled
in frequency (CTYPE3='FREQ'). The section describes how one can define
alternative FITS headers to deal with different velocity definitions.
We want to examine this exercise in more detail than provided in the
article to illustrate how a FITS header can be modified and serve as an
alternate header.
The topocentric spectral properties in the FITS header from the paper are::
CTYPE3= 'FREQ'
CRVAL3= 1.37835117405e9
CDELT3= 9.765625e4
CRPIX3= 32
CUNIT3= 'Hz'
RESTFRQ= 1.420405752e+9
SPECSYS='TOPOCENT'
.. note::
For a pixel coordinate :math:`N`, reference pixel :math:`N_{ref}` with reference
world coordinate :math:`W_{ref}` and a step size in
world coordinates :math:`\Delta W`, the world coordinate :math:`W` is calculated with:
.. math::
:label: eq2
W(N) = W_{ref} + (N - N_{ref}) \times \Delta W
If *CTYPE* contains a code for a non linear conversion algorithm
(as in CTYPE='VOPT-F2W') then this relation cannot be applied.
As stated in the note above, code for a conversion algorithm is important.
The statements can be verified with the following script:
.. code-block:: python
#!/usr/bin/env python
from kapteyn import wcs
Z0 = 9120000 # Barycentric optical reference velocity
dZ0 = -2.1882651e+4 # Increment in barycentric optical velocity
N = 32 # Pixel coordinate of reference pixel
header = { 'NAXIS' : 1,
'RESTWAV' : 0.211061140507, # [m]
'CTYPE1' : 'VOPT',
'CRVAL1' : Z0, # [m/s]
'CDELT1' : dZ0, # [m/s]
'CRPIX1' : N,
'CUNIT1' : 'm/s'
}
spec = wcs.Projection(header)
print "From VOPT: Pixel, velocity wcs, velocity linear (%s)" % spec.units
pixels = range(30,35)
Vwcs = spec.toworld1d(pixels)
for p,v in zip(pixels, Vwcs):
print p, v/1000.0, (Z0 + (p-N)*dZ0)/1000.0
header = { 'NAXIS' : 1,
'CNAME1' : 'Barycentric optical velocity',
'RESTWAV' : 0.211061140507, # [m]
'CTYPE1' : 'VOPT-F2W',
'CRVAL1' : Z0, # [m/s]
'CDELT1' : dZ0, # [m/s]
'CRPIX1' : N,
'CUNIT1' : 'm/s'
}
spec = wcs.Projection(header)
print "From VOPT-F2W: Pixel, velocity wcs, velocity linear (%s)" % spec.units
pixels = range(30,35)
Vwcs = spec.toworld1d(pixels)
for p,v in zip(pixels, Vwcs):
print p, v/1000.0, (Z0 + (p-N)*dZ0)/1000.0
# Output:
#
# From VOPT: Pixel, velocity wcs, velocity linear (m/s)
# Conversion is linear; no differences
# 30 9163.765302 9163.765302
# 31 9141.882651 9141.882651
# 32 9120.0 9120.0
# 33 9098.117349 9098.117349
# 34 9076.234698 9076.234698
# From VOPT-F2W: Pixel, velocity wcs, velocity linear (m/s)
# Conversion is not linear
# 30 9163.77150335 9163.765302
# 31 9141.88420123 9141.882651
# 32 9120.0 9120.0
# 33 9098.11889901 9098.117349
# 34 9076.24089759 9076.234698
Relation optical velocity and barycentric/lsrk reference frequency
--------------------------------------------------------------------
Let's start to find the alternate header information for the header from article in
[Ref3]_ .
The extra information about the velocity there is that we have an optical barycentric
velocity of 9120 km/s (as required by an observer) stored as an alternate FITS keyword CRVAL3Z.::
CTYPE3Z= 'VOPT-F2W'
CRVAL3Z= 9.120e+6 / [m/s]
The relation between frequency and optical velocity requires a rest frequency (RESTFRQ=).
The relation is:
.. math::
:label: eq5
Z = c\ \bigl(\frac{\lambda - \lambda_0}{\lambda_0}\bigr) = c\ \bigl(\frac{\nu_0 - \nu}{\nu}\bigr)
We adopted variable Z for velocities following the optical definition.
The header tells us that equal steps in pixel coordinates are equal steps in frequency
and the formula above shows that these steps in terms of optical velocity
depends on the frequency in a non-linear way. Therefore we set the conversion algorithm
to **F2W** which indicates that there is a non linear conversion from frequency to wavelength
(optical velocities are associated with wavelength, see [Ref3]_ .). Note that we can use wildcards
for the non linear conversion algorithm, so *CTYPE3Z='VOPT-???'* is also allowed in
our programs.
We can rewrite equation 1 into:
.. math::
:label: eq10
\nu = \frac{\nu_0}{(1+Z/c)}
If we enter the numbers we get a **barycentric** HI reference frequency:
.. math::
:label: eq20
\nu_b = \frac{1.420405752\times 10^9}{(1+9120000/299792458.0)} = 1378471216.43\ Hz
and we have part of a new alternate header::
CTYPE3F= 'FREQ'
CRVAL3F= 1.37847121643e+9 / [Hz]
So given an optical velocity in a reference system (in our case the barycentric system),
we can calculate which barycentric frequency we can use as a reference frequency.
For a conversion between a barycentric frequency and a barycentric velocity we
also need to know what the barycentric frequency increment is.
Barycentric/lsrk frequency increments
--------------------------------------
.. image:: topocentriccorrection.png
:width: 400
:align: center
*fig.1 Overview of velocities and frequencies of barycenter (B) and Earth (E) w.r.t. source.
The arrows represent velocities. The object and the Earth are moving. The longest arrow represents the
(relativistic) addition of two velocities*
Let's use index *b* for variables bound to the barycentric system and *e*
for the topocentric system.
This frequency, :math:`\nu_b` =1.37847121643 GHz is greater than the reference frequency
:math:`\nu_e` at the observatory (FITS keyword `CRVAL3=` 1.37835117405 GHz).
**The difference between frequencies in the topocentric and barycentric system
is caused by the difference between the velocities of reference frames B and E at
the time of observation.**
This velocity is a *true* velocity. It is called the *topocentric correction*.
Let's try to find an expression for this topocentric correction in terms of frequencies.
The relation between a true velocity and
a shift in frequency is given by the formula
.. math::
:label: eq30
\nu = \nu_0\sqrt{\frac{1-v/c}{1+v/c}} = \nu_0\sqrt{\frac{c-v}{c+v}} =
\nu_0 \frac{c-v}{\sqrt{c^2-v^2}}
If we want to express the apparent radial velocity in terms of frequencies, then this can be written as:
.. math::
:label: eq40
v = c\ \frac{\nu_0^2-\nu^2}{\nu_0^2+\nu^2}
For the apparent radial velocities :math:`v_b` and :math:`v_e` we have:
.. math::
:label: eq50
v_b = c\ \frac{\nu_0^2-\nu_b^2}{\nu_0^2+\nu_b^2}=299792458.0 \ \frac{1420405752.0^2-1378471216.43^2}{1420405752.0^2+1378471216.43^2} = 8981342.29811\ m/s
and:
.. math::
:label: eq60
v_e = c\ \frac{\nu_0^2-\nu_e^2}{\nu_0^2+\nu_e^2}=299792458.0 \ \frac{1420405752.0^2-1378351174.05^2}{1420405752.0^2+1378351174.05^2} = 9007426.97201\ m/s
The relativistic addition of velocities in fig. 1. requires:
.. math::
:label: eq70
v_e = \frac{v_b + v _{t}}{1 + \frac{v_b v_{t}}{c^2}}
which gives the topocentric correction as:
.. math::
:label: eq80
v_t = \frac{v_e - v_{b}}{1 - \frac{v_b v_{e}}{c^2}}
With the numbers inserted we find:
.. math::
:label: eq90
v_t = \frac{9007426.97201 - 8981342.29811}{1 - \frac{8981342.29811\times 9007426.97201}{299792458.0^2}} = 26108.1743997\ m/s
If the FITS header has keywords with the position of the source, the time of observation and
the location of the observatory then one can calculate the topocentric correction by hand.
This information was needed at the observatory to set a frequency for a given barycentric
velocity. However many FITS files do not have enough information to calculate the topocentric correction.
Also it is not needed if one knows the shifted frequencies :math:`\nu_e` and :math:`\nu_b` , then
we can calculate the topocentric velocity without calculating the apparent radial velocities.
This can be shown if we insert the expressions for velocities :math:`v_e` and :math:`v_b`
in the expression for :math:`v_t` . Then after some rearranging one finds:
.. math::
:label: eq100
v_t = c\ \frac{\nu_b^2-\nu_e^2}{\nu_b^2+\nu_e^2}
and with the numbers:
.. math::
:label: eq110
v_t = 299792458.0\ \frac{1378471216.43^2-1378351174.05^2}{1378471216.43^2+1378351174.05^2} = 26108.1743998\ m/s
which is consistent with :eq:`eq90`.
::
VELOSYSZ=26108 / [m/s]
With a given topocentric correction and the reference frequency in the barycenter
we can reconstruct the reference frequency at the
observatory with :eq:`eq100` written as:
.. math::
:label: eq120
\nu_e =\nu_b\sqrt{\frac{c-v_t}{c+v_t}}
.. note::
1) It is important to realize that the reference frequency at E is smaller
than the reference frequency at B because w.r.t. the source E moves faster than B.
So if there is a change in the velocity of the source, the frequencies in B and E will change,
but
the topocentric correction keeps the same value and therefore the relation between
the frequencies :math:`\nu_e` and :math:`\nu_b` remains the same (eq. :eq:`eq120`).
.. note::
2) If we forget about the source and we have an *event on E* with a certain frequency
then an *observer* in barycenter *B* will observe a *lower* frequency.
This is because on the line that connects the source and B, the observatory at E moves away
from B which decreases the remote frequency.
So if we change a frequency on E by tuning the receiver at the observatory at frequency
:math:`\nu_e + \Delta \nu_e` ,
then the observer at B would observe a smaller frequency
:math:`\nu_b + \Delta \nu_b` .
The amount of the decrease is related to the topocentric correction as follows:
.. math::
:label: eq130
\nu_b+\Delta \nu_b = (\nu_e+\Delta \nu_e) \sqrt{\frac{c-v_t}{c+v_t}}
and therefore we can write for the frequency bandwidth in B:
.. math::
:label: eq140
\Delta \nu_b =\Delta \nu_e\sqrt{\frac{c-v_t}{c+v_t}}
At first it seems that this contradicts eq. :eq:`eq120`
(where the indices seem to be swapped), but this is not true
because we changed the frame of the observer from Earth to the barycenter.
The event was in E and it is observed in B.
.. math::
:label: eq150
\Delta \nu_b = 97656.25\ \frac{\sqrt{299792458.0-26108.1743998}}{\sqrt{299792458.0+26108.1743998}} = 97647.745732\ Hz
The increment in frequency therefore becomes 97.64775 kHz::
CDELT3F= 9.764775e+4 / [Hz]
So if we change CRVAL1 and CDELT1 in our demonstration script to the barycentric values,
we get the barycentric optical convention velocities for the pixels. As a check we listed
the script and the value for pixel 32 which is exactly 9120 (km/s):
.. code-block:: python
#!/usr/bin/env python
from kapteyn import wcs
header = { 'NAXIS' : 1,
'CTYPE1' : 'FREQ',
'CRVAL1' : 1378471216.4292786,
'CRPIX1' : 32,
'CUNIT1' : 'Hz',
'CDELT1' : 97647.745732,
'RESTFRQ': 1.420405752e+9
}
spec = wcs.Projection(header).spectra('VOPT-F2W')
pixels = range(30,35)
Vwcs = spec.toworld1d(pixels)
print "Pixel, velocity (%s)" % spec.units
for p,v in zip(pixels, Vwcs):
print p, v/1000.0
print "Pixel at velocity 9120 km/s: ", spec.topixel1d(9120000)
# Output
# Pixel, velocity (m/s)
# 30 9163.77150423
# 31 9141.88420167
# 32 9120.0
# 33 9098.11889856
# 34 9076.2408967
# Pixel at velocity 9120 km/s: 32.0
Note: A closure test is added with method `topixel1d()`
.. note::
In the previous two sections we started with a topocentric frequency and
a topocentric frequency increment and derived values for a barycentric frequency
and a barycentric frequency increment. These values can be used
to set an alternate header (barycentric frequency system 'F') for which we
can convert between frequency and optical velocity.
For GIPSY legacy headers these steps are used to convert between
topocentric frequencies and velocities in another reference system,
See :ref:`spectral_gipsy`
Increment in barycentric/lsrk optical velocity
-----------------------------------------------
The optical velocity was given by:
.. math::
:label: eq160
Z = c\ \bigl(\frac{\nu_0 - \nu}{\nu}\bigr) = c\ \bigl(\frac{\nu_0}{\nu} - 1\bigr)
Its derivative is:
.. math::
:label: eq170
\frac{dZ}{d\nu} = \frac{-c \nu_0}{\nu^2}
But for :math:`\nu` we have the expression:
.. math::
:label: eq180
\nu = \frac{\nu_0}{(1+\frac{Z}{c})}
so we end up with:
.. math::
:label: eq190
dZ = \frac{-c}{\nu_0}\ {\bigl(1+\frac{Z}{c}\bigr)}^2\ d\nu
With :math:`d\nu = \Delta \nu_b` and the given barycentric velocity
:math:`Z_b` = 9120000 m/s,
this gives an increment in optical velocity of:
.. math::
:label: eq200
dZ_b = \frac{-299792458.0}{1420405752.0}\ {\bigl(1+\frac{9120000.0}{299792458.0}\bigr)}^2\ 97647.745732 =-21882.651\ m/s
With these values we explained some other alternate header
keywords in the basic spectral-line example::
CDELT3Z= -2.1882651e+4 / [m/s]
SPECSYSZ= 'BARYCENT' / Velocities w.r.t. barycenter
SSYSOBSZ= 'TOPOCENT' / Observation was made from the 'TOPOCENT' frame
Barycentric/lsrk radio velocity
--------------------------------
For radio velocities one needs to apply the definition:
.. math::
:label: eq210
V_{radio} = V = c\ \bigl(\frac{\nu_0 - \nu}{\nu_0}\bigr)
and for the shifted frequency we derive from this equation:
.. math::
:label: eq220
\nu = \nu_0\ \bigl(1-\frac{V}{c}\bigr)
and the spectral translation code becomes:
`proj.spectra('VRAD')`
In the next code example we demonstrate for a barycentric radio velocity
*V* = 8850.750904 km/s how to calculate the barycentric velocities at arbitrary pixels.
This velocity is derived from the optical example in a way that shifted frequency and
topocentric correction are the same. One can use the formula
.. math::
:label: eq225
\frac{V_b}{Z_b} = \frac{\nu_b}{\nu_0}
to find the value of :math:`V_b = 1.37847121643*9120/1.420405752 = 8850.750904` km/s (with the frequencies in GHz and the velocity in km/s).
In a next section we will derive this value in another way; see :eq:`eq230` and :eq:`eq240`
.. code-block:: python
#!/usr/bin/env python
from kapteyn import wcs
import numpy as n
c = 299792458.0 # Speed of light (m/s)
f = 1.37835117405e9 # Topocentric reference frequency (Hz)
df = 9.765625e4 # Topocentric frequency increment (Hz)
f0 = 1.420405752e+9 # Rest frequency (Hz)
V = 8850750.904 # Barycentric radio velocity (m/s)
fb = f0*(1-V/c)
print "Barycentric freq.: ", fb
v = c * ((fb*fb-f*f)/(fb*fb+f*f))
print "VELOSYSR= Topocentric correction:", v, "m/s"
dfb = df*(c-v)/n.sqrt(c*c-v*v)
print "CDELT3F= Delta in frequency in the barycentric frame eq.4): ", dfb
header = { 'NAXIS' : 1,
'CTYPE1' : 'FREQ',
'CRVAL1' : fb,
'CRPIX1' : 32,
'CUNIT1' : 'Hz',
'CDELT1' : dfb,
'RESTFRQ': 1.420405752e+9
}
line = wcs.Projection(header).spectra('VRAD')
pixels = range(30,35)
Vwcs = line.toworld1d(pixels)
for p,v in zip(pixels, Vwcs):
print p, v/1000
# Output:
# Barycentric freq.: 1378471216.43
# VELOSYSR= Topocentric correction: 26108.1745986 m/s
# CDELT3F= Delta in frequency in the barycentric frame eq.4): 97647.745732
#
# Output Radio velocities (km/s)
# 30 8891.97019316
# 31 8871.36054858
# 32 8850.750904
# 33 8830.14125942
# 34 8809.53161484
Frequency to Radio velocity
-----------------------------
From the definition of radio velocity:
.. math::
:label: eq230
V = c\ \bigl(\frac{\nu_0 - \nu}{\nu_0}\bigr)
we can find a radio velocity that corresponds to the value of the optical
velocity. This (barycentric) optical velocity (9120 Km/s) caused a shift of the rest frequency.
The new frequency became :math:`\nu_b = 1.37847122 \times 10^9 Hz`.
If we insert this number in the equation above we find:
.. math::
:label: eq240
V_b = c\ \bigl(\frac{1420405752.0 - 1378471216.43}{1420405752.0}\bigr) = 8850750.90419\ m/s
The formula for a direct conversion from optical to radio velocity can be derived by
inserting the formula for the frequency shift corresponding to optical velocity, into
the expression for the radio velocity:
.. math::
:label: eq250
V = c\ \bigl(1 - \frac{1}{1+\frac{Z}{c}}\bigr)
With eq. :eq:`eq230` it is easy to find the increment of the velocity if the increment in frequency
at the reference frequency is given:
.. math::
:label: eq260
dV = \frac{-c}{\nu_0}\ d\nu
Note that this increment in frequency is the **increment in the barycentric system**!
Inserting the numbers with :math:`d\nu = \Delta \nu_b` we find:
.. math::
:label: eq270
dV_b = \frac{-299792458.0}{1420405752.0}\times 97647.7457312 = -20609.644582\ m/s
This gives us another two values for the alternate header keywords::
CTYPE3R= 'VRAD'
CRVAL3R= 8.85075090419e+6 / [m/s]
CDELT3R= -2.0609645e+4 / [m/s]
Note that *CTYPE3R= 'VRAD'* indicates that the conversion between frequency and radio velocity
is linear.
The next script shows how we can use these new header values to get a list of
radio velocities as function of pixel. We commented out the rest frequency. Its value
is not necessary because we can rewrite the formulas for the velocity in terms of
:math:`\nu/\nu_0` and :math:`\Delta \nu/\nu_0`
.. code-block:: python
#!/usr/bin/env python
from kapteyn import wcs
header = { 'NAXIS' : 1,
'CTYPE1' : 'VRAD',
'CRVAL1' : 8850750.904193053,
'CRPIX1' : 32,
'CUNIT1' : 'm/s',
'CDELT1' : -20609.644582145629,
# 'RESTFRQ': 1.420405752e+9
}
line = wcs.Projection(header)
pixels = range(30,35)
Vwcs = line.toworld1d(pixels)
for p,v in zip(pixels, Vwcs):
print p, v/1000
#
# Output barycentric radio velocity in km/s:
# 30 8891.97019336
# 31 8871.36054878
# 32 8850.75090419
# 33 8830.14125961
# 34 8809.53161503
Alternatively use the spectral translation method `spectra()`
with the values of the barycentric frequency and frequency increment
as follows to get (exactly) the same output:
.. code-block:: python
#!/usr/bin/env python
from kapteyn import wcs
header = { 'NAXIS' : 1,
'CTYPE1' : 'FREQ',
'CRVAL1' : 1378471216.4292786,
'CRPIX1' : 32,
'CUNIT1' : 'Hz',
'CDELT1' : 97647.745732,
'RESTFRQ': 1.420405752e+9
}
line = wcs.Projection(header).spectra('VRAD')
pixels = range(30,35)
Vwcs = line.toworld1d(pixels)
for p,v in zip(pixels, Vwcs):
print p, v/1000
#
# Output barycentric radio velocity in km/s:
# 30 8891.97019336
# 31 8871.36054878
# 32 8850.75090419
# 33 8830.14125961
# 34 8809.53161503
Frequency to Apparent radial velocity
-------------------------------------
As written before, the relation between a true velocity and a shifted frequency is:
.. math::
:label: eq300
v = c\ \frac{\nu_0^2-\nu^2}{\nu_0^2+\nu^2}
Observed from the barycenter the source has an apparent radial velocity:
.. math::
:label: eq310
v_b = 299792458.0 \ \frac{1420405752.0^2-1378471216.42927^2}{1420405752.0^2+1378471216.42927^2} = 8981342.29811\ m/s
::
CTYPE3V= 'VELO-F2V'
CRVAL3V= 8.98134229811e+6 / [m/s]
Note that *CTYPE3V= 'VELO-F2V'* indicates that we derived these velocities from a system in which
the frequency is linear with the pixel value.
For the increment of the apparent radial velocity we need to find the derivative of eq. :eq:`eq40`
.. math::
:label: eq320
\frac{dv}{d\nu} = c(\nu_0^2-\nu^2)\frac{d}{d\nu}{(\nu_0^2+\nu^2)}^{-1} + c{(\nu_0^2+\nu^2)}^{-1}\frac{d}{d\nu}(\nu_0^2-\nu^2)
This works out as:
.. math::
:label: eq330
dv = \frac{-4 c \nu \nu_0^2}{{(\nu_0^2+\nu^2)}^2}\ d\nu
and with the appropriate numbers inserted for :math:`d\nu = \Delta \nu_b`
and :math:`\nu = \nu_b`:
.. math::
:label: eq340
dv_b = \frac{-4 \times 299792458.0\times 1378471216.4292786\times 1420405752.0^2}{{(1420405752.0^2+1378471216.4292786^2)}^2}\ 97647.745732 = -21217.55136
which reveals the value of another keyword from the header in the article's example::
CDELT3V= -2.1217551e+4 / [m/s]
Sometimes you might encounter an alternative formula that doesn't list the frequency.
It uses eq. :eq:`eq30` to express the frequency in terms of the apparent radial velocity and the rest
frequency.
.. math::
:label: eq350
\nu = \nu_0\sqrt{\frac{1-v/c}{1+v/c}}
If you insert this into:
.. math::
:label: eq360
dv = \frac{-4 c \nu \nu_0^2}{{(\nu_0^2+\nu^2)}^2}\ d\nu
then after some rearrangements you end up with the expression:
.. math::
:label: eq370
dv = \frac{-c}{\nu_0}\ \sqrt{(1-\frac{v}{c})}\ {(1+\frac{v}{c})}^{\frac{3}{2}}\ d\nu
If you insert v = 8981342.29811 (m/s) in this expression you will get exactly the same
apparent radial velocity increment (-2.1217551e+4 m/s).
We found an apparent radial velocity and
calculated the increment for this radial velocity. With a short script and
a minimal header we demonstrate how to use WCSLIB to get
an apparent radial velocity for an arbitrary pixel:
.. code-block:: python
#!/usr/bin/env python
from kapteyn import wcs
header = { 'NAXIS' : 1,
'CTYPE1' : 'VELO-F2V',
'CRVAL1' : 8981342.2981121931,
'CRPIX1' : 32,
'CUNIT1' : 'm/s',
'CDELT1' : -21217.5513673598,
'RESTFRQ': 1.420405752e+9
}
line = wcs.Projection(header)
pixels = range(30,35)
Vwcs = line.toworld1d(pixels)
for p,v in zip(pixels, Vwcs):
print p, v/1000
# Output:
# 30 9023.78022672
# 31 9002.56055595
# 32 8981.34229811
# 33 8960.12545322
# 34 8938.9100213
How can this work?
From eq. :eq:`eq350` and eq. :eq:`eq360` it is obvious that WCSLIB can calculate the reference frequency from
the reference apparent radial velocity. For this reference frequency and the increment
in apparent radial velocity it can calculate the increment in frequency at this reference frequency.
Then we have all the information to use eq. :eq:`eq350` to calculate radial
velocities for different frequencies (i.e. different pixels). Note that the step in
frequency is linear and the step in radial velocity is **not** (which explains the
extension 'F2V' in the CTYPE keyword).
Next script and header is an alternative to get exactly the same results. The header lists the barycentric
frequency and frequency increment. We need a spectral translation with
method `spectra()` to tell WCSLIB to calculate apparent radial velocities:
.. code-block:: python
#!/usr/bin/env python
from kapteyn import wcs
header = { 'NAXIS' : 1,
'CTYPE1' : 'FREQ',
'CRVAL1' : 1378471216.4292786,
'CRPIX1' : 32,
'CUNIT1' : 'Hz',
'CDELT1' : 97647.745732,
'RESTFRQ': 1.420405752e+9
}
line = wcs.Projection(header).spectra('VELO-F2V')
pixels = range(30,35)
Vwcs = line.toworld1d(pixels)
for p,v in zip(pixels, Vwcs):
print p, v/1000
# Output:
# 30 9023.78022672
# 31 9002.56055595
# 32 8981.34229811
# 33 8960.12545322
# 34 8938.9100213
Frequency to Wavelength
------------------------
The rest wavelength is given by the relation:
.. math::
:label: eq380
\lambda_0 = \frac{c}{\nu_0}
Inserting the right numbers we find:
.. math::
:label: eq390
\lambda_0 = \frac{299792458.0}{1420405752.0} = 0.211061140507\ m
For the barycentric wavelength we need to insert the barycentric frequency.
.. math::
:label: eq400
\lambda = \frac{299792458.0}{1378471216.43} = 0.217481841062\ m
The increment in wavelength as function of the increment in
(barycentric) frequency is:
.. math::
:label: eq410
d\lambda = \frac{-c}{\nu^2} d\nu
With the right numbers:
.. math::
:label: eq420
d\lambda = \frac{-299792458.0}{1378471216.43^2}\ 97647.745732 = -1.54059158176\times 10^{-5}\ m
This gives us the alternate header keywords::
RESTWAVZ= 0.211061140507 / [m]
::
CTYPE3W= 'WAVE-F2W'
CRVAL3W= 0.217481841062 / [m]
CDELT3W= -1.5405916e-05 / [m]
CUNIT3W= 'm'
RESTWAVW= 0.211061140507 / [m]
Note that CTYPE indicates that there is a non linear conversion from frequency
to wavelength.
From the standard definition of optical velocity:
.. math::
:label: eq430
Z = c\ \frac{\lambda-\lambda_0}{\lambda_0}
it follows that the increment in optical velocity as function of increment of wavelength
is given by:
.. math::
:label: eq440
dZ = \frac{c}{\lambda_0}\ d\lambda
Then with the numbers we find:
.. math::
:label: eq450
dZ_b = \frac{299792458.0}{0.211061140507}\times -1.54059158176\times 10^{-5}= -21882.6514422\ m/s
which is the increment in optical velocity earlier given for CDELT3Z.
This is one of the possible conversions between wavelength and velocity. Others are listed
in `scs.pdf `_
table 3 of E.W. Greisen et al. page 750.
Conclusions
-----------
* Note that the inertial system is set by a (FITS) header using a special keyword
(e.g. VELREF=) or it is coded in the CTYPEn keyword. It doesn't change anything in
the calculations above. Conversions between inertial reference systems is not possible because
headers do (usually) not contain the relevant information to calculate the topocentric
correction w.r.t. that system (one needs time of observation, position of observatory
and position of the observed source).
* From a header with CTYPEn='FREQ' we can derive optical, radio and apparent radial velocities
with method *spectra()*:
* *proj = wcs.Projection(header).spectra('VOPT-F2W')*
* *proj = wcs.Projection(header).spectra('VRAD')*
* *proj = wcs.Projection(header).spectra('VELO-F2V')*
This applies also to alternate axis descriptions. So if CTYPE1='VRAD' one can derive
one of the other velocity definitions by adding the *spectra()* method with
the appropriate argument.
Here is an example:
.. code-block:: python
#!/usr/bin/env python
from kapteyn import wcs
wcs.debug = True
header = { 'NAXIS' : 1,
'CTYPE1' : 'VRAD',
'CRVAL1' : 8850750.904193053,
'CRPIX1' : 32,
'CUNIT1' : 'm/s',
'CDELT1' : -20609.644582145629,
'RESTFRQ': 1.420405752e+9
}
line = wcs.Projection(header).spectra('VOPT-F2W')
pixels = range(30,35)
Vwcs = line.toworld1d(pixels)
for p,v in zip(pixels, Vwcs):
print p, v/1000
# Output:
# Velocities in km/s converted from 'VRAD' to 'VOPT-F2W'
# 30 9163.77150423
# 31 9141.88420167
# 32 9120.0
# 33 9098.11889856
# 34 9076.2408967
Note that the rest frequency is required now.
Note also that we added statement *wcs.debug = True* to get some debug
information from WCSLIB.
* Axis types 'FREQ-HEL' and 'FREQ-LSR' (AIPS definitions) are recognized by WCSLIB
and are treated as 'FREQ'. No conversions are done. Internally the keyword *SPECSYS=*
gets a value.
The complete alternate axis descriptions
-----------------------------------------
In this section we summarize the alternate axis descriptions and we add
a small script that proves that these descriptions are consistent::
CNAME= 'Topocentric Frequency. Basic header'
CTYPE3= 'FREQ'
CRVAL3= 1.37835117405e9
CDELT3= 9.765625e4
CRPIX3= 32
CUNIT3= 'Hz'
RESTFRQ= 1.420405752e+9
SPECSYS='TOPOCENT'
CNAME3Z= 'Barycentric optical velocity'
RESTWAVZ= 0.211061140507 / [m]
CTYPE3Z= 'VOPT-F2W'
CRVAL3Z= 9.120e+6 / [m/s]
CDELT3Z= -2.1882651e+4 / [m/s]
CRPIX3Z= 32
CUNIT3Z= 'm/s'
SPECSYSZ='BARYCENT' / Velocities w.r.t. barycenter
SSYSOBSZ='TOPOCENT' / Observation was made from the 'TOPOCENT' frame
VELOSYSZ= 26108 / [m/s]
CNAME3F= 'Barycentric frequency'
CTYPE3F= 'FREQ'
CRVAL3F= 1.37847121643e+9 / [Hz]
CDELT3F= 9.764775e+4 / [Hz]
CRPIX3F= 32
CUNIT3F= 'Hz'
RESTFRQF= 1.420405752e+9
SPECSYSF='BARYCENT'
SSYSOBSF='TOPOCENT'
VELOSYSF= 26108 / [m/s]
CNAME3R= 'Barycentric radio velocity'
CTYPE3R= 'VRAD'
CRVAL3R= 8.85075090419e+6 / [m/s]
CDELT3R= -2.0609645e+4 / [m/s]
CRPIX3R= 32
CUNIT3R= 'm/s'
RESTFRQR= 1.420405752e+9
SPECSYSR='BARYCENT'
SSYSOBSR='TOPOCENT'
VELOSYSR= 26108 / [m/s]
CNAME3V= 'Barycentric apparent radial velocity'
RESTFRQV= 1.420405752e+9 / [Hz]
CTYPE3V= 'VELO-F2V'
CRVAL3V= 8.98134229811e+6 / [m/s]
CDELT3V= -2.1217551e+4 / [m/s]
CRPIX3V= 32
CUNIT3V= 'm/s'
SPECSYSV='BARYCENT'
SSYSOBSV='TOPOCENT'
VELOSYSV= 26108 / [m/s]
CNAME3W= 'Barycentric wavelength'
CTYPE3W= 'WAVE-F2W'
CRVAL3W= 0.217481841062 / [m]
CDELT3W= -1.5405916e-05 / [m]
CRPIX3W= 32
CUNIT3W= 'm'
RESTWAVW= 0.211061140507 / [m]
SPECSYSW='BARYCENT'
SSYSOBSW='TOPOCENT'
VELOSYSW= 26108 / [m/s]
To check the validity and completeness of these alternate axis descriptions,
we wrote a small script that loops over all the mnemonic letter codes in a header that
is composed from the header fragments above. We only changed axisnumber 3 to 1.
The output is the same within the boundaries of the given precision of the numbers.
To change the axis description in a header we use the *alter* parameter
when we create the projection object.
Parameter *alter* is an optional
letter from 'A' through 'Z', indicating an alternative WCS axis description:
.. code-block:: python
#!/usr/bin/env python
from kapteyn import wcs
header = { 'NAXIS' : 1,
'CTYPE1' : 'FREQ',
'CRVAL1' : 1378471216.4292786,
'CRPIX1' : 32,
'CUNIT1' : 'Hz',
'CDELT1' : 97647.745732,
'RESTFRQ' : 1.420405752e+9,
'CNAME1Z' : 'Barycentric optical velocity',
'RESTWAVZ' : 0.211061140507, # [m]
'CTYPE1Z' : 'VOPT-F2W',
'CRVAL1Z' : 9.120e+6, # [m/s]
'CDELT1Z' : -2.1882651e+4, # [m/s]
'CRPIX1Z' : 32,
'CUNIT1Z' : 'm/s',
'SPECSYSZ' : 'BARYCENT', # Velocities w.r.t. barycenter,
'SSYSOBSZ' : 'TOPOCENT', # Observation was made from the 'TOPOCENT' frame,
'VELOSYSZ' : 26108, # [m/s]
'CNAME1F' : 'Barycentric frequency',
'CTYPE1F' : 'FREQ',
'CRVAL1F' : 1.37847121643e+9, # [Hz]
'CDELT1F' : 9.764775e+4, # [Hz]
'CRPIX1F' : 32,
'CUNIT1F' : 'Hz',
'RESTFRQF' : 1.420405752e+9,
'SPECSYSF' : 'BARYCENT',
'SSYSOBSF' : 'TOPOCENT',
'VELOSYSF' : 26108, # [m/s]
'CNAME1W' : 'Barycentric wavelength',
'CTYPE1W' : 'WAVE-F2W',
'CRVAL1W' : 0.217481841062, # [m]
'CDELT1W' : -1.5405916e-05, # [m]
'CRPIX1W' : 32,
'CUNIT1W' : 'm',
'RESTWAVW' : 0.211061140507, # [m]
'SPECSYSW' : 'BARYCENT',
'SSYSOBSW' : 'TOPOCENT',
'VELOSYSW' : 26108, # [m/s]
'CNAME1R' : 'Barycentric radio velocity',
'CTYPE1R' : 'VRAD',
'CRVAL1R' : 8.85075090419e+6, # [m/s]
'CDELT1R' : -2.0609645e+4, # [m/s]
'CRPIX1R' : 32,
'CUNIT1R' : 'm/s',
'RESTFRQR' : 1.420405752e+9,
'SPECSYSR' : 'BARYCENT',
'SSYSOBSR' : 'TOPOCENT',
'VELOSYSR' : 26108, # [m/s]
'CNAME1V' : 'Barycentric apparent radial velocity',
'CTYPE1V' : 'VELO-F2V',
'CRVAL1V' : 8.98134229811e+6, # [m/s]
'CDELT1V' : -2.1217551e+4, # [m/s]
'CRPIX1V' : 32,
'CUNIT1V' : 'm/s',
'RESTFRQV' : 1.420405752e+9, # [Hz]
'SPECSYSV' : 'BARYCENT',
'SSYSOBSV' : 'TOPOCENT',
'VELOSYSV' : 26108 # [m/s]
}
# Loop over all the alternative headers
for alt in ['F', 'Z', 'W', 'R', 'V']:
spec = wcs.Projection(header, alter=alt).spectra('VOPT-F2W')
pixels = range(30,35)
Vwcs = spec.toworld1d(pixels)
cname = header['CNAME1'+alt] # Just a header text
print "VOPT-F2W from %s" % (cname,)
print "Pixel, velocity (%s)" % spec.units
for p,v in zip(pixels, Vwcs):
print p, v/1000.0
# Output
# VOPT-F2W from Barycentric frequency
# Pixel, velocity (m/s)
# 30 9163.77150598
# 31 9141.88420246
# 32 9119.99999984
# 33 9098.11889745
# 34 9076.24089463
# VOPT-F2W from Barycentric optical velocity
# Pixel, velocity (m/s)
# 30 9163.77150335
# 31 9141.88420123
# 32 9120.0
# 33 9098.11889901
# 34 9076.24089759
# VOPT-F2W from Barycentric wavelength
# Pixel, velocity (m/s)
# 30 9163.77150495
# 31 9141.88420213
# 32 9120.0000002
# 33 9098.1188985
# 34 9076.24089638
# VOPT-F2W from Barycentric radio velocity
# Pixel, velocity (m/s)
# 30 9163.77150512
# 31 9141.88420211
# 32 9120.0
# 33 9098.11889812
# 34 9076.24089581
# VOPT-F2W from Barycentric apparent radial velocity
# Pixel, velocity (m/s)
# 30 9163.77150347
# 31 9141.88420129
# 32 9120.0
# 33 9098.11889894
# 34 9076.24089746
Alternative conversions
+++++++++++++++++++++++++
Conversion between radio and optical velocity
-----------------------------------------------
In the next two sections we give some formula's that could be handy if you want to verify
numbers. They are not used in WCSLIB.
With the definitions for radio and optical velocity it is easy to derive:
.. math::
:label: eq500
\frac{V}{Z} = \frac{\nu}{\nu_0}
This can be verified with:
* Z = 9120000.00000 m/s
* V = 8850750.90419 m/s
* :math:`\nu_0` = 1420405752.00 Hz
* :math:`\nu_b` = 1378471216.43 Hz
Both ratios are equal to 1.030421045482.
Conversion between apparent radial velocity and optical/radio velocity
-------------------------------------------------------------------------
It is possible to find a relation between the true velocity and the optical velocity
using eq. :eq:`eq10` and eq. :eq:`eq50`.
The apparent radial velocity can be written as:
.. math::
:label: eq510
\frac{v}{c} = \frac{\frac{\nu_0^2}{\nu^2}-1}{\frac{\nu_0^2}{\nu^2}+1}
The frequency shift for an optical velocity is:
.. math::
:label: eq520
\frac{\nu_0}{\nu} = \bigl(1+\frac{Z}{c}\bigr)
Then:
.. math::
:label: eq530
\frac{v}{c} = \frac{{(1+Z/c)}^2-1}{{(1+Z/c)}^2+1} = \frac{Z^2+2cZ}{Z^2+2cZ+2c^2}
This equation is used in AIPS memo 27 [Aipsmemo]_ to relate an optical velocity to an
apparent radial velocity.
If we insert :math:`Z_b` = 9120000 (m/s) then we find :math:`v_b` = 8981342.29811 (m/s) as expected
(eq. :eq:`eq50`, :eq:`eq310`)
For radio velocities we find in a similar way:
.. math::
:label: eq540
\frac{\nu_0}{\nu} = \frac{1}{\bigl(1-\frac{V}{c}\bigr)}
which gives the relation between apparent radial velocity and radio velocity:
.. math::
:label: eq550
\frac{v}{c} = \frac{2cV-V^2}{V^2-2cV+2c^2}
If we substitute the calculated barycentric radio velocity :math:`V_b` = 8850750.90419 (m/s)
then one finds again: :math:`v_b` = 8981342.29811 (m/s) (see also (eq. :eq:`eq50`, :eq:`eq310`)
Note that the last formula is equation 4 in AIPS memo 27 [Aipsmemo]_
Non-Linear Coordinate Systems in AIPS. However that formula lacks a minus sign
in the nominator and therefore does not give a correct result.
Legacy headers
++++++++++++++
.. _spectral_gipsy:
A recipe for modification of Nmap/GIPSY FITS data
--------------------------------------------------
For FITS headers produced by Nmap/GIPSY we don't have an increment
in velocity available so we cannot use them as input for WCSLIB (otherwise we
would treat them like the FELO axis recognized by AIPS). The Python interface to
WCSLIB applies a conversion for these headers before they are
processed by WCSLIB. From the previous steps we can summarize how
the data in the Nmap/GIPSY FITS header is changed:
* The extension in CTYPEn is '-OHEL', '-OLSR', '-RHEL' or '-RLSR'
* The velocity is retrieved from FITS keyword VELR= (always in m/s) or DRVALn= (in units of DUNITn)
* Convert reference frequency to a frequency in Hz.
* Calculate the reference frequency in the barycentric system using eq. :eq:`eq10`
if the velocity is optical and eq. :eq:`eq220` if the velocity is a radio velocity.
* Calculate the topocentric velocity using eq. :eq:`eq100`
* Convert frequency increment to an increment in Hz
* Calculate the increment in frequency in the selected reference system (HEL, LSR) using
eq. :eq:`eq140`.
* Change CRVALn and CDELTn to the barycentric values
* Change CTYPEn to 'FREQ'
* Create a projection object with spectral translation, e.g. **proj.spectra('VOPT-F2W')**
In the following script we show:
* the (invisible) conversion to the heliocentric system
* how to get the same output by applying the appropriate formulas
* the approximation that GIPSY uses
.. literalinclude:: EXAMPLES/gipsyspectralheader.py
:language: python
Output::
VELR is the reference velocity given in the velocity frame
coded in CTYPE (e.g. HEL, LSR)
The velocity is either an optical or a radio velocity. This
is also coded in CTYPE (e.g. 'O', 'R')
VOPT-F2W with spectral translation:
29 1000.194731
30 1016.794655
31 1033.396411
32 1050.000000
33 1066.605422
34 1083.212677
VOPT calculated:
29 1000.194731
30 1016.794655
31 1033.396411
32 1050.000000
33 1066.605422
34 1083.212677
VOPT with native GIPSY formula, which is an approximation:
29 1000.191559
30 1016.792540
31 1033.395354
32 1050.000000
33 1066.606480
34 1083.214793
The Python interface allows for an easy implementation for these special exceptions.
Here is a script that uses this facility. The conversion here is triggered by the CTYPE
extension **OHEL**. So as long this is unique to GIPSY spectral axes, you are safe to
use it. Note that we converted the frequencies to optical, radio and apparent radial velocities.
This is added value to the existing GIPSY implementation where these conversions are not
possible. These WCSLIB conversions are explained in previous sections:
.. code-block:: python
#!/usr/bin/env python
from kapteyn import wcs
header = { 'NAXIS' : 1,
'CTYPE1' : 'FREQ-OHEL',
'CRVAL1' : 1.37835117405e9,
'CRPIX1' : 32,
'CUNIT1' : 'Hz',
'CDELT1' : 9.765625e4,
'RESTFRQ': 1.420405752e+9,
'DRVAL1' : 9120000.0,
# 'VELR' : 9120000.0
'DUNIT1' : 'm/s'
}
proj = wcs.Projection(header)
pixels = range(30,35)
voptical = proj.spectra('VOPT-F2W')
Vwcs = voptical.toworld1d(pixels)
print "\nPixel, optical velocity (%s)" % voptical.units
for p,v in zip(pixels, Vwcs):
print p, v/1000.0
vradio = proj.spectra('VRAD')
Vwcs = vradio.toworld1d(pixels)
print "\nPixel, radio velocity (%s)" % vradio.units
for p,v in zip(pixels, Vwcs):
print p, v/1000.0
vradial = proj.spectra('VELO-F2V')
Vwcs = vradial.toworld1d(pixels)
print "\nPixel, apparent radial velocity (%s)" % vradial.units
for p,v in zip(pixels, Vwcs):
print p, v/1000.0
# Output:
# Pixel, optical velocity (m/s)
# 30 9163.77150423
# 31 9141.88420167
# 32 9120.0
# 33 9098.11889856
# 34 9076.2408967
#
# Pixel, radio velocity (m/s)
# 30 8891.97019336
# 31 8871.36054878
# 32 8850.75090419
# 33 8830.14125961
# 34 8809.53161503
#
# Pixel, apparent radial velocity (m/s)
# 30 9023.78022672
# 31 9002.56055595
# 32 8981.34229811
# 33 8960.12545322
# 34 8938.9100213
.. note::
Note that changing DRVAL1 to VELR gives the same output. Both are recognized as keywords
that store a velocity. The value in VELR should always be in m/s.
Note also how we created different sub-projections (one for each type of velocity)
from the same main projection. All these objects can coexist.
AIPS axis type FELO
--------------------
Next script and output shows that with the optical reference velocity and the corresponding
increment in velocity (CDELT3Z), we can get velocities without spectral translation.
WCSLIB recognizes the axis type 'FELO' which is regularly
gridded in frequency but expressed in velocity units in the optical convention. It is therefore not a
surprise that the output is the same as the list with optical velocities derived from
the spectral translation 'VOPT-F2W'.
We can prove this if we calculate the barycentric reference frequency and its increment. If
*Zr* is the optical reference velocity then we find the barycentric reference
frequency with:
.. math::
:label: eq552
\nu_r = \frac{\nu_0}{\bigl(1+\frac{Z_r}{c}\bigr)}
and from
.. math::
:label: eq553
dZ = \frac{-c}{\nu_0}\ {\bigl(1+\frac{Z_r}{c}\bigr)}^2\ d\nu
we derive:
.. math::
:label: eq554
d\nu = \frac{-\nu_0}{{\bigl(1+\frac{Z_r}{c}\bigr)}^2}\ dZ
which we rewrite in:
.. math::
:label: eq555
d\nu = \frac{-\nu_0 c} {{(c+Z_r)}^2}\ dZ
So if we have a barycentric reference velocity and a barycentric velocity increment, then
according to the formulas above it is easy to retrieve the values for the barycentric
reference frequency and the barycentric frequency increment. The script below proves that
indeed with these values the optical velocities are derived from a linear frequency axis
and not from a linear velocity axis (see the last option in this script):
.. literalinclude:: EXAMPLES/aipsfelo.py
:language: python
Output::
Pixel, velocity (km/s) with native header with FELO-HEL
30 9163.77150423
31 9141.88420167
32 9120.0
33 9098.11889857
34 9076.24089671
Calculated a reference frequency: 1378471216.43
Calculated a frequency increment: 97647.7457311
Pixel, velocity (km/s) with barycentric reference frequency and increment:
30 9163.77150423
31 9141.88420167
32 9120.0
33 9098.11889857
34 9076.24089671
Pixel, velocity (km/s) with spectral translation VOPT-F2W
30 9163.77150423
31 9141.88420167
32 9120.0
33 9098.11889857
34 9076.24089671
Pixel, velocity (km/s) with CUNIT='FELO', which is unrecognized
and therefore linear. This deviates from the previous output.
The second velocity is calculated manually.
30 9163.76530288 9163.76530288
31 9141.88265144 9141.88265144
32 9120.0 9120.0
33 9098.11734856 9098.11734856
34 9076.23469712 9076.23469712
So in this script we demonstrated the use of a special velocity axis type which
originates from a classic AIPS data FITS file.
It is called 'FELO'. WCSLIB (and not our Python interface) recognizes this type as an
**optical velocity** and performs
the necessary internal conversions as we can see in the source code:
.. code-block:: python
if (strcmp(wcs->ctype[i], "FELO") == 0) {
strcpy(wcs->ctype[i], "VOPT-F2W");
The source code also reveals that the extensions in CUNITn are translated into values
for FITS keyword *SPECSYS*:
.. code-block:: python
if (strcmp(scode, "-LSR") == 0) {
strcpy(wcs->specsys, "LSRK");
} else if (strcmp(scode, "-HEL") == 0) {
strcpy(wcs->specsys, "BARYCENT");
} else if (strcmp(scode, "-OBS") == 0) {
strcpy(wcs->specsys, "TOPOCENT");
**Conclusions**
* The extension HEL or LSR after FELO in *CTYPE1* is not used in the calculations.
But when you omit
a valid extension the axis will be treated as a linear axis.
* In the example above one can replace *FELO-HEL* in *CTYPE1*
by FITS standard *VOPT-F2W* showing that for WCSLIB *FELO-HEL*
is in fact the same as *VOPT-F2W*.
AIPS axis type VELO
-------------------
In this section we want to address the question what WCSLIB does if it
encounters an AIPS VELO-XXX axis as in *CTYPE1='VELO-HEL'* or *'VELO-LSR'*.
From the AIPS documentation we learn that VELO is regularly gridded in
velocity (m/s) in the optical convention, unless overridden by use of
the *VELREF* keyword. VELREF is an integer. From the documentation
of WCSLIB we learn that for Classic Aips:
1. LSR kinematic, originally described simply as "LSR" without distinction
between the kinematic and dynamic definitions.
2. Barycentric, originally described as "HEL" meaning heliocentric.
3. Topocentric, originally described as "OBS" meaning geocentric
but widely interpreted as topocentric.
And for AIPS++ extensions to VELREF which are also recognized:
4. LSR dynamic.
5. Geocentric.
6. Source rest frame.
7. Galactocentric.
.. note::
From the WCSLIB documentation:
For an AIPS 'VELO' axis, a radio convention velocity is denoted by adding 256 to VELREF,
otherwise an optical velocity is indicated (not applicable to 'FELO' axes).
Unrecognized values of VELREF are simply ignored.
VELREF takes precedence over CTYPEia in defining the Doppler frame.
.. note::
Only WCSLIB (versions >= 4.5.1) do recognize keyword *VELREF*.
We show the use of *VELREF* with the following script:
.. literalinclude:: EXAMPLES/aipsvelo.py
:language: python
Output::
The velocity increment is constant and equal to 5.000000 (km/s):
Allowed spectral translations [('FREQ', 'Hz'), ('ENER', 'J'), ('WAVN', '/m'),
('VOPT-F2W', 'm/s'), ('VRAD', 'm/s'), ('VELO-F2V', 'm/s'), ('WAVE-F2W', 'm'),
('ZOPT-F2W', ''), ('AWAV-F2A', 'm'), ('BETA-F2V', '')]
T1. Radio velocity directly from header and optical velocity
from spectral translation. VELO is a radio velocity here because
VELREF > 256
Pixel Vradio in (km/s) and Voptical (km/s)
30 -253.000000 -252.786669
31 -248.000000 -247.795014
32 -243.000000 -242.803193
33 -238.000000 -237.811206
34 -233.000000 -232.819052
T2. Now insert CTYPE1='VRAD' in the header and convert to VOPT-F2W
with a spectral translation (Z1) and with a calculation (Z2)
This should give the same results as in table T1.
With CTYPE='RAD' and spec.trans 'VOPT-F2W': Pixel , Vrad, Z1 (km/s), Z2 (km/s)
30 -253.0 -252.786668992 -252.786668992
31 -248.0 -247.795014311 -247.795014311
32 -243.0 -242.803193261 -242.803193261
33 -238.0 -237.811205834 -237.811205834
34 -233.0 -232.819052022 -232.819052022
T3. We set CTYPE1 to VELO-HEL and VELREF to 2 (Helio) and
derive optical and radio velocities from it. Compare these with
the relativistic velocity in Table T4.
Allowed spectral translations for VELO as optical velocity [('FREQ-W2F', 'Hz'),
('ENER-W2F', 'J'), ('WAVN-W2F', '/m'), ('VOPT', 'm/s'), ('VRAD-W2F', 'm/s'),
('VELO-W2V', 'm/s'), ('WAVE', 'm'), ('ZOPT', ''), ('AWAV-W2A', 'm'),
('BETA-W2V', '')]
Pixel Voptical in (km/s) and Vradio (km/s)
30 -253.000000 -253.213691
31 -248.000000 -248.205325
32 -243.000000 -243.197126
33 -238.000000 -238.189094
34 -233.000000 -233.181229
T4. Next a list with optical velocities calculated from relativistic
velocity with constant increment.
If these values are different from the previous optical velocity then
obviously the velocities derived from the header are not relativistic
as in pre 4.5.1 versions of WCSLIB.
30 -252.893335
31 -247.897507
32 -242.901597
33 -237.905603
34 -232.909526
We used eq. :eq:`eq30` to calculate a frequency for a given apparent radial velocity.
This frequency is used in eq. :eq:`eq5` to calculate the optical velocity.
The script proves:
* Axis VELO-HEL is processed as an optical velocity and if keyword *VELREF*
is present and its value is greater than 256, then VELO-HEL is processed
as a radio velocity. In versions of WCSLIB < 4.5.1, the VELO-XXX axis was
processed as VELO i.e. a relativistic velocity.
.. note::
From the WCSLIB API documentation:
AIPS-convention celestial projection types, NCP and GLS, and spectral types,
'{FREQ,FELO,VELO}-{OBS,HEL,LSR}' as in
'FREQ-LSR', 'FELO-HEL', etc., set in CTYPEia are translated on-the-fly by
wcsset() but without modifying the relevant ctype[], pv[] or specsys members
of the wcsprm struct. That is, only the information extracted from ctype[]
is translated when wcsset() fills in wcsprm::cel (celprm struct) or wcsprm::spc (spcprm struct).
On the other hand, these routines do change the values of wcsprm::ctype[],
wcsprm::pv[], wcsprm::specsys and other wcsprm struct members as appropriate to
produce the same result as if the FITS header itself had been translated.
Definitions and formulas from AIPS and GIPSY
--------------------------------------------
AIPS
*****
A radio velocity is defined by:
.. math::
:label: eq600
V = c \bigl( \frac{\nu_0 - \nu^{'}}{\nu_0} \bigr)
where :math:`\nu` is the Doppler shifted rest frequency, given by:
.. math::
:label: eq610
\nu' = \nu_0\sqrt{(\frac{c-v}{c+v})}
Equivalent to the relativistic addition of apparent radial velocities we can derive a relation for
radio velocities if the velocities in given in different reference systems.
The addition of apparent radial velocities is given in
AIPS memo 27 [Aipsmemo]_
Non-Linear Coordinate Systems in AIPS (Eric W. Greisen, NRAO) Greisen, is
.. math::
:label: eq620
v = \frac{v_s + v _{obs}}{1 + \frac{v_s v_{obs}}{c^2}}
To stay close to our previous examples and definitions we set :math:`v_s`
which is the apparent radial velocity of an object w.r.t. an inertial system, to
be equal to :math:`v_b` (our inertial system in this case is barycentric).
The other velocity, :math:`v_{obs}` is equal to the topocentric correction: :math:`v_t`
and the result :math:`v = v_e`, the apparent radial velocity of the object as we would observe
it on earth.
Then we get the familiar formula (eq. :eq:`eq70`):
.. math::
:label: eq630
v_e = \frac{v_b + v_t}{1 + \frac{v_b v_t}{c^2}}
With the relation between V and v and the relativistic addition of velocities we find that
the radio velocities in different systems are related according to the equation:
.. math::
:label: eq640
V_e = V_b + V_t - V_b V_t/ c
(see also AIPS memo 27 [Aipsmemo]_ ).
The barycentric radio velocity was calculated in a previous section.
Its value was :math:`V_b` = 8850750.90404 m/s. With the topocentric reference frequency
1378351174.05 Hz we find :math:`V_e` = 8876087.18567 m/s. We know from fig. 1 that the
topocentric correction is positive. To calculate the corresponding radio velocity :math:`V_t`
we use:
.. math::
:label: eq650
V_t = c(\frac{\nu_b - \nu_e}{\nu_b}) = 299792458.0\times\frac{(1378471216.43-1378351174.05)}{1378471216.43}=26107.03781\ m/s
With these values for :math:`V_b` and :math:`V_t` you can verify that the
expression for :math:`V_e` is valid.
.. math::
:label: eq660
V_e = 8850750.90404 +26107.03781 - \frac{8850750.90404 \times 26107.03781}{299792458.0} = 8876087.18567\ m/s
which is the value of :math:`V_e` that we found before using the topocentric reference frequency, so we can have
confidence in the relation for radio velocities as found in the AIPS memo [Aipsmemo]_ .
But this radio velocity :math:`V_e` (w.r.t. observer on Earth) for a pixel N is also given by the relation:
.. math::
:label: eq670
V_e(N) = -\frac{c}{\nu_0}(\nu_e(N)-\nu_0) = -\frac{c}{\nu_0}(\nu_e+\delta_{\nu}(N-N_{\nu})-\nu_0)
It is important to emphasize the meaning of the variables:
* :math:`\nu_e` = topocentric reference frequency).
* :math:`\delta_\nu` = the increment in frequency per pixel in the topocentric system
* :math:`N_\nu` = the frequency reference pixel
* :math:`N` = the pixel
If we use the previous formulas we can also write:
.. math::
:label: eq680
V_e(N_V) = V'_b + V_t - V'_b V_t/ c
.. math::
:label: eq690
V_e(N_V) = -\frac{c}{\nu_0}(\nu_e+\delta_{\nu}(N_V-N_{\nu})-\nu_0)
The velocity :math:`V^{'}_b` is the barycentric reference velocity
at velocity reference pixel :math:`N_V`.
From these relations we observe:
.. math::
:label: eq700
V_b(N) = \frac{V_e(N)-V_t}{1-\frac{V_t}{c}}
and from eq. :eq:`eq680` with :math:`V^{'}_b = V_b(N_V)`:
.. math::
:label: eq710
V_t = \frac{V_e(N_V)-V_b(N_V)}{1-\frac{V_b(N_V)}{c}}
Using also the equations with the frequencies, we can derive the following expression
for :math:`V_b(N)`:
.. math::
:label: eq720
V_b(N) = V_b(N_V) - \frac{\delta_\nu \bigl(c-V_b(N_V)\bigr)(N-N_V)}{\nu_e + \delta_\nu (N_V-N_\nu)}
or in an alternative notation:
.. math::
:label: eq730
V_b(N) = V_b(N_V) + \delta_V (N-N_V)
Note that in AIPS memo 27 [Aipsmemo]_ the variable :math:`V_R` is used for :math:`V_b(N_V)` and
:math:`V_R` and :math:`N_V` are stored in AIPS headers as alternative
reference information (if frequency is in the main axis description).
The difference between the velocity and frequency reference pixel can be expressed in terms of
the radio velocities :math:`V_b(N_V)` and :math:`V_b(N_\nu)`.
It follows from eq. :eq:`eq720`) that for :math:`N = N_{\nu}` and a little rearranging:
.. math::
:label: eq740
N_V - N_{\nu} = \frac{\nu_e\ \bigl[V_b(N_{\nu}) - V_b(N_V)\bigr]}{\delta_\nu\ \bigl[c-V_b(N_{\nu})\bigr]}
We conclude that either one calculates (barycentric) radio velocities using the reference frequency
and the frequency increment from the header, or one calculates these velocities using
a reference velocity and a velocity increment from the header.
Note that we assumed that the frequency increment in the barycentric system
is the same as in the the system of the observer, which is not correct.
However the differences are small (less than 0.01% for 100 pixels from the reference pixel
for typical observations as in our examples).
For optical velocities Greisen derives:
.. math::
:label: eq750
Z_e = Z_b + Z_{t} + Z_b Z_{t} / c
and:
.. math::
:label: eq760
Z_b(N) = Z_b(N_V) - \frac{\delta_\nu\ \bigl(c+Z_b(N_V)\bigr)\ (N-N_Z)}{\nu_e + \delta_\nu (N-N_\nu)}
The difference between the velocity and frequency reference pixels in terms of
optical velocity is:
.. math::
:label: eq770
N_Z - N_{\nu} = \frac{\nu_e\ \bigl[Z_b(N_{\nu}) - Z_b(N_Z)\bigr]}{\delta_\nu\ \bigl[c+Z_b(N_{\nu})\bigr]}
Next script demonstrates how we reconstruct the topocentric optical velocity and the reference pixel for
that velocity as it is used in the AIPS formula. Then we compare the output of the WCSLIB method
and the AIPS formula::
#!/usr/bin/env python
from kapteyn import wcs
import numpy
c = 299792458.0 # m/s From literature
f0 = 1.42040575200e+9 # Rest frequency HI (Hz)
fR = 1.37835117405e+9 # Topocentric reference frequency (Hz)
dfR = 9.765625e+4 # Increment in topocentric frequency (Hz)
fb = 1.3784712164292786e+9 # Barycentric reference frequency (Hz)
dfb = 97647.745732 # Increment in barycentric frequency (Hz)
Zb = 9120.0e+3 # Barycentric optical velocity in m/s
Nf = 32 # Reference pixel for frequency
header = { 'NAXIS' : 1,
'CTYPE1' : 'FREQ',
'CRVAL1' : fb,
'CRPIX1' : Nf,
'CUNIT1' : 'Hz',
'CDELT1' : dfR,
'RESTFRQ': f0
}
line = wcs.Projection(header).spectra('VOPT-F2W')
pixels = numpy.array(range(30,35))
Vwcs = line.toworld1d(pixels) / 1000
print """Optical velocities from WCSLIB with spectral
translation and with barycentric ref. freq. (km/s):"""
for p,v in zip(pixels, Vwcs):
print p, v
# Select an arbitrary velocity reference pixel
Nz = 44.0
# then calculate corresponding velocity
Zb2 = (fR*Zb-dfR*c*(Nz-Nf))/(fR+dfR*(Nz-Nf))
print "Zb(Nz) =", Zb2
dN = fR*(Zb-Zb2)/(dfR*(c+Zb2))
Nz = dN + Nf
print "Closure test for selected reference pixel: Nz=", Nz
print "\nOptical velocities using AIPS formula (km/s):"
Zs = Zb2 - dfR*(c+Zb2)*(pixels-Nz)/(fR+dfR*(pixels-Nf))
Zs /= 1000
for p,z in zip(pixels, Zs):
print p, z
fx = fR + dfR*(Nz-Nf)
dZ = -dfR*(c+Zb2) / fx
print "Velocity increment: ", dZ
header = { 'NAXIS' : 1,
'CTYPE1' : 'VOPT-F2W',
'CRVAL1' : Zb2,
'CRPIX1' : Nz,
'CUNIT1' : 'm/s',
'CDELT1' : dZ,
'RESTFRQ': f0
}
line2 = wcs.Projection(header)
Vwcs = line2.toworld1d(pixels) / 1000
print """\nOptical velocities from WCSLIB without spectral
translation with barycentric Z (km/s):"""
for p,v in zip(pixels, Vwcs):
print p, v
# Output:
# Optical velocities from WCSLIB with spectral
# translation and with barycentric ref. freq. (km/s):
# 30 9163.77531689
# 31 9141.88610773
# 32 9120.0
# 33 9098.11699305
# 34 9076.23708621
# Zb(Nz) = 8857585.54671
# Closure test for selected reference pixel: Nz= 44.0
#
# Optical velocities using AIPS formula (km/s):
# 30 9163.77912988
# 31 9141.88801395
# 32 9120.0
# 33 9098.11508736
# 34 9076.23327538
# Velocity increment: -21849.2948239
#
# Optical velocities from WCSLIB without spectral
# translation with barycentric Z (km/s):
# 30 9163.77912988
# 31 9141.88801395
# 32 9120.0
# 33 9098.11508736
# 34 9076.23327538
Note that we used the topocentric frequency increment in the WCSLIB call
for a better comparison with the AIPS formula. The output of velocities with the AIPS formula
is exactly the same as WCSLIB with optical velocities using the velocity increment
calculated with the AIPS method (as to be expected). And these velocities
are very close to the velocities calculates with WCSLIB using the barycentric
frequency that corresponds to the given optical velocity. The differences can be explained
by the fact that the different methods are used to calculate a velocity increment.
What did we prove with this script? We selected an arbitrary pixel as reference pixel for
the velocity. This velocity has a relation with the initial optical velocity (9120 km/s)
through the difference in reference pixels. We calculated that velocity and showed that
the AIPS formula generates results that are almost equal to WCSLIB with the barycentric
reference frequency. If we use the AIPS formulas to calculate a velocity increment,
we can use the values in WCSLIB if we set CTYPE to 'VOPT-F2W'. This generates exactly
the same results as with the AIPS formula for velocities. So in frequency mode
WCSLIB calculates topocentric frequencies (and *topocentric* velocities
if we use a spectral translation method)
and in velocity mode it calculates barycentric velocities. AIPS axis type FELO can be
used as input for WCSLIB without modification.
**Conclusions**
* In AIPS the reference pixel for the reference velocity differs from the
frequency reference pixel. There is a relation between this reference velocity and
the barycentric velocity and these reference pixels. To us it is not clear
what this reference velocity represents and why it is not changed to a velocity at
the same reference pixel as the frequency.
* In the AIPS approach it is assumed that the increment in frequency is the same in
different reference systems. This assumption is not correct, but the deviations are
usually very small.
GIPSY
*****
The formulas used in GIPSY to convert frequencies to velocities are described in section:
`spectral coordinates `_
in the GIPSY programmers guide.
There is a formula for optical velocities and one for radio velocities.
Both formulas are derived from the standard formulas for velocities but the result is
split into a reference velocity and a part that is a non linear function of the increment in
frequency.
Optical
^^^^^^^^^
For optical velocities we use symbol *Z*.
The conversion from frequencies to **optical** velocities is not linear. One can try to
approximate a constant step in velocity, and to apply the standard linear
transformation :math:`Z(N) = Z_r + (N-crpix) \times dZ`, but this approximation can
deviate significantly in certain circumstances.
Therefore most reduction and analysis packages provide
functionality to calculate velocities also for the non-linear cases. Like Classic AIPS,
GIPSY provides a system for these transformations (e.g. function ``velpro.c``), but
it turns out that these transformations are also approximations because
where a barycentric or lsrk frequency should be used, GIPSY uses values from the
FITS header and for FITS files made by Newstar/Nmap for data observed before 2006-07-03, these
frequencies are topocentric.
In this section we show how GIPSY transforms frequencies to optical velocities. Also we
derive formulas for a linear transformation (i.e. for a constant velocity increment)
which can be used if one wants to compose a modified header for a linear transformation
:math:`Z(N) = Z_r + (N-crpix) \times dZ`
Given a barycentric (or lsrk) frequency one calculates an optical velocity *Z*
in that system with:
.. math::
:label: eq850
Z = -c (\frac{\nu_b-\nu_0}{\nu_b})
Assume for channel :math:`N`:
.. math::
:label: eq852
\nu(N) = \nu_{br} + (N-N_{ref}) \delta_{\nu_b} = \nu_{br} + {\bf n} \delta_{\nu_b}
For :math:`(N-N_{ref})` we wrote :math:`\bf n`.
The frequencies are related to the barycentric (or lrsk) reference system.
:math:`N_{ref}` is the reference pixel (*CRPIX*) given in a FITS header,
:math:`\nu_{br}` is the reference frequency in this barycentric system and
:math:`\delta_{\nu_b}` is the barycentric frequency increment.
Inserting :eq:`eq852` into :eq:`eq850` gives:
.. math::
:label: eq854
Z(N) = -c\bigl(\frac{\nu_{br}+\bf n\delta_{\nu_b}-\nu_0}{\nu_{br}+\bf n\delta_{\nu_b}}\bigr) = -c \bigl(\frac{\nu_{br}-\nu_0}{\nu_{br}}\bigr) + {\bf n}dZ = Z_r + {\bf n}dZ
:math:`Z_r` is the given reference velocity in the barycentric/lsrk reference system.
Solve this equation for :math:`{\bf n}dZ` to get an expression for the increment:
.. math::
:label: eq856
{\bf n}dZ = {\bf n}\ \frac{-c\nu_0 \delta_{\nu_b}}{(\nu_{br}+{\bf n}\delta_{\nu_b})\nu_{br}} = c \nu_0 \bigl(\frac{1}{(\nu_{br}+{\bf n}\delta_{\nu_b})} - \frac{1}{\nu_{br}}\bigr)
The formula to calculate optical velocities then becomes:
.. math::
:label: gipsynonlinear
Z(N) = Z_r + c \nu_0 \bigl(\frac{1}{(\nu_{br}+{\bf n}\delta_{\nu_b})} - \frac{1}{\nu_{br}}\bigr)
with:
* :math:`Z(N)` is the barycentric optical velocity for pixel :math:`N`
* :math:`\nu_{br}` is the barycentric reference frequency
* :math:`\delta_{\nu_b}` is the increment in barycentric frequency
**This is the formula that GIPSY uses to calculate optical velocities. However, GIPSY
uses the topocentric reference frequency and the topocentric frequency increment.**
If we want to express the optical velocity at pixel N as a function
of the reference velocity and a **constant** velocity increment as in
:math:`Z(N) = Z_r + {\bf n}dZ`, then we need to find an expression for *dZ* which does not depend
on *n*. Rewrite *ndZ* into:
.. math::
:label: eq862
{\bf n}dZ = {\bf n} \frac{-c \nu_0 \delta _{\nu_b}}{(\nu_{br}+n\delta_{\nu_b}) \nu_{br}}
Then, with the observation that :math:`{\bf n}\delta_{\nu_b} << \nu_{br}`:
.. math::
:label: eq864
{\bf n}dZ \approx {\bf n} \frac{-c \nu_0 \delta_{\nu_b}}{{\nu_{br}}^2}
and thereby:
.. math::
:label: eq866
dZ \approx \frac{-c\nu_0\delta_{\nu_b}}{{\nu_{br}}^2}
This is the formula that is documented in the programmers manual to get a
value for GIPSY's keyword *DDELT* (one of the alternative keywords from the list
DRVAL, DDELT, DRPIX, DUNIT which describe an alternative coordinate system with a higher
priority than the system described by the corresponding keywords that start with 'C').
However the formula is never used in GIPSY to explicitly set the value of DDELT.
Only when DDELT is given in a header, it is used as an increment.
So the formula to calculate optical velocities, without the use of
the rest frequency, is:
.. math::
:label: gipsylinearwithf0
Z(N) = Z_r + {\bf n}\frac{-c\nu_0\delta_{\nu_b}}{{\nu_{br}}^2}
In the formulas above we included the rest frequency. But it is not necessary to
know its value because we can express this rest frequency in terms of optical
velocity:
.. math::
:label: eq868a
Z = -c (\frac{\nu_b-\nu_0}{\nu_b}) \rightarrow \nu_0 = \nu_{br} \bigl(1+\frac{Z_r}{c} \bigr)
Then:
.. math::
:label: eq868b
Z(N) = Z_r + c \nu_{br} \bigl(1+\frac{Z_r}{c}\bigr) \bigl(\frac{1}{(\nu_{br}+{\bf n}\delta_{\nu_b})} - \frac{1}{\nu_{br}}\bigr)
from which we derive in a straightforward way:
.. math::
:label: gipsynonlinearwithoutf0
Z(N) = \frac{Z_r\nu_{br} - c {\bf n}\delta_{\nu_b}}{\nu_{br} + {\bf n} \delta_{\nu_b}}
**The formula above is the method used by GIPSY's function velpro.c to get
velocities if the rest frequency is unknown.**
And again, if we want to express the optical velocity at pixel N as a function
of the reference velocity and a **constant** velocity increment as in
:math:`Z(N) = Z_r + {\bf n}dZ` then we need to find an expression for *dZ* which does not depend
on *n*. Note that :math:`{\bf n}\delta \nu_b << \nu_{br}`, then
.. math::
:label: gipsylinearwithoutf0
Z(N) \approx \frac{Z_r\nu_{br} - {\bf n} c \delta_{\nu_b}}{\nu_{br}}
= Z_r + {\bf n} \bigl(-c\frac{\delta_{\nu_b}}{\nu_{br}}\bigr)
Next script implements these formulas and show the deviations. The first three columns show the
correct result.
.. literalinclude:: EXAMPLES/gipsy2vels.py
:language: python
Output::
Topocentric correction (km/s): 9.57140206387
Barycentric frequency and increment (Hz): 1418966870.14 -9765.3132202
pix WCSLIB GIP+bary GIP+bary-f0
61.9940 299.869536 299.869536 299.869536
62.9940 301.934754 301.934754 301.934754
63.9940 304.000000 304.000000 304.000000
64.9940 306.065274 306.065274 306.065274
65.9940 308.130577 308.130577 308.130577
GIP+topo Linear+f0 Linear-f0
299.869141 299.869479 299.873664
301.934556 301.934740 301.936832
304.000000 304.000000 304.000000
306.065472 306.065260 306.063168
308.130973 308.130521 308.126336
The columns in the output are:
1. *pix*: The (non integer) pixel value at which a velocity is calculated.
2. *WCSLIB*: The optical velocity (km/s) as calculated by WCSLIB. The extension in
CTYPE is recognized and the frequencies are replaced by their
barycentric counterparts according to the recipe in :ref:`spectral_gipsy`.
3. *GIP+bary*: The optical velocity (km/s) calculated with GIPSY formula in eq. :eq:`gipsynonlinear`
using barycentric reference frequency and barycentric frequency increment.
4. *GIP+bary-f0*: The optical velocity (km/s) calculated with GIPSY formula
without the rest frequency as in eq. :eq:`gipsynonlinearwithoutf0`
using barycentric reference frequency and barycentric frequency increment.
5. *GIP+topo*: The optical velocity (km/s) calculated with GIPSY formula in eq. :eq:`gipsynonlinear`
using topocentric/geocentric reference frequency and frequency increment.
6. *Linear+f0*: The optical velocity (km/s) calculated with GIPSY formula in eq. :eq:`gipsylinearwithf0`
using a rest frequency.
7. *Linear-f0*: The optical velocity (km/s) calculated with GIPSY formula in eq. :eq:`gipsylinearwithoutf0`
without a rest frequency.
If you do some experiments with the values in this script,
you will observe that the GIPSY formula with topocentric instead of the barycentric/lsrk
values is not
a bad approximation although it is sensitive to the channel number (*p*).
The linear approximations are worse and should be avoided if high precision is required.
What remains is the question how good GIPSY's approximation is.
With :eq:`gipsynonlinear` we write:
.. math::
:label: gipsyapprox0
Z_{\nu_{b}}(N) - Z_{\nu_{t}}(N) = c \nu_0 \Bigl(
\frac{1} {\nu_{br}+{\bf n}\delta_{\nu_{b}}} - \frac{1}{\nu_{br}}-
\bigl(\frac{1} {\nu_{tr}+{\bf n}\delta_{\nu_{t}}} - \frac{1}{\nu_{tr}} \bigr)\Bigr)
With the parameters:
* :math:`Z_{\nu_{t}}(N)` the optical velocity at pixel *N* using topocentric values
* :math:`\nu_{tr}` the topocentric frequency at the reference pixel
* :math:`\delta_{\nu_{t}}` the topocentric frequency increment
Rewrite this in:
.. math::
:label: gipsyapprox1
Z_{\nu_{b}}(N) - Z_{\nu_{t}}(N) = -{\bf n} c \nu_0 \bigl(
\frac{\delta_{\nu_{b}}} {\nu_{br}(\nu_{br}+{\bf n}\delta_{\nu_{b}})} -
\frac{\delta_{\nu_{t}}} {\nu_{tr}(\nu_{tr}+{\bf n}\delta_{\nu_{t}})} \bigr)
Note that :math:`{\bf n}\delta \nu_b << \nu_{br}` and :math:`{\bf n}\delta \nu_t << \nu_{tr}`. Then write
the difference in increment as function of *N* as:
.. math::
:label: gipsyapprox2
Z_{\nu_{b}}(N) - Z_{\nu_{t}}(N) \approx -{\bf n} c \nu_0 \bigl( \frac{\delta_{\nu_{b}}} {\nu^2_{br}} -
\frac{\delta_{\nu_{t}}} {\nu^2_{tr}} \bigr)
This expression explains the different values in the output of our previous script and it shows
that the differences depend on **n**.
Use :eq:`eq120` to write:
.. math::
:label: gipsyapprox2a
\nu_{tr} = \nu_{br} \sqrt{\frac{c-v_{tc}}{c+v_{tc}}}
and from :eq:`eq140`
.. math::
:label: gipsyapprox2b
\delta_{\nu_b} = \delta_{\nu_t} \sqrt{\frac{c-v_{tc}}{c+v_{tc}}}
Define :math:`q = \sqrt{\frac{c-v_{tc}}{c+v_{tc}}}` then :math:`\nu_{br} = q / \nu_{tr}`
and :math:`\delta_{\nu_b} = q * \delta_{\nu_t}`
Insert this in :eq:`gipsyapprox2` to obtain:
.. math::
:label: gipsyapprox2c
Z_{\nu_{b}}(N) - Z_{\nu_{t}}(N) \approx -{\bf n} c \nu_0 \frac{\delta_{\nu_t}}{\nu^2_{tr}} (q^3-1)
The topocentric correction :math:`v_{tc}` has a range between -30 Km/s and 30 Km/s.
For :math:`v_{tc} = 30000` *m/s* this corresponds to a maximum of q:
:math:`q = 0.99989993577786473`. With this maximum for *q* we find for
:eq:`gipsyapprox2c` approximately 0.62 m/s
Note that the difference is a function of **n**, so
after 64 channels the deviation is almost 40 m/s.
In our example, the channel separation is approximately 2 km/s and the deviations
are therefore small (2%).
For the example at the start of this chapter, the reference velocity was 9120 km/s.
The channel separation (*CDELT3Z*) is approximately 20 km/s. For the listed topocentric
frequency and the calculated barycentric frequency we find with :eq:`gipsyapprox2c` an error
of approximately 6.6 m/s. After 64 channels the deviation is approximately 420 m/s
(2%).
With :eq:`eq866` we get an relative error:
.. math::
:label: gipsyapprox2d
\frac{ Z_{\nu_{b}}(N) - Z_{\nu_{t}}(N)}{dZ} = {\bf n} (q^3-1) \frac{\delta_{\nu_t}}{\nu^2_{tr}} \frac{\nu^2_{br}}{\delta_{\nu_t}} \approx {\bf n}(q^3-1)
With the maximum value of *q* we find a maximum percentage of 0.03% for 1 channel.
After 64 channels the
deviation is almost 2%. After 512 channels it is more than 15%.
**Conclusions**
* The GIPSY formulas assume constant frequency increments in the system of the reference
system. When these are topocentric, there are small deviations from the result with WCSLIB
which assume the frequencies in the same reference system as the given velocity.
* The formula that GIPSY routines use to calculate optical velocities
is an approximation. The deviations are small but
depend on the pixel i.e. :math:`(N-N_{\nu})`. This approximation is not necessary because
when the optical velocity in the barycenter is given, then
one can calculate the barycentric reference frequency (see eq. :eq:`eq10`)
and use that frequency in the GIPSY formula to get the exact result.
* The deviation is more sensitive to the topocentric correction
(velocity between observatory on earth and barycenter/lsrk) than the
reference frequency and the frequency increment. Also there is a maximum
value for the topocentric velocity which results in a maximum deviation of
0.03% for one channel.
For the data in the previous script, we used the code below (which should be added to the previous script)
to calculate the percentages:
.. code-block:: python
q = sqrt((c-Vtopo)/(c+Vtopo))
delta = -c*f0*df/fr/fr * (q*q*q-1)
d = (p-crpix) * delta
# Now change the topocentric correction to its maximum.
Vtopo = 30000.0
qmax = sqrt((c-Vtopo)/(c+Vtopo))
deltamax = -c*f0*df/fr/fr * (qmax*qmax*qmax-1)
dmax = (p-crpix) * deltamax
perc = abs(100*deltamax/dZ)
print "dZ, deltamax:", dZ, deltamax
print "Percentage deviation for 1 channel: ", perc
print "Approximate percentage: ", abs(100 * (qmax*qmax*qmax-1))
print "Percentage deviation for 64 channel: ", 64*perc
print "Approximate percentage: ", abs(100 * 64*(qmax*qmax*qmax-1))
print "Percentage deviation for 64 channel: ", 512*perc
print "Approximate percentage: ", abs(100 * 512*(qmax*qmax*qmax-1))
print "\nThe approximate difference and the real difference"
print "between topocentric nd barycentric increments"
for pixel, d1,d2,d3 in zip(pixrange, d, Z2-Z4, dmax):
print "%10.4f %14f %14f %14f" % (pixel, d1/1000, d2/1000, d3/1000)
Output::
dZ, deltamax: -21236.6115174 6.57007047211
Percentage deviation for 1 channel: 0.0309374707295
Approximate percentage: 0.0300162628862
Percentage deviation for 64 channel: 1.97999812669
Approximate percentage: 1.92104082472
Percentage deviation for 64 channel: 15.8399850135
Approximate percentage: 15.3683265977
The approximate difference and the real difference
between topocentric and barycentric increments and
the maximum deviation as function of the pixel:
61.9940 -0.011436 -0.011438 -0.013140
62.9940 -0.005718 -0.005719 -0.006570
63.9940 0.000000 0.000000 0.000000
64.9940 0.005718 0.005717 0.006570
65.9940 0.011436 0.011433 0.013140
61.9940 -0.011436 -0.011438 -0.013140
62.9940 -0.005718 -0.005719 -0.006570
63.9940 0.000000 0.000000 0.000000
64.9940 0.005718 0.005717 0.006570
65.9940 0.011436 0.011433 0.013140
Radio
^^^^^^^
Given a frequency, a radio velocity is calculated with the formula:
.. math::
:label: gipradio10
V = -c (\frac{\nu'-\nu_0}{\nu_0})
Assume for channel :math:`N`:
.. math::
:label: gipradio20
\nu(N) = \nu_{br} + (N-N_{ref}) \delta_{\nu_b} = \nu_{br} + {\bf n} \delta_{\nu_b}
For :math:`(N-N_{ref})` we wrote :math:`\bf n`.
The frequencies are related to the barycentric (or lrsk) reference system.
:math:`N_{ref}` is the reference pixel (*CRPIX*) given in a FITS header,
:math:`\nu_{br}` is the reference frequency in this barycentric system and
:math:`\delta_{\nu_b}` is the barycentric frequency increment.
Inserting :eq:`gipradio10` into :eq:`gipradio20` gives:
.. math::
:label: gipradio30
V_b(N) = -c\bigl( \frac{\nu_{br}+{\bf n} \delta_{\nu_b} - \nu_0}{\nu_0} \bigr) = V_r + {\bf n}\frac{-c\delta_{\nu_b}}{\nu_0}
with:
* :math:`V_b(N)` is the barycentric radio velocity for pixel *N* using barycentric
frequency increments
* :math:`\nu_{br}` is the barycentric reference frequency
* :math:`\delta_{\nu_b}` is the increment in barycentric frequency
This increment in radio velocity was also derived in eq. :eq:`eq260`.
The increment in radio velocity is a linear function of the increment in frequency.
The frequencies in the FITS and GIPSY headers for pre July, 2006 WSRT/Nmap FITS files
are the topocentric frequencies.
We show the difference between the velocities derived from the barycentric/lsrk values and
the velocities derived from the topocentric values.
.. literalinclude:: EXAMPLES/gipsyradiovels.py
:language: python
Output::
Topocentric correction (km/s): 9.26313531147
Barycentric frequency and increment (Hz): 1418965411.07 -9765.32326156
pix WCSLIB GIP+bary GIP+bary-f0 GIP+topo
61.9940 299.877839 299.877839 299.877839 299.877712
62.9940 301.938920 301.938920 301.938920 301.938856
63.9940 304.000000 304.000000 304.000000 304.000000
64.9940 306.061080 306.061080 306.061080 306.061144
65.9940 308.122161 308.122161 308.122161 308.122288
The second, third and fourth column represent :math:`V_b` and the last column is :math:`V_t`.
The difference between the exact and approximate velocities as function of **n** is given by:
.. math::
:label: gipradio40
V_t(N) - V_b(N) = -{\bf n} \frac{c}{\nu_0}(\delta_{\nu_t}-\delta_{\nu_b})
With the parameters:
* :math:`V_t(N)` the barycentric radio velocity at pixel *N* using topocentric
frequency increments
* :math:`\delta_{\nu_{t}}` the topocentric frequency increment
The topocentric correction :math:`v_{tc}` has a range between -30 Km/s and 30 Km/s.
Rewrite :eq:`eq140` into:
.. math::
:label: gipradio50
\frac{\delta_{\nu_b}} {\delta_{\nu_t}} = \sqrt{\frac{c-v_{tc}}{c+v_{tc}}}
For :math:`v_{tc} = 30000` *m/s* this corresponds to a maximum
:math:`q = {\delta_{\nu_b} / \delta_{\nu_t}} = 0.99989993577786473`
which is equivalent to:
.. math::
:label: gipradio60
\frac{c}{\nu_0}(1-q)\delta_{\nu_t} \approx 0.2\ m/s
Note that the difference is a function of **n**, so
after 64 channels the deviation is more than 12 m/s.
In our example, the channel separation is approximately 2 km/s and the deviations
are therefore small.
Header items in a (legacy) WSRT FITS file
-----------------------------------------
Program *nmap* (part of NEWSTAR which is a package developed
to process WSRT and ATCA data)
is/was used to create FITS files with WSRT line data.
We investigated the meaning or interpretation of the various FITS header items.
The program generates it own descriptors related to velocities and frequencies.
For example:
* VEL: Velocity (m/s)
* VELC: Velocity code
* 0=continuum,
* 1=heliocentric radio
* 2=LSR radio
* 3=heliocentric optical
* 4=LSR optical
* VELR: Velocity at reference frequency (FRQC)
* INST: Instrument code (0=WSRT, 1=ATCA)
* FRQ0: Rest frequency for line (MHz)
* FRQV: Real frequency for line (MHz)
* FRQC: Centre frequency for line (MHz)
One of functions in *nmap* is called *nmawfh.for*.
It writes a FITS header using the values in the *nmap* descriptors.
The value of *CRVAL3* is set to *FRQV* if the velocity code is one of
combinations of optical and radio velocity with heliocentric or local standard of rest
reference systems (i.e. RHEL, RLSR, OHEL, OLSR).
The value of *CRPIX3* is equal to *FRQV* -lowest frequency divided by the
channel separation. 'lowest frequency' is the frequency of the input channel with
the lowest frequency.
* The value for FITS keyword VEL= is equal to *nmap* descriptor VEL, the centre velocity in m/s
* The value for FITS keyword VELR= is equal to *nmap* descriptor VELR, the Reference velocity
* The value for FITS keyword FREQR= is equal to *nmap* descriptor FRQC, the Reference frequency (Hz)
* The value for FITS keyword FREQ0= is equal to *nmap* descriptor FRQ0, the Rest frequency (Hz)
::
VEL !CENTRE VELOCITY (M/S)
VELCODE !VELOCITY CODE
VELR !REFERENCE VELOCITY (M/S)
FREQR !REFERENCE FREQUENCY (HERTZ)
FREQ0 !REST FREQUENCY (HERTZ)
WCSLIB in a GIPSY task
+++++++++++++++++++++++
GIPSY (`Groningen Image Processing SYstem `_ )
is one of the oldest image processing and data analysis systems. Python can be used
to create GIPSY tasks. The Kapteyn Package is integrated in GIPSY. Here we
give a small example how to use both.
Assuming you have a data set with three axes and the last axis is the spectral axis,
the next script is a very small GIPSY program that asks the user for the name of this set
and then calculates the optical velocities for a number of pixels in the neighborhood of
the reference pixel (CRPIX3).
GIPSY data have a descriptor which contains FITS header
items (e.g. CRVAL1=) and GIPSY specific keywords but not only attached to the
set but also to subsets (slices) of the data. Not only planes or lines can have
their own header but even pixels can. The script below reads its information from
top level (which hosts the global description of the data cube itself):
.. code-block:: python
#!/usr/bin/env python
from gipsy import *
from kapteyn import wcs
init()
while True:
try:
set = Set(usertext('INSET=', 'Input set'))
break
except:
reject('INSET=', 'Cannot open set')
proj = wcs.Projection(set).sub((3,))
s = "Ref. freq at that pixel: %f Hz" % (set['CRVAL3'],)
anyout(s)
s = "Velocity: %f m/s" % (set['DRVAL3'],)
anyout(s)
crpix = set['CRPIX3']
proj2 = proj.spectra('VOPT-F2W')
for i in range(-2,+3):
world = proj2.toworld((crpix+i,))[0]/1000.0 # to world coordinates
anyout(str(world)+' km/s')
finis()
This little GIPSY task simulates the functionality of GIPSY task *COORDS*
which lists world coordinates for data slices. The two most important differences between
this task and *COORDS* are:
* With WCSLIB it is simple to change the
output velocity to radio or apparent radial by changing the spectral translation.
* The Python interface to WCSLIB prepares
the GIPSY header information to give correct barycentric or lsrk velocities (i.e. it
also converts the frequency increment to the barycentric or lsrk system).
Read more about GIPSY tasks written in Python in:
`Python recipes for GIPSY `_
References
----------
.. [Aipsmemo] `AIPS memo 27 `_
Non-Linear Coordinate Systems in AIPS (Eric W. Greisen, NRAO)