Module Celestial

This document describes functions from the Python module celestial (celestial.py) which provides a programmer with a basic set of routines to transform a world coordinate in a given sky system into a world coordinate of another system assuming zero proper motion, parallax, and recessional velocity.

The most important function builds a matrix for conversions of positions between sky systems, celestial reference systems and epochs of the equinox. This function is called skymatrix() and it can be used in the following contexts:

  • Implicit, in module wcs, using the Transformation class as in:

    world_eq = (192.25, 27.4)   # FK4 coordinates of galactic pole
    tran = wcs.Transformation("equatorial fk4_no_e B1950.0", "galactic")
    print(tran(world_eq))
    
  • As stand alone utility in scripts or in an interactive Python session. Usually one uses function sky2sky() to transform longitudes and latitudes:

    M = celestial.sky2sky( (celestial.eq, celestial.fk5), celestial.gal,
                           (0,0,1.0), (10,20,20) )
    
  • Hidden in the topixel() and toworld() methods in module wcs. There the sky system is read from a (FITS) header and the sky system for which we want the transformed coordinates is set with attribute skyout of the projection object.

See also

Tutorial material:

Sky definitions

A sky definition can consist of a sky system, a reference system, an equinox and an epoch of observation. It is either a string or it is a tuple with one or more elements. It can also be a single element. The elements in a tuple representing a sky- or reference system are symbols from the table below. For a string, the parts of the string representing a sky- or reference system are minimal matched against the strings in the table below. The match is case insensitive.

Sky systems

Symbol

String

Description

eq, equatorial

EQUATORIAL

Equatorial coordinates (\(\alpha\) , \(\delta\) ) See also next table with reference systems

ecl, ecliptic

ECLIPTIC

Ecliptic coordinates (\(\lambda\) , \(\beta\) ) referred to the ecliptic and mean equinox

gal, galactic

GALACTIC

Galactic coordinates (lII, bII)

sgal, supergalactic

SUPERGALACTIC

De Vaucouleurs Supergalactic coordinates (sgl, sgb)

Reference systems

Symbol

String

Description

fk4

FK4

Mean place pre-IAU 1976 system. FK4 is the old barycentric (i.e. w.r.t. the common center of mass) equatorial coordinate system, which should be qualified by an Equinox value. For accurate work FK4 coordinate systems should also be qualified by an Epoch value. This is the epoch of observation.

fk4_no_e

FK4_NO_E, FK4-NO-E

The old FK4 (barycentric) equatorial system but without the E-terms of aberration. This coordinate system should also be qualified by both an Equinox and an Epoch value.

fk5

FK5

Mean place post IAU 1976 system. Also a barycentric equatorial coordinate system. This should be qualified by an Equinox value (only).

icrs

ICRS

The International Celestial Reference System, for optical data realized through the Hipparcos catalog. By definition, ICRS is not an equatorial system, but it is very close to the FK5 (J2000) system. No Equinox value is required.

j2000, dynj2000

DYNJ2000

This is an equatorial coordinate system based on the mean dynamical equator and equinox at epoch J2000. The dynamical equator and equinox differ slightly compared to the equator and equinox of FK5 at J2000 and the ICRS system. This system need not be qualified by an Equinox value

Note

Reference systems are stored in FITS headers under keyword RADESYS=.

Note

Standard in FITS: RADESYS defaults to IRCS unless EQUINOX is given alone, in which case it defaults to FK4 prior to 1984 and FK5 after 1984.

EQUINOX defaults to 2000 unless RADESYS is FK4, in which case it defaults to 1950.

Note

In routines dealing with sky definitions tne names are minimal matched against a list with full names.

Epochs for the equinox and epoch of observation

An epoch can be set in various ways. The options are distinguished by a prefix. Only the ‘B’ and ‘J’ epochs can be negative.

Prefix

Epoch

B

Besselian epoch. Example: 'B 1950', 'b1950', 'B1983.5', '-B1100'

J

Julian epoch. Example: 'j2000.7', 'J 2000', '-j100.0'

JD

Julian date. This number of days (with decimals) that have elapsed since the initial epoch defined as noon Universal Time (UT) Monday, January 1, 4713 BC in the proleptic Julian calendar Example: 'JD2450123.7'

MJD

The Modified Julian Day (MJD) is the number of days that have elapsed since midnight at the beginning of Wednesday November 17, 1858. In terms of the Julian day: MJD = JD - 2400000.5 Example: 'mJD 24034', 'MJD50123.2'

RJD

The Reduced Julian Day (RJD): Julian date counted from nearly the same day as the MJD, but lacks the additional offset of 12 hours that MJD has. It therefore starts from the previous noon UT or TT, on Tuesday November 16, 1858. It is defined as: RJD = JD - 2400000 Example: 'rJD50123.2', 'Rjd 23433'

F

Various FITS formats:

  • DD/MM/YY Old FITS format. Example: 'F29/11/57'

  • YYYY-MM-DD FITS format. Example: 'F2000-01-01'

  • YYYY-MM-DDTHH:MM:SS FITS format with date and time. Example: 'F2002-04-04T09:42:42.1'

Epoch of observation.

Reference system FK4 is not an inertial system. It is slowly rotating and positions are further away from the true mean places if the date of observation is greater than B1950. FK5 is an inertial system. If we convert coordinates from FK4 to FK5, the accuracy of the FK5 position can be improved if we know the date of the observation. So in all transformations where a conversion between FK4 and FK5 is involved, an epoch of observation can be part of the sky definition. Note that this also involves a conversion between galactic coordinates and equatorial, FK5 coordinates because that conversion is done in steps and one step involves FK4.

To be able to distinguish an equinox from an epoch of observation, an epoch of observation is followed by an underscore character and some arbitrary characters to indicate that it is a special epoch (e.q. “B1960_OBS”). Only the underscore is obligatory.

Note

If a sky definition is entered as a string, there cannot be a space between the prefix and the epoch, because a space is a separator for the parser in celestial.skyparser().

Note

An epoch of observation is either the second epoch in your input or or the epoch string has a suffix ‘_’ which may be followed by arbitrary characters (e.g. “B1963.5_OBS”).

Input Examples

Input string

Description

Remarks

“eq”

Equatorial, ICRS

ICRS because no reference system and no equinox is given.

“Eclip”

Ecliptic, ICRS

Ecliptic coordinates

“ecl fk5”

Ecliptic, FK5

Ecliptic coordinates with a non default reference system

“GALACtic”

Galactic II

Minimal match is case insensitive

“s”

Supergalactic

Shortest string to identify system.

“fk4”

Equatorial, FK4

Only a reference system is entered. Sky system is assumed to be equatorial

“B1960”

Equatorial, FK4

Only an equinox is given. This is a date before 1984 so FK4 is assumed. Therefore the sky system is equatorial

“EQ, fk4_no_e, B1960”

Equatorial, FK4 no e-terms

Sky system, reference system, and an equinox

“EQ, fk4-no-e, B1960”

Equatorial, FK4 no e-terms

Same as above but underscores replaced by hyphens.

“fk4,J1983.5_OBS”

Equatorial, FK4 + epobs

FK4 with an epoch of observation. Note that only the underscore is important.

“J1983.5_OBS”

Equatorial, FK4 + epobs

Only a date of observation. Then reference system FK4 is assumed.

“EQ,fk4,B1960, B1983.5_O”

Equatorial, FK4 + epobs

A complete description of an equatorial system.

“B1983.5_O fk4 B1960,eq”

Equatorial, FK4 + epobs

The same as above, showing that the order of the elements are unimportant.

Code examples

To show that one can use both the tuple and the string representation of a system, we use both for the same system and compare a transformed position. The result should be 0 for both coordinates.

>>> world_eq = numpy.array([192.25, 27.4])     # FK4 coordinates of galactic pole
>>> tran1 = wcs.Transformation("equatorial fk4_no_e B1950.0", "galactic")
>>> tran2 = wcs.Transformation((wcs.equatorial, wcs.fk4_no_e, 'B1950.0'), wcs.galactic)
>>> print(tran1(world_eq)-tran2(world_eq))
[ 0.  0.]

Module level data

skyrefsystems

An object from class skyrefset which is a container with a list with systems and two dictionaries with systems.

>>> for s in skyrefsystems.skyrefs_list:
>>>    print(s.fullname, s.description, s.idnum)

For programmers who need to access the id’s of the sky and reference systems: External modules can set their own variables. Here are some examples how one can do this.

Example with copy of celestial’s variables:

  • eq = celestial.eq

  • ec = celestial.ecl

  • ga = celestial.gal etc.

Example with minimal match:

  • eq = celestial.skyrefsystems.minmatch2skyref('EQUA')[0].idnum

  • ec = celestial.skyrefsystems.minmatch2skyref('ecli')[0].idnum

Read this as: get the object for which a minimal match is found. Item [0] is the object (the other is the number of times a match is found). The ‘idnum’ is the integer for which we can identify a system.

Or use the equivalent with method skyrefset.minmatch2id():

  • eq = celestial.skyrefsystems.minmatch2id('EQUA')

  • ec = celestial.skyrefsystems.minmatch2id('ecli')

Example with full name (case sensitive!):

  • eq = celestial.skyrefsystems.fullname2id('EQUATORIAL')

  • ec = celestial.skyrefsystems.fullname2id('ECLIPTIC')

Classes

class kapteyn.celestial.skyrefsys(fullname, idnum, description, refsystem)[source]

Class creates an object that describes a sky- or reference system. This module initializes a set of systems. They are accessible through methods in class celestial.skyrefset

Parameters:
  • fullname (String) – Complete name to identify the system, e.g. “EQUATORIAL”

  • idnum (Integer) – A unique integer to identify the system

  • description (String) – A short description of the system

  • refsystem (Boolean) – Is this system a reference system?

Attributes:

fullname

A string to identify a system, e.g. “EQUATORIAL”.

idnum

A unique integer to identify the system.

description

A string to describe the system.

refsystem

If True then this system is a reference system. Else it is a sky system.

class kapteyn.celestial.skyrefset[source]

A container with sky- and reference system objects from class skyrefsys. It is used to initialize variables that can be used as identifiers for sky- or reference systems. Applications can use its methods to retrieve information given an integer identifier or (part of) a string.

For example when we want a list with all the supported systems then type:

>>> for s in skyrefsystems.skyrefs_list:
>>>    print(s.fullname, s.description, s.idnum)

Attributes:

skyrefs_list

The list with systems

skyrefs_id

A dictionary with the systems and with id’s as keys

skyrefs_fullname

A dictionary with the systems and with full names as keys

Examples:

Next short script shows how to get a list with sky systems and how to use methods of this class to get data for a system if an (integer) id is found:

from kapteyn.celestial import skyrefsystems

for s in skyrefsystems.skyrefs_list:
   print(s.fullname, s.description, s.idnum)
   i = s.idnum
   print("Full name using id2fullname:", skyrefsystems.id2fullname(i))
   print("Description using id2description:", skyrefsystems.id2description(i))
   print("id of %s with minimal match: %d" % \
         (s.fullname[:3], skyrefsystems.minmatch2skyref(s.fullname[:3])[0].idnum))
   print("id of %s with minimal match, alternative: %d" % \
         (s.fullname[:3], skyrefsystems.minmatch2id(s.fullname[:3])))
   print("id of %s with full name: %d" % \
         (s.fullname[:3], skyrefsystems.fullname2id(s.fullname)))
append(skyrefsys)[source]

Append sky reference systems

Parameters:

skyrefsys (Instance of class skyrefsys) – Append this system to the list with supported systems

Returns:

A unique integer id which can be used to identify a system.

fullname2id(fullname)[source]

This is the fastest method to get an integer id from a string which represents a sky system or a reference system. Note that the routine is case sensitive because it uses the full names as keys in a dictionary. The parameter fullname therefore must be in in capitals!

Parameters:

fullname (String) – The full descriptive name of a system e.g. “EQUATORIAL”

Returns:

Integer id of the found system or None if nothing was found.

id2description(idnum)[source]

Given an integer id of a system, return the description of the corresponding system.

Parameters:

idnum (Integer) – Integer id of a system

Returns:

A short description of the corresponding system or an empty string if nothing was found.

id2fullname(idnum)[source]

Given an integer id of a system, return the full name of the corresponding system.

Parameters:

idnum (Integer) – Integer id of a system

Returns:

Full name (e.g. “EQUATORIAL”) of the corresponding system or an empty string if nothing was found.

id2skyref(idnum)[source]

Given an integer id of a system, return the corresponding system as an instance of class skyrefsys. Usually the calling environment will deal with the attributes of this object, for instance to write a short description of the system.

Parameters:

idnum (Integer) – Integer id of a system

Returns:

Instance of class skyrefsys or None if there was not a corresponding system.

minmatch2id(s)[source]

From the found skyrefsys object corresponding to string s, return the idnum attribute. Case insensitive minimal match is used to find the sky- or reference system. Return None if there was no match or more than one match.

Parameters:

s (String) – Part of the string name of a system

Returns:

Instance of class skyrefsys or None if there was not a match or more than one match.

minmatch2skyref(s)[source]

Return the relevant skyrefsys object with the number of times it is matched or return None if nothing was found.

Parameters:

s (String) – Part of the string name of a system

Returns:

Instance of class skyrefsys and the number of times that the input string gives a match.

qqqqmembers: append, minmatch2skyref, minmatch2id, fullname2id, id2skyref, id2fullname, id2description

Core Functions

kapteyn.celestial.skyparser(skyin)[source]

Parse a string, tuple or single integer that represents a sky definition. A sky definition can consist of a sky system, a reference system, an equinox and an epoch of observation. See also the description at Sky definitions. The elements in the string are separated by a comma or a space. The order of the elements is not important. The string is converted to a tuple by parseskydef().

The parser is used in function skymatrix() and sky2sky(). External applications can use this function to check whether user input is valid.

Definitions in strings are usually used to define output sky definitions in prompts or on command lines. Applications can use integer id’s for the sky- and reference systems. These integer id’s are global constants See also Sky systems and Reference systems.

The sky system and reference system strings are minimal matched (case INsensitive) with the strings in the table in the documentation at Sky systems and Reference systems.

For the epoch syntax read the documentation at Epochs for the equinox and epoch of observation. Note that an epoch of observation is either a second epoch in the string (the first is always the equinox) or the epoch string has a suffix ‘_’ which may be follwed by arbitrary characters.

Parameters:

skyin (String, tuple or integer) – Represents a sky definition. See examples.

Returns:

A tuple with the ‘coded’ system where strings for sky- and reference systems are replaced by integer id’s. Missing values are filled in with defaults.

If an error occurred then an exception will be raised.

Raises:
ValueError

From parseskydef():

  • Empty string!

  • Too many items for sky definition!

  • … is ambiguous sky or reference system!

  • … is not a valid epoch or sky/ref system!

From this function:

  • Sky definition is not a string nor a tuple!

  • Too many elements in sky definition (max. 4)!

  • Two sky systems given!

  • Two reference systems given!

  • Invalid number for sky- or reference system!

  • Cannot determine the sky system!

  • Input contains an element that is not an integer or a string!

Examples:
>>> print(celestial.skyparser("B1983.5_O fk4 B1960,eq"))
(0, 4, 1960.0, 1983.5)
>>> print(celestial.skyparser("su"))
(3, None, None, None)
>>> print(celestial.skyparser("supergal"))
(3, None, None, None)
Notes:

This is the parser for a sky definition. In this definition one can specify the sky system, the reference system, an equinox and an epoch of observation if the reference system is fk4. The order of these elements is not important.

The rules for the defaults are:

  • What if the sky system is not defined? If there is a reference system then we assume it is equatorial (could have been ecliptic).

  • If there no sky system and no reference system but there is an equinox, assume sky system is equatorial (could have been ecliptic).

  • If there no sky system and no reference system and no equinox but there is an epoch of observation, assume sky system is equatorial.

  • Assume we have a sky system. What if there is no reference system? Standard in FITS: RADESYS (i.e our reference system) defaults to IRCS unless EQUINOX is given alone, in which case it defaults to FK4 prior to 1984 and FK5 after 1984.

  • Assume we have a sky system and a reference system and the sky system was ecliptic or equatorial. What if we don’t have an equinox? Standard in FITS: EQUINOX defaults to 2000 unless RADESYS is FK4, in which case it defaults to 1950.

  • We have one item to address and that is the epoch of observation. This epoch of observation only applies to the reference systems FK4 and FK4_NO_E. In ‘Representations of celestial coordinates in FITS’ (Calabretta & Greisen) we read that all reference systems are allowed for both equatorial- and ecliptic coordinates, except FK4-NO-E, which is only allowed for equatorial coordinates. If FK4-NO-E is given in combination with an ecliptic sky system then silently FK4 is assumed.

kapteyn.celestial.skymatrix(skyin, skyout)[source]

Create a transformation matrix to be used to transform a position from one sky system to another (including epoch transformations). For a description of the sky definitions see Sky definitions.

Parameters:
  • skyin (Integer or tuple with one to four elements) – One of the supported sky systems or a tuple for equatorial systems which are identified with an equinox an reference system. This is the sky system from which you want to transform to another sky system (skyout).

  • skyout – The destination sky system

Returns:

Three elements:

  • The transformation matrix M for the transformation of positions in (x,y,z) as in XYZskyout = M * XYZskyin

  • followed by ‘None’ or a tuple with the e-term vector belonging input epoch.

  • followed by None or a tuple with the e-term vector belonging to the output epoch.

See also notes below.

Notes:

The reference systems FK4 and FK4_NO_E are special. We consider FK4 as a catalog position where the e-terms are included. So besides a transformation matrix, this function should also return a flag for the addition or removal of e-terms. This flag is either None or the e-term vector which depends on the epoch.

The structure of the output then is as follows: M, (A1,A2,A3), (A4,A5,A6) where:

  • M: The 3x3 transformation matrix

  • (A1,A2,A3) or None: for adding or removing e-terms in the input sky system using this e-term vector (A1,A2,A3).

  • (A4,A5,A6) or None: for adding or removing e-terms in the output sky system using this e-term vector (A4,A5,A6).

This function is the main function of this module. It calls function skyparser() for the parsing of the input and rotmatrix() to get the rotation matrix. There utility function sky2sky() transforms a sequence of longitudes and latitudes from one sky system to another. It is a valuable tool for experiments in an interactive Python session.

Examples:

Some examples of transformations between sky systems using either strings or tuples. We advise to use strings which is more safe then using variables from celestial (which can be accidentally replaced by other values). Note that for transformations where FK4 is involved, the matrix is followed by a vector with e-terms.

>>> from kapteyn import celestial
>>> print(skymatrix(celestial.gal,(celestial.eq,"j2000",celestial.fk5)))
(matrix([[-0.05487554,  0.49410945, -0.86766614],
         [-0.8734371 , -0.44482959, -0.19807639],
         [-0.48383499,  0.74698225,  0.45598379]]),
      None,
      None)
>>> print(skymatrix(celestial.fk4, celestial.fk5))
(matrix([[  9.99925679e-01,  -1.11814832e-02,  -4.85900382e-03],
         [  1.11814832e-02,   9.99937485e-01,  -2.71625947e-05],
         [  4.85900377e-03,  -2.71702937e-05,   9.99988195e-01]]),
      (-1.6255503575995309e-06,
         -3.1918587795578522e-07,
         -1.3842701121066153e-07), None)
>>> print(skymatrix("eq,B1950.0,fk4_no_e","eq,B1950.0,fk4"))
(matrix([[ 1.,  0.,  0.],
         [ 0.,  1.,  0.],
         [ 0.,  0.,  1.]]),
      None,
      (-1.6255503575995309e-06,
         -3.1918587795578522e-07,
         -1.3842701121066153e-07))
>>> print(skymatrix("eq b1950 fk4 j1983.5", "eq J2000 fk5"))
(matrix([[  9.99925679e-01,  -1.11818698e-02,  -4.85829658e-03],
         [  1.11818699e-02,   9.99937481e-01,  -2.71546879e-05],
         [  4.85829648e-03,  -2.71721706e-05,   9.99988198e-01]]),
      (-1.6255503575995309e-06,
         -3.1918587795578522e-07,
         -1.3842701121066153e-07),
      None)
>>> print(skymatrix("eq J2000 fk4 F1984-1-1T0:30", "eq J2000 fk5"))
(matrix([[  1.00000000e+00,  -5.45185721e-06,  -3.39404820e-07],
         [  5.45185723e-06,   1.00000000e+00,   2.24950276e-08],
         [  3.39404701e-07,  -2.24971595e-08,   1.00000000e+00]]),
      (-1.6181121582090453e-06,
         -3.4112123324131958e-07,
         -1.4789407828956555e-07),
      None)

See Epochs for the equinox and epoch of observation for the possible epoch formats.

kapteyn.celestial.sky2sky(skyin, skyout, lons, lats)[source]

Utility function to facilitate command line use of skymatrix.

Parameters:
  • skyin (See function skymatrix()) – The input sky definition

  • skyout (See function skymatrix()) – The output sky definition

  • lons (Floating point number(s), scalar, list or tuple) – Input longitude(s)

  • lats (Floating point number(s), scalar, list or tuple) – Input latitude(s)

Returns:

Matrix. One position per row. See example below how to extract rows, columns and elements from this matrix.

Example:

Interactive Python session:

>>> from kapteyn import celestial
>>> M = celestial.sky2sky( (celestial.eq, celestial.fk5), celestial.gal,
                            (0,0,1.0), (10,20,20) )
>>> M
matrix([[ 102.6262244 ,  -50.83256452],
        [ 106.78021643,  -41.25289649],
        [ 107.9914125 ,  -41.49143448]])
>>> M[2,0]
107.99141249678289
>>> M[0]         # Extract first transformed long, lat
matrix([[ 102.6262244 ,  -50.83256452]])
>>> M[:,1]       # Extract second column with latitudes
matrix([[-50.83256452],
        [-41.25289649],
        [-41.49143448]])
Notes:

This function illustrates the core use of module celestial. First it converts the input of world coordinates into a matrix. This matrix is converted to spatial positions (X,Y,Z) with function longlat2xyz(). The function dotrans() transforms these positions (X,Y,Z) to positions (X2,Y2,Z2) in the output sky system. Then the function xyz2longlat() converts these positions into longitudes and latitudes and finally a matrix with these values is returned:

lonlat = n.array( [(lons,lats)] )
xyz = longlat2xyz(lonlat)
xyz2 = dotrans(skymatrix(skyin, skyout), xyz)
newlonlats = xyz2longlat(xyz2)
return newlonlats
kapteyn.celestial.epochs(spec)[source]

Flexible epoch parser. The functions in this module have different input parameters (Julian epoch, Besselian epochs, Julian dates) because the algorithms came from different sources. What we needed was a routine that could convert a string which represents a date in various formats, to values for a Julian epoch, Besselian epochs and a Julian date. This function returns these value for any valid input date.

For the epoch syntax read the documentation at Epochs for the equinox and epoch of observation. Note that an epoch of observation is either a second epoch in the string (the first is always the equinox) or the epoch string has a suffix ‘_’ which may be follwed by arbitrary characters.

Parameters:

spec (String) – An epoch specification (see below)

Returns:

Calculated corresponding Besselian epoch, Julian epoch and Julian date. Return in order: B, J, JD

Reference:

Various sources listing Julian dates.

Notes:

Examples:

Some checks:

>>> celestial.epochs('F2008-03-31T8:09')  # should return:
    (2008.2474210134737, 2008.2459673739454, 2454556.8395833336)
>>> celestial.epochs('F2007-01-14T13:18:59.9')
    (2007.0378545262108, 2007.0364267212976, 2454115.0548599539)
>>> celestial.epochs("j2007.0364267212976")
    (2007.0378545262108, 2007.0364267212976, 2454115.0548599539)
>>> celestial.epochs("b2007.0378545262108")
    (2007.0378545262108, 2007.0364267212976, 2454115.0548599539)

Utility functions

kapteyn.celestial.JD(year, month, day)[source]

Calculate Julian day number (Julian date)

Parameters:
  • year (Integer) – Year (nnnn)

  • month (Integer) – Month (nn)

  • day (Floating point number) – Day (nn.n…)

Returns:

Julian day number jd.

Reference:

Meeus, Astronomical formula for Calculators, 2nd ed, 1982

Notes:

Months start at 1. Days start at 1. The Julian day begins at Greenwich mean noon, i.e. at 12h. So Jan 1, 1984 at 0h is entered as JD(1984,1,1) and Jan 1, 1984 at 12h is entered as JD(1984,1,1.5)

There is a jump at JD(1582,10,15) caused by a change of calendars. For dates after 1582-10-15 one enters a date from the Julian calendar and before this date you enter a date from the Gregorian calendar.

Examples:
  • Julian date of JD reference: print(celestial.JD(-4712,1,1.5)) ==> 0.0

  • The first day of 1 B.C.: print(celestial.JD(0,1,1)) ==> 1721057.5

  • Last day before Gregorian reform: print(celestial.JD(1582,10,4)) ==> 2299159.5

  • First day of Gregorian reform: print(celestial.JD(1582,10,15)) ==> 2299170.5

  • Half a day later: print(celestial.JD(1582,10,15.5)) ==> 2299161.0

  • Unix reference: print(celestial.JD(1970,1,1)) ==> 2440587.5

kapteyn.celestial.lon2hms(a, prec=1, delta=None, tex=False)[source]

Convert an angle in degrees to hours, minutes, seconds format.

Parameters:
  • a (Floating point number) – Angle (in degrees) for which we want to create a formatted text label.

  • prec (Integer) – The required number of decimals in the seconds part of output. If a value is omitted, then the default is 1.

  • delta (None or a floating point number) – If one labels world coordinates along an axis then the default labels are in hours, minutes and seconds with some decimal number. This is probably not want you want if the step size between subsequent positions is for example an integer number of degrees or minutes. Then you want labels showing only hours or hours and minutes. This function tries to find out whether this is the case (given a value for delta) or not. If so, a minimum length label is returned.

  • tex (Boolean) – The default is False. If set to True, the string is formatted in LaTeX. Such labels can be plotted in, for example, Matplotlib.

Returns:

Formatted string representing the input angle.

Notes:

Longitudes are forced into the range, 360 deg. and then converted to hours, minutes and seconds.

Examples:

Format a position in hms and dms:

>>> ra = 359.9999
>>> dec = 0.0000123
>>> print(celestial.lon2hms(ra),  celestial.lat2dms(dec))
    00h 00m  0.0s +00d 00m  0.0s
>>> print(celestial.lon2hms(ra, 2),  celestial.lat2dms(dec, 2))
    23h 59m 59.98s +00d 00m  0.04s
>>> print(celestial.lon2hms(ra, 4),  celestial.lat2dms(dec, 4))
    23h 59m 59.9760s +00d 00m  0.0443s
kapteyn.celestial.lat2dms(a, prec=1, delta=None, tex=False)[source]

Convert an angle in degrees into the degrees, minutes, seconds format assuming it was a latitude. Its value should be in the range -90 to 90 degrees

Parameters:
  • a (Floating point number) – Angle (in degrees) for which we want to create a formatted text label.

  • prec (Integer) – The required number of decimals in the seconds part of output. If a value is omitted, then the default is 1.

  • delta (None or a floating point number) – If one labels world coordinates along an axis then the default labels are in degrees, minutes and seconds with some decimal number. This is probably not want you want if the step size between subsequent positions is for example an integer number of degrees or minutes. Then you want labels showing only degrees or degrees and minutes. This function tries to find out whether this is the case (given a value for delta) or not. If so, a minimum length label is returned.

  • tex (Boolean) – The default is False. If set to True, the string is formatted in LaTeX. Such labels can be plotted in, for example, Matplotlib.

Returns:

Formatted string representing the input angle or a string with ‘#’ characters indicating that the input was out of range.

Notes:

The HMS and DMS format should be treated differently because their ranges in world coordinates are different. Longitudes should be in range of (0,360) degrees. So -10 deg is in fact 350 deg. and 370 deg is in fact 10 deg. Latitudes range from -90 to 90 degrees. Then 91 degrees is in fact 89 degrees but at a longitude that is separated 180 deg. from the stated longitude. But we don’t have control over the longitudes here so the only thing we can do is reject the value and return a dummy string.

kapteyn.celestial.lon2dms(a, prec=1, delta=None, tex=False)[source]

Convert an angle in degrees to degrees, minutes, seconds format, assuming the input is a longitude but not associated with an equatorial system.

Parameters:
  • a (Floating point number) – Angle (in degrees) for which we want to create a formatted text label

  • prec (Integer) – The required number of decimals in the seconds part of output If a value is omitted, then the default is 1.

  • delta (None or a floating point number) – If one labels world coordinates along an axis then the default labels are in hours, minutes and seconds with some decimal number. This is probably not want you want if the step size between subsequent positions is for example an integer number of degrees or minutes. Then you want labels showing only degrees or degrees and minutes. This function tries to find out whether this is the case (given a value for delta) or not. If so, a minimum length label is returned.

  • tex (Boolean) – The default is False. If set to True, the string is formatted in LaTeX. Such labels can be plotted in, for example, Matplotlib.

Returns:

Formatted string representing the input angle.

Notes:

Longitudes are forced into the range 0, 360 deg. and then converted to hours, minutes and seconds.

Examples:

Format a longitude to dms:

>>> print(celestial.lon2dms(167.342, 4))
   167d 20m 31.2000s
>>> print(celestial.lon2dms(-10, 4))
   350d  0m  0.0000s
kapteyn.celestial.JD2epochBessel(JD)[source]

Convert a Julian date to a Besselian epoch.

Parameters:

JD (Floating point number) – Julian date (e.g. 2445700.5)

Returns:

Besselian epoch (e.g. 1983.9)

Reference:

Standards Of Fundamental Astronomy,

http://www.iau-sofa.rl.ac.uk/2003_0429/sofa/epb.html

Notes:

e.g. 2445700.5 -> 1983.99956681

One Tropical Year is 365.242198781 days and JD(1900) = 2415020.31352

If we know the JD then the Besselian epoch can be calculated with:

BE = B[1900 + (JD - 2415020.31352)/365.242198781]

Expression corresponds to the IAU SOFA expression in the reference with: 2451545-36524.68648 = 2415020.31352

kapteyn.celestial.epochBessel2JD(Bepoch)[source]

Convert a Besselian epoch to a Julian date

Parameters:

Bepoch (Floating point number) – Besselian epoch in format nnnn.nn

Returns:

Julian date

Reference:

See: JD2epochBessel()

Notes:

e.g. 1983.99956681 converts into 2445700.5 It’s the inverse of JD2epochBessel()

kapteyn.celestial.JD2epochJulian(JD)[source]

Convert a Julian date to a Julian epoch

Parameters:

JD (Floating point number) – Julian date

Returns:

Julian epoch

Reference:

Standards Of Fundamental Astronomy,

http://www.iau-sofa.rl.ac.uk/2003_0429/sofa/epj.html

Notes:

e.g. 2445700.5 converts into 1983.99863107 Assuming years of exactly 365.25 days, we can calculate a Julian epoch from a Julian date. Expression corresponds to IAU SOFA routine ‘epj’

kapteyn.celestial.epochJulian2JD(Jepoch)[source]

Convert a Julian epoch to a Julian date

Parameters:

Jepoch (Floating point number) – Julian epoch (in format nnnn.nn)

Returns:

Julian date

Reference:

See JD2epochJulian()

Notes:

e.g. 1983.99863107 converts into 2445700.5 It’s the inverse of function JD2epochJulian

kapteyn.celestial.obliquity1980(jd)[source]

What is the obliquity of the ecliptic at this Julian date? (IAU 1980 model)

Parameters:

jd (Floating point number) – Julian date

Returns:

Mean obliquity in degrees

Reference:

Explanatory Supplement to the Astronomical Almanac, P. Kenneth Seidelmann (ed), University Science Books (1992), Expression 3.222-1 (p114).

Notes:

The epoch is entered in Julian date and the time is calculated w.r.t. J2000.

The obliquity is the angle between the mean equator and ecliptic, or, between the ecliptic pole and mean celestial pole of date

kapteyn.celestial.obliquity2000(jd)[source]

What is the obliquity of the ecliptic at this Julian date? (IAU model 2000)

Parameters:

jd (Floating point number) – Julian date

Returns:

Mean obliquity in degrees

Reference:

Fukushima, T. 2003, AJ, 126,1 Kaplan, H., 2005, The IAU Resolutions on Astronomical Reference Systems, Time Scales, and Earth Rotation Models, United States Naval Observatory circular no. 179, http://aa.usno.navy.mil/publications/docs/Circular_179.pdf (page 44)

Notes:

The epoch is entered in Julian date and the time is calculated w.r.t. J2000.

The obliquity is the angle between the mean equator and ecliptic, or, between the ecliptic pole and mean celestial pole of date.

kapteyn.celestial.IAU2006precangles(epoch)[source]

Calculate IAU 2000 precession angles for precession from input epoch to J2000.

Parameters:

epoch (Floating point number) – Julian epoch of observation.

Returns:

Angles \(\zeta\) (zeta), z, \(\theta\) (theta) in degrees to setup a rotation matrix to transform from J2000 to input epoch.

Reference:

Capitaine N. et al., IAU 2000 precession A&A 412, 567-586 (2003)

Notes:

Input are Julian epochs! T = (jd-2451545.0)/36525.0 Combined with jd = Jepoch-2000.0)*365.25 + 2451545.0 gives: (see module code at function epochJulian2JD(epoch)) T = (epoch-2000.0)/100.0

This function should be updated as soon as there are IAU2006 adopted angles to replace the angles used in this function.

kapteyn.celestial.Lieskeprecangles(jd1, jd2)[source]

Calculate IAU 1976 precession angles for a precession of epoch corresponding to Julian date jd1 to epoch corresponds to Julian date jd2.

Parameters:
  • jd1 (Floating point number) – Julian date for start epoch

  • jd2 (Floating point number) – Julian date for end epoch

Returns:

Angles \(\zeta\) (zeta), z, \(\theta\) (theta) degrees

Reference:

Lieske,J.H., 1979. Astron.Astrophys.,73,282. equations (6) & (7), p283.

Notes:

The ES (Explanatory Supplement to the Astronomical Almanac) lists for a IAU1976 precession from 1984, January 1d0h to J2000 the angles in arcsec: xi_a=368.9985, ze_a=369.0188 and th_a=320.7279 Using the functions in this module, this can be calculated by applying:

>>> jd1 = celestial.JD(1984,1,1)
>>> jd2 = celestial.JD(2000,1,1.5)
>>> print(celestial.Lieskeprecangles(jd1, jd2))
   (0.10249958598931658, 0.10250522534285664, 0.089091092843880629)
>>> print([a*3600 for a in angles])
    [368.99850956153966, 369.01881123428387, 320.72793423797026]

The function returns values in degrees, while literature values often are listed in seconds of arc.

Lieske’s fit belongs to the so called Quasi-Linear Types Below a table with the precision (according to IAU SOFA):

  • 1960AD to 2040AD: < 0.1”

  • 1640AD to 2360AD: < 1”

  • 500BC to 3000AD: < 3”

  • 1200BC to 3900AD: > 10”

  • < 4200BC or > 5600AD: > 100”

  • < 6800BC or > 8200AD: > 1000”

kapteyn.celestial.Newcombprecangles(epoch1, epoch2)[source]

Calculate precession angles for a precession in FK4, using Newcomb’s method (Woolard and Clemence angles)

Parameters:
  • epoch1 (Floating point number) – Besselian start epoch

  • epoch2 (Floating point number) – Besselian end epoch

Returns:

Angles \(\zeta\) (zeta), z, \(\theta\) (theta) degrees

Reference:

ES 3.214 p.106

Notes:

Newcomb’s precession angles for old catalogs (FK4), see ES 3.214 p.106. Input are Besselian epochs! Adopted accumulated precession angles from equator and equinox at B1950 to 1984 January 1d 0h according to ES (table 3.214.1, p 107) are: zeta=783.7092, z=783.8009 and theta=681.3883 The Woolard and Clemence angles (derived in this routine) are: zeta=783.70925, z=783.80093 and theta=681.38830 (see same ES table as above).

This routine found (in seconds of arc): zeta,z,theta =  783.709246271 783.800934641 681.388298284 for t1 = 0.1 and t2 = 0.133999566814 using the lines in the next example.

Examples:

From an interactive Python session:

>>> b1 = 1950.0
>>> b2 = celestial.epochs("F1984-01-01")[0]
>>> print([x*3600 for x in celestial.Newcombprecangles(be1, be2)])
    [783.70924627097793, 783.80093464073127, 681.38829828393466]

Rotation matrices

kapteyn.celestial.MatrixEqB19502Gal()[source]

Create matrix to convert equatorial fk4 coordinates (without e-terms) to IAU 1958 lII,bII system of galactic coordinates.

Parameters:

None

Results:

3x3 Matrix M as in XYZgal = M * XYZb1950

Reference:
  1. Blaauw, A., Gum C.S., Pawsey, J.L., Westerhout, G.: 1958,

  2. Monthly Notices Roy. Astron. Soc. 121, 123,

  3. Blaauw, A., 2007. Private communications.

Notes:

Original definitions from 1.:

  • The new north galactic pole lies in the direction alpha = 12h49m (192.25 deg), delta=27.4 deg (equinox 1950.0).

  • The new zero of longitude is the great semicircle originating at the new north galactic pole at the position angle theta = 123 deg with respect to the equatorial pole for 1950.0.

  • Longitude increases from 0 to 360 deg. The sense is such that, on the galactic equator increasing galactic longitude corresponds to increasing Right Ascension. Latitude increases from -90 deg through 0 deg to 90 deg at the new galactic pole.

Given the RA and Dec of the galactic pole, and using the Euler angles scheme:

M = rotZ(a3).rotY(a2).rotZ(a1)

We first rotate the spin vector of the XY plane about an angle a1 = ra_pole and then rotate the spin vector in the XZ plane (i.e. around the Y axis) with an angle a2=90-dec_pole to point it in the right declination.

Now think of a circle with the galactic pole as its center. The radius is equal to the distance between this center and the equatorial pole. The zero point now is on the circle and opposite to this pole.

We need to rotate along this circle (i.e. a rotation around the new Z-axis) in a way that the angle between the zero point and the equatorial pole is equal to 123 deg. So first we need to compensate for the 180 deg of the current zero longitude, opposite to the pole. Then we need to rotate about an angle 123 deg but in a way that increasing galactic longitude corresponds to increasing Right Ascension which is opposite to the standard rotation of this circle (note that we rotated the original X axis about 192.25 deg). The last rotation angle therefore is a3=+180-123:

M = rotZ(180-123.0)*rotY(90-27.4)*rotZ(192.25)

The composed rotation matrix is the same as in Slalib’s ‘ge50.f’ and the matrix in eq. (32) of Murray (1989).

kapteyn.celestial.MatrixGal2Sgal()[source]

Transform galactic to supergalactic coordinates

Parameters:

None

Returns:

Matrix M as in XYZsgal = M * XYZgal

Reference:

Lahav, O., The supergalactic plane revisited with the Optical Redshift Survey Mon. Not. R. Astron. Soc. 312, 166-176 (2000)

Notes:

The Supergalactic equator is conceptually defined by the plane of the local (Virgo-Hydra-Centaurus) supercluster, and the origin of supergalactic longitude is at the intersection of the supergalactic and galactic planes. (de Vaucouleurs)

North SG pole at l=47.37 deg, b=6.32 deg. Node at l=137.37, sgl=0 (inclination 83.68 deg).

Older references give for he position of the SG node 137.29 which differs from 137.37 deg in the official definition.

For the rotation matrix we chose the scheme Rz.Ry.Rz Then first we rotate about 47.37 degrees along the Z-axis followed by a rotation about 90-6.32 degrees is needed to set the pole to the right declination. The new plane intersects the old one at two positions. One of them is l=137.37, b=0 (in galactic coordinates). If we want this to be sgl=0 we have to rotate this plane along the new Z-axis about an angle of 90 degrees. So the composed rotation matrix is:

M = Rotz(90)*Roty(90-6.32)*Rotz(47.37)
kapteyn.celestial.MatrixEq2Ecl(epoch, S1)[source]

Calculate a rotation matrix to convert equatorial coordinates to ecliptical coordinates

Parameters:
  • epoch (Floating point number) – Epoch of the equator and equinox of date

  • S1 (Integer) – equatorial system to determine if one entered epoch in B or J coordinates.

Returns:

3x3 Matrix M as in XYZecl = M * XYZeq

Reference:

Representations of celestial coordinates in FITS, Calabretta. M.R., & Greisen, E.W., (2002) Astronomy & Astrophysics, 395, 1077-1122. http://www.atnf.csiro.au/people/mcalabre/WCS/ccs.pdf

Notes:
  1. The origin for ecliptic longitude is the vernal equinox. Therefore the coordinates of a fixed object is subject to shifts due to precession. The rotation matrix uses the obliquity to do the conversion to the wanted ecliptic coordinates. So we always need to enter an epoch. Usually this is J2000, but it can also be the epoch of date. The additional reference system indicates whether we need a Besselian or a Julian epoch.

  2. In the FITS paper of Calabretta and Greisen (2002), one observes the following relations to FITS:

    -Keyword RADESYSa sets the catalog system FK4, FK4-NO-E or FK5 This applies to equatorial and ecliptical coordinates with the exception of FK4-NO-E.

    -FK4 coordinates are not strictly spherical since they include a contribution from the elliptic terms of aberration, the so-called e-terms which amount to max. 343 milliarcsec. FITS paper: ‘Strictly speaking, therefore, a map obtained from, say, a radio synthesis telescope, should be regarded as FK4-NO-E unless it has been appropriately re-sampled or a distortion correction provided. In common usage, however, CRVALia for such maps is usually given in FK4 coordinates. In doing so, the e-terms are effectively corrected to first order only.’. (See also ES, eq. 3.531-1 page 170.

    -Keyword EQUINOX sets the epoch of the mean equator and equinox.

    -Keyword EPOCH is often used in older FITS files. It is a deprecated keyword and should be replaced by EQUINOX. It does not require keyword RADESYS. From its value we derive whether the reference system is FK4 or FK5 (the marker value is 1984.0)

    -Ecliptic coordinates require the epoch of the equator and equinox of date. This will be taken as the time of observation rather than EQUINOX.

    FITS paper: ‘The time of observation may also be required for other astrometric purposes in addition to the usual astrophysical uses, for example, to specify when the mean place was correct in accounting for proper motion, including “fictitious” proper motions in the conversion between the FK4 and FK5 systems. The old *DATE-OBS keyword may be used for this purpose. However, to provide a more convenient specification we here introduce the new keyword MJD-OBS’.*

    So MJD-OBS is the modified Julian Date (JD - 2400000.5) of the start of the observation.

  3. Equatorial to ecliptic transformations use the time dependent obliquity of the equator (also known as the obliquity of the ecliptic). Again, start with:

    M = rotZ(0).rotX(eps).rotZ(0) = E.rotX(eps).E = rotX(eps)
    

    In fact this is only a rotation around the X axis

kapteyn.celestial.FK42FK5Matrix(t=None)[source]

Create a matrix to precess from B1950 in FK4 to J2000 in FK5 following to Murray’s (1989) procedure.

Parameters:

t (Floating point number) – Besselian epoch as epoch of observation.

Returns:

3x3 matrix M as in XYZfk5 = M * XYZfk4

Reference:
  • Murray, C.A. The Transformation of coordinates between the systems B1950.0 and J2000.0, and the principal galactic axis referred to J2000.0, Astronomy and Astrophysics (ISSN 0004-6361), vol. 218, no. 1-2, July 1989, p. 325-329.

  • Poppe P.C.R.,, Martin, V.A.F., Sobre as Bases de Referencia Celeste SitientibusSerie Ciencias Fisicas

Notes:

Murray precesses from B1950 to J2000 using a precession matrix by Lieske. Then applies the equinox correction and ends up with a transformation matrix X(0) as given in this function.

In Murray’s article it is proven that using the procedure as described in the article, r_fk5 = X(0).r_fk4 for extra galactic sources where we assumed that the proper motion in FK5 is zero. This procedure is independent of the epoch of observation. Note that the matrix is not a rotation matrix.

FK4 is not an inertial coordinate frame (because of the error in precession and the motion of the equinox. This has consequences for the proper motions. e.g. a source with zero proper motion in FK5 has a fictitious proper motion in FK4. This affects the actual positions in a way that the correction is bigger if the epoch of observation is further away from 1950.0 The focus of this library is on data of which we do not have information about the proper motions. So for positions of which we allow non zero proper motion in FK5 one needs to supply the epoch of observation.

Examples:

Print the difference between the rotation matrix for 1970 and 1980:

>>> M1 = celestial.FK42FK5Matrix(1970)
>>> M2 = celestial.FK42FK5Matrix(1980)
>>> M2 - M1
matrix([[ -2.64546940e-10,  -1.15396722e-07,   2.11108953e-07],
        [  1.15403817e-07,  -1.29040234e-09,   2.36016437e-09],
        [ -2.11125281e-07,  -5.60232514e-10,   1.02585540e-09]])
kapteyn.celestial.ICRS2FK5Matrix()[source]

Create a rotation matrix to convert a position from ICRS to fk5, J2000

Parameters:

None

Returns:

3x3 rotation matrix M as in XYZfk5 = M * XYZicrs

Reference:

Kaplan G.H., The IAU Resolutions on Astronomical Reference systems, Time scales, and Earth Rotation Models, US Naval Observatory, Circular No. 179

Notes:

Return a matrix that converts a position vector in ICRS to FK5, J2000. We do not use the first or second order approximations given in the reference, but use the three rotation matrices from the same paper to obtain the exact result:

M =  rotX(-eta0)*rotY(xi0)*rotZ(da0)

eta0 = -19.9 mas, xi0 = 9.1 mas and da0 = -22.9 mas

kapteyn.celestial.ICRS2J2000Matrix()[source]

Return a rotation matrix for conversion of a position in the ICRS to the dynamical reference system based on the dynamical mean equator and equinox of J2000.0 (called the dynamical J2000 system)

Parameters:

None

Returns:

Rotation matrix to transform positions from ICRS to dyn J2000

Reference:
  • Hilton and Hohenkerk (2004), Astronomy and Astrophysics 413, 765-770

  • Kaplan G.H., The IAU Resolutions on Astronomical Reference systems, Time scales, and Earth Rotation Models, US Naval Observatory, Circular No. 179

Notes:

Return a matrix that converts a position vector in ICRS to Dyn. J2000. We do not use the first or second order approximations given in the reference, but use the three rotation matrices to obtain the exact result:

M = rotX(-eta0)*rotY(xi0)*rotZ(da0)

eta0 = -6.8192 mas, xi0 = -16.617 mas and da0 = -14.6 mas

kapteyn.celestial.JMatrixEpoch12Epoch2(Jepoch1, Jepoch2)[source]

Precession from one epoch to another in the fk5 system. It uses Lieskeprecangles() to calculate the precession angles.

Parameters:
  • Jepoch1 (Floating point number) – Julian start epoch

  • Jepoch2 (Floating point number) – Julian epoch to precess to.

Returns:

3x3 rotation matrix M as in XYZepoch2 = M * XYZepoch1

Reference:

Seidelman, P.K., 1992. Explanatory Supplement to the Astronomical Almanac. University Science Books, Mill Valley. 3.214 p 106

Notes:

The precession matrix is:

M = rotZ(-z).rotY(+theta).rotZ(-zeta)
kapteyn.celestial.BMatrixEpoch12Epoch2(Bepoch1, Bepoch2)[source]

Precession from one epoch to another in the fk4 system. It uses Newcombprecangles() to calculate the precession angles.

Parameters:
  • Bepoch1 (Floating point number) – Besselian start epoch

  • Bepoch2 (Floating point number) – Besselian epoch to precess to.

Returns:

3x3 rotation matrix M as in XYZepoch2 = M * XYZepoch1

Reference:

Seidelman, P.K., 1992. Explanatory Supplement to the Astronomical Almanac. University Science Books, Mill Valley. 3.214 p 106

Notes:

The precession matrix is:

M = rotZ(-z).rotY(+theta).rotZ(-zeta)
kapteyn.celestial.IAU2006MatrixEpoch12Epoch2(epoch1, epoch2)[source]

Create a rotation matrix for a precession based on IAU 2000/2006 expressions, see function IAU2006precangles()

Parameters:
  • epoch1 (Floating point number) – Julian start epoch

  • epoch2 (Floating point number) – Julian epoch to precess to.

Returns:

Matrix to transform equatorial coordinates from epoch1 to epoch2 as in XYZepoch2 = M * XYZepoch1

Reference:

Capitaine N. et al.: IAU 2000 precession A&A 412, 567-586 (2003)

Notes:

Note that we apply this precession only to equatorial coordinates in the system of dynamical J2000 coordinates. When converting from ICRS coordinates this means applying a frame bias. Therefore the angles differ from the precession Fukushima-Williams angles (IAU 2006)

The precession matrix is:

M = rotZ(-z).rotY(+theta).rotZ(-zeta)
kapteyn.celestial.MatrixEpoch12Epoch2(epoch1, epoch2, S1, S2, epobs=None)[source]

Helper function for skymatrix(). It handles precession and the transformation between equatorial systems. This function includes also conversions between reference systems.

Parameters:
  • epoch1 (Floating point number) – Epoch belonging to system S1 depending on the reference system either Besselian or Julian.

  • epoch2 – Epoch belonging to system S2 depending on the reference system either Besselian or Julian.

  • S1 (Integer) – Input reference system

  • S2 (Integer) – Output rreferencesystem

  • epobs (Floating point number) – Epoch of observation. Only valid for conversions between FK4 and FK5.

Returns:

Rotation matrix to transform a position in one of the reference systems S1 with epoch1 to an equatorial system with equator and equinox at epoch2 in reference system S2.

Notes:

Return matrix to transform equatorial coordinates from epoch1 to epoch2 in either reference system FK4 or FK5. Or transform from epoch, FK4 or FK5 to ICRS or J2000 vice versa. Note that each transformation between FK4 and one of the other reference systems involves a conversion to FK5 and therefore the epoch of observation will be involved.

Note that if no systems are entered and the one epoch is > 1984 and the other < 1984, then the transformation involves both sky reference systems FK4 and FK5.

Examples:

Calculate rotation matrix for a conversion between FK4, epoch 1940 to FK5, epoch 1960, while the date of observation was 1950.

>>> from kapteyn import celestial
>>> celestial.MatrixEpoch12Epoch2(1940, 1960, celestial.fk4, celestial.fk5, 1950)
matrix([[  9.99988107e-01,  -4.47301372e-03,  -1.94362889e-03],
        [  4.47301372e-03,   9.99989996e-01,  -4.34712255e-06],
        [  1.94362889e-03,  -4.34680782e-06,   9.99998111e-01]])