sky1.h - astronomical coordinate conversion routines, IAU 1980 More...

Data Structures
struct	Sky1_Prec1976
	Precession angles. More...

struct	Sky1_Nut1980
	Nutation angles and obliquity. More...

Functions
void	sky1_frameBiasFK5 (V3D_Matrix *biasM)
	Create the frame bias matrix from the IAU 2000 precession-nutation model. More...

void	sky1_precessionIAU1976 (double t0, double t1, Sky1_Prec1976 *terms)
	This procedure calculates the equatorial precession parameters ζ, z, and θ which represent the rotation required to transform the FK5 equatorial reference system of Julian epoch t0 to that of Julian epoch t1, according to the IAU 1976 precession constants. More...

void	sky1_createPrec1976Matrix (const Sky1_Prec1976 terms, V3D_Matrix precM)
	This routine calculates the precession matrix, based on angles ζ, z, and θ. More...

void	sky1_nutationIAU1980 (double t_cy, int precision, Sky1_Nut1980 *nut)
	Calculates the nutation in longitude and obliquity, according to the IAU 1980 Nutation Theory. More...

void	sky1_epsilon1980 (double t_cy, Sky1_Nut1980 *nut)
	Calculate the obliquity of the ecliptic and the equation of the equinoxes. More...

void	sky1_createNut1980Matrix (const Sky1_Nut1980 nut, V3D_Matrix nutM)
	This routine calculates the Nutation matrix, using the nutation in longitude, the nutation in obliquity, and the mean obliquity of the equator. More...

void	sky1_createNPmatrix (double t0_cy, double t1_cy, int precision, V3D_Matrix *npM)
	Call the various precession and nutation routines in this module to create a combined precession and nutation rotation matrix, suitable for converting a celestial object's coordinates from mean place at epoch t0_cy to apparent coordinates at epoch t1_cy. More...

double	sky1_gmSiderealTimeIAU1982 (double du)
	Calculate the Greenwich mean sidereal time. More...

void	sky1_appToTirs (const V3D_Vector appV, double j2kUT1_d, double eqEq_rad, V3D_Vector terInterV)
	Convert a position in geocentric apparent coordinates to geocentric coordinates in the Terrestrial Intermediate Reference System. More...

Detailed Description

sky1.h - astronomical coordinate conversion routines, IAU 1980

Author: David Hoadley

This is one of three alternative modules: sky0.h / sky0.c, sky1.h / sky1.c and sky2.h / sky2.c. They contain routines for transforming astronomical positions from frame to another: precession, nutation, sidereal time etc. and they reflect changes in the International Astronomical Union's precession and nutation theory. The differences are:

sky0.h / sky0.c: nutation, obliquity and sidereal time routines from the NREL Solar Position Algorithm document. These are based on the IAU 1980 nutation theory.
sky1.h / sky1.c: precession, nutation and their associated rotation matrices, obliquity and sidereal time. These are the IAU 1980 precession and nutation theory.
sky2.h / sky2.c: precession, nutation and their associated rotation matrices, obliquity and sidereal time. These are the newer IAU 2000 precession and nutation theory.

This module (sky1.h / sky1.c) contains precession and nutation routines. The precession routine uses the IAU 1976 algorithm. The nutation routine uses the IAU 1980 algorithm.

Definition in file sky1.h.

Function Documentation

◆ sky1_frameBiasFK5()

void sky1_frameBiasFK5 ( V3D_Matrix * biasM )

Create the frame bias matrix from the IAU 2000 precession-nutation model.

This matrix is used to convert coordinates from the Geocentric Celestial Reference System (GCRS) to FK5 reference system. (Although anything to do with the ICRS/GCRS is related to the IAU2000+ precession-nutation theory, this routine is included here to allow ICRS star catalogue positions to be converted to apparent coordinates using the IAU1980 precession-nutation theory, which is quite good enough for tracking. It is not good enough for astrometry, but it will get us to better than within an arcsecond of our desired apparent position).

Parameters

[out] biasM Frame bias matrix B

Reference:: Astronomical Almanac 2007, page B28

Definition at line 339 of file sky1.c.

◆ sky1_precessionIAU1976()

void sky1_precessionIAU1976	(	double	t0_cy,
		double	t1_cy,
		Sky1_Prec1976 *	terms
	)

This procedure calculates the equatorial precession parameters ζ, z, and θ which represent the rotation required to transform the FK5 equatorial reference system of Julian epoch t0 to that of Julian epoch t1, according to the IAU 1976 precession constants.

Parameters

[in]	t0_cy	Julian centuries since J2000.0 of initial epoch, TT timescale
[in]	t1_cy	Julian centuries since J2000.0 of final epoch, TT timescale
[out]	terms	The three precession angles ζ, z, θ.

References:: Lieske,J.H., Astron. Astrophys, Vol 73, pp282-284, 1979
Supplement to The Astronomical Almanac 1984, HMNAO, ppS18-19

When to call this function: It is quite likely that you will not need to call this function directly. If you are tracking a celestial object, using star_getApparent() or star_getTopocentric() or calling planet_getApparent() or planet_getTopocentric(), they will call this routine for you. Likewise, if you are tracking the object using the skyfast module, the call to skyfast_init() will call either star_getApparent() or planet_getApparent(), and therefore call this routine for you.

: The values calculated by this routine change only slowly. So if you are calling it yourself, you can call it infrequently. Intervals of up to an hour between calls will not introduce much error. This function is also called by sky1_createNPmatrix().

Definition at line 375 of file sky1.c.

◆ sky1_createPrec1976Matrix()

void sky1_createPrec1976Matrix	(	const Sky1_Prec1976 *	terms,
		V3D_Matrix *	precM
	)

This routine calculates the precession matrix, based on angles ζ, z, and θ.

precM is the combined orthogonal rotation matrix P, required for rigorous precession transformations using rectangular coordinates and matrix methods: V1 = P * V0 It appears to be calculated from
P = R3(-z) × R2(θ) × R3(-ζ)

Parameters

[in]	terms	The three precession angles ζ, z, θ, as returned by sky1_precessionIAU1976()
[out]	precM	Precession rotation matrix

References:: Lieske,J.H., Astron. Astrophys, Vol 73, pp282-284, 1979
Supplement to The Astronomical Almanac 1984, HMNAO, ppS18-19

When to call this function: Call this function after new values of ζ, z and θ have been calculated by routine sky1_precessionIAU1976(). So, again, it need only be called infrequently. But as mentioned for that function, it is quite likely that you will not need to call this function directly. Any of the functions sky1_createNPmatrix(), star_getApparent(), star_getTopocentric(), planet_getApparent() or planet_getTopocentric() will call this routine for you.

Definition at line 437 of file sky1.c.

◆ sky1_nutationIAU1980()

void sky1_nutationIAU1980	(	double	t_cy,
		int	precision,
		Sky1_Nut1980 *	nut
	)

Calculates the nutation in longitude and obliquity, according to the IAU 1980 Nutation Theory.

Calculates first the fundamental nutation arguments, and then a series of terms. There are 106 terms in the series, but fewer may be selected for faster execution.

Parameters

[in]	t_cy	Julian centuries since J2000.0, TT timescale
[in]	precision	How much precision do you want? Valid range [0, 4]. Values outside this range will be clamped to the range. 0 = full precision, use full 106-term series 1 = ignore terms < 0.2 milliarcseconds. 72-term series 2 = ignore terms < 0.3 milliarcseconds. 63-term series 3 = ignore terms < 0.5 milliarcseconds. 49-term series 4 = ignore terms < 1.0 milliarcseconds. 35-term series
[out]	nut	field nut->dPsi_rad - Nutation in longitude Δψ (radian) field nut->dEps_rad - Nutation in obliquity Δε (radian)

References:: Final report of the IAU Working Group on Nutation, chairman P.K.Seidelmann, 1980, published in Celestial Mechanics, Vol 27, pp79-106, 1982.
Kaplan,G.H. USNO circular no. 163, pp A3-6, 1981.
Supplement to the Astronomical Almanac 1984.
The form of the data tables was inspired by, and extended from, Reda, I. and Andreas, A. (2003), "Solar Position Algorithm for Solar Radiation Applications", National Renewable Energy Laboratory, NREL publication no. NREL/TP-560-34302.

Note

With precision = 0, exact agreement with tabulated values in the Astronomical Almanac 1987 pages B24 - B31 was found (to 3 decimal places, as tabulated in the Almanac). Likewise for the 1990 edition pages B24 - B31. But this IAU 1980 theory has been superseded by the more accurate IAU 2000 theory. This is evident when comparing outputs with the tabulated values in the 2007 Almanac, pages B34 - B41 (to 4 decimal places). There are differences of several milliarcseconds.
With precision = 1, differences up to 1.0 milliarcsecond were seen compared to the precision = 0 values during testing.
With precision = 2, differences up to 1.2 milliarcseconds were seen
With precision = 3, differences up to 1.5 milliarcseconds were seen
With precision = 4, differences of up to 5 milliarcseconds were seen.

For tracking purposes, these are still very small. The testing was not exhaustive, so take this as a rough guide only.

When to call this function: It is quite likely that you will not need to call this function directly. If you are tracking a celestial object, using star_getApparent() or star_getTopocentric() or calling planet_getApparent() or planet_getTopocentric(), they will call this routine for you. Likewise, if you are tracking the object using the skyfast module, the call to skyfast_init() will call either star_getApparent() or planet_getApparent(), and therefore call this routine for you.

: The values calculated by this routine change only slowly. So if you are calling it yourself, you can call it infrequently. Intervals of up to an hour between calls will not introduce much error. This function is also called by sky1_createNPmatrix().

Definition at line 481 of file sky1.c.

◆ sky1_epsilon1980()

void sky1_epsilon1980	(	double	t_cy,
		Sky1_Nut1980 *	nut
	)

Calculate the obliquity of the ecliptic and the equation of the equinoxes.

Parameters

[in]	t_cy	Julian centuries since J2000.0, TT timescale
[in,out]	nut	[in] field nut->dPsi_rad - Nutation in longitude Δψ, as returned by function sky1_nutationIAU1980() (radian) [in] field nut->dEps_rad - Nutation in obliquity Δε, as returned by function sky1_nutationIAU1980() (radian) [out] field nut->eps0_rad - Mean obliquity of the ecliptic ε0 (radian) [out] field nut->eqEq_rad - Equation of the equinoxes = Δψ * cos(ε0 + Δε) (radian) Note: not seconds

References: Astronomical Almanac 1987 page B18
International Astronomical Union, Standards of Fundamental Astronomy software collection, release 2017-04-20. http://www.iausofa.org routine iauObl80()

When to call this function: Call this function after new values of Δψ, and Δε have been calculated by routine sky1_nutationIAU1980(). So, again, it need only be called infrequently. But as mentioned for that function, it is quite likely that you will not need to call this function directly. Any of the functions sky1_createNPmatrix(), star_getApparent(), star_getTopocentric(), planet_getApparent() or planet_getTopocentric() will call this routine for you.

Definition at line 631 of file sky1.c.

◆ sky1_createNut1980Matrix()

void sky1_createNut1980Matrix	(	const Sky1_Nut1980 *	nut,
		V3D_Matrix *	nutM
	)

This routine calculates the Nutation matrix, using the nutation in longitude, the nutation in obliquity, and the mean obliquity of the equator.

We calculate the Nutation matrix in a completely rigorous manner from the product of three rotation matrices:
N = R1(-ε) × R3(-Δψ) × R1(ε0) where ε = ε0 + Δε (true obliquity of the equator)

Parameters

[in]	nut	Angles Δψ, Δε and ε0, as returned by functions sky1_nutationIAU1980() and sky1_epsilon1980()
[out]	nutM	Nutation rotation matrix N

References:: Astronomical Almanac 2007, page B31
also Mueller, p75

When to call this function: Call it after you have called sky1_nutationIAU1980() and sky1_epsilon1980(). As mentioned for the other functions above, this need only be done infrequently. But also as mentioned, it is quite likely that you will not need to call this function directly. Any of the functions sky1_createNPmatrix(), star_getApparent(), star_getTopocentric(), planet_getApparent() or planet_getTopocentric() will call this routine for you.

Definition at line 676 of file sky1.c.

◆ sky1_createNPmatrix()

void sky1_createNPmatrix	(	double	t0_cy,
		double	t1_cy,
		int	precision,
		V3D_Matrix *	npM
	)

Call the various precession and nutation routines in this module to create a combined precession and nutation rotation matrix, suitable for converting a celestial object's coordinates from mean place at epoch t0_cy to apparent coordinates at epoch t1_cy.

Parameters

[in]	t0_cy	Julian centuries since J2000.0 of initial epoch, TT timescale
[in]	t1_cy	Julian centuries since J2000.0 of final epoch, TT timescale
[in]	precision	Precision specification to be passed to routine sky1_nutationIAU1980(). See that routine for explanation
[out]	npM	Combined precession and nutation matrix NP = N × P

This function calls routines sky1_precessionIAU1976(), sky1_createPrec1976Matrix(), sky1_nutationIAU1980(), sky1_epsilon1980() and sky1_createNut1980Matrix() to calculate the precession and nutation matrices.

When to call this function: Call this function if you are tracking non-solar system objects. The resulting matrix changes only very slowly, so this function needs to be called only once per hour.

Definition at line 721 of file sky1.c.

◆ sky1_gmSiderealTimeIAU1982()

double sky1_gmSiderealTimeIAU1982 ( double du )

Calculate the Greenwich mean sidereal time.

Returns: Greenwich Mean Sidereal Time (radian)

Parameters

[in] du Days since J2000.0, UT1 timescale

Reference:: Astronomical Almanac 1987, page B6

When to call this function: If you are tracking a celestial object, you need not call this function directly. So long as you call sky1_appToTirs() every time around your control loop, it will call this function for you.

Definition at line 766 of file sky1.c.

◆ sky1_appToTirs()

void sky1_appToTirs	(	const V3D_Vector *	appV,
		double	j2kUT1_d,
		double	eqEq_rad,
		V3D_Vector *	terInterV
	)

Convert a position in geocentric apparent coordinates to geocentric coordinates in the Terrestrial Intermediate Reference System.

This is the first stage of converting apparent coordinates to topocentric coordinates. The resulting vector depends upon the current rotational position of the Earth. (For the second stage, to obtain topocentric coordinates, call routine sky_siteTirsToTopo()).

Parameters

[in]	appV	Position vector of apparent place (unit vector in equatorial coordinates)
[in]	j2kUT1_d	days since J2000.0, UT1 timescale, as returned by function sky_updateTimes() in the j2kUT1_d field of the Sky_Times struct.
[in]	eqEq_rad	Equation of the equinoxes (radian), as returned by function sky1_epsilon1980() in the eqEq_rad field of the Sky1_Nut1980 struct.
[out]	terInterV	Position vector in Terrestrial Intermediate Ref System

When to call this function: When you have the position of a celestial object expressed in Apparent coordinates, use this function to convert it to Terrestrial coordinates at the relevant rotational position of the earth. (Effectively you are converting from Right Ascension and Declination to the negative of the Greenwich Hour Angle and Declination.) If you are running a control loop to enable continuous tracking of this object, you will need to call this function (once) every time around your control loop.

: Follow this function with a call to sky_siteTirsToTopo() to obtain the object's position in topocentric coordinates at the observing site.

Definition at line 810 of file sky1.c.

Data Structures

Functions

Detailed Description

Function Documentation

◆ sky1_frameBiasFK5()

◆ sky1_precessionIAU1976()

◆ sky1_createPrec1976Matrix()

◆ sky1_nutationIAU1980()

◆ sky1_epsilon1980()

◆ sky1_createNut1980Matrix()

◆ sky1_createNPmatrix()

◆ sky1_gmSiderealTimeIAU1982()

◆ sky1_appToTirs()