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:
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.
Symbol  String  Description 

eq, equatorial  EQUATORIAL  Equatorial coordinates (α, δ), See also next table with reference systems 
ecl, ecliptic  ECLIPTIC  Ecliptic coordinates (λ, β) referred to the ecliptic and mean equinox 
gal, galactic  GALACTIC  Galactic coordinates (lII, bII) 
sgal, supergalactic  SUPERGALACTIC  De Vaucouleurs Supergalactic coordinates (sgl, sgb) 
Symbol  String  Description 

fk4  FK4  Mean place preIAU 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, FK4NOE  The old FK4 (barycentric) equatorial system but without the Eterms 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.
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:

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 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 eterms  Sky system, reference system, and an equinox 
“EQ, fk4noe, B1960”  Equatorial, FK4 no eterms  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. 
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.]
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')
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: 


Attributes:
A string to identify a system, e.g. “EQUATORIAL”.
A unique integer to identify the system.
A string to describe the system.
If True then this system is a reference system. Else it is a sky system.
A container with sky and reference system objects from class celestial.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
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. 
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. 
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. 
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. 
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. 
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. 
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. 
Attributes:
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))


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 celestial.parseskydefs().
The parser is used in function celestial.skymatrix() and celestial.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: 

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:

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: 


Returns:  Three elements:
See also notes below. 
Notes:  The reference systems FK4 and FK4_NO_E are special. We consider FK4 as a catalog position where the eterms are included. So besides a transformation matrix, this function should also return a flag for the addition or removal of eterms. This flag is either None or the eterm vector which depends on the epoch. The structure of the output then is as follows: M, (A1,A2,A3), (A4,A5,A6) where:
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 eterms. >>> 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.99925679e01, 1.11814832e02, 4.85900382e03],
[ 1.11814832e02, 9.99937485e01, 2.71625947e05],
[ 4.85900377e03, 2.71702937e05, 9.99988195e01]]),
(1.6255503575995309e06,
3.1918587795578522e07,
1.3842701121066153e07), 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.6255503575995309e06,
3.1918587795578522e07,
1.3842701121066153e07))
>>> print skymatrix("eq b1950 fk4 j1983.5", "eq J2000 fk5")
(matrix([[ 9.99925679e01, 1.11818698e02, 4.85829658e03],
[ 1.11818699e02, 9.99937481e01, 2.71546879e05],
[ 4.85829648e03, 2.71721706e05, 9.99988198e01]]),
(1.6255503575995309e06,
3.1918587795578522e07,
1.3842701121066153e07),
None)
>>> print skymatrix("eq J2000 fk4 F198411T0:30", "eq J2000 fk5")
(matrix([[ 1.00000000e+00, 5.45185721e06, 3.39404820e07],
[ 5.45185723e06, 1.00000000e+00, 2.24950276e08],
[ 3.39404701e07, 2.24971595e08, 1.00000000e+00]]),
(1.6181121582090453e06,
3.4112123324131958e07,
1.4789407828956555e07),
None)
See Epochs for the equinox and epoch of observation for the possible epoch formats. 
Utility function to facilitate command line use of skymatrix.
Parameters: 


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

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('F20080331T8:09') # should return:
(2008.2474210134737, 2008.2459673739454, 2454556.8395833336)
>>> celestial.epochs('F20070114T13: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)

Calculate Julian day number (Julian date)
Parameters: 


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 15821015 one enters a date from the Julian calendar and before this date you enter a date from the Gregorian calendar. 
Examples: 

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


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

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: 


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. 
Convert an angle in degrees to degrees, minutes, seconds format, assuming the input is a longitude but not associated with an equatorial system.
Parameters: 


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

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, 
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: 245154536524.68648 = 2415020.31352 
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() 
Convert a Julian date to a Julian epoch
Parameters:  JD (Floating point number) – Julian date 

Returns:  Julian epoch 
Reference:  Standards Of Fundamental Astronomy, 
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’ 
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 
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.2221 (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 
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. 
Calculate IAU 2000 precession angles for precession from input epoch to J2000.
Parameters:  epoch (Floating point number) – Julian epoch of observation. 

Returns:  Angles ζ (zeta), z, θ (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, 567586 (2003) 
Notes:  Input are Julian epochs! T = (jd2451545.0)/36525.0 Combined with jd = Jepoch2000.0)*365.25 + 2451545.0 gives: (see module code at function epochJulian2JD(epoch)) T = (epoch2000.0)/100.0 This function should be updated as soon as there are IAU2006 adopted angles to replace the angles used in this function. 
Calculate IAU 1976 precession angles for a precession of epoch corresponding to Julian date jd1 to epoch corresponds to Julian date jd2.
Parameters: 


Returns:  Angles ζ (zeta), z, θ (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 QuasiLinear Types Below a table with the precision (according to IAU SOFA):

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


Returns:  Angles ζ (zeta), z, θ (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("F19840101")[0]
>>> print [x*3600 for x in celestial.Newcombprecangles(be1, be2)]
[783.70924627097793, 783.80093464073127, 681.38829828393466]

Create matrix to convert equatorial fk4 coordinates (without eterms) to IAU 1958 lII,bII system of galactic coordinates.
Parameters:  None 

Results:  3x3 Matrix M as in XYZgal = M * XYZb1950 
Reference: 

Notes:  Original definitions from 1.:
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=90dec_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 Zaxis) 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=+180123: M = rotZ(180123.0)*rotY(9027.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). 
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, 166176 (2000) 
Notes:  The Supergalactic equator is conceptually defined by the plane of the local (VirgoHydraCentaurus) 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 Zaxis followed by a rotation about 906.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 Zaxis about an angle of 90 degrees. So the composed rotation matrix is: M = Rotz(90)*Roty(906.32)*Rotz(47.37)

Calculate a rotation matrix to convert equatorial coordinates to ecliptical coordinates
Parameters: 


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, 10771122. http://www.atnf.csiro.au/people/mcalabre/WCS/ccs.pdf 
Notes: 

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: 

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.64546940e10, 1.15396722e07, 2.11108953e07],
[ 1.15403817e07, 1.29040234e09, 2.36016437e09],
[ 2.11125281e07, 5.60232514e10, 1.02585540e09]])

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 
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: 

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 
Precession from one epoch to another in the fk5 system. It uses Lieskeprecangles() to calculate the precession angles.
Parameters: 


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)

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


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)

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


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, 567586 (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 FukushimaWilliams angles (IAU 2006) The precession matrix is: M = rotZ(z).rotY(+theta).rotZ(zeta)

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


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.99988107e01, 4.47301372e03, 1.94362889e03],
[ 4.47301372e03, 9.99989996e01, 4.34712255e06],
[ 1.94362889e03, 4.34680782e06, 9.99998111e01]])
