Refractive index of air and wave length conversion

Refractive index of air

The Python module ref_index has functions for calculating the refractive index of air. It also has functions for converting between wave length of light in vacuum and that in air.

The equations used in this module are taken from the documentation for the NIST online refractive index of air calculator.

Examples

In the following examples, wave lengths are in nano-meters, pressure in Pascal, temperature in degree Celsius, and relative humidity is a number between 0 and 100.

Converting wave length of light from vacuum to that in air

>>> ref_index.vac2air(633.0)
632.82500476826874

>>> ref_index.vac2air(np.array([633.0, 550.0, 400.0]))
array([ 632.82500477,  549.84723175,  399.88692724])

Default calculation is carried out at a temperature of 15 degree Celsius, pressure of 101325Pa, relative humidity of 0, and a CO2 concentration of 450micro-mole/mole. These can be changed by passing appropriate parameters to the function.

>>> ref_index.vac2air(wave=633.0, t=20, p=101325, rh=50)
632.8282676547044

Converting wave length of light from air to that in vacuum

The equations need wave length in vacuum, but I use the wave lengths in air as an approximation. The full conversion from vacuum to air and back agrees to ~1e-5nm, as shown below.

>>> x = ref_index.vac2air(np.array([633.0, 550.0, 400.0]))

>>> ref_index.air2vac(x)
array([ 633.0000014 ,  550.00000164,  400.00000243])

Refractive index

Refractive index can be calculated using two different equations: one due to Edlen and another due to Ciddor. The wave length conversion functions, vac2air() and air2vac() both use the Ciddor formula.

>>> ref_index.ciddor(wave=633.0, t=20, p=101325, rh=20)
1.0002716285340578
>>> ref_index.edlen(wave=633.0, t=20, p=101325, rh=20)

1.0002716291691649
>>> ref_index.edlen(wave=633.0, t=20, p=101325, rh=80)
1.0002711197635226
>>> ref_index.ciddor(wave=633.0, t=20, p=101325, rh=80)

1.0002711183472626
>>> ref_index.edlen(wave=633.0, t=60, p=101325, rh=80)
1.0002339748542823
>>> ref_index.ciddor(wave=633.0, t=60, p=101325, rh=80)

1.0002340241754055

Installation

The module can be installed using pip and easy_install.

$ pip install ref_index

or,

$ easy_install ref_index

The code repository is at Github: http://github.com/phn/ref_index.

Details

The following comments are based on the discussions presented in the NIST documentation. It is intended as a brief overview. See http://emtoolbox.nist.gov/Wavelength/Documentation.asp, for detailed discussions.

Refractive index of air can be calculated using two different algorithms: one due to Edlen (updated by Birch and Down), and one due to Ciddor. The latter has been adopted by the International Association of Geodesy (IAG) as the reference equation for calculating refractive index of air. Functions for calculating refractive index using either of these are defined in this module.

The vacuum to air and air to vacuum wave length conversion functions in this module use the Ciddor equation, in the form presented in the NIST documentation.

Uncertainties in refractive index, and hence in wave length conversions, due to uncertainties in measured values of temperature, pressure, and humidity exceeds that due to the intrinsic uncertainty in the equations used.

An uncertainty of 1e-6 in refractive index can result from a combination of:

  • an error of 1 degree Celsius (1.8 degree F) in air temperature
  • an error of 0.4kPa (3mm of Hg) in air pressure
  • an error of 50% in relative humidity at sufficiently high air
    temperatures (near 35 degree Celsius)

Valid range for input parameters for the refractive index calculations
are presented below. The online calculator issues a warning if input
parameters are outside a smaller interval within the maximum
range. Functions in this module do not raise a warning by default. But
they accept a keyword warn, which when set to True will result
in warnings, when the input parameters are outside the accepted range.

  • Wavelength [300nm – 1700nm]

    Warning is issued if value is outside [350nm – 1600nm].

  • Pressure [10kPa – 140kPa]

    Warning is issued if value is outside [60kPa – 120kPa].

  • Temperature [-40 degree C – 100 degree C].

    Warning is issued if value is outside [0 degree C – 40 degree C].

  • Humidity [0 – 100]

    Can be given as relative humidity, dew point, frost point or partial pressure of water vapour. A warning is given if the mole fraction of water vapour exceeds 20% or, equivalently, relative humidity exceeds 85%. A warning is issued if relative humidity is less than 1%.

  • CO2 concentration [0micro-mole/mole – 2000micro-mole/mole]

    The common value to use is 450. Outdoor values are rarely below 300 and indoor can be as high as 600. A difference of 150 will lead to a difference of only ~ 2e-8 in index of refraction.

    A warning is issued if a value other than 450 is used.

Most of these options are probably not relevant for astronomy applications. The defaults should be enough. See the next section for comparison with IDLASTRO IDL code.

Comparison with IDLASTRO code

In astronomy, the convention is to use the refraction correction for wave length greater than 200nm, even though the equations are not strictly valid at wave lengths shorter than 300nm. For example, the popular IDLASTRO IDL code vactoair.pro and airtovac.pro will accept any wave length greater than 2000 angstrom.

To accommodate this type of usage, instead of limiting the possible input wave lengths, functions in this module will accept any wave length value. It is up to the user to decide if a particular wave length is to be used as an input to the equations.

Comparison with the IDLASTRO vactoair.pro and airtovac.pro algorithms show that the equivalent functions in this module, vac2air and air2vac, give results that agree to within 1e-4nm, over a range of wavelengths from 200nm to 1700nm. This uncertainty translates to a velocity difference of 150m/s to 17m/s, over the wave length range 1700nm to 200nm.

The IDLASTRO code uses a fixed value of temperature and humidity which is not documented in the IDL code. The above comparison was carried out at a temperature of 15 degree C and a relative humidity of 0.

The IDL code used for testing was downloaded on 2011/10/07. The revision history indicates that the IDL code in vactoair.pro and airtovac.pro were last modified in March 2011.

See the functions with names starting with the string _test for details on how the tests were performed.

Posted in Astronomy, Python | Tagged , | Leave a comment

Line id plot

[Update 2012/04/22: Following suggestions from Jane Rigby (see comments below), I have created a couple of functions for changing properties of labels. These exists as separate code in the following Github Gist: https://gist.github.com/2326396 .]

Manually placing labels in a plot, especially if it is crowded, is cumbersome.

The Python module lineid_plot, has functions for automatically placing labels in a plot, in such a way that the labels do not overlap with each other. This is very useful, for example, in creating plots that have labels identifying features in a spectrum.

This Python module is an adaptation of the IDL module lineid_plot.pro, in the IDLAstro library.

The code that adjusts the label positions is a direct translation of the appropriate section of the IDL code. Other features are implemented using the facilities provided by Matplotlib.

The module can be installed using pip or easy_install:

$ easy_install lineid_plot

or

$ pip install lineid_plot

The module can also be downloaded from http://pypi.python.org/pypi/lineid_plot/. Documentation is available at http://packages.python.org/lineid_plot/.

The source code repository for the module is available on Github: http://github.com/phn/lineid_plot.

An example is show below. The plot can be customized in several ways, as shown in the documentation.

import numpy as np

from matplotlib import pyplot as plt
import lineid_plot

wave = 1240 + np.arange(300) * 0.1

flux = np.random.normal(size=300)

line_wave = [1242.80, 1260.42, 1264.74, 1265.00, 1265.2, 1265.3, 1265.35]
line_label1 = ['N V', 'Si II', 'Si II', 'Si II', 'Si II', 'Si II', 'Si II']

lineid_plot.plot_line_ids(wave, flux, line_wave, line_label1)

plt.show()
LineID plot example.

Example showing automatic placement of labels in a plot.

Posted in Astronomy, Python | Tagged , | 13 Comments

Julian dates for proleptic Gregorian and Julian calendars

Julian dates for proleptic Gregorian and Julian calendars

The Python module jdcal provides a few functions for calculating Julian dates for Gregorian and Julian calendar dates.  Julian dates are stored as two floating point numbers for maximum precision.

Installation

The module can be installed using easy_install or pip.

$ pip install jdcal

or,

$ easy_install jdcal

The module can be downloaded from the pypi page for the module.

The source code repository is at http://github.com/phn/jdcal .

The module is available under the BSD license (http://www.opensource.org/licenses/bsd-license.php).

Functions

The function gcal2jd accepts the year, month and day of a Gregorian calendar date and returns corresponding Julian date. The function jcal2jd accepts the year, month and day of a Julian calendar date and returns the corresponding Julian date. Functions jd2gcal and jd2jcal take Julian dates and returns Gregorian and Julian calendar dates, respectively.

Different regions of the world switched to Gregorian calendar from Julian calendar on different dates. Having separate functions for Julian and Gregorian calendars allow maximum flexibility in choosing the relevant calendar.

All the above functions are “proleptic”. This means that they work for dates on which the concerned calendar is not valid. For example, Gregorian calendar was not used prior to around October 1582.

Julian dates are stored in two floating point numbers (double). Julian dates, and Modified Julian dates, are large numbers. If only one number is used, then the precision of the time stored is limited. Using two numbers, time can be split in a manner that will allow maximum precision. For example, the first number could be the Julian date for the beginning of a day and the second number could be the fractional day. Calculations that need the latter part can now work with maximum precision.

The function is_leap can be used to test whether a given Gregorian calendar year is a leap year or not.

Zero point of Modified Julian Date (MJD) and the MJD of 2000/1/1 12:00:00 are defined as MJD_0 and MJD_J2000 respectively.

This module is based on the TPM C library, by Jeffery W. Percival.

Examples

Gregorian calendar to Julian dates

Function gcal2jd returns the Julian date for midnight of the given Calendar date.

>>> gcal2jd(2000,1,1)
(2400000.5, 51544.0)

>>> 2400000.5 + 51544.0 + 0.5
2451545.0

Negative months and days are valid. For example, 2000/-2/-4 => 1999/+12-2/-4 => 1999/10/-4 => 1999/9/30-4 => 1999/9/26.

>>> gcal2jd(2000, -2, -4)
(2400000.5, 51447.0)
>>> gcal2jd(1999, 9, 26)
(2400000.5, 51447.0)

>>> gcal2jd(2000, 2, -1)
(2400000.5, 51573.0)
>>> gcal2jd(2000, 1, 30)
(2400000.5, 51573.0)

>>> gcal2jd(2000, 3, -1)
(2400000.5, 51602.0)
>>> gcal2jd(2000, 2, 28)
(2400000.5, 51602.0)

Month 0 becomes previous month.

>>> gcal2jd(2000, 0, 1)
(2400000.5, 51513.0)
>>> gcal2jd(1999, 12, 1)
(2400000.5, 51513.0)

Day number 0 becomes last day of previous month.

>>> gcal2jd(2000, 3, 0)
(2400000.5, 51603.0)
>>> gcal2jd(2000, 2, 29)
(2400000.5, 51603.0)

If day is greater than the number of days in month, then it gets carried over to the next month.

>>> gcal2jd(2000,2,30)
(2400000.5, 51604.0)

>>> gcal2jd(2000,3,1)
(2400000.5, 51604.0)

>>> gcal2jd(2001,2,30)
(2400000.5, 51970.0)

>>> gcal2jd(2001,3,2)
(2400000.5, 51970.0)

Function jd2gcal converts Julian dates to Gregorian calendar dates. The values returned are year, month, day and fraction of the day.

>>> jd2gcal(*gcal2jd(2000,1,1))
(2000, 1, 1, 0.0)

>>> jd2gcal(*gcal2jd(1950,1,1))
(1950, 1, 1, 0.0)

>>> gcal2jd(2000,1,1)
(2400000.5, 51544.0)
>>> jd2gcal(2400000.5, 51544.0)
(2000, 1, 1, 0.0)

>>> jd2gcal(2400000.5, 51544.5)
(2000, 1, 1, 0.5)
>>> jd2gcal(2400000.5, 51544.245)
(2000, 1, 1, 0.24500000000261934)

>>> jd2gcal(2400000.5, 51544.1)
(2000, 1, 1, 0.099999999998544808)
>>> jd2gcal(2400000.5, 51544.75)
(2000, 1, 1, 0.75)

Out of range months and days are carried over to the next/previous year or next/previous month.

>>> jd2gcal(*gcal2jd(1999,10,12))
(1999, 10, 12, 0.0)

>>> jd2gcal(*gcal2jd(2000,2,30))
(2000, 3, 1, 0.0)

>>> jd2gcal(*gcal2jd(-1999,10,12))
(-1999, 10, 12, 0.0)

>>> jd2gcal(*gcal2jd(2000, -2, -4))
(1999, 9, 26, 0.0)

Julian calendar to Julian dates

Julian calendar dates can be converted to Julian dates using jcal2jd. Unlike gcal2jd, negative months and days can result in incorrect Julian dates.

>>> jcal2jd(2000, 1, 1)
 (2400000.5, 51557.0)

The function jd2jcal converts Julian dates to Julian calendar dates.

>>> jd2jcal(*jcal2jd(2000, 1, 1))
(2000, 1, 1, 0.0)

>>> jd2jcal(*jcal2jd(-4000, 10, 11))
(-4000, 10, 11, 0.0)

>>> jcal2jd(2000, 1, 1)
(2400000.5, 51557.0)
>>> jd2jcal(2400000.5, 51557.0)
(2000, 1, 1, 0.0)

>>> jd2jcal(2400000.5, 51557.5)
(2000, 1, 1, 0.5)
>>> jd2jcal(2400000.5, 51557.245)
(2000, 1, 1, 0.24500000000261934)

>>> jd2jcal(2400000.5, 51557.1)
(2000, 1, 1, 0.099999999998544808)
>>> jd2jcal(2400000.5, 51557.75)
(2000, 1, 1, 0.75)

Leap years

Use the function is_leap to determine whether a Gregorian calendar year is a leap year or not.

>>> is_leap(2000)
True
>>> is_leap(2100)

False

Constants

The zero point of Modified Julian Date (MJD) is stored in MJD_0. The MJD for 2000/1/1 12:00:00 is stored in MJD_J2000.

>>> print MJD_0

2400000.5
>>> print MJD_JD2000
51544.5
Posted in Astronomy, Python | Tagged , , , | 6 Comments