mirror/astronomy
1
0
Fork 0
mirror of https://github.com/cosinekitty/astronomy.git synced 2026-01-03 13:10:43 -05:00
Code Issues Packages Projects Releases Wiki Activity
Files
python_precision_problem
astronomy/generate/novas/novas.c
Don Cross 30a9b0b757 Fixed warnings from gcc 12.2. 
I'm taking gcc 12.2 for a test drive today.
It reports a few warnings that slipped through earlier versions.
None of the warnings concern me for actual code safety,
but I went ahead and resolved them to keep the build clean.

Also provide a hook for a CPP environment variable to override
the C++ compiler to use, instead of forcing g++.
2022-12-08 13:41:12 -05:00

10199 lines
258 KiB
C
Raw Permalink Blame History

/*
Naval Observatory Vector Astrometry Software (NOVAS)
C Edition, Version 3.1
novas.c: Main library
U. S. Naval Observatory
Astronomical Applications Dept.
Washington, DC
http://www.usno.navy.mil/USNO/astronomical-applications
*/
#ifndef _NOVAS_
#include "novas.h"
#endif
#include <math.h>
/*
Global variables.
'PSI_COR' and 'EPS_COR' are celestial pole offsets for high-
precision applications. See function 'cel_pole' for more details.
*/
static double PSI_COR = 0.0;
static double EPS_COR = 0.0;
/********app_star */
short int app_star (double jd_tt, cat_entry *star, short int accuracy,
double *ra, double *dec)
/*
------------------------------------------------------------------------
PURPOSE:
Computes the apparent place of a star at date 'jd_tt', given its
catalog mean place, proper motion, parallax, and radial velocity.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Explanatory Supplement to the Astronomical Almanac (1992),
Chapter 3.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date for apparent place.
*star (struct cat_entry)
Pointer to catalog entry structure containing catalog data for
the object in the ICRS (defined in novas.h).
accuracy (short int)
Code specifying the relative accuracy of the output position.
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*ra (double)
Apparent right ascension in hours, referred to true equator
and equinox of date 'jd_tt'.
*dec (double)
Apparent declination in degrees, referred to true equator
and equinox of date 'jd_tt'.
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
> 10 ... Error code from function 'make_object'.
> 20 ... Error code from function 'place'.
GLOBALS
USED:
SIZE_OF_OBJ_NAME novas.h
FUNCTIONS
CALLED:
make_object novas.c
place novas.c
strcpy string.h
VER./DATE/
PROGRAMMER:
V1.0/12-05/JAB (USNO/AA)
V1.1/10-06/JAB (USNO/AA): Incorporate 'output' structure.
V1.2/10-08/JAB (USNO/AA): Add 'accuracy' option to input.
V1.3/05-10/JAB (USNO/AA): Fix bug in set-up of 'obj_name'.
NOTES:
1. This function is the "C" version of Fortran NOVAS routine
'apstar'.
2. SIZE_OF_OBJ_NAME is defined in novas.h
------------------------------------------------------------------------
*/
{
char obj_name[SIZE_OF_OBJ_NAME];
short int error = 0;
short int type, number, coord_sys;
double delta_t = 0.0;
object cel_obj;
observer location;
sky_pos output;
*ra = *dec = 0.0; /* [Don Cross]: fix warnings */
/*
Set 'obj_name' equal to 'starname' in the 'star' structure. Length
will be checked in 'make_object'.
*/
strcpy (obj_name, star->starname);
/*
Set up a structure of type 'object' containing the star data.
*/
type = 2; /* Object located outside the solar system */
number = 0; /* Set to zero */
if ((error = make_object (type,number,obj_name,star, &cel_obj)) != 0)
{
error += 10;
return (error);
}
/*
Compute the apparent place with a call to function 'place'.
*/
location.where = 0; /* Geocenter */
coord_sys = 1; /* True equator and equinox of date */
if ((error = place (jd_tt,&cel_obj,&location,delta_t,coord_sys,
accuracy, &output)) != 0)
{
*ra = 0.0;
*dec = 0.0;
error += 20;
}
else
{
*ra = output.ra;
*dec = output.dec;
}
return (error);
}
/********virtual_star */
short int virtual_star (double jd_tt, cat_entry *star,
short int accuracy,
double *ra, double *dec)
/*
------------------------------------------------------------------------
PURPOSE:
Computes the virtual place of a star at date 'jd_tt', given its
catalog mean place, proper motion, parallax, and radial velocity.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Explanatory Supplement to the Astronomical Almanac (1992),
Chapter 3.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date for virtual place.
*star (struct cat_entry)
Pointer to catalog entry structure containing catalog data for
the object in the ICRS (defined in novas.h).
accuracy (short int)
Code specifying the relative accuracy of the output position.
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*ra (double)
Virtual right ascension in hours, referred to the GCRS.
*dec (double)
Virtual declination in degrees, referred to the GCRS.
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
> 10 ... Error code from function 'make_object'.
> 20 ... Error code from function 'place'.
GLOBALS
USED:
SIZE_OF_OBJ_NAME novas.h
FUNCTIONS
CALLED:
make_object novas.c
place novas.c
strcpy string.h
VER./DATE/
PROGRAMMER:
V1.0/12-05/JAB (USNO/AA)
V1.1/10-06/JAB (USNO/AA): Incorporate 'output' structure.
V1.2/10-08/JAB (USNO/AA): Add 'accuracy' option to input.
V1.3/05-10/JAB (USNO/AA): Fix bug in set-up of 'obj_name'.
NOTES:
1. This function is the "C" version of Fortran NOVAS routine
'vpstar'.
2. SIZE_OF_OBJ_NAME is defined in novas.h
------------------------------------------------------------------------
*/
{
char obj_name[SIZE_OF_OBJ_NAME];
short int error = 0;
short int type, number, coord_sys;
double delta_t = 0.0;
object cel_obj;
observer location;
sky_pos output;
/*
Set 'obj_name' equal to 'starname' in the 'star' structure. Length
will be checked in 'make_object'.
*/
strcpy (obj_name, star->starname);
/*
Set up a structure of type 'object' containing the star data.
*/
type = 2; /* Object located outside the solar system */
number = 0; /* Set to zero */
if ((error = make_object (type,number,obj_name,star, &cel_obj)) != 0)
{
error += 10;
return (error);
}
/*
Compute the virtual place with a call to function 'place'.
*/
location.where = 0; /* Geocenter */
coord_sys = 0; /* GCRS */
if ((error = place (jd_tt,&cel_obj,&location,delta_t,coord_sys,
accuracy, &output)) != 0)
{
*ra = 0.0;
*dec = 0.0;
error += 20;
}
else
{
*ra = output.ra;
*dec = output.dec;
}
return (error);
}
/********astro_star */
short int astro_star (double jd_tt, cat_entry *star, short int accuracy,
double *ra, double *dec)
/*
------------------------------------------------------------------------
PURPOSE:
Computes the astrometric place of a star at date 'jd_tt', given
its catalog mean place, proper motion, parallax, and radial
velocity.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Explanatory Supplement to the Astronomical Almanac (1992),
Chapter 3.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date for astrometric place.
*star (struct cat_entry)
Pointer to catalog entry structure containing catalog data for
the object in the ICRS (defined in novas.h).
accuracy (short int)
Code specifying the relative accuracy of the output position.
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*ra (double)
Astrometric right ascension in hours (referred to the ICRS,
without light deflection or aberration).
*dec (double)
Astrometric declination in degrees (referred to the ICRS,
without light deflection or aberration).
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
> 10 ... Error code from function 'make_object'.
> 20 ... Error code from function 'place'.
GLOBALS
USED:
SIZE_OF_OBJ_NAME novas.h
FUNCTIONS
CALLED:
make_object novas.c
place novas.c
strcpy string.h
VER./DATE/
PROGRAMMER:
V1.0/12-05/JAB (USNO/AA)
V1.1/10-06/JAB (USNO/AA): Incorporate 'output' structure.
V1.2/10-08/JAB (USNO/AA): Add 'accuracy' option to input.
V1.3/05-10/JAB (USNO/AA): Fix bug in set-up of 'obj_name'.
NOTES:
1. This function is the "C" version of Fortran NOVAS routine
'asstar'.
2. SIZE_OF_OBJ_NAME is defined in novas.h
------------------------------------------------------------------------
*/
{
char obj_name[SIZE_OF_OBJ_NAME];
short int error = 0;
short int type, number, coord_sys;
double delta_t = 0.0;
object cel_obj;
observer location;
sky_pos output;
/*
Set 'obj_name' equal to 'starname' in the 'star' structure. Length
will be checked in 'make_object'.
*/
strcpy (obj_name, star->starname);
/*
Set up a structure of type 'object' containing the star data.
*/
type = 2; /* Object located outside the solar system */
number = 0; /* Set to zero */
if ((error = make_object (type,number,obj_name,star, &cel_obj)) != 0)
{
error += 10;
return (error);
}
/*
Compute the astrometric place with a call to function 'place'.
*/
location.where = 0; /* Geocenter */
coord_sys = 3; /* ICRS astrometric coordinates */
if ((error = place (jd_tt,&cel_obj,&location,delta_t,coord_sys,
accuracy, &output)) != 0)
{
*ra = 0.0;
*dec = 0.0;
error += 20;
}
else
{
*ra = output.ra;
*dec = output.dec;
}
return (error);
}
/********app_planet */
short int app_planet (double jd_tt, object *ss_body, short int accuracy,
double *ra, double *dec, double *dis)
/*
------------------------------------------------------------------------
PURPOSE:
Compute the apparent place of a solar system body.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Explanatory Supplement to the Astronomical Almanac (1992),
Chapter 3.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date for apparent place.
*ss_body (struct object)
Pointer to structure containing the body designation for the
solar system body (defined in novas.h).
accuracy (short int)
Code specifying the relative accuracy of the output position.
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*ra (double)
Apparent right ascension in hours, referred to true equator
and equinox of date.
*dec (double)
Apparent declination in degrees, referred to true equator
and equinox of date.
*dis (double)
True distance from Earth to the body at 'jd_tt' in AU.
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
= 1 ... Invalid value of 'type' in structure 'ss_body'.
> 10 ... Error code from function 'place'.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
place novas.c
VER./DATE/
PROGRAMMER:
V1.0/12-05/JAB (USNO/AA)
V1.1/10-06/JAB (USNO/AA): Incorporate 'output' structure.
V1.2/10-08/JAB (USNO/AA): Add 'accuracy' option to input.
NOTES:
1. This function is the "C" version of Fortran NOVAS routine
'applan'.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int coord_sys;
double delta_t = 0.0;
observer location;
sky_pos output;
/*
Check for a valid value of 'type' in structure 'ss_body'.
*/
if ((ss_body->type < 0) || (ss_body->type > 1))
{
error = 1;
return (error);
}
/*
Compute the apparent place with a call to function 'place'.
*/
location.where = 0; /* Geocenter */
coord_sys = 1; /* True equator and equinox of date */
if ((error = place (jd_tt,ss_body,&location,delta_t,coord_sys,
accuracy, &output)) != 0)
{
*ra = 0.0;
*dec = 0.0;
*dis = 0.0;
error += 10;
}
else
{
*ra = output.ra;
*dec = output.dec;
*dis = output.dis;
}
return (error);
}
/********virtual_planet */
short int virtual_planet (double jd_tt, object *ss_body,
short int accuracy,
double *ra, double *dec, double *dis)
/*
------------------------------------------------------------------------
PURPOSE:
Compute the virtual place of a solar system body.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Explanatory Supplement to the Astronomical Almanac (1992),
Chapter 3.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date for virtual place.
*ss_body (struct object)
Pointer to structure containing the body designation for the
solar system body (defined in novas.h).
accuracy (short int)
Code specifying the relative accuracy of the output position.
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*ra (double)
Virtual right ascension in hours, referred to the GCRS.
*dec (double)
Virtual declination in degrees, referred to the GCRS.
*dis (double)
True distance from Earth to the body in AU.
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
= 1 ... Invalid value of 'type' in structure 'ss_body'.
> 10 ... Error code from function 'place'.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
place novas.c
VER./DATE/
PROGRAMMER:
V1.0/12-05/JAB (USNO/AA)
V1.1/10-06/JAB (USNO/AA): Incorporate 'output' structure.
V1.2/10-08/JAB (USNO/AA): Add 'accuracy' option to input.
NOTES:
1. This function is the "C" version of Fortran NOVAS routine
'vpplan'.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int coord_sys;
double delta_t = 0.0;
observer location;
sky_pos output;
/*
Check for a valid value of 'type' in structure 'ss_body'.
*/
if ((ss_body->type < 0) || (ss_body->type > 1))
{
error = 1;
return (error);
}
/*
Compute the virtual place with a call to function 'place'.
*/
location.where = 0; /* Geocenter */
coord_sys = 0; /* GCRS */
if ((error = place (jd_tt,ss_body,&location,delta_t,coord_sys,
accuracy, &output)) != 0)
{
*ra = 0.0;
*dec = 0.0;
*dis = 0.0;
error += 10;
}
else
{
*ra = output.ra;
*dec = output.dec;
*dis = output.dis;
}
return (error);
}
/********astro_planet */
short int astro_planet (double jd_tt, object *ss_body,
short int accuracy,
double *ra, double *dec, double *dis)
/*
------------------------------------------------------------------------
PURPOSE:
Compute the astrometric place of a solar system body.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Explanatory Supplement to the Astronomical Almanac (1992),
Chapter 3.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date for astrometric place.
*ss_body (struct object)
Pointer to structure containing the body designation for the
solar system body (defined in novas.h).
accuracy (short int)
Code specifying the relative accuracy of the output position.
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*ra (double)
Astrometric right ascension in hours (referred to the ICRS,
without light deflection or aberration).
*dec (double)
Astrometric declination in degrees (referred to the ICRS,
without light deflection or aberration).
*dis (double)
True distance from Earth to the body in AU.
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
= 1 ... Invalid value of 'type' in structure 'ss_body'.
> 10 ... Error code from function 'place'.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
place novas.c
VER./DATE/
PROGRAMMER:
V1.0/12-05/JAB (USNO/AA)
V1.1/10-06/JAB (USNO/AA): Incorporate 'output' structure.
V1.2/10-08/JAB (USNO/AA): Add 'accuracy' option to input.
NOTES:
1. This function is the "C" version of Fortran NOVAS routine
'asplan'.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int coord_sys;
double delta_t = 0.0;
observer location;
sky_pos output;
/*
Check for a valid value of 'type' in structure 'ss_body'.
*/
if ((ss_body->type < 0) || (ss_body->type > 1))
{
error = 1;
return (error);
}
/*
Compute the astrometric place with a call to function 'place'.
*/
location.where = 0; /* Geocenter */
coord_sys = 3; /* ICRS astrometric coordinates */
if ((error = place (jd_tt,ss_body,&location,delta_t,coord_sys,
accuracy, &output)) != 0)
{
*ra = 0.0;
*dec = 0.0;
*dis = 0.0;
error += 10;
}
else
{
*ra = output.ra;
*dec = output.dec;
*dis = output.dis;
}
return (error);
}
/********topo_star */
short int topo_star (double jd_tt, double delta_t, cat_entry *star,
on_surface *position, short int accuracy,
double *ra, double *dec)
/*
------------------------------------------------------------------------
PURPOSE:
Computes the topocentric place of a star at date 'jd_tt', given its
catalog mean place, proper motion, parallax, and radial velocity.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Explanatory Supplement to the Astronomical Almanac (1992),
Chapter 3.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date for topocentric place.
delta_t (double)
Difference TT-UT1 at 'jd_tt', in seconds of time.
*star (struct cat_entry)
Pointer to catalog entry structure containing catalog data for
the object in the ICRS (defined in novas.h).
*position (struct on_surface)
Specifies the position of the observer (structure defined in
novas.h).
accuracy (short int)
Code specifying the relative accuracy of the output position.
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*ra (double)
Topocentric right ascension in hours, referred to true equator
and equinox of date 'jd_tt'.
*dec (double)
Topocentric declination in degrees, referred to true equator
and equinox of date 'jd_tt'.
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
= 1 ... Invalid value of 'where' in structure 'location'.
> 10 ... Error code from function 'make_object'.
> 20 ... Error code from function 'place'.
GLOBALS
USED:
SIZE_OF_OBJ_NAME novas.h
FUNCTIONS
CALLED:
make_observer novas.c
make_object novas.c
place novas.c
strcpy string.h
VER./DATE/
PROGRAMMER:
V1.0/01-06/JAB (USNO/AA)
V1.1/01-06/WKP (USNO/AA): Fixed minor syntax problem.
V1.2/10-06/JAB (USNO/AA): Incorporate 'output' structure.
V1.3/10-08/JAB (USNO/AA): Add 'accuracy' option to input.
V1.4/05-10/JAB (USNO/AA): Fix bug in set-up of 'obj_name'.
NOTES:
1. This function is the "C" version of Fortran NOVAS routine
'tpstar'.
2. SIZE_OF_OBJ_NAME is defined in novas.h
------------------------------------------------------------------------
*/
{
char obj_name[SIZE_OF_OBJ_NAME];
short int error = 0;
short int type, number, coord_sys;
in_space dummy;
observer location;
object cel_obj;
sky_pos output;
/*
Initialize the 'dummy' structure.
*/
dummy.sc_pos[0] = dummy.sc_pos[1] = dummy.sc_pos[2] = 0.0;
dummy.sc_vel[0] = dummy.sc_vel[1] = dummy.sc_vel[2] = 0.0;
/*
Set up a structure of type 'observer' containing the position
of the observer.
*/
if ((error = make_observer (1,position,&dummy, &location)) != 0)
{
error = 1;
return (error);
}
/*
Set 'obj_name' equal to 'starname' in the 'star' structure. Length
will be checked in 'make_object'.
*/
strcpy (obj_name, star->starname);
/*
Set up a structure of type 'object' containing the star data.
*/
type = 2; /* Object located outside the solar system */
number = 0; /* Set to zero */
if ((error = make_object (type,number,obj_name,star, &cel_obj)) != 0)
{
error += 10;
return (error);
}
/*
Compute the topocentric place with a call to function 'place'.
*/
coord_sys = 1; /* True equator and equinox of date */
if ((error = place (jd_tt,&cel_obj,&location,delta_t,coord_sys,
accuracy, &output)) != 0)
{
*ra = 0.0;
*dec = 0.0;
error += 20;
}
else
{
*ra = output.ra;
*dec = output.dec;
}
return (error);
}
/********local_star */
short int local_star (double jd_tt, double delta_t, cat_entry *star,
on_surface *position, short int accuracy,
double *ra, double *dec)
/*
------------------------------------------------------------------------
PURPOSE:
Computes the local place of a star at date 'jd_tt', given its
catalog mean place, proper motion, parallax, and radial velocity.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Explanatory Supplement to the Astronomical Almanac (1992),
Chapter 3.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date for local place.
delta_t (double)
Difference TT-UT1 at 'jd_tt', in seconds of time.
*star (struct cat_entry)
Pointer to catalog entry structure containing catalog data for
the object in the ICRS (defined in novas.h).
*position (struct on_surface)
Specifies the position of the observer (structure defined in
novas.h).
accuracy (short int)
Code specifying the relative accuracy of the output position.
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*ra (double)
Local right ascension in hours, referred to the 'local GCRS'.
*dec (double)
Local declination in degrees, referred to the 'local GCRS'.
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
= 1 ... Invalid value of 'where' in structure 'location'.
> 10 ... Error code from function 'make_object'.
> 20 ... Error code from function 'place'.
GLOBALS
USED:
SIZE_OF_OBJ_NAME novas.h
FUNCTIONS
CALLED:
make_observer novas.c
make_object novas.c
place novas.c
strcpy string.h
VER./DATE/
PROGRAMMER:
V1.0/01-06/JAB (USNO/AA)
V1.1/10-06/JAB (USNO/AA): Incorporate 'output' structure.
V1.2/10-08/JAB (USNO/AA): Add 'accuracy' option to input.
V1.3/05-10/JAB (USNO/AA): Fix bug in set-up of 'obj_name'.
NOTES:
1. This function is the "C" version of Fortran NOVAS routine
'lpstar'.
2. SIZE_OF_OBJ_NAME is defined in novas.h
------------------------------------------------------------------------
*/
{
char obj_name[SIZE_OF_OBJ_NAME];
short int error = 0;
short int type, number, coord_sys;
in_space dummy;
observer location;
object cel_obj;
sky_pos output;
/*
Initialize the 'dummy' structure.
*/
dummy.sc_pos[0] = dummy.sc_pos[1] = dummy.sc_pos[2] = 0.0;
dummy.sc_vel[0] = dummy.sc_vel[1] = dummy.sc_vel[2] = 0.0;
/*
Set up a structure of type 'observer' containing the position
of the observer.
*/
if ((error = make_observer (1,position,&dummy, &location)) != 0)
{
error = 1;
return (error);
}
/*
Set 'obj_name' equal to 'starname' in the 'star' structure. Length
will be checked in 'make_object'.
*/
strcpy (obj_name, star->starname);
/*
Set up a structure of type 'object' containing the star data.
*/
type = 2; /* Object located outside the solar system */
number = 0; /* Set to zero */
if ((error = make_object (type,number,obj_name,star, &cel_obj)) != 0)
{
error += 10;
return (error);
}
/*
Compute the local place with a call to function 'place'.
*/
coord_sys = 0; /* "Local GCRS" */
if ((error = place (jd_tt,&cel_obj,&location,delta_t,coord_sys,
accuracy, &output)) != 0)
{
*ra = 0.0;
*dec = 0.0;
error += 20;
}
else
{
*ra = output.ra;
*dec = output.dec;
}
return (error);
}
/********topo_planet */
short int topo_planet (double jd_tt, object *ss_body, double delta_t,
on_surface *position, short int accuracy,
double *ra, double *dec, double *dis)
/*
------------------------------------------------------------------------
PURPOSE:
Computes the topocentric place of a solar system body.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Explanatory Supplement to the Astronomical Almanac (1992),
Chapter 3.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date for topocentric place.
*ss_body (struct object)
Pointer to structure containing the body designation for the
solar system body (defined in novas.h).
delta_t (double)
Difference TT-UT1 at 'jd_tt', in seconds of time.
*position (struct on_surface)
Specifies the position of the observer (structure defined in
novas.h).
accuracy (short int)
Code specifying the relative accuracy of the output position.
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*ra (double)
Topocentric right ascension in hours, referred to true equator
and equinox of date.
*dec (double)
Topocentric declination in degrees, referred to true equator
and equinox of date.
*dis (double)
True distance from Earth to the body at 'jd_tt' in AU.
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
= 1 ... Invalid value of 'where' in structure 'location'.
> 10 ... Error code from function 'place'.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
make_observer novas.c
place novas.c
VER./DATE/
PROGRAMMER:
V1.0/01-06/JAB (USNO/AA)
V1.1/01-06/WKP (USNO/AA): Fixed minor syntax problem.
V1.2/10-06/JAB (USNO/AA): Incorporate 'output' structure.
V1.3/10-08/JAB (USNO/AA): Add 'accuracy' option to input.
V1.4/07-10/JLB (USNO/AA): Fixed minor documentation error in prolog
NOTES:
1. This function is the "C" version of Fortran NOVAS routine
'tpplan'.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int coord_sys;
in_space dummy;
observer location;
sky_pos output;
/*
Initialize the 'dummy' structure.
*/
dummy.sc_pos[0] = dummy.sc_pos[1] = dummy.sc_pos[2] = 0.0;
dummy.sc_vel[0] = dummy.sc_vel[1] = dummy.sc_vel[2] = 0.0;
/*
Set up a structure of type 'observer' containing the position
of the observer.
*/
if ((error = make_observer (1,position,&dummy, &location)) != 0)
{
error = 1;
return (error);
}
/*
Compute the topocentric place with a call to function 'place'.
*/
coord_sys = 1; /* True equator and equinox of date */
if ((error = place (jd_tt,ss_body,&location,delta_t,coord_sys,
accuracy, &output)) != 0)
{
*ra = 0.0;
*dec = 0.0;
*dis = 0.0;
error += 10;
}
else
{
*ra = output.ra;
*dec = output.dec;
*dis = output.dis;
}
return (error);
}
/********local_planet */
short int local_planet (double jd_tt, object *ss_body,
double delta_t, on_surface *position,
short int accuracy,
double *ra, double *dec, double *dis)
/*
------------------------------------------------------------------------
PURPOSE:
Computes the local place of a solar system body.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Explanatory Supplement to the Astronomical Almanac (1992),
Chapter 3.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date for local place.
*ss_body (struct object)
Pointer to structure containing the body designation for the
solar system body (defined in novas.h).
delta_t (double)
Difference TT-UT1 at 'jd_tt', in seconds of time.
*position (struct on_surface)
Specifies the position of the observer (structure defined in
novas.h).
accuracy (short int)
Code specifying the relative accuracy of the output position.
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*ra (double)
Local right ascension in hours, referred to the 'local GCRS'.
*dec (double)
Local declination in degrees, referred to the 'local GCRS'.
*dis (double)
True distance from Earth to the body in AU.
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
= 1 ... Invalid value of 'where' in structure 'location'.
> 10 ... Error code from function 'place'.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
make_observer novas.c
place novas.c
VER./DATE/
PROGRAMMER:
V1.0/01-06/JAB (USNO/AA)
V1.1/01-06/WKP (USNO/AA): Fixed minor syntax problem.
V1.2/10-06/JAB (USNO/AA): Incorporate 'output' structure.
V1.3/10-08/JAB (USNO/AA): Add 'accuracy' option to input.
NOTES:
1. This function is the "C" version of Fortran NOVAS routine
'lpplan'.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int coord_sys;
in_space dummy;
observer location;
sky_pos output;
/*
Initialize the 'dummy' structure.
*/
dummy.sc_pos[0] = dummy.sc_pos[1] = dummy.sc_pos[2] = 0.0;
dummy.sc_vel[0] = dummy.sc_vel[1] = dummy.sc_vel[2] = 0.0;
/*
Set up a structure of type 'observer' containing the position
of the observer.
*/
if ((error = make_observer (1,position,&dummy, &location)) != 0)
{
error = 1;
return (error);
}
/*
Compute the local place with a call to function 'place'.
*/
coord_sys = 0; /* "Local GCRS" */
if ((error = place (jd_tt,ss_body,&location,delta_t,coord_sys,
accuracy, &output)) != 0)
{
*ra = 0.0;
*dec = 0.0;
*dis = 0.0;
error += 10;
}
else
{
*ra = output.ra;
*dec = output.dec;
*dis = output.dis;
}
return (error);
}
/********mean_star */
short int mean_star (double jd_tt, double ra, double dec,
short int accuracy,
double *ira, double *idec)
/*
------------------------------------------------------------------------
PURPOSE:
Computes the ICRS position of a star, given its apparent place
at date 'jd_tt'. Proper motion, parallax and radial velocity
are assumed to be zero.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Explanatory Supplement to the Astronomical Almanac (1992),
Chapter 3.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date of apparent place.
ra (double)
Apparent right ascension in hours, referred to true equator
and equinox of date.
dec (double)
Apparent declination in degrees, referred to true equator
and equinox of date.
accuracy (short int)
Code specifying the relative accuracy of the output position.
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*ira (double)
ICRS right ascension in hours.
*idec (double)
ICRS declination in degrees.
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
= 1 ... Iterative process did not converge after 30 iterations.
= 2 ... length of dummy star name out of bounds.
= 3 ... length of dummy catalog name out of bounds.
> 10 ... Error from function 'vector2radec'.
> 20 ... Error from function 'app_star'.
GLOBALS
USED:
T0 novascon.c
FUNCTIONS
CALLED:
make_cat_entry novas.c
starvectors novas.c
precession novas.c
app_star novas.c
vector2radec novas.c
fabs math.h
VER./DATE/
PROGRAMMER:
V1.0/12-05/JAB (USNO/AA)
V1.1/03-08/WKP (USNO/AA): Fixed subtle bug in exit condition of
iteration and updated variable names.
V1.2/10-08/JAB (USNO/AA): Add 'accuracy' option to input.
V1.3/05-10/JAB (USNO/AA): Use 'make_cat_entry' to set up
local structure 'tempstar'.
NOTES:
1. This function is the "C" version of Fortran NOVAS routine
'mpstar'.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int iter = 0;
double pos[3], dum[3], pos2[3], newira, newidec, oldira, oldidec, ra2,
dec2, deltara, deltadec;
cat_entry tempstar;
/*
Set up the 'tempstar' structure, then use it to create a position
vector based on the apparent RA and declination of the star.
*/
if ((error = make_cat_entry ("dummy","CAT",0,ra,dec,0.0,0.0,0.0,0.0,
&tempstar)) != 0)
{
return (error + 1);
}
starvectors (&tempstar, pos,dum);
/*
Get initial approximation by precessing star position at 'jd_tt'
to its position at J2000.0.
*/
precession (jd_tt,pos,T0, pos2);
if ((error = vector2radec (pos2, &newira,&newidec)) != 0)
{
return (error + 10);
}
/*
Iteratively find ICRS coordinates that produce input apparent place
of star at date 'jd_tt'.
*/
do
{
oldira = newira;
oldidec = newidec;
tempstar.ra = oldira;
tempstar.dec = oldidec;
if ((error = app_star (jd_tt,&tempstar,accuracy, &ra2,&dec2)) != 0)
{
*ira = 0.0;
*idec = 0.0;
return (error + 20);
}
deltara = ra2 - oldira;
deltadec = dec2 - oldidec;
if (deltara < -12.0)
deltara += 24.0;
if (deltara > 12.0)
deltara -= 24.0;
newira = ra - deltara;
newidec = dec - deltadec;
if (iter >= 30)
{
*ira = 0.0;
*idec = 0.0;
return (error = 1);
}
else
{
iter++;
}
} while (!(fabs (newira - oldira) <= 1.0e-12) ||
!(fabs (newidec - oldidec) <= 1.0e-11));
*ira = newira;
*idec = newidec;
if (*ira < 0.0)
*ira += 24.0;
if (*ira >= 24.0)
*ira -= 24.0;
return (error);
}
/********place */
short int place (double jd_tt, object *cel_object,
observer *location, double delta_t,
short int coord_sys, short int accuracy,
sky_pos *output)
/*
------------------------------------------------------------------------
PURPOSE:
This function computes the apparent direction of a star or solar
system body at a specified time and in a specified coordinate
system.
REFERENCES:
Kaplan, G. et al. (1989), Astronomical Journal 97, 1197-1210.
Klioner, S. (2003), Astronomical Journal 125, 1580-1597.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date for place.
*cel_object (struct object)
Specifies the celestial object of interest (structure defined
in novas.h).
*location (struct observer)
Specifies the location of the observer (structure defined in
novas.h).
delta_t (double)
Difference TT-UT1 at 'jd_tt', in seconds of time.
coord_sys (short int)
Code specifying coordinate system of the output position.
= 0 ... GCRS or "local GCRS"
= 1 ... true equator and equinox of date
= 2 ... true equator and CIO of date
= 3 ... astrometric coordinates, i.e., without light
deflection or aberration.
accuracy (short int)
Code specifying the relative accuracy of the output position.
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*output (struct sky_pos)
Output data specifying object's place on the sky at time
'jd_tt', with respect to the specified output coordinate system
(struct defined in novas.h).
RETURNED
VALUE:
= 0 ... No problems.
= 1 ... invalid value of 'coord_sys'
= 2 ... invalid value of 'accuracy'
= 3 ... Earth is the observed object, and the observer is
either at the geocenter or on the Earth's surface
(not permitted)
> 10, < 40 ... 10 + error from function 'ephemeris'
> 40, < 50 ... 40 + error from function 'geo_posvel'
> 50, < 70 ... 50 + error from function 'light_time'
> 70, < 80 ... 70 + error from function 'grav_def'
> 80, < 90 ... 80 + error from function 'cio_location'
> 90, < 100 ... 90 + error from function 'cio_basis'
GLOBALS
USED:
T0, C_AUDAY novascon.c
FUNCTIONS
CALLED:
make_cat_entry novas.c
make_object novas.c
tdb2tt novas.c
ephemeris novas.c
geo_posvel novas.c
starvectors novas.c
d_light novas.c
proper_motion novas.c
bary2obs novas.c
light_time novas.c
limb_angle novas.c
grav_def novas.c
aberration novas.c
frame_tie novas.c
precession novas.c
nutation novas.c
cio_location novas.c
cio_basis novas.c
rad_vel novas.c
vector2radec novas.c
fabs math.h
sqrt math.h
VER./DATE/
PROGRAMMER:
V1.0/01-05/JAB (USNO/AA)
V1.1/09-05/WKP (USNO/AA) Removed 'star' from the input argument
list since it is nested in the 'object'
structure.
V1.2/01-06/WKP (USNO/AA) Updated error codes and replaced 'mode'
with 'accuracy' where appropriate.
V1.3/01-06/JAB (USNO/AA) Added more static variables.
V1.4/07-06/JAB (USNO/AA) Implement 'cio_location' construct.
V1.5/10-06/JAB (USNO/AA) Add radial velocity calculation and
implement 'sky_pos' structure.
V1.6/01-07/JAB (USNO/AA) Update to accommodate new radial velocity
algorithm.
V1.7/10-08/JAB (USNO/AA) Modify calls to 'ephemeris' to support
two-part input Julian date.
V1.8/07-10/JLB (USNO/AA) Corrected citation to Kaplan et al.
NOTES:
1. Values of 'location->where' and 'coord_sys' dictate the various
standard kinds of place:
location->where = 0 and coord_sys = 1: apparent place
location->where = 1 and coord_sys = 1: topocentric place
location->where = 0 and coord_sys = 0: virtual place
location->where = 1 and coord_sys = 0: local place
location->where = 0 and coord_sys = 3: astrometric place
location->where = 1 and coord_sys = 3: topocentric astrometric
place
2. Input value of 'delta_t' is used only when 'location->where'
equals 1 or 2 (observer is on surface of Earth or in a
near-Earth satellite).
------------------------------------------------------------------------
*/
{
static short int first_time = 1;
short int error = 0;
short int loc, rs, i;
static double tlast1 = 0.0;
static double tlast2 = 0.0;
static double jd_tdb, peb[3], veb[3], psb[3], vsb[3], px[3], py[3],
pz[3];
double x, secdif, jd[2], pog[3], vog[3], pob[3], vob[3], pos1[3],
vel1[3], dt, pos2[3], pos3[3], t_light, t_light0, pos4[3], frlimb,
pos5[3], pos6[3], pos7[3], pos8[3], r_cio, d_obs_geo, d_obs_sun,
d_obj_sun;
cat_entry null_star;
static object earth, sun;
/*
Check for invalid value of 'coord_sys' or 'accuracy'.
*/
if ((coord_sys < 0) || (coord_sys > 3))
return (error = 1);
if ((accuracy < 0) || (accuracy > 1))
return (error = 2);
/*
Create a null star 'cat_entry' and Earth and Sun 'object's.
*/
if (first_time)
{
make_cat_entry ("NULL_STAR"," ",0L,0.0,0.0,0.0,0.0,0.0,0.0,
&null_star);
make_object (0,3,"Earth",&null_star, &earth);
make_object (0,10,"Sun",&null_star, &sun);
first_time = 0;
}
/*
---------------------------------------------------------------------
Check on Earth as an observed object. Earth can only be an observed
object when 'location' is a near-Earth satellite.
---------------------------------------------------------------------
*/
if ((cel_object->type == 0) && (cel_object->number == 3) &&
(location->where != 2))
return (error = 3);
/*
---------------------------------------------------------------------
Get position and velocity of Earth (geocenter) and Sun.
---------------------------------------------------------------------
*/
if (fabs (jd_tt - tlast1) > 1.0e-8)
{
/*
Compute 'jd_tdb', the TDB Julian date corresponding to 'jd_tt'.
*/
jd_tdb = jd_tt;
tdb2tt (jd_tdb, &x,&secdif);
jd_tdb = jd_tt + secdif / 86400.0;
/*
Get position and velocity of Earth wrt barycenter of solar system,
in ICRS.
*/
jd[0] = jd_tdb;
jd[1] = 0.0;
if ((error = ephemeris (jd,&earth,0,accuracy, peb,veb)) != 0)
return (error += 10);
/*
Get position and velocity of Sun wrt barycenter of solar system,
in ICRS.
*/
if ((error = ephemeris (jd,&sun,0,accuracy, psb,vsb)) != 0)
return (error += 10);
tlast1 = jd_tt;
}
/*
---------------------------------------------------------------------
Get position and velocity of observer.
---------------------------------------------------------------------
*/
if ((location->where == 1) || (location->where == 2))
{
/*
For topocentric place, get geocentric position and velocity vectors
of observer (observer is on surface of Earth or in a near-Earth
satellite).
*/
if ((error = geo_posvel (jd_tt,delta_t,accuracy,location, pog,vog))
!= 0)
return (error += 40);
loc = 1;
}
else
{
/*
For geocentric place, there is nothing to do (observer is at
geocenter).
*/
for (i = 0; i < 3; i++)
{
pog[i] = 0.0;
vog[i] = 0.0;
}
loc = 0;
}
/*
Compute position and velocity of observer wrt barycenter of
solar system (Galilean transformation).
*/
for (i = 0; i < 3; i++)
{
pob[i] = peb[i] + pog[i];
vob[i] = veb[i] + vog[i];
}
/*
---------------------------------------------------------------------
Find geometric position of observed object.
---------------------------------------------------------------------
*/
if (cel_object->type == 2) /* Observed object is star. */
{
/*
Get position of star updated for its space motion.
*/
starvectors (&cel_object->star, pos1,vel1);
dt = d_light (pos1,pob);
proper_motion (T0,pos1,vel1,(jd_tdb + dt), pos2);
/*
Get position of star wrt observer (corrected for parallax).
*/
bary2obs (pos2,pob, pos3,&t_light);
output->dis = 0.0;
}
else /* Observed object is solar system body. */
{
/*
Get position of body wrt barycenter of solar system.
*/
jd[0] = jd_tdb;
jd[1] = 0.0;
if ((error = ephemeris (jd,cel_object,0,accuracy, pos1,vel1)) != 0)
return (error += 10);
/*
Get position of body wrt observer, and true (Euclidian) distance.
*/
bary2obs (pos1,pob, pos2,&t_light0);
output->dis = t_light0 * C_AUDAY;
/*
Get position of body wrt observer, antedated for light-time.
*/
if ((error = light_time (jd_tdb,cel_object,pob,t_light0,accuracy,
pos3,&t_light)) != 0)
return (error += 50);
}
/*
---------------------------------------------------------------------
Apply gravitational deflection of light and aberration.
---------------------------------------------------------------------
*/
if (coord_sys == 3)
{
/*
These calculations are skipped for astrometric place.
*/
for (i = 0; i < 3; i++)
{
pos5[i] = pos3[i];
}
}
else
{
/*
Variable 'loc' determines whether Earth deflection is included.
*/
if (loc == 1)
{
limb_angle (pos3,pog, &x,&frlimb);
if (frlimb < 0.8)
loc = 0;
}
/*
Compute gravitational deflection and aberration.
*/
if ((error = grav_def (jd_tdb,loc,accuracy,pos3,pob, pos4)) != 0)
return (error += 70);
aberration (pos4,vob,t_light, pos5);
}
/*
---------------------------------------------------------------------
Transform, if necessary, to output coordinate system.
---------------------------------------------------------------------
*/
switch (coord_sys)
{
case (1): /* Transform to equator and equinox of date. */
frame_tie (pos5,1, pos6);
precession (T0,pos6,jd_tdb, pos7);
nutation (jd_tdb,0,accuracy,pos7, pos8);
break;
case (2): /* Transform to equator and CIO of date. */
if (fabs (jd_tdb - tlast2) > 1.0e-8 )
{
/*
Obtain the basis vectors, in the GCRS, of the celestial intermediate
system.
*/
if ((error = cio_location (jd_tdb,accuracy, &r_cio,
&rs)) != 0)
return (error += 80);
if ((error = cio_basis (jd_tdb,r_cio,rs,accuracy,
px,py,pz)) != 0)
return (error += 90);
tlast2 = jd_tdb;
}
/*
Transform position vector to celestial intermediate system.
*/
pos8[0] = px[0] * pos5[0] + px[1] * pos5[1] + px[2] * pos5[2];
pos8[1] = py[0] * pos5[0] + py[1] * pos5[1] + py[2] * pos5[2];
pos8[2] = pz[0] * pos5[0] + pz[1] * pos5[1] + pz[2] * pos5[2];
break;
default: /* No transformation -- keep coordinates in GCRS, */
/* or ICRS for astrometric coordinates. */
for (i = 0; i < 3; i++)
{
pos8[i] = pos5[i];
}
}
/*
---------------------------------------------------------------------
Compute radial velocity.
---------------------------------------------------------------------
*/
/*
Compute distances: observer-geocenter, observer-Sun, object-Sun.
*/
d_obs_geo = sqrt ((pob[0] - peb[0]) * (pob[0] - peb[0]) +
( pob[1] - peb[1]) * (pob[1] - peb[1]) +
( pob[2] - peb[2]) * (pob[2] - peb[2]));
d_obs_sun = sqrt ((pob[0] - psb[0]) * (pob[0] - psb[0]) +
( pob[1] - psb[1]) * (pob[1] - psb[1]) +
( pob[2] - psb[2]) * (pob[2] - psb[2]));
d_obj_sun = sqrt ((pos1[0] - psb[0]) * (pos1[0] - psb[0]) +
( pos1[1] - psb[1]) * (pos1[1] - psb[1]) +
( pos1[2] - psb[2]) * (pos1[2] - psb[2]));
rad_vel (cel_object,pos3,vel1,vob,d_obs_geo,d_obs_sun,d_obj_sun,
&output->rv);
/*
---------------------------------------------------------------------
Finish up.
---------------------------------------------------------------------
*/
vector2radec (pos8, &output->ra,&output->dec);
x = sqrt (pos8[0] * pos8[0] + pos8[1] * pos8[1] + pos8[2] * pos8[2]);
for (i = 0; i < 3; i++)
{
output->r_hat[i] = pos8[i] / x;
}
return (error);
}
/********equ2gal */
void equ2gal (double rai, double deci,
double *glon, double *glat)
/*
------------------------------------------------------------------------
PURPOSE:
To convert ICRS right ascension and declination to galactic
longitude and latitude.
REFERENCES:
Hipparcos and Tycho Catalogues, Vol. 1, Section 1.5.3.
INPUT
ARGUMENTS:
rai (double)
ICRS right ascension in hours.
deci (double)
ICRS declination in degrees.
OUTPUT
ARGUMENTS:
*glon (double)
Galactic longitude in degrees.
*glat (double)
Galactic latitude in degrees.
RETURNED
VALUE:
None.
GLOBALS
USED:
DEG2RAD, RAD2DEG novascon.c
FUNCTIONS
CALLED:
sin math.h
cos math.h
sqrt math.h
atan2 math.h
VER./DATE/
PROGRAMMER:
V1.0/11-03/JAB (USNO/AA)
V1.1/03-06/JAB (USNO/AA): Fixed initialization of 'ag'.
V1.2/03-11/WKP (USNO/AA): Added braces to 2-D array initialization
to quiet gcc warnings.
NOTES:
1. This function uses the coordinate transformation specified
in the reference.
2. This function is the C version of NOVAS Fortran routine
'eqgal'.
------------------------------------------------------------------------
*/
{
double r, d, pos1[3], pos2[3], xyproj, g;
/*
Rotation matrix A_g from Hipparcos documentation eq. 1.5.11.
*/
double ag[3][3] = {
{-0.0548755604, +0.4941094279, -0.8676661490},
{-0.8734370902, -0.4448296300, -0.1980763734},
{-0.4838350155, +0.7469822445, +0.4559837762}};
/*
Form position vector in equatorial system from input coordinates.
*/
r = rai * 15.0 * DEG2RAD;
d = deci * DEG2RAD;
pos1[0] = cos (d) * cos (r);
pos1[1] = cos (d) * sin (r);
pos1[2] = sin (d);
/*
Rotate position vector to galactic system, using Hipparcos
documentation eq. 1.5.13.
*/
pos2[0] = ag[0][0] * pos1[0] + ag[1][0] * pos1[1] +
ag[2][0] * pos1[2];
pos2[1] = ag[0][1] * pos1[0] + ag[1][1] * pos1[1] +
ag[2][1] * pos1[2];
pos2[2] = ag[0][2] * pos1[0] + ag[1][2] * pos1[1] +
ag[2][2] * pos1[2];
/*
Decompose galactic vector into longitude and latitude.
*/
xyproj = sqrt (pos2[0] * pos2[0] + pos2[1] * pos2[1]);
if (xyproj > 0.0)
g = atan2 (pos2[1], pos2[0]);
else
g = 0.0;
*glon = g * RAD2DEG;
if (*glon < 0.0)
*glon += 360.0;
g = atan2 (pos2[2], xyproj);
*glat = g * RAD2DEG;
return;
}
/********equ2ecl */
short int equ2ecl (double jd_tt, short int coord_sys,
short int accuracy, double ra, double dec,
double *elon, double *elat)
/*
------------------------------------------------------------------------
PURPOSE:
To convert right ascension and declination to ecliptic longitude
and latitude.
REFERENCES:
None.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date of equator, equinox, and ecliptic used for
coordinates.
coord_sys (short int)
Coordinate system selection.
= 0 ... mean equator and equinox of date 'jd_tt'
= 1 ... true equator and equinox of date 'jd_tt'
= 2 ... ICRS
(ecliptic is always the mean plane)
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
ra (double)
Right ascension in hours, referred to specified equator and
equinox of date.
dec (double)
Declination in degrees, referred to specified equator and
equinox of date.
OUTPUT
ARGUMENTS:
*elon (double)
Ecliptic longitude in degrees, referred to specified ecliptic
and equinox of date.
*elat (double)
Ecliptic latitude in degrees, referred to specified ecliptic
and equinox of date.
RETURNED
VALUE:
(short int)
= 0 ... everything OK
= 1 ... invalid value of 'coord_sys'
GLOBALS
USED:
DEG2RAD, RAD2DEG novascon.c
FUNCTIONS
CALLED:
equ2ecl_vec novas.c
sin math.h
cos math.h
sqrt math.h
atan2 math.h
VER./DATE/
PROGRAMMER:
V1.0/11-03/JAB (USNO/AA)
V1.1/01-06/WKP (USNO/AA) Changed 'mode' to 'accuracy'.
V1.2/05-06/JAB (USNO/AA) Use vector transformation function.
V1.3/05-08/WKP (USNO/AA) Changed values of coord_sys to be
more consistent with gcrs2equ.
NOTES:
1. To convert ICRS RA and dec to ecliptic coordinates (mean
ecliptic and equinox of J2000.0), set 'coord_sys' = 2; the value
of 'jd_tt' can be set to anything, since J2000.0 is assumed.
Except for the input to this case, all input coordinates
are dynamical.
2. This function is the C version of NOVAS Fortran routine
'eqecl'.
------------------------------------------------------------------------
*/
{
short int error = 0;
double r, d, pos1[3], pos2[3], xyproj, e;
/*
Form position vector in equatorial system from input coordinates.
*/
r = ra * 15.0 * DEG2RAD;
d = dec * DEG2RAD;
pos1[0] = cos (d) * cos (r);
pos1[1] = cos (d) * sin (r);
pos1[2] = sin (d);
/*
Convert the vector from equatorial to ecliptic system.
*/
if ((error = equ2ecl_vec (jd_tt,coord_sys,accuracy,pos1, pos2)) != 0)
return (error);
/*
Decompose ecliptic vector into ecliptic longitude and latitude.
*/
xyproj = sqrt (pos2[0] * pos2[0] + pos2[1] * pos2[1]);
if (xyproj > 0.0)
e = atan2 (pos2[1], pos2[0]);
else
e = 0.0;
*elon = e * RAD2DEG;
if (*elon < 0.0)
*elon += 360.0;
e = atan2 (pos2[2], xyproj);
*elat = e * RAD2DEG;
return (error);
}
/********equ2ecl_vec */
short int equ2ecl_vec (double jd_tt, short int coord_sys,
short int accuracy, double *pos1,
double *pos2)
/*
------------------------------------------------------------------------
PURPOSE:
Converts an equatorial position vector to an ecliptic position
vector.
REFERENCES:
None.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date of equator, equinox, and ecliptic used for
coordinates.
coord_sys (short int)
Coordinate system selection.
= 0 ... mean equator and equinox of date 'jd_tt'
= 1 ... true equator and equinox of date 'jd_tt'
= 2 ... ICRS
(ecliptic is always the mean plane)
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
pos1[3] (double)
Position vector, referred to specified equator and equinox of
date.
OUTPUT
ARGUMENTS:
pos2[3] (double)
Position vector, referred to specified ecliptic and equinox
of date.
RETURNED
VALUE:
(short int)
= 0 ... everything OK
= 1 ... invalid value of 'coord_sys'
GLOBALS
USED:
T0, DEG2RAD novascon.c
FUNCTIONS
CALLED:
tdb2tt novas.c
e_tilt novas.c
frame_tie novas.c
fabs math.h
sin math.h
cos math.h
VER./DATE/
PROGRAMMER:
V1.0/05-06/JAB (USNO/AA)
V1.1/05-08/WKP (USNO/AA) Changed values of coord_sys to be
more consistent with gcrs2equ.
NOTES:
1. To convert an ICRS vector to an ecliptic vector (mean ecliptic
and equinox of J2000.0 only), set 'coord_sys' = 2; the value
of 'jd_tt' can be set to anything, since J2000.0 is assumed.
Except for the input to this case, all vectors are assumed
to be with respect to a dynamical system.
2. This function is the C version of NOVAS Fortran routine
'eqec'.
------------------------------------------------------------------------
*/
{
short int error = 0;
static double t_last = 0.0;
static double ob2000 = 0.0;
static double oblm, oblt;
double t, secdiff, jd_tdb, pos0[3], w, x, y, z, obl;
/*
'jd_tdb' is the TDB Julian date corresponding to 'jd_tt'.
*/
tdb2tt (jd_tt, &t,&secdiff);
jd_tdb = jd_tt + secdiff / 86400.0;
/*
Get obliquity, depending upon the "system" of the input coordinates.
*/
switch (coord_sys)
{
case 0: /* Input: mean equator and equinox of date */
case 1: /* Input: true equator and equinox of date */
pos0[0] = pos1[0];
pos0[1] = pos1[1];
pos0[2] = pos1[2];
if (fabs (jd_tt - t_last) > 1.0e-8)
{
e_tilt (jd_tdb,accuracy, &oblm,&oblt,&x,&y,&z);
t_last = jd_tt;
}
obl = DEG2RAD * ((coord_sys==0) ? oblm : oblt);
break;
case 2: /* Input: ICRS */
frame_tie (pos1,1, pos0);
if (ob2000 == 0.0)
{
e_tilt (T0,accuracy, &oblm,&w,&x,&y,&z);
ob2000 = oblm;
}
obl = ob2000 * DEG2RAD;
break;
default:
return (error = 1);
}
/*
Rotate position vector to ecliptic system.
*/
pos2[0] = pos0[0];
pos2[1] = pos0[1] * cos (obl) + pos0[2] * sin (obl);
pos2[2] = -pos0[1] * sin (obl) + pos0[2] * cos (obl);
return (error);
}
/********ecl2equ_vec */
short int ecl2equ_vec (double jd_tt, short int coord_sys,
short int accuracy, double *pos1,
double *pos2)
/*
------------------------------------------------------------------------
PURPOSE:
Converts an ecliptic position vector to an equatorial position
vector.
REFERENCES:
None.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date of equator, equinox, and ecliptic used for
coordinates.
coord_sys (short int)
Coordinate system selection.
= 0 ... mean equator and equinox of date
= 1 ... true equator and equinox of date
= 2 ... ICRS
(ecliptic is always the mean plane)
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
pos1[3] (double)
Position vector, referred to specified ecliptic and equinox of
date. If 'coord_sys' = 2, 'pos1' must be on mean ecliptic and
equinox of J2000.0; see Note 1 below.
OUTPUT
ARGUMENTS:
POS2[3] (double)
Position vector, referred to specified equator and equinox
of date.
RETURNED
VALUE:
(short int)
= 0 ... everything OK
= 1 ... invalid value of 'coord_sys'
GLOBALS
USED:
T0, DEG2RAD novascon.c
FUNCTIONS
CALLED:
tdb2tt novas.c
e_tilt novas.c
frame_tie novas.c
fabs math.h
sin math.h
cos math.h
VER./DATE/
PROGRAMMER:
V1.0/05-06/JAB (USNO/AA)
V1.1/05-08/WKP (USNO/AA) Changed values of coord_sys to be
more consistent with gcrs2equ.
V1.2/09-10/WKP (USNO/AA) Initialized 'obl' variable to silence
compiler warning.
NOTES:
1. To convert an ecliptic vector (mean ecliptic and equinox of
J2000.0 only) to an ICRS vector, set 'coord_sys' = 2; the value
of 'jd_tt' can be set to anything, since J2000.0 is assumed.
Except for the output from this case, all vectors are assumed to
be with respect to a dynamical system.
2. This function is the C version of NOVAS Fortran routine
'eceq'.
------------------------------------------------------------------------
*/
{
short int error = 0;
static double t_last = 0.0;
static double ob2000 = 0.0;
static double oblm, oblt;
double t, secdiff, jd_tdb, pos0[3], w, x, y, z, obl = 0.0;
/*
'jd_tdb' is the TDB Julian date.
*/
tdb2tt (jd_tt, &t,&secdiff);
jd_tdb = jd_tt + secdiff / 86400.0;
/*
Get obliquity, depending upon the "system" of the input coordinates.
*/
switch (coord_sys)
{
case 0: /* Output: mean equator and equinox of date */
case 1: /* Output: true equator and equinox of date */
if (fabs (jd_tt - t_last) > 1.0e-8)
{
e_tilt (jd_tdb,accuracy, &oblm,&oblt,&x,&y,&z);
t_last = jd_tt;
}
switch (coord_sys)
{
case 0: /* Use mean obliquity of date */
obl = oblm * DEG2RAD;
break;
case 1: /* Use true obliquity of date */
obl = oblt * DEG2RAD;
break;
}
break;
case 2: /* Output: ICRS */
if (ob2000 == 0.0)
{
e_tilt (T0,accuracy, &oblm,&w,&x,&y,&z);
ob2000 = oblm;
}
obl = ob2000 * DEG2RAD;
break;
default:
return (error = 1);
}
/*
Rotate position vector to ecliptic system.
*/
pos2[0] = pos1[0];
pos2[1] = pos1[1] * cos (obl) - pos1[2] * sin (obl);
pos2[2] = pos1[1] * sin (obl) + pos1[2] * cos (obl);
/*
Case where output vector is to be in ICRS, rotate from dynamical
system to ICRS.
*/
if (coord_sys == 2)
{
pos0[0] = pos2[0];
pos0[1] = pos2[1];
pos0[2] = pos2[2];
frame_tie (pos0,-1, pos2);
}
return (error);
}
/********equ2hor */
void equ2hor (double jd_ut1, double delta_t, short int accuracy,
double xp, double yp, on_surface *location, double ra,
double dec, short int ref_option,
double *zd, double *az, double *rar, double *decr)
/*
------------------------------------------------------------------------
PURPOSE:
This function transforms topocentric right ascension and
declination to zenith distance and azimuth. It uses a method
that properly accounts for polar motion, which is significant at
the sub-arcsecond level. This function can also adjust
coordinates for atmospheric refraction.
REFERENCES:
Kaplan, G. (2008). USNO/AA Technical Note of 28 Apr 2008,
"Refraction as a Vector."
INPUT
ARGUMENTS:
jd_ut1 (double)
UT1 Julian date.
delta_t (double)
Difference TT-UT1 at 'jd_ut1', in seconds.
accuracy (short int)
Selection for method and accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
xp (double)
Conventionally-defined x coordinate of celestial intermediate
pole with respect to ITRS reference pole, in arcseconds.
yp (double)
Conventionally-defined y coordinate of celestial intermediate
pole with respect to ITRS reference pole, in arcseconds.
*location (struct on_surface)
Pointer to structure containing observer's location (defined
in novas.h).
ra (double)
Topocentric right ascension of object of interest, in hours,
referred to true equator and equinox of date.
dec (double)
Topocentric declination of object of interest, in degrees,
referred to true equator and equinox of date.
ref_option (short int)
= 0 ... no refraction
= 1 ... include refraction, using 'standard' atmospheric
conditions.
= 2 ... include refraction, using atmospheric parameters
input in the 'location' structure.
OUTPUT
ARGUMENTS:
*zd (double)
Topocentric zenith distance in degrees, affected by
refraction if 'ref_option' is non-zero.
*az (double)
Topocentric azimuth (measured east from north) in degrees.
*rar (double)
Topocentric right ascension of object of interest, in hours,
referred to true equator and equinox of date, affected by
refraction if 'ref_option' is non-zero.
*decr (double)
Topocentric declination of object of interest, in degrees,
referred to true equator and equinox of date, affected by
refraction if 'ref_option' is non-zero.
RETURNED
VALUE:
None.
GLOBALS
USED:
DEG2RAD, RAD2DEG novascon.c
FUNCTIONS
CALLED:
ter2cel novas.c
refract novas.c
sin math.h
cos math.h
sqrt math.h
atan2 math.h
fabs math.h
VER./DATE/
PROGRAMMER:
V1.0/06-98/JAB (USNO/AA)
V2.0/01-06/JAB (USNO/AA)
V2.1/01-06/WKP (USNO/AA): Changed 'mode' to 'accuracy'.
V2.2/11-07/JAB (USNO/AA): Cosmetic changes.
V2.2/05-08/WKP (USNO/AA): Updated the refraction algorithm.
V2.3/06-08/WKP (USNO/AA): Tweaked convergence criteria.
V2.4/03-09/JAB (USNO/AA): Conformed input variables to IERS
conventions.
NOTES:
1. 'xp' and 'yp' can be set to zero if sub-arcsecond accuracy is
not needed. 'ra' and 'dec' can be obtained from functions
'tpstar' or 'tpplan'.
2. The directions 'zd'= 0 (zenith) and 'az'= 0 (north) are here
considered fixed in the terrestrial system. Specifically, the
zenith is along the geodetic normal, and north is toward
the ITRS pole.
3. If 'ref_option'= 0, then 'rar'='ra' and 'decr'='dec'.
4. This function is the C version of NOVAS Fortran routine
'zdaz'.
------------------------------------------------------------------------
*/
{
short int j;
double sinlat, coslat, sinlon, coslon, sindc, cosdc, sinra, cosra,
uze[3], une[3], uwe[3], uz[3], un[3], uw[3], p[3], pz, pn, pw,
proj, zd0, zd1, refr, sinzd, coszd, sinzd0, coszd0, pr[3];
/*
Preliminaries.
*/
*rar = ra;
*decr = dec;
sinlat = sin (location->latitude * DEG2RAD);
coslat = cos (location->latitude * DEG2RAD);
sinlon = sin (location->longitude * DEG2RAD);
coslon = cos (location->longitude * DEG2RAD);
sindc = sin (dec * DEG2RAD);
cosdc = cos (dec * DEG2RAD);
sinra = sin (ra * 15.0 * DEG2RAD);
cosra = cos (ra * 15.0 * DEG2RAD);
/*
Set up orthonormal basis vectors in local Earth-fixed system.
Define vector toward local zenith in Earth-fixed system (z axis).
*/
uze[0] = coslat * coslon;
uze[1] = coslat * sinlon;
uze[2] = sinlat;
/*
Define vector toward local north in Earth-fixed system (x axis).
*/
une[0] = -sinlat * coslon;
une[1] = -sinlat * sinlon;
une[2] = coslat;
/*
Define vector toward local west in Earth-fixed system (y axis).
*/
uwe[0] = sinlon;
uwe[1] = -coslon;
uwe[2] = 0.0;
/*
Obtain vectors in celestial system.
Rotate Earth-fixed orthonormal basis vectors to celestial system
(wrt equator and equinox of date).
*/
ter2cel (jd_ut1,0.0,delta_t,1,accuracy,1,xp,yp,uze, uz);
ter2cel (jd_ut1,0.0,delta_t,1,accuracy,1,xp,yp,une, un);
ter2cel (jd_ut1,0.0,delta_t,1,accuracy,1,xp,yp,uwe, uw);
/*
Define unit vector 'p' toward object in celestial system
(wrt equator and equinox of date).
*/
p[0] = cosdc * cosra;
p[1] = cosdc * sinra;
p[2] = sindc;
/*
Compute coordinates of object wrt orthonormal basis.
Compute components of 'p' - projections of 'p' onto rotated
Earth-fixed basis vectors.
*/
pz = p[0] * uz[0] + p[1] * uz[1] + p[2] * uz[2];
pn = p[0] * un[0] + p[1] * un[1] + p[2] * un[2];
pw = p[0] * uw[0] + p[1] * uw[1] + p[2] * uw[2];
/*
Compute azimuth and zenith distance.
*/
proj = sqrt (pn * pn + pw * pw);
if (proj > 0.0)
*az = -atan2 (pw, pn) * RAD2DEG;
if (*az < 0.0)
*az += 360.0;
if (*az >= 360.0)
*az -= 360.0;
*zd = atan2 (proj, pz) * RAD2DEG;
/*
Apply atmospheric refraction if requested.
*/
if (ref_option != 0)
{
/*
Get refraction in zenith distance.
Iterative process is required because refraction algorithms are
always a function of observed (not computed) zenith distance.
Require convergence to 0.1 arcsec (actual accuracy less).
*/
zd0 = *zd;
do
{
zd1 = *zd;
refr = refract (location,ref_option,*zd);
*zd = zd0 - refr;
} while (fabs (*zd - zd1) > 3.0e-5);
/*
Apply refraction to celestial coordinates of object.
*/
if ((refr > 0.0) && (*zd > 3.0e-4))
{
/*
Shift position vector of object in celestial system to account
for refraction (see USNO/AA Technical Note 1998-09).
*/
sinzd = sin (*zd * DEG2RAD);
coszd = cos (*zd * DEG2RAD);
sinzd0 = sin (zd0 * DEG2RAD);
coszd0 = cos (zd0 * DEG2RAD);
/*
Compute refracted position vector.
*/
for (j = 0; j < 3; j++)
pr[j] = ((p[j] - coszd0 * uz[j]) / sinzd0) * sinzd + uz[j] *
coszd;
/*
Compute refracted right ascension and declination.
*/
proj = sqrt (pr[0] * pr[0] + pr[1] * pr[1]);
if (proj > 0.0)
*rar = atan2 (pr[1],pr[0]) * RAD2DEG / 15.0;
if (*rar < 0.0)
*rar += 24.0;
if (*rar >= 24.0)
*rar -= 24.0;
*decr = atan2 (pr[2],proj) * RAD2DEG;
}
}
return;
}
/********gcrs2equ */
short int gcrs2equ (double jd_tt, short int coord_sys,
short int accuracy, double rag, double decg,
double *ra, double *dec)
/*
------------------------------------------------------------------------
PURPOSE:
This function converts GCRS right ascension and declination
to coordinates with respect to the equator of date (mean or true).
For coordinates with respect to the true equator of date, the
origin of right ascension can be either the true equinox or the
celestial intermediate origin (CIO).
REFERENCES:
None.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date of equator to be used for output coordinates.
coord_sys (short int)
Coordinate system selection for output coordinates.
= 0 ... mean equator and equinox of date
= 1 ... true equator and equinox of date
= 2 ... true equator and CIO of date
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
rag (double)
GCRS right ascension in hours.
decg (double)
GCRS declination in degrees.
OUTPUT
ARGUMENTS:
ra (double)
Right ascension in hours, referred to specified equator and
right ascension origin of date.
dec (double)
Declination in degrees, referred to specified equator of date.
RETURNED
VALUE:
(short int)
= 0 ... everything OK.
< 0 ... error from function 'vector2radec'
> 10 ... 10 + error from function 'cio_location'
> 20 ... 20 + error from function 'cio_basis'
GLOBALS
USED:
DEG2RAD, T0 novascon.c
FUNCTIONS
CALLED:
tdb2tt novas.c
frame_tie novas.c
precession novas.c
nutation novas.c
cio_location novas.c
cio_basis novas.c
vector2radec novas.c
sin math.h
cos math.h
VER./DATE/
PROGRAMMER:
V1.0/09-04/JAB (USNO/AA)
V1.1/01-06/WKP (USNO/AA) Changed 'mode' to 'accuracy'.
V1.2/07-06/JAB (USNO/AA) Implement 'cio_location' construct.
V1.3/04-08/WKP (USNO/AA) Updated variable names.
NOTES:
1. Set input value of 'accuracy' equal to any short int if
'coord_sys' equals 0 or 1. It is not used in these cases.
2. This function only supports the CIO-based method.
3. This function is the C version of NOVAS Fortran routine
'gcrseq'.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int rs;
double t1, t, secdiff, r, d, pos1[3], pos2[3], pos3[3], pos4[3],
r_cio, x[3], y[3], z[3];
/*
't1' is the TDB Julian date.
*/
tdb2tt (jd_tt, &t,&secdiff);
t1 = jd_tt + secdiff / 86400.0;
/*
Form position vector in equatorial system from input coordinates.
*/
r = rag * 15.0 * DEG2RAD;
d = decg * DEG2RAD;
pos1[0] = cos (d) * cos (r);
pos1[1] = cos (d) * sin (r);
pos1[2] = sin (d);
/*
Transform the position vector based on the value of 'coord_sys'.
*/
if (coord_sys <= 1 )
{
/*
Transform the position vector from GCRS to mean equator and equinox
of date.
*/
frame_tie (pos1,1, pos2);
precession (T0,pos2,t1, pos3);
/*
If requested, transform further to true equator and equinox of date.
*/
if (coord_sys == 1)
{
nutation (t1,0,accuracy,pos3, pos4);
}
else
{
pos4[0] = pos3[0];
pos4[1] = pos3[1];
pos4[2] = pos3[2];
}
}
else
{
/*
Obtain the basis vectors, in the GCRS, of the celestial intermediate
system.
*/
if ((error = cio_location (t1,accuracy, &r_cio,&rs)) != 0)
return (error += 10);
if ((error = cio_basis (t1,r_cio,rs,accuracy, x,y,z)) != 0)
return (error += 20);
/*
Transform position vector to the celestial intermediate system
(which has the CIO as its origin of right ascension).
*/
pos4[0] = x[0] * pos1[0] + x[1] * pos1[1] + x[2] * pos1[2];
pos4[1] = y[0] * pos1[0] + y[1] * pos1[1] + y[2] * pos1[2];
pos4[2] = z[0] * pos1[0] + z[1] * pos1[1] + z[2] * pos1[2];
}
/*
Convert the position vector to equatorial spherical coordinates.
*/
if ((error = vector2radec (pos4, ra,dec)) != 0)
{
*ra = 0.0;
*dec = 0.0;
return (-1 * error);
}
return (error);
}
/********sidereal_time */
short int sidereal_time (double jd_high, double jd_low,
double delta_t,short int gst_type,
short int method, short int accuracy,
double *gst)
/*
------------------------------------------------------------------------
PURPOSE:
Computes the Greenwich sidereal time, either mean or apparent, at
Julian date 'jd_high' + 'jd_low'.
REFERENCES:
Kaplan, G. (2005), US Naval Observatory Circular 179.
INPUT
ARGUMENTS:
jd_high (double)
High-order part of UT1 Julian date.
jd_low (double)
Low-order part of UT1 Julian date.
delta_t (double)
Difference TT-UT1 at 'jd_high'+'jd_low', in seconds
of time.
gst_type (short int)
= 0 ... compute Greenwich mean sidereal time
= 1 ... compute Greenwich apparent sidereal time
method (short int)
Selection for method
= 0 ... CIO-based method
= 1 ... equinox-based method
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*gst (double)
Greenwich (mean or apparent) sidereal time, in hours.
RETURNED
VALUE:
(short int)
= 0 ... everything OK
= 1 ... invalid value of 'accuracy'
= 2 ... invalid value of 'method'
> 10, < 30 ... 10 + error from function 'cio_rai'
GLOBALS
USED:
T0, RAD2DEG novascon.c
FUNCTIONS
CALLED:
tdb2tt novas.c
era novas.c
e_tilt novas.c
cio_location novas.c
cio_basis novas.c
nutation novas.c
precession novas.c
frame_tie novas.c
fabs math.h
atan2 math.h
fmod math.h
VER./DATE/
PROGRAMMER:
V1.0/06-92/TKB (USNO/NRL Optical Interfer.) Translate Fortran.
V1.1/08-93/WTH (USNO/AA) Update to C programing standards.
V1.2/03-98/JAB (USNO/AA) Expand documentation.
V1.3/08-98/JAB (USNO/AA) Match flow of the Fortran counterpart.
V2.0/09-03/JAB (USNO/AA) Incorporate the 2nd-reference changes.
V2.1/08-04/JAB (USNO/AA) Incorporate the 1st-reference changes.
V2.2/12-05/WKP (USNO/AA) Updated error handling.
V2.3/01-06/WKP (USNO/AA) Changed 'mode' to 'method' and 'accuracy'.
V2.4/04-06/JAB (USNO/AA) Use precession-in-RA terms in mean
sidereal time from third reference.
V2.5/07-06/JAB (USNO/AA) Implement 'cio_location' construct.
V2.6/06-08/WKP (USNO/AA) Changed value of direction argument in
call to 'nutation' from 1 to -1 for
consistency.
V2.7/03-11/WKP (USNO/AA) Updated prolog description to clarify
this function computes either mean or
apparent sidereal time, and removed
Note 1 for consistency with Fortran.
NOTES:
1. The Julian date may be split at any point, but for highest
precision, set 'jd_high' to be the integral part of the Julian
date, and set 'jd_low' to be the fractional part.
2. This function is the C version of NOVAS Fortran routine
'sidtim'.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int ref_sys;
static double ee;
static double jd_last = -99.0;
double unitx[3] = {1.0, 0.0, 0.0};
double jd_ut, jd_tt, jd_tdb, tt_temp, t, theta, a, b, c, d,
ra_cio, x[3], y[3], z[3], w1[3], w2[3], eq[3], ha_eq, st,
secdiff, eqeq;
/*
Invalid value of 'accuracy'.
*/
if ((accuracy < 0) || (accuracy > 1))
return (error = 1);
/*
Time argument for precession and nutation components of sidereal
time is TDB. First approximation is TDB = TT, then refine.
*/
jd_ut = jd_high + jd_low;
jd_tt = jd_ut + (delta_t / 86400.0);
jd_tdb = jd_tt;
tdb2tt (jd_tdb, &tt_temp,&secdiff);
jd_tdb = jd_tt + (secdiff / 86400.0);
t = (jd_tdb - T0) / 36525.0;
/*
Compute the Earth Rotation Angle. Time argument is UT1.
*/
theta = era (jd_high, jd_low);
/*
Compute the equation of the equinoxes if needed, depending upon the
input values of 'gst_type' and 'method'. If not needed, set to zero.
*/
if (((gst_type == 0) && (method == 0)) || /* GMST; CIO-TIO */
((gst_type == 1) && (method == 1))) /* GAST; equinox */
{
if (fabs (jd_tdb - jd_last) > 1.0e-8)
{
e_tilt (jd_tdb,accuracy, &a,&b,&ee,&c,&d);
jd_last = jd_tdb;
}
eqeq = ee * 15.0;
}
else
{
eqeq = 0.0;
}
/*
Compute Greenwich sidereal time depending upon input values of
'method' and 'gst_type'.
*/
switch (method)
{
case (0):
/*
Use 'CIO-TIO-theta' method. See Circular 179, Section 6.5.4.
*/
/*
Obtain the basis vectors, in the GCRS, of the celestial intermediate
system.
*/
if ((error = cio_location (jd_tdb,accuracy, &ra_cio,
&ref_sys)) != 0)
{
*gst = 99.0;
return (error += 10);
}
cio_basis (jd_tdb,ra_cio,ref_sys,accuracy, x,y,z);
/*
Compute the direction of the true equinox in the GCRS.
*/
nutation (jd_tdb,-1,accuracy,unitx, w1);
precession (jd_tdb,w1,T0, w2);
frame_tie (w2,-1, eq);
/*
Compute the hour angle of the equinox wrt the TIO meridian
(near Greenwich, but passes through the CIP and TIO).
*/
ha_eq = theta - atan2 ((eq[0] * y[0] + eq[1] * y[1] +
eq[2] * y[2]), (eq[0] * x[0] + eq[1] * x[1] +
eq[2] * x[2])) * RAD2DEG;
/*
For mean sidereal time, subtract the equation of the equinoxes.
*/
ha_eq -= (eqeq / 240.0);
ha_eq = fmod (ha_eq, 360.0) / 15.0;
if (ha_eq < 0.0)
ha_eq += 24.0;
*gst = ha_eq;
break;
case (1):
/*
Use equinox method. See Circular 179, Section 2.6.2.
*/
/*
Precession-in-RA terms in mean sidereal time taken from third
reference, eq. (42), with coefficients in arcseconds.
*/
st = eqeq + 0.014506 +
(((( - 0.0000000368 * t
- 0.000029956 ) * t
- 0.00000044 ) * t
+ 1.3915817 ) * t
+ 4612.156534 ) * t;
/*
Form the Greenwich sidereal time.
*/
*gst = fmod ((st / 3600.0 + theta), 360.0) / 15.0;
if (*gst < 0.0)
*gst += 24.0;
break;
/*
Invalid value of 'method'.
*/
default:
*gst = 99.0;
error = 2;
break;
}
return (error);
}
/********era */
double era (double jd_high, double jd_low)
/*
------------------------------------------------------------------------
PURPOSE:
This function returns the value of the Earth Rotation Angle
(theta) for a given UT1 Julian date. The expression used is
taken from the note to IAU Resolution B1.8 of 2000.
REFERENCES:
IAU Resolution B1.8, adopted at the 2000 IAU General Assembly,
Manchester, UK.
Kaplan, G. (2005), US Naval Observatory Circular 179.
INPUT
ARGUMENTS:
jd_high (double)
High-order part of UT1 Julian date.
jd_low (double)
Low-order part of UT1 Julian date.
OUTPUT
ARGUMENTS:
None.
RETURNED
VALUE:
(double)
The Earth Rotation Angle (theta) in degrees.
GLOBALS
USED:
T0 novascon.c
FUNCTIONS
CALLED:
fmod math.h
VER./DATE/
PROGRAMMER:
V1.0/09-03/JAB (USNO/AA)
NOTES:
1. The algorithm used here is equivalent to the canonical
theta = 0.7790572732640 + 1.00273781191135448 * t,
where t is the time in days from J2000 (t = jd_high +
jd_low - T0), but it avoids many two-PI 'wraps' that decrease
precision (adopted from SOFA Fortran routine iau_era00; see
also expression at top of page 35 of IERS Conventions (1996)).
2. This function is the C version of NOVAS Fortran routine
'erot'.
------------------------------------------------------------------------
*/
{
double theta, thet1, thet2, thet3;
thet1 = 0.7790572732640 + 0.00273781191135448 * (jd_high - T0);
thet2 = 0.00273781191135448 * jd_low;
thet3 = fmod (jd_high, 1.0) + fmod (jd_low, 1.0);
theta = fmod (thet1 + thet2 + thet3, 1.0) * 360.0;
if (theta < 0.0)
theta += 360.0;
return theta;
}
/********ter2cel */
short int ter2cel (double jd_ut_high, double jd_ut_low, double delta_t,
short int method, short int accuracy, short int option,
double xp, double yp, double *vec1,
double *vec2)
/*
------------------------------------------------------------------------
PURPOSE:
This function rotates a vector from the terrestrial to the
celestial system. Specifically, it transforms a vector in the
ITRS (rotating earth-fixed system) to the GCRS (a local space-
fixed system) by applying rotations for polar motion, Earth
rotation, nutation, precession, and the dynamical-to-GCRS
frame tie.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Kaplan, G. H. (2003), 'Another Look at Non-Rotating Origins',
Proceedings of IAU XXV Joint Discussion 16.
INPUT
ARGUMENTS:
jd_ut_high (double)
High-order part of UT1 Julian date.
jd_ut_low (double)
Low-order part of UT1 Julian date.
delta_t (double)
Value of Delta T (= TT - UT1) at the input UT1 Julian date.
method (short int)
Selection for method
= 0 ... CIO-based method
= 1 ... equinox-based method
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
option (short int)
= 0 ... The output vector is referred to GCRS axes.
= 1 ... The output vector is produced with respect to the
equator and equinox of date.
See Note 2 below.
xp (double)
Conventionally-defined X coordinate of celestial intermediate
pole with respect to ITRS pole, in arcseconds.
yp (double)
Conventionally-defined Y coordinate of celestial intermediate
pole with respect to ITRS pole, in arcseconds.
vec1[3] (double)
Position vector, geocentric equatorial rectangular coordinates,
referred to ITRS axes (terrestrial system) in the normal case
where 'option' = 0.
OUTPUT
ARGUMENTS:
vec2[3] (double)
Position vector, geocentric equatorial rectangular coordinates,
referred to GCRS axes (celestial system) or with respect to
the equator and equinox of date, depending on 'option'.
RETURNED
VALUE:
= 0 ... everything is ok.
= 1 ... invalid value of 'accuracy'
= 2 ... invalid value of 'method'
> 10 ... 10 + error from function 'cio_location'
> 20 ... 20 + error from function 'cio_basis'
GLOBALS
USED:
T0 novascon.c
FUNCTIONS
CALLED:
tdb2tt novas.c
wobble novas.c
cio_location novas.c
cio_basis novas.c
era novas.c
spin novas.c
sidereal_time novas.c
nutation novas.c
precession novas.c
frame_tie novas.c
VER./DATE/
PROGRAMMER:
V1.0/08-93/WTH (USNO/AA) Translate Fortran.
V2.0/10-03/JAB (USNO/AA) Formerly function 'pnsw.' Update
for IAU 2000 resolutions.
V2.1/01-05/JAB (USNO/AA) Update to include 'mode' flexibility.
V2.2/12-05/WKP (USNO/AA) Check error codes.
V2.3/01-06/WKP (USNO/AA) Changed 'mode' to 'method' and 'accuracy'.
V2.4/07-06/JAB (USNO/AA) Implement 'cio_location' construct.
V2.5/06-08/WKP (USNO/AA) Corrected sign of direction argument in
'frame_tie' call and changed value of
direction argument in call to 'nutation'
for consistency.
V2.6/10-08/WKP (USNO/AA) Renamed input JD variables to include
time system.
V2.7/03-09/JAB (USNO/AA) Conformed input variables to IERS
conventions.
V2.8/12-10/JAB (USNO/AA) Cosmetic changes for consistency with
Fortran.
NOTES:
1. 'xp' = 'yp' = 0 means no polar motion transformation.
2. The 'option' flag only works for the equinox-based method.
3. This function is the C version of NOVAS Fortran routine
'tercel'.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int rs, j;
double jd_ut1, jd_tt, dummy, secdiff, jd_tdb, gast, r_cio, theta,
v1[3], v2[3], v3[3], v4[3], x[3], y[3], z[3];
vec2[0] = vec2[1] = vec2[2] = 0.0; /* [Don Cross]: fix build warning. */
/*
Invalid value of 'accuracy'.
*/
if ((accuracy < 0) || (accuracy > 1))
return (error = 1);
/*
Compute the TT Julian date corresponding to the input UT1 Julian
date.
*/
jd_ut1 = jd_ut_high + jd_ut_low;
jd_tt = jd_ut1 + (delta_t / 86400.0);
/*
Compute the TDB Julian date corresponding to the input UT1 Julian
date.
*/
jd_tdb = jd_tt;
tdb2tt (jd_tdb, &dummy,&secdiff);
jd_tdb = jd_tt + secdiff / 86400.0;
switch (method)
{
case (0):
/*
'CIO-TIO-THETA' method.
See second reference, eq. (3) and (4).
Apply polar motion, transforming the vector to the terrestrial
intermediate system.
*/
if ((xp == 0.0) && (yp == 0.0))
{
v1[0] = vec1[0];
v1[1] = vec1[1];
v1[2] = vec1[2];
}
else
wobble (jd_tdb,0,xp,yp,vec1, v1);
/*
Obtain the basis vectors, in the GCRS, of the celestial intermediate
system.
*/
if ((error = cio_location (jd_tdb,accuracy, &r_cio,&rs)) != 0)
return (error += 10);
if ((error = cio_basis (jd_tdb,r_cio,rs,accuracy, x,y,z)) != 0)
return (error += 20);
/*
Compute and apply the Earth rotation angle, 'theta', transforming the
vector to the celestial intermediate system.
*/
theta = era (jd_ut_high,jd_ut_low);
spin (-theta,v1, v2);
/*
Transform the vector from the celestial intermediate system to the
GCRS.
*/
vec2[0] = x[0] * v2[0] + y[0] * v2[1] + z[0] * v2[2];
vec2[1] = x[1] * v2[0] + y[1] * v2[1] + z[1] * v2[2];
vec2[2] = x[2] * v2[0] + y[2] * v2[1] + z[2] * v2[2];
break;
case (1):
/*
Equinox mode.
Apply polar motion.
*/
if ((xp == 0.0) && (yp == 0.0))
{
for (j = 0; j < 3; j++)
{
v1[j] = vec1[j];
}
}
else
wobble (jd_tdb,0,xp,yp,vec1, v1);
/*
Apply Earth rotation.
*/
sidereal_time (jd_ut_high,jd_ut_low,delta_t,1,1,accuracy, &gast);
spin (-gast * 15.0,v1, v2);
/*
'option' = 1 skips remaining transformations.
*/
if (option == 1)
{
vec2[0] = v2[0];
vec2[1] = v2[1];
vec2[2] = v2[2];
}
else
{
/*
Apply precession, nutation, and frame tie.
*/
nutation (jd_tdb,-1,accuracy,v2, v3);
precession (jd_tdb,v3,T0, v4);
frame_tie (v4,-1, vec2);
}
break;
/*
Invalid value of 'method'.
*/
default:
error = 2;
break;
}
return (error);
}
/********cel2ter */
short int cel2ter (double jd_ut_high, double jd_ut_low, double delta_t,
short int method, short int accuracy, short int option,
double xp, double yp, double *vec1,
double *vec2)
/*
------------------------------------------------------------------------
PURPOSE:
This function rotates a vector from the celestial to the
terrestrial system. Specifically, it transforms a vector in the
GCRS (a local space-fixed system) to the ITRS (a rotating earth-
fixed system) by applying rotations for the GCRS-to-dynamical
frame tie, precession, nutation, Earth rotation, and polar motion.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Kaplan, G. H. (2003), 'Another Look at Non-Rotating Origins',
Proceedings of IAU XXV Joint Discussion 16.
INPUT
ARGUMENTS:
jd_ut_high (double)
High-order part of UT1 Julian date.
jd_ut_low (double)
Low-order part of UT1 Julian date.
delta_t (double)
Value of Delta T (= TT - UT1) at the input UT1 Julian date.
method (short int)
Selection for method
= 0 ... CIO-based method
= 1 ... equinox-based method
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
option (short int)
= 0 ... The input vector is referred to GCRS axes.
= 1 ... The input vector is produced with respect to the
equator and equinox of date.
See Note 2 below.
xp (double)
Conventionally-defined X coordinate of celestial intermediate
pole with respect to ITRS pole, in arcseconds.
yp (double)
Conventionally-defined Y coordinate of celestial intermediate
pole with respect to ITRS pole, in arcseconds.
vec1[3] (double)
Position vector, geocentric equatorial rectangular coordinates,
referred to GCRS axes (celestial system) or with respect to
the equator and equinox of date, depending on 'option'.
OUTPUT
ARGUMENTS:
vec2[3] (double)
Position vector, geocentric equatorial rectangular coordinates,
referred to ITRS axes (terrestrial system).
RETURNED
VALUE:
= 0 ... everything is ok.
= 1 ... invalid value of 'accuracy'
= 2 ... invalid value of 'method'
> 10 ... 10 + error from function 'cio_location'
> 20 ... 20 + error from function 'cio_basis'
GLOBALS
USED:
T0 novascon.c
FUNCTIONS
CALLED:
tdb2tt novas.c
wobble novas.c
cio_location novas.c
cio_basis novas.c
era novas.c
spin novas.c
sidereal_time novas.c
nutation novas.c
precession novas.c
frame_tie novas.c
VER./DATE/
PROGRAMMER:
V1.0/03-09/JAB (USNO/AA)
V1.1/12-10/JAB (USNO/AA) Fix confusion with input/output vectors;
cosmetic changes to better match Fortran.
NOTES:
1. 'xp' = 'yp' = 0 means no polar motion transformation.
2. 'option' = 1 only works for the equinox-based method
('method' = 1).
3. This function is the C version of NOVAS Fortran routine
'celter'.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int rs, j;
double jd_ut1, jd_tt, dummy, secdiff, jd_tdb, gast, r_cio, theta,
v1[3], v2[3], v3[3], v4[3], x[3], y[3], z[3];
/*
Invalid value of 'accuracy'.
*/
if ((accuracy < 0) || (accuracy > 1))
return (error = 1);
/*
Compute the TT Julian date corresponding to the input UT1 Julian
date.
*/
jd_ut1 = jd_ut_high + jd_ut_low;
jd_tt = jd_ut1 + (delta_t / 86400.0);
/*
Compute the TDB Julian date corresponding to the input UT1 Julian
date.
*/
jd_tdb = jd_tt;
tdb2tt (jd_tdb, &dummy,&secdiff);
jd_tdb = jd_tt + secdiff / 86400.0;
switch (method)
{
case (0):
/*
'CIO-TIO-THETA' method.
See second reference, eq. (3) and (4).
Obtain the basis vectors, in the GCRS, of the celestial intermediate
system.
*/
if ((error = cio_location (jd_tdb,accuracy, &r_cio,&rs)) != 0)
return (error += 10);
if ((error = cio_basis (jd_tdb,r_cio,rs,accuracy, x,y,z)) != 0)
return (error += 20);
/*
Transform the vector from the GCRS to the celestial intermediate
system.
*/
v1[0] = x[0] * vec1[0] + x[1] * vec1[1] + x[2] * vec1[2];
v1[1] = y[0] * vec1[0] + y[1] * vec1[1] + y[2] * vec1[2];
v1[2] = z[0] * vec1[0] + z[1] * vec1[1] + z[2] * vec1[2];
/*
Compute and apply the Earth rotation angle, 'theta', transforming the
vector to the terrestrial intermediate system.
*/
theta = era (jd_ut_high,jd_ut_low);
spin (theta,v1, v2);
/*
Apply polar motion, transforming the vector to the ITRS.
*/
if ((xp == 0.0) && (yp == 0.0))
{
vec2[0] = v2[0];
vec2[1] = v2[1];
vec2[2] = v2[2];
}
else
wobble (jd_tdb,1,xp,yp,v2, vec2);
break;
case (1):
/*
Equinox mode.
'option' = 1 skips initial transformations.
*/
if (option == 1)
{
v3[0] = vec1[0];
v3[1] = vec1[1];
v3[2] = vec1[2];
}
else
{
/*
Apply frame tie, nutation, and precession.
*/
frame_tie (vec1,1, v1);
precession (T0,v1,jd_tdb, v2);
nutation (jd_tdb,0,accuracy,v2, v3);
}
/*
Apply Earth rotation.
*/
sidereal_time (jd_ut_high,jd_ut_low,delta_t,1,1,accuracy, &gast);
spin (gast * 15.0,v3, v4);
/*
Apply polar motion.
*/
if ((xp == 0.0) && (yp == 0.0))
{
for (j = 0; j < 3; j++)
{
vec2[j] = v4[j];
}
}
else
wobble (jd_tdb,1,xp,yp,v4, vec2);
break;
/*
Invalid value of 'method'.
*/
default:
error = 2;
break;
}
return (error);
}
/********spin */
void spin (double angle, double *pos1,
double *pos2)
/*
------------------------------------------------------------------------
PURPOSE:
This function transforms a vector from one coordinate system
to another with same origin and axes rotated about the z-axis.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
INPUT
ARGUMENTS:
angle (double)
Angle of coordinate system rotation, positive counterclockwise
when viewed from +z, in degrees.
pos1[3] (double)
Position vector.
OUTPUT
ARGUMENTS:
pos2[3] (double)
Position vector expressed in new coordinate system rotated
about z by 'angle'.
RETURNED
VALUE:
None.
GLOBALS
USED:
DEG2RAD novascon.c
FUNCTIONS
CALLED:
sin math.h
cos math.h
VER./DATE/
PROGRAMMER:
V1.0/08-93/WTH (USNO/AA) Translate Fortran.
V2.0/10-03/JAB (USNO/AA) Update for IAU 2000 resolutions.
V2.1/01-05/JAB (USNO/AA) Generalize the function.
NOTES:
1. This function is the C version of NOVAS Fortran routine 'spin'.
------------------------------------------------------------------------
*/
{
static double ang_last = -999.0;
static double xx, yx, zx, xy, yy, zy, xz, yz, zz;
double angr, cosang, sinang;
if (fabs (angle - ang_last) >= 1.0e-12)
{
angr = angle * DEG2RAD;
cosang = cos (angr);
sinang = sin (angr);
/*
Rotation matrix follows.
*/
xx = cosang;
yx = sinang;
zx = 0.0;
xy = -sinang;
yy = cosang;
zy = 0.0;
xz = 0.0;
yz = 0.0;
zz = 1.0;
ang_last = angle;
}
/*
Perform rotation.
*/
pos2[0] = xx * pos1[0] + yx * pos1[1] + zx * pos1[2];
pos2[1] = xy * pos1[0] + yy * pos1[1] + zy * pos1[2];
pos2[2] = xz * pos1[0] + yz * pos1[1] + zz * pos1[2];
return;
}
/********wobble */
void wobble (double tjd, short int direction, double xp, double yp,
double *pos1,
double *pos2)
/*
------------------------------------------------------------------------
PURPOSE:
Corrects a vector in the ITRS (rotating Earth-fixed system)
for polar motion, and also corrects the longitude origin
(by a tiny amount) to the Terrestrial Intermediate Origin
(TIO). The ITRS vector is thereby transformed to the terrestrial
intermediate system, based on the true (rotational) equator and
TIO. Because the true equator is the plane orthogonal to the
direction of the Celestial Intermediate Pole (CIP), the components
of the output vector are referred to z and x axes toward the CIP
and TIO, respectively.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
Lambert & Bizouard (2002), Astronomy and Astrophysics 394,
317-321.
INPUT
ARGUMENTS:
tjd (double)
TT or UT1 Julian date.
direction (short int)
Flag determining 'direction' of transformation;
direction = 0 transformation applied, ITRS to terrestrial
intermediate system
direction != 0 inverse transformation applied, terrestrial
intermediate system to ITRS
xp (double)
Conventionally-defined X coordinate of Celestial Intermediate
Pole with respect to ITRS pole, in arcseconds.
yp (double)
Conventionally-defined Y coordinate of Celestial Intermediate
Pole with respect to ITRS pole, in arcseconds.
pos1[3] (double)
Position vector, geocentric equatorial rectangular coordinates,
referred to ITRS axes.
OUTPUT
ARGUMENTS:
pos2[3] (double)
Position vector, geocentric equatorial rectangular coordinates,
referred to true equator and TIO.
RETURNED
VALUE:
None.
GLOBALS
USED:
T0, ASEC2RAD novascon.c
FUNCTIONS
CALLED:
sin math.h
cos math.h
VER./DATE/
PROGRAMMER:
V1.0/08-93/WTH (USNO/AA) Translate Fortran.
V2.0/10-03/JAB (USNO/AA) Update for IAU 2000 resolutions.
V2.1/01-05/JAB (USNO/AA) Include higher-order terms in the
rotation.
V2.2/06-05/JAB (USNO/AA) Corrected sign of first term in
expression for 'yy' (should be -).
V2.3/03-09/JAB (USNO/AA) Add 'direction' flag and inverse
transformation; conformed input
variables to IERS conventions.
NOTES:
1. This function is the C version of NOVAS Fortran routine
'wobble'.
------------------------------------------------------------------------
*/
{
double xpole, ypole, t, sprime, tiolon, sinx, cosx, siny, cosy, sinl,
cosl, xx, yx, zx, xy, yy, zy, xz, yz, zz;
xpole = xp * ASEC2RAD;
ypole = yp * ASEC2RAD;
/*
Compute approximate longitude of TIO, using eq. (10) of the second
reference.
Note that 'tiolon,' the longitude correction, is negligible for
most astronomical purposes.
*/
t = (tjd - T0) / 36525.0;
sprime = -47.0e-6 * t;
tiolon = -sprime * ASEC2RAD;
/*
Compute elements of rotation matrix.
Equivalent to R3(-s')R2(x)R1(y) as per IERS Conventions (2003).
*/
sinx = sin (xpole);
cosx = cos (xpole);
siny = sin (ypole);
cosy = cos (ypole);
sinl = sin (tiolon);
cosl = cos (tiolon);
xx = cosx * cosl;
yx = sinx * siny * cosl + cosy * sinl;
zx = -sinx * cosy * cosl + siny * sinl;
xy = -cosx * sinl;
yy = -sinx * siny * sinl + cosy * cosl;
zy = sinx * cosy * sinl + siny * cosl;
xz = sinx;
yz = -cosx * siny;
zz = cosx * cosy;
/*
Perform rotation from ITRS to terrestrial intermediate system.
*/
if (direction == 0)
{
pos2[0] = xx * pos1[0] + yx * pos1[1] + zx * pos1[2];
pos2[1] = xy * pos1[0] + yy * pos1[1] + zy * pos1[2];
pos2[2] = xz * pos1[0] + yz * pos1[1] + zz * pos1[2];
}
else
/*
Perform rotation from terrestrial intermediate system to ITRS.
*/
{
pos2[0] = xx * pos1[0] + xy * pos1[1] + xz * pos1[2];
pos2[1] = yx * pos1[0] + yy * pos1[1] + yz * pos1[2];
pos2[2] = zx * pos1[0] + zy * pos1[1] + zz * pos1[2];
}
return;
}
/********terra */
void terra (on_surface *location, double st,
double *pos, double *vel)
/*
------------------------------------------------------------------------
PURPOSE:
Computes the position and velocity vectors of a terrestrial
observer with respect to the center of the Earth.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
INPUT
ARGUMENTS:
*location (struct on_surface)
Pointer to structure containing observer's location (defined
in novas.h).
st (double)
Local apparent sidereal time at reference meridian in hours.
OUTPUT
ARGUMENTS:
pos[3] (double)
Position vector of observer with respect to center of Earth,
equatorial rectangular coordinates, referred to true equator
and equinox of date, components in AU.
vel[3] (double)
Velocity vector of observer with respect to center of Earth,
equatorial rectangular coordinates, referred to true equator
and equinox of date, components in AU/day.
RETURNED
VALUE:
None.
GLOBALS
USED:
AU_KM, ERAD, F novascon.c
ANGVEL, DEG2RAD novascon.c
FUNCTIONS
CALLED:
sin math.h
cos math.h
sqrt math.h
VER./DATE/
PROGRAMMER:
V1.0/04-93/WTH (USNO/AA): Translate Fortran.
V1.1/06-98/JAB (USNO/AA): Move constants 'f' and 'omega' to
file 'novascon.c'.
V1.2/10-03/JAB (USNO/AA): Updates Notes; removed call to 'pow'.
V1.3/12-04/JAB (USNO/AA): Update to use 'on_surface" structure.
V1.4/09-09/WKP (USNO/AA): Moved ht_km calculation from first_entry
block.
NOTES:
1. If reference meridian is Greenwich and st=0, 'pos' is
effectively referred to equator and Greenwich.
2. This function ignores polar motion, unless the
observer's longitude and latitude have been corrected for it,
and variation in the length of day (angular velocity of earth).
3. The true equator and equinox of date do not form an
inertial system. Therefore, with respect to an inertial system,
the very small velocity component (several meters/day) due to
the precession and nutation of the Earth's axis is not accounted
for here.
4. This function is the C version of NOVAS Fortran routine
'terra'.
------------------------------------------------------------------------
*/
{
static short int first_entry = 1;
short int j;
static double erad_km, ht_km;
double df, df2, phi, sinphi, cosphi, c, s, ach, ash, stlocl, sinst,
cosst;
if (first_entry)
{
erad_km = ERAD / 1000.0;
first_entry = 0;
}
/*
Compute parameters relating to geodetic to geocentric conversion.
*/
df = 1.0 - F;
df2 = df * df;
phi = location->latitude * DEG2RAD;
sinphi = sin (phi);
cosphi = cos (phi);
c = 1.0 / sqrt (cosphi * cosphi + df2 * sinphi * sinphi);
s = df2 * c;
ht_km = location->height / 1000.0;
ach = erad_km * c + ht_km;
ash = erad_km * s + ht_km;
/*
Compute local sidereal time factors at the observer's longitude.
*/
stlocl = (st * 15.0 + location->longitude) * DEG2RAD;
sinst = sin (stlocl);
cosst = cos (stlocl);
/*
Compute position vector components in kilometers.
*/
pos[0] = ach * cosphi * cosst;
pos[1] = ach * cosphi * sinst;
pos[2] = ash * sinphi;
/*
Compute velocity vector components in kilometers/sec.
*/
vel[0] = -ANGVEL * ach * cosphi * sinst;
vel[1] = ANGVEL * ach * cosphi * cosst;
vel[2] = 0.0;
/*
Convert position and velocity components to AU and AU/DAY.
*/
for (j = 0; j < 3; j++)
{
pos[j] /= AU_KM;
vel[j] /= AU_KM;
vel[j] *= 86400.0;
}
return;
}
/********e_tilt */
void e_tilt (double jd_tdb, short int accuracy,
double *mobl, double *tobl, double *ee, double *dpsi,
double *deps)
/*
------------------------------------------------------------------------
PURPOSE:
Computes quantities related to the orientation of the Earth's
rotation axis at Julian date 'jd_tdb'.
REFERENCES:
None.
INPUT
ARGUMENTS:
jd_tdb (double)
TDB Julian Date.
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*mobl (double)
Mean obliquity of the ecliptic in degrees at 'jd_tdb'.
*tobl (double)
True obliquity of the ecliptic in degrees at 'jd_tdb'.
*ee (double)
Equation of the equinoxes in seconds of time at 'jd_tdb'.
*dpsi (double)
Nutation in longitude in arcseconds at 'jd_tdb'.
*deps (double)
Nutation in obliquity in arcseconds at 'jd_tdb'.
RETURNED
VALUE:
None.
GLOBALS
USED:
PSI_COR, EPS_COR novas.c
T0, ASEC2RAD novascon.c
DEG2RAD novascon.c
FUNCTIONS
CALLED:
nutation_angles novas.c
ee_ct novas.c
mean_obliq novas.c
fabs math.h
cos math.h
VER./DATE/
PROGRAMMER:
V1.0/08-93/WTH (USNO/AA) Translate Fortran.
V1.1/06-97/JAB (USNO/AA) Incorporate IAU (1994) and IERS (1996)
adjustment to the "equation of the
equinoxes".
V1.2/10-97/JAB (USNO/AA) Implement function that computes
arguments of the nutation series.
V1.3/07-98/JAB (USNO/AA) Use global variables 'PSI_COR' and
'EPS_COR' to apply celestial pole offsets
for high-precision applications.
V2.0/10-03/JAB (USNO/AA) Update function for IAU 2000 resolutions.
V2.1/12-04/JAB (USNO/AA) Add 'mode' argument.
V2.2/01-06/WKP (USNO/AA) Changed 'mode' to 'accuracy'.
NOTES:
1. Values of the celestial pole offsets 'PSI_COR' and 'EPS_COR'
are set using function 'cel_pole', if desired. See the prolog
of 'cel_pole' for details.
2. This function is the C version of NOVAS Fortran routine
'etilt'.
------------------------------------------------------------------------
*/
{
static short int accuracy_last = 0;
short int acc_diff;
static double jd_last = 0.0;
static double dp, de, c_terms;
double t, d_psi, d_eps, mean_ob, true_ob, eq_eq;
/*
Compute time in Julian centuries from epoch J2000.0.
*/
t = (jd_tdb - T0) / 36525.0;
/*
Check for difference in accuracy mode from last call.
*/
acc_diff = accuracy - accuracy_last;
/*
Compute the nutation angles (arcseconds) if the input Julian date
is significantly different from the last Julian date, or the
accuracy mode has changed from the last call.
*/
if (((fabs (jd_tdb - jd_last)) > 1.0e-8) || (acc_diff != 0))
{
nutation_angles (t,accuracy, &dp,&de);
/*
Obtain complementary terms for equation of the equinoxes in
arcseconds.
*/
c_terms = ee_ct (jd_tdb,0.0,accuracy) / ASEC2RAD;
/*
Reset the values of the last Julian date and last mode.
*/
jd_last = jd_tdb;
accuracy_last = accuracy;
}
/*
Apply observed celestial pole offsets.
*/
d_psi = dp + PSI_COR;
d_eps = de + EPS_COR;
/*
Compute mean obliquity of the ecliptic in arcseconds.
*/
mean_ob = mean_obliq (jd_tdb);
/*
Compute true obliquity of the ecliptic in arcseconds.
*/
true_ob = mean_ob + d_eps;
/*
Convert obliquity values to degrees.
*/
mean_ob /= 3600.0;
true_ob /= 3600.0;
/*
Compute equation of the equinoxes in seconds of time.
*/
eq_eq = d_psi * cos (mean_ob * DEG2RAD) + c_terms;
eq_eq /= 15.0;
/*
Set output values.
*/
*dpsi = d_psi;
*deps = d_eps;
*ee = eq_eq;
*mobl = mean_ob;
*tobl = true_ob;
return;
}
/********cel_pole */
short int cel_pole (double tjd, short int type, double dpole1,
double dpole2)
/*
------------------------------------------------------------------------
PURPOSE:
This function allows for the specification of celestial pole
offsets for high-precision applications. Each set of offsets is
a correction to the modeled position of the pole for a specific
date, derived from observations and published by the IERS.
REFERENCES:
Kaplan, G. (2005), US Naval Observatory Circular 179.
Kaplan, G. (2003), USNO/AA Technical Note 2003-03.
INPUT
ARGUMENTS:
tjd (double)
TDB or TT Julian date for pole offsets.
type (short int)
Type of pole offset
= 1 for corrections to angular coordinates of modeled pole
referred to mean ecliptic of date, that is,
delta-delta-psi and delta-delta-epsilon.
= 2 for corrections to components of modeled pole unit
vector referred to GCRS axes, that is, dx and dy.
dpole1 (double)
Value of celestial pole offset in first coordinate,
(delta-delta-psi or dx) in milliarcseconds.
dpole2 (double)
Value of celestial pole offset in second coordinate,
(delta-delta-epsilon or dy) in milliarcseconds.
OUTPUT
ARGUMENTS:
None.
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
= 1 ... Invalid value of 'type'.
GLOBALS
USED:
PSI_COR, EPS_COR novas.c
T0, ASEC2RAD novascon.c
FUNCTIONS
CALLED:
frame_tie novas.c
precession novas.c
mean_obliq novas.c
sin math.h
VER./DATE/
PROGRAMMER:
V1.0/12-04/JAB (USNO/AA)
NOTES:
1. This function sets the values of global variables 'PSI_COR'
and 'EPS_COR' declared at the top of file 'novas.c'. These
global variables are used only in NOVAS function 'e_tilt'.
2. This function, if used, should be called before any other
NOVAS functions for a given date. Values of the pole offsets
specified via a call to this function will be used until
explicitly changed.
3. 'tjd' is used only for 'type' = 2, to transform dx and dy to
the equivalent delta-delta-psi and delta-delta-epsilon values.
4. For 'type' = 2, dx and dy are unit vector component
corrections, but are expressed in milliarcseconds simply by
multiplying by 206264806, the number of milliarcseconds in one
radian.
5. This function is the C version of NOVAS Fortran routine
'celpol'.
------------------------------------------------------------------------
*/
{
short int error = 0;
double dx, dy, t, mean_ob, sin_e, x, dz, dp1[3], dp2[3], dp3[3];
switch (type)
{
case (1):
/*
Angular coordinates of modeled pole referred to mean ecliptic of
date, that is,delta-delta-psi and delta-delta-epsilon.
*/
PSI_COR = dpole1 * 1.0e-3;
EPS_COR = dpole2 * 1.0e-3;
break;
case (2):
/*
Components of modeled pole unit vector referred to GCRS axes, that
is, dx and dy.
*/
dx = dpole1;
dy = dpole2;
t = (tjd - T0) / 36525.0;
/*
Compute sin_e of mean obliquity of date.
*/
mean_ob = mean_obliq (tjd);
sin_e = sin (mean_ob * ASEC2RAD);
/*
The following algorithm, to transform dx and dy to
delta-delta-psi and delta-delta-epsilon, is from eqs. (7)-(9) of the
second reference.
Trivial model of pole trajectory in GCRS allows computation of dz.
*/
x = (2004.190 * t) * ASEC2RAD;
dz = - (x + 0.5 * x * x * x) * dx;
/*
Form pole offset vector (observed - modeled) in GCRS.
*/
dp1[0] = dx * 1.0e-3 * ASEC2RAD;
dp1[1] = dy * 1.0e-3 * ASEC2RAD;
dp1[2] = dz * 1.0e-3 * ASEC2RAD;
/*
Precess pole offset vector to mean equator and equinox of date.
*/
frame_tie (dp1,1, dp2);
precession (T0,dp2,tjd, dp3);
/*
Compute delta-delta-psi and delta-delta-epsilon in arcseconds.
*/
PSI_COR = (dp3[0] / sin_e) / ASEC2RAD;
EPS_COR = dp3[1] / ASEC2RAD;
break;
default:
/*
Invalid value of 'type'.
*/
error = 1;
break;
}
return (error);
}
/********ee_ct */
double ee_ct (double jd_high, double jd_low, short int accuracy)
/*
------------------------------------------------------------------------
PURPOSE:
To compute the "complementary terms" of the equation of the
equinoxes.
REFERENCES:
Capitaine, N., Wallace, P.T., and McCarthy, D.D. (2003). Astron. &
Astrophys. 406, p. 1135-1149. Table 3.
IERS Conventions (2010), Chapter 5, p. 60, Table 5.2e.
(Table 5.2e presented in the printed publication is a truncated
series. The full series, which is used in NOVAS, is available
on the IERS Conventions Center website in file tab5.2e.txt.)
ftp://tai.bipm.org/iers/conv2010/chapter5/
INPUT
ARGUMENTS:
jd_high (double)
High-order part of TT Julian date.
jd_low (double)
Low-order part of TT Julian date.
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
None
RETURNED
VALUE:
(double)
Complementary terms, in radians.
GLOBALS
USED:
T0, ASEC2RAD novascon.c
TWOPI novascon.c
FUNCTIONS
CALLED:
norm_ang novas.c
fund_args novas.c
fmod math.h
sin math.h
cos math.h
VER./DATE/
PROGRAMMER:
V1.0/09-03/JAB (USNO/AA)
V1.1/12-04/JAB (USNO/AA) Added low-accuracy formula.
V1.2/01-06/WKP (USNO/AA) Changed 'mode' to 'accuracy'.
V1.3/11-10/JAB (USNO/AA) Updated reference and notes.
V1.4/03-11/WKP (USNO/AA) Added braces to 2-D array initialization
to quiet gcc warnings.
NOTES:
1. The series used in this function was derived from the first
reference. This same series was also adopted for use in the IAU's
Standards of Fundamental Astronomy (SOFA) software (i.e.,
subroutine eect00.for and function eect00.c).
2. The low-accuracy series used in this function is a simple
implementation derived from the first reference, in which terms
smaller than 2 microarcseconds have been omitted.
3. This function is based on NOVAS Fortran routine 'eect2000',
with the low-accuracy formula taken from NOVAS Fortran routine
'etilt'.
------------------------------------------------------------------------
*/
{
short int i, j;
double t, fa[14], fa2[5], s0, s1, a, c_terms;
/*
Argument coefficients for t^0.
*/
const short int ke0_t[33][14] = {
{0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 2, -2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 2, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 2, -2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 2, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 1, 2, -2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 1, 2, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 4, -4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 1, -1, 1, 0, -8, 12, 0, 0, 0, 0, 0, 0},
{0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 0, 2, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 1, -2, 2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 1, -2, 2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 8,-13, 0, 0, 0, 0, 0, -1},
{0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{2, 0, -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 1, 2, -2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 4, -2, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 2, -2, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 0, -2, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 0, -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
/*
Argument coefficients for t^1.
*/
const short int ke1[14] =
{0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0};
/*
Sine and cosine coefficients for t^0.
*/
const double se0_t[33][2] = {
{+2640.96e-6, -0.39e-6},
{ +63.52e-6, -0.02e-6},
{ +11.75e-6, +0.01e-6},
{ +11.21e-6, +0.01e-6},
{ -4.55e-6, +0.00e-6},
{ +2.02e-6, +0.00e-6},
{ +1.98e-6, +0.00e-6},
{ -1.72e-6, +0.00e-6},
{ -1.41e-6, -0.01e-6},
{ -1.26e-6, -0.01e-6},
{ -0.63e-6, +0.00e-6},
{ -0.63e-6, +0.00e-6},
{ +0.46e-6, +0.00e-6},
{ +0.45e-6, +0.00e-6},
{ +0.36e-6, +0.00e-6},
{ -0.24e-6, -0.12e-6},
{ +0.32e-6, +0.00e-6},
{ +0.28e-6, +0.00e-6},
{ +0.27e-6, +0.00e-6},
{ +0.26e-6, +0.00e-6},
{ -0.21e-6, +0.00e-6},
{ +0.19e-6, +0.00e-6},
{ +0.18e-6, +0.00e-6},
{ -0.10e-6, +0.05e-6},
{ +0.15e-6, +0.00e-6},
{ -0.14e-6, +0.00e-6},
{ +0.14e-6, +0.00e-6},
{ -0.14e-6, +0.00e-6},
{ +0.14e-6, +0.00e-6},
{ +0.13e-6, +0.00e-6},
{ -0.11e-6, +0.00e-6},
{ +0.11e-6, +0.00e-6},
{ +0.11e-6, +0.00e-6}};
/*
Sine and cosine coefficients for t^1.
*/
const double se1[2] =
{ -0.87e-6, +0.00e-6};
/*
Interval between fundamental epoch J2000.0 and current date.
*/
t = ((jd_high - T0) + jd_low) / 36525.0;
/*
High accuracy mode.
*/
if (accuracy == 0)
{
/*
Fundamental Arguments.
Mean Anomaly of the Moon.
*/
fa[0] = norm_ang ((485868.249036 +
(715923.2178 +
( 31.8792 +
( 0.051635 +
( -0.00024470)
* t) * t) * t) * t) * ASEC2RAD
+ fmod (1325.0*t, 1.0) * TWOPI);
/*
Mean Anomaly of the Sun.
*/
fa[1] = norm_ang ((1287104.793048 +
(1292581.0481 +
( -0.5532 +
( +0.000136 +
( -0.00001149)
* t) * t) * t) * t) * ASEC2RAD
+ fmod (99.0*t, 1.0) * TWOPI);
/*
Mean Longitude of the Moon minus Mean Longitude of the Ascending
Node of the Moon.
*/
fa[2] = norm_ang (( 335779.526232 +
( 295262.8478 +
( -12.7512 +
( -0.001037 +
( 0.00000417)
* t) * t) * t) * t) * ASEC2RAD
+ fmod (1342.0*t, 1.0) * TWOPI);
/*
Mean Elongation of the Moon from the Sun.
*/
fa[3] = norm_ang ((1072260.703692 +
(1105601.2090 +
( -6.3706 +
( 0.006593 +
( -0.00003169)
* t) * t) * t) * t) * ASEC2RAD
+ fmod (1236.0*t, 1.0) * TWOPI);
/*
Mean Longitude of the Ascending Node of the Moon.
*/
fa[4] = norm_ang (( 450160.398036 +
(-482890.5431 +
( 7.4722 +
( 0.007702 +
( -0.00005939)
* t) * t) * t) * t) * ASEC2RAD
+ fmod (-5.0*t, 1.0) * TWOPI);
fa[ 5] = norm_ang (4.402608842 + 2608.7903141574 * t);
fa[ 6] = norm_ang (3.176146697 + 1021.3285546211 * t);
fa[ 7] = norm_ang (1.753470314 + 628.3075849991 * t);
fa[ 8] = norm_ang (6.203480913 + 334.0612426700 * t);
fa[ 9] = norm_ang (0.599546497 + 52.9690962641 * t);
fa[10] = norm_ang (0.874016757 + 21.3299104960 * t);
fa[11] = norm_ang (5.481293872 + 7.4781598567 * t);
fa[12] = norm_ang (5.311886287 + 3.8133035638 * t);
fa[13] = (0.024381750 + 0.00000538691 * t) * t;
/*
Evaluate the complementary terms.
*/
s0 = 0.0;
s1 = 0.0;
for (i = 32; i >= 0; i--)
{
a = 0.0;
for (j = 0; j < 14; j++)
{
a += (double) ke0_t[i][j] * fa[j];
}
s0 += (se0_t[i][0] * sin (a) + se0_t[i][1] * cos (a));
}
a = 0.0;
for (j = 0; j < 14; j++)
{
a += (double) (ke1[j]) * fa[j];
}
s1 += (se1[0] * sin (a) + se1[1] * cos (a));
c_terms = (s0 + s1 * t);
}
else
/*
Low accuracy mode: Terms smaller than 2 microarcseconds omitted.
*/
{
fund_args (t, fa2);
c_terms =
2640.96e-6 * sin (fa2[4])
+ 63.52e-6 * sin (2.0 * fa2[4])
+ 11.75e-6 * sin (2.0 * fa2[2] - 2.0 * fa2[3] + 3.0 * fa2[4])
+ 11.21e-6 * sin (2.0 * fa2[2] - 2.0 * fa2[3] + fa2[4])
- 4.55e-6 * sin (2.0 * fa2[2] - 2.0 * fa2[3] + 2.0 * fa2[4])
+ 2.02e-6 * sin (2.0 * fa2[2] + 3.0 * fa2[4])
+ 1.98e-6 * sin (2.0 * fa2[2] + fa2[4])
- 1.72e-6 * sin (3.0 * fa2[4])
- 0.87e-6 * t * sin (fa2[4]);
}
return (c_terms *= ASEC2RAD);
}
/********frame_tie */
void frame_tie (double *pos1, short int direction,
double *pos2)
/*
------------------------------------------------------------------------
PURPOSE:
To transform a vector from the dynamical reference system to the
International Celestial Reference System (ICRS), or vice versa.
The dynamical reference system is based on the dynamical mean
equator and equinox of J2000.0. The ICRS is based on the space-
fixed ICRS axes defined by the radio catalog positions of several
hundred extragalactic objects.
REFERENCES:
Hilton, J. and Hohenkerk, C. (2004), Astronomy and Astrophysics
413, 765-770, eq. (6) and (8).
IERS (2003) Conventions, Chapter 5.
INPUT
ARGUMENTS:
pos1[3] (double)
Position vector, equatorial rectangular coordinates.
direction (short int)
Set 'direction' < 0 for dynamical to ICRS transformation.
Set 'direction' >= 0 for ICRS to dynamical transformation.
OUTPUT
ARGUMENTS:
pos2[3] (double)
Position vector, equatorial rectangular coordinates.
RETURNED
VALUE:
None.
GLOBALS
USED:
ASEC2RAD novascon.c
FUNCTIONS
CALLED:
None.
VER./DATE/
PROGRAMMER:
V1.0/09-03/JAB (USNO/AA)
V1.1/02-06/WKP (USNO/AA) Added second-order corrections to diagonal
elements.
NOTES:
1. For geocentric coordinates, the same transformation is used
between the dynamical reference system and the GCRS.
2. This function is the C version of NOVAS Fortran routine 'frame'.
------------------------------------------------------------------------
*/
{
static short int compute_matrix = 1;
/*
'xi0', 'eta0', and 'da0' are ICRS frame biases in arcseconds taken
from IERS (2003) Conventions, Chapter 5.
*/
const double xi0 = -0.0166170;
const double eta0 = -0.0068192;
const double da0 = -0.01460;
static double xx, yx, zx, xy, yy, zy, xz, yz, zz;
/*
Compute elements of rotation matrix to first order the first time
this function is called. Elements will be saved for future use and
not recomputed.
*/
if (compute_matrix == 1)
{
xx = 1.0;
yx = -da0 * ASEC2RAD;
zx = xi0 * ASEC2RAD;
xy = da0 * ASEC2RAD;
yy = 1.0;
zy = eta0 * ASEC2RAD;
xz = -xi0 * ASEC2RAD;
yz = -eta0 * ASEC2RAD;
zz = 1.0;
/*
Include second-order corrections to diagonal elements.
*/
xx = 1.0 - 0.5 * (yx * yx + zx * zx);
yy = 1.0 - 0.5 * (yx * yx + zy * zy);
zz = 1.0 - 0.5 * (zy * zy + zx * zx);
compute_matrix = 0;
}
/*
Perform the rotation in the sense specified by 'direction'.
*/
if (direction < 0)
{
/*
Perform rotation from dynamical system to ICRS.
*/
pos2[0] = xx * pos1[0] + yx * pos1[1] + zx * pos1[2];
pos2[1] = xy * pos1[0] + yy * pos1[1] + zy * pos1[2];
pos2[2] = xz * pos1[0] + yz * pos1[1] + zz * pos1[2];
}
else
{
/*
Perform rotation from ICRS to dynamical system.
*/
pos2[0] = xx * pos1[0] + xy * pos1[1] + xz * pos1[2];
pos2[1] = yx * pos1[0] + yy * pos1[1] + yz * pos1[2];
pos2[2] = zx * pos1[0] + zy * pos1[1] + zz * pos1[2];
}
return;
}
/********proper_motion */
void proper_motion (double jd_tdb1, double *pos, double *vel,
double jd_tdb2,
double *pos2)
/*
------------------------------------------------------------------------
PURPOSE:
Applies proper motion, including foreshortening effects, to a
star's position.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
INPUT
ARGUMENTS:
jd_tdb1 (double)
TDB Julian date of first epoch.
pos[3] (double)
Position vector at first epoch.
vel[3] (double)
Velocity vector at first epoch.
jd_tdb2 (double)
TDB Julian date of second epoch.
OUTPUT
ARGUMENTS:
pos2[3] (double)
Position vector at second epoch.
RETURNED
VALUE:
None.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
None.
VER./DATE/
PROGRAMMER:
V1.0/01-93/TKB (USNO/NRL Optical Interfer.) Translate Fortran.
V1.1/08-93/WTH (USNO/AA) Updated to C programming standards.
NOTES:
1. This function is the C version of NOVAS Fortran routine
'propmo'.
------------------------------------------------------------------------
*/
{
short int j;
for (j = 0; j < 3; j++)
{
pos2[j] = pos[j] + (vel[j] * (jd_tdb2 - jd_tdb1));
}
return;
}
/********bary2obs */
void bary2obs (double *pos, double *pos_obs,
double *pos2, double *lighttime)
/*
------------------------------------------------------------------------
PURPOSE:
This function moves the origin of coordinates from the barycenter
of the solar system to the observer (or the geocenter); i.e.,
this function accounts for parallax (annual+geocentric or just
annual).
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
INPUT
ARGUMENTS:
pos1[3] (double)
Position vector, referred to origin at solar system barycenter,
components in AU.
pos_obs[3] (double)
Position vector of observer (or the geocenter), with respect
to origin at solar system barycenter, components in AU.
OUTPUT
ARGUMENTS:
pos2[3] (double)
Position vector, referred to origin at center of mass of the
Earth, components in AU.
*lighttime (double)
Light time from object to Earth in days.
RETURNED
VALUE:
None.
GLOBALS
USED:
C_AUDAY novascon.c
FUNCTIONS
CALLED:
sqrt math.h
VER./DATE/
PROGRAMMER:
V1.0/01-93/TKB (USNO/NRL Optical Interfer.) Translate Fortran.
V1.1/08-93/WTH (USNO/AA) Update to C Standards.
V1.2/10-03/JAB (USNO/AA) Remove calls to 'pow'.
V2.0/01-05/JAB (USNO/AA) Update documentation; change
'earthvector' to 'pos_obs'.
NOTES:
1. This function is the C version of NOVAS Fortran routine
'geocen'.
------------------------------------------------------------------------
*/
{
short int j;
/*
Translate vector to geocentric coordinates.
*/
for (j = 0; j < 3; j++)
{
pos2[j] = pos[j] - pos_obs[j];
}
/*
Calculate length of vector in terms of light time.
*/
*lighttime = sqrt (pos2[0] * pos2[0] + pos2[1] * pos2[1] +
pos2[2] * pos2[2]) / C_AUDAY;
return;
}
/********geo_posvel */
short int geo_posvel (double jd_tt, double delta_t, short int accuracy,
observer *obs,
double *pos, double *vel)
/*
------------------------------------------------------------------------
PURPOSE:
This function computes the geocentric position and velocity
of an observer on the surface of the earth or on a near-earth
spacecraft. The final vectors are expressed in the GCRS.
REFERENCES:
None.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date.
delta_t (double)
Value of Delta T (= TT - UT1) at 'jd_tt'.
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
*obs (struct observer)
Data specifying the location of the observer (struct defined
in novas.h).
OUTPUT
ARGUMENTS:
*pos (double)
Position vector of observer, with respect to origin at
geocenter, referred to GCRS axes, components in AU.
*vel (double)
Velocity vector of observer, with respect to origin at
geocenter, referred to GCRS axes, components in AU/day.
RETURNED
VALUE:
(short int)
= 0 ... everything OK.
= 1 ... invalid value of 'accuracy'.
GLOBALS
USED:
AU_KM, T0 novascon.c
FUNCTIONS
CALLED:
tdb2tt novas.c
sidereal_time novas.c
e_tilt novas.c
terra novas.c
nutation novas.c
precession novas.c
frame_tie novas.c
fabs math.h
VER./DATE/
PROGRAMMER:
V1.0/12-04/JAB (USNO/AA)
V2.2/12-05/WKP (USNO/AA) Updated error handling.
V2.3/01-06/WKP (USNO/AA) Changed 'mode' to 'accuracy'.
V2.4/02-06/WKP (USNO/AA) Changed time argument in sidereal_time
call from 'jd_tt' to 'jd_ut1'.
V2.5/04-06/WKP (USNO/AA) Corrected calculation of 'fac' and
changed it to a static.
V2.6/10-06/JAB (USNO/AA) Remove 'jd_tdb' input argument; TDB
is approximated by TT.
V2.7/02-07/JAB (USNO/AA) Compute 'jd_tdb' corresponding to input
'jd_tt'.
NOTES:
1. This function is the C version of NOVAS Fortran routine
'geopos'. This C function differs from its Fortran counterpart
in that it does not use the input TT time as a substitute for
the TDB time.
------------------------------------------------------------------------
*/
{
static double t_last = 0;
static double gast, fac;
static short int first_time = 1;
double x, secdif, gmst, x1, x2, x3, x4, eqeq, pos1[3], vel1[3],
pos2[3], vel2[3], pos3[3], vel3[3], jd_tdb, jd_ut1;
short int error = 0;
/*
Invalid value of 'accuracy'.
*/
if ((accuracy < 0) || (accuracy > 1))
return (error = 1);
/*
Compute 'jd_tdb', the TDB Julian date corresponding to 'jd_tt'.
*/
jd_tdb = jd_tt;
tdb2tt (jd_tdb, &x,&secdif);
jd_tdb = jd_tt + secdif / 86400.0;
switch (obs->where)
{
/*
Observer at geocenter. Trivial case.
*/
case (0):
pos[0] = 0.0;
pos[1] = 0.0;
pos[2] = 0.0;
vel[0] = 0.0;
vel[1] = 0.0;
vel[2] = 0.0;
return (error = 0);
break;
/*
Other two cases: Get geocentric position and velocity vectors of
observer wrt equator and equinox of date.
Observer on surface of Earth.
*/
case (1):
/*
Compute UT1 and sidereal time.
*/
jd_ut1 = jd_tt - (delta_t / 86400.0);
if (fabs (jd_ut1 - t_last) > 1.0e-8 )
{
sidereal_time (jd_ut1,0.0,delta_t,0,1,accuracy, &gmst);
e_tilt (jd_tdb,accuracy, &x1,&x2,&eqeq,&x3,&x4);
gast = gmst + eqeq / 3600.0;
t_last = jd_ut1;
}
/*
Function 'terra' does the hard work, given sidereal time.
*/
terra (&obs->on_surf,gast, pos1,vel1);
break;
/*
Observer on near-earth spacecraft.
*/
case (2):
/*
Convert units to AU and AU/day.
*/
if (first_time)
{
fac = AU_KM / 86400.0;
first_time = 0;
}
pos1[0] = obs->near_earth.sc_pos[0] / AU_KM;
pos1[1] = obs->near_earth.sc_pos[1] / AU_KM;
pos1[2] = obs->near_earth.sc_pos[2] / AU_KM;
vel1[0] = obs->near_earth.sc_vel[0] / fac;
vel1[1] = obs->near_earth.sc_vel[1] / fac;
vel1[2] = obs->near_earth.sc_vel[2] / fac;
break;
}
/*
Transform geocentric position vector of observer to GCRS.
*/
nutation (jd_tdb,-1,accuracy,pos1, pos2);
precession (jd_tdb,pos2,T0, pos3);
frame_tie (pos3,-1, pos);
/*
Transform geocentric velocity vector of observer to GCRS.
*/
nutation (jd_tdb,-1,accuracy,vel1, vel2);
precession (jd_tdb,vel2,T0, vel3);
frame_tie (vel3,-1, vel);
return (error = 0);
}
/********light_time */
short int light_time (double jd_tdb, object *ss_object, double pos_obs[3],
double tlight0, short int accuracy,
double pos[3], double *tlight)
/*
------------------------------------------------------------------------
PURPOSE:
This function computes the geocentric position of a solar system
body, as antedated for light-time.
REFERENCES:
None.
INPUT
ARGUMENTS:
jd_tdb (double)
TDB Julian date of observation.
*ss_object (struct object)
Pointer to structure containing the designation for the
solar system body (defined in novas.h).
pos_obs[3] (double)
Position vector of observer (or the geocenter), with respect to
origin at solar system barycenter, referred to ICRS axes,
components in AU.
tlight0 (double)
First approximation to light-time, in days (can be set to 0.0
if unknown).
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
pos[3] (double)
Position vector of body, with respect to origin at observer (or
the geocenter), referred to ICRS axes, components in AU.
tlight (double)
Final light-time, in days.
RETURNED
VALUE:
(short int)
= 0 ... everything OK.
= 1 ... algorithm failed to converge after 10 iterations.
> 10 ... error is 10 + error from function 'solarsystem'.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
ephemeris novas.c
bary2obs novas.c
fabs math.h
VER./DATE/
PROGRAMMER:
V1.0/09-04/JAB (USNO/AA)
V1.1/06-05/JAB (USNO/AA): Changed the convergence tolerance from
1.0e-8 to 1.0e-10
V1.2/04-06/JAB (USNO/AA): Changed the convergence tolerance from
1.0e-10 to 1.0e-9
V2.0/12-06/JAB (USNO/AA): Support use of split Julian date to
improve accuracy when the full accuracy
option is selected.
V2.1/11-07/JAB (USNO/AA): Changed high-accuracy convergence
tolerance from 1.0e-11 to 1.0e-12.
V2.2/10-08/JAB (USNO/AA): Use function 'ephemeris' instead of
direct calls to the solar system
functions.
V2.3/09-10/WKP (USNO/AA): Initialized 't3' variable to silence
compiler warning.
NOTES:
1. This function is the C version of NOVAS Fortran routine
'littim'.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int iter = 0;
double tol, jd[2], t1, t2, t3 = 0.0, pos1[3], vel1[3];
/*
Set light-time convergence tolerance. If full-accuracy option has
been selected, split the Julian date into whole days + fraction of
day.
*/
if (accuracy == 0)
{
tol = 1.0e-12;
jd[0] = (double) ((long int) jd_tdb);
t1 = jd_tdb - jd[0];
t2 = t1 - tlight0;
}
else
{
tol = 1.0e-9;
jd[0] = 0.0;
jd[1] = 0.0;
t1 = jd_tdb;
t2 = jd_tdb - tlight0;
}
/*
Iterate to obtain correct light-time (usually converges rapidly).
*/
do
{
if (iter > 10)
{
error = 1;
*tlight = 0.0;
break;
}
if (iter > 0)
{
t2 = t3;
}
jd[1] = t2;
error = ephemeris (jd,ss_object,0,accuracy, pos1, vel1);
if (error != 0)
{
error += 10;
*tlight = 0.0;
break;
}
bary2obs (pos1,pos_obs, pos,tlight);
t3 = t1 - *tlight;
iter++;
} while (fabs (t3 - t2) > tol);
return (error);
}
/********d_light */
double d_light (double *pos1, double *pos_obs)
/*
------------------------------------------------------------------------
PURPOSE:
This function returns the difference in light-time, for a star,
between the barycenter of the solar system and the observer (or
the geocenter).
REFERENCES:
None.
INPUT
ARGUMENTS:
pos1[3] (double)
Position vector of star, with respect to origin at solar
system barycenter.
pos_obs[3] (double)
Position vector of observer (or the geocenter), with respect
to origin at solar system barycenter, components in AU.
OUTPUT
ARGUMENTS:
None.
RETURNED
VALUE:
(double)
Difference in light time, in the sense star to barycenter minus
star to earth, in days.
GLOBALS
USED:
C_AUDAY novascon.c
FUNCTIONS
CALLED:
sqrt math.h
VER./DATE/
PROGRAMMER:
V1.0/12-04/JAB (USNO/AA)
NOTES:
1. Alternatively, this function returns the light-time from the
observer (or the geocenter) to a point on a light ray that is
closest to a specific solar system body. For this purpose, 'pos1'
is the position vector toward observed object, with respect to
origin at observer (or the geocenter); 'pos_obs' is the position
vector of solar system body, with respect to origin at observer
(or the geocenter), components in AU; and the returned value is
the light time to point on line defined by 'pos1' that is closest
to solar system body (positive if light passes body before hitting
observer, i.e., if 'pos1' is within 90 degrees of 'pos_obs').
2. This function is the C version of NOVAS Fortran routine
'dlight'.
------------------------------------------------------------------------
*/
{
double dis, u1[3], diflt;
/*
From 'pos1', form unit vector 'u1' in direction of star or light
source.
*/
dis = sqrt (pos1[0] * pos1[0] + pos1[1] * pos1[1] +
pos1[2] * pos1[2]);
u1[0] = pos1[0] / dis;
u1[1] = pos1[1] / dis;
u1[2] = pos1[2] / dis;
/*
Light-time returned is the projection of vector 'pos_obs' onto the
unit vector 'u1' (formed from 'pos1'), divided by the speed of light.
*/
diflt = (pos_obs[0] * u1[0] + pos_obs[1] * u1[1] +
pos_obs[2] * u1[2]) / C_AUDAY;
return (diflt);
}
/********grav_def */
short int grav_def (double jd_tdb, short int loc_code,
short int accuracy, double *pos1, double *pos_obs,
double *pos2)
/*
------------------------------------------------------------------------
PURPOSE:
This function computes the total gravitational deflection of
light for the observed object due to the major gravitating bodies
in the solar system. This function valid for an observed body
within the solar system as well as for a star.
REFERENCES:
Klioner, S. (2003), Astronomical Journal 125, 1580-1597,
Section 6.
INPUT
ARGUMENTS:
jd_tdb (double)
TDB Julian date of observation.
loc_code (short int)
Code for location of observer, determining whether the
gravitational deflection due to the earth itself is applied.
= 0 ... No earth deflection (normally means observer
is at geocenter)
= 1 ... Add in earth deflection (normally means
observer is on or above surface of earth, including earth
orbit)
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
*pos1 (double)
Position vector of observed object, with respect to origin at
observer (or the geocenter), referred to ICRS axes, components
in AU.
*pos_obs (double)
Position vector of observer (or the geocenter), with respect to
origin at solar system barycenter, referred to ICRS axes,
components in AU.
OUTPUT
ARGUMENTS:
*pos2 (double)
Position vector of observed object, with respect to origin at
observer (or the geocenter), referred to ICRS axes, corrected
for gravitational deflection, components in AU.
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
< 30 ... Error from function 'ephemeris'.
> 30 ... Error from function 'make_object'.
GLOBALS
USED:
C_AUDAY novascon.c
RMASS[12] novascon.c
FUNCTIONS
CALLED:
make_cat_entry novas.c
make_object novas.c
ephemeris novas.c
bary2obs novas.c
d_light novas.c
grav_vec novas.c
sqrt math.h
VER./DATE/
PROGRAMMER:
V1.0/01-05/JAB (USNO/AA)
V1.1/04-06/WKP (USNO/AA): Added second note (superceded).
V1.2/10-08/JAB (USNO/AA): Substituted calls to 'ephemeris' for
calls to 'solarsystem'; added
actions based on 'accuracy' input.
NOTES:
1. This function is the C version of NOVAS Fortran routine
'grvdef'.
2. If 'accuracy' is set to zero (full accuracy), three bodies
(Sun, Jupiter, and Saturn) are used in the calculation. If
the reduced-accuracy option is set, only the Sun is used in the
calculation. In both cases, if the observer is not at the
geocenter, the deflection due to the Earth is included.
3. The number of bodies used at full and reduced accuracy can be
set by making a change to the code in this function as indicated
in the comments.
------------------------------------------------------------------------
*/
{
/*
The following list of body numbers and corresponding names identifies
which gravitating bodies (aside from the Earth) are potentially
used -- list is taken from Klioner's table 1, the order based on area
of sky affected (col 2). Order is Sun, Jupiter, Saturn, Moon, Venus,
Uranus, Neptune.
*/
char body_name[7][8] = {"Sun", "Jupiter", "Saturn", "Moon",
"Venus", "Uranus", "Neptune"};
const short int body_num[7] = {10, 5, 6, 11, 2, 7, 8};
static short int first_time = 1;
static short int nbodies_last = 0;
short int error = 0;
short int nbodies, i;
double jd[2], tlt, pbody[3], vbody[3], pbodyo[3], x, dlt, tclose;
cat_entry dummy_star;
static object body[7], earth;
jd[1] = 0.0;
/*
Initialize output vector of observed object to equal input vector.
*/
for (i = 0; i < 3; i++)
{
pos2[i] = pos1[i];
}
/*
Set the number of bodies -- and hence the bodies used -- based on the
value of the 'accuracy' flag.
Change value of 'nbodies' to include or exclude gravitating bodies
('nbodies' <= 0 means no deflection calculated, 'nbodies' = 1 means
Sun only, 'nbodies' = 2 means Sun + Jupiter, etc.) Default is
'nbodies' = 3: Sun + Jupiter + Saturn.
*/
if (accuracy == 0)
{
nbodies = 3;
}
else
{
nbodies = 1;
}
/*
Set up the structures of type 'object' containing the body
information.
*/
if ((first_time == 1) || (nbodies != nbodies_last))
{
for (i = 0; i < nbodies; i++)
{
if (i == 0)
{
make_cat_entry ("dummy"," ",0,0.0,0.0,0.0,0.0,0.0,0.0,
&dummy_star);
make_object (0,3,"Earth",&dummy_star, &earth);
}
if ((error = make_object (0,body_num[i],body_name[i],
&dummy_star, &body[i])) != 0)
{
return (error += 30);
}
}
first_time = 0;
nbodies_last = nbodies;
}
/*
Compute light-time to observed object.
*/
tlt = sqrt (pos1[0] * pos1[0] + pos1[1] * pos1[1] +
pos1[2] *pos1[2]) / C_AUDAY;
/*
Cycle through gravitating bodies.
*/
for (i = 0; i < nbodies; i++)
{
/*
Get position of gravitating body wrt ss barycenter at time 'jd_tdb'.
*/
jd[0] = jd_tdb;
if ((error = ephemeris (jd,&body[i],0,accuracy, pbody,vbody))
!= 0)
{
return (error);
}
/*
Get position of gravitating body wrt observer at time 'jd_tdb'.
*/
bary2obs (pbody,pos_obs, pbodyo,&x);
/*
Compute light-time from point on incoming light ray that is closest
to gravitating body.
*/
dlt = d_light (pos2,pbodyo);
/*
Get position of gravitating body wrt ss barycenter at time when
incoming photons were closest to it.
*/
tclose = jd_tdb;
if (dlt > 0.0)
tclose = jd_tdb - dlt;
if (tlt < dlt)
tclose = jd_tdb - tlt;
jd[0] = tclose;
if ((error = ephemeris (jd,&body[i],0,accuracy, pbody,vbody))
!= 0)
{
return (error);
}
/*
Compute deflection due to gravitating body.
*/
grav_vec (pos2,pos_obs,pbody,RMASS[body_num[i]], pos2);
}
/*
If observer is not at geocenter, add in deflection due to Earth.
*/
if (loc_code != 0 )
{
/*
Get position of Earth wrt solar system barycenter at time 'jd_tdb'.
*/
jd[0] = jd_tdb;
if ((error = ephemeris (jd,&earth,0,accuracy, pbody,vbody))
!= 0)
{
return (error);
}
/*
Compute deflection due to Earth.
*/
grav_vec (pos2,pos_obs,pbody,RMASS[3], pos2);
}
return (error);
}
/********grav_vec */
void grav_vec (double *pos1, double *pos_obs, double *pos_body,
double rmass,
double *pos2)
/*
------------------------------------------------------------------------
PURPOSE:
This function corrects position vector for the deflection
of light in the gravitational field of an arbitrary body. This
function valid for an observed body within the solar system as
well as for a star.
REFERENCES:
Murray, C.A. (1981) Mon. Notices Royal Ast. Society 195, 639-648.
See also formulae in Section B of the Astronomical Almanac, or
Kaplan, G. et al. (1989) Astronomical Journal 97, 1197-1210,
section iii f.
INPUT
ARGUMENTS:
pos1[3] (double)
Position vector of observed object, with respect to origin at
observer (or the geocenter), components in AU.
pos_obs [3] (double)
Position vector of observer (or the geocenter), with respect
to origin at solar system barycenter, components in AU.
pos_body[3] (double)
Position vector of gravitating body, with respect to origin
at solar system barycenter, components in AU.
rmass (double)
Reciprocal mass of gravitating body in solar mass units,
that is, Sun mass / body mass.
OUTPUT
ARGUMENTS:
pos2[3] (double)
Position vector of observed object, with respect to origin at
observer (or the geocenter), corrected for gravitational
deflection, components in AU.
RETURNED
VALUE:
None.
GLOBALS
USED:
C, AU, GS novascon.c
FUNCTIONS
CALLED:
sqrt math.h
fabs math.h
VER./DATE/
PROGRAMMER:
V1.0/12-04/JAB (USNO/AA)
NOTES:
1. This function is the C version of NOVAS Fortran routine
'grvd'.
------------------------------------------------------------------------
*/
{
short int i;
double pq[3], pe[3], pmag, emag, qmag, phat[3], ehat[3], qhat[3],
pdotq, edotp, qdote, fac1, fac2, p2i;
/*
Construct vector 'pq' from gravitating body to observed object and
construct vector 'pe' from gravitating body to observer.
*/
for (i = 0; i < 3; i++)
{
pq[i] = pos_obs[i] + pos1[i] - pos_body[i];
pe[i] = pos_obs[i] - pos_body[i];
}
/*
Compute vector magnitudes and unit vectors.
*/
pmag = sqrt (pos1[0] * pos1[0] + pos1[1] * pos1[1] +
pos1[2] * pos1[2]);
emag = sqrt (pe[0] * pe[0] + pe[1] * pe[1] + pe[2] * pe[2]);
qmag = sqrt (pq[0] * pq[0] + pq[1] * pq[1] + pq[2] * pq[2]);
for (i = 0; i < 3; i++)
{
phat[i] = pos1[i] / pmag;
ehat[i] = pe[i] / emag;
qhat[i] = pq[i] / qmag;
}
/*
Compute dot products of vectors.
*/
pdotq = phat[0] * qhat[0] + phat[1] * qhat[1] + phat[2] * qhat[2];
edotp = ehat[0] * phat[0] + ehat[1] * phat[1] + ehat[2] * phat[2];
qdote = qhat[0] * ehat[0] + qhat[1] * ehat[1] + qhat[2] * ehat[2];
/*
If gravitating body is observed object, or is on a straight line
toward or away from observed object to within 1 arcsec, deflection
is set to zero; set 'pos2' equal to 'pos1'.
*/
if (fabs (edotp) > 0.99999999999)
{
for (i = 0; i < 3; i++)
{
pos2[i] = pos1[i];
}
}
else
{
/*
Compute scalar factors.
*/
fac1 = 2.0 * GS / (C * C * emag * AU * rmass);
fac2 = 1.0 + qdote;
/*
Construct corrected position vector 'pos2'.
*/
for (i = 0; i < 3; i++)
{
p2i = phat[i] + fac1 * (pdotq * ehat[i] - edotp * qhat[i]) /
fac2;
pos2[i] = p2i * pmag;
}
}
return;
}
/********aberration */
void aberration (double *pos, double *ve, double lighttime,
double *pos2)
/*
------------------------------------------------------------------------
PURPOSE:
Corrects position vector for aberration of light. Algorithm
includes relativistic terms.
REFERENCES:
Murray, C. A. (1981) Mon. Notices Royal Ast. Society 195, 639-648.
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
INPUT
ARGUMENTS:
pos[3] (double)
Position vector, referred to origin at center of mass of the
Earth, components in AU.
ve[3] (double)
Velocity vector of center of mass of the Earth, referred to
origin at solar system barycenter, components in AU/day.
lighttime (double)
Light time from object to Earth in days.
OUTPUT
ARGUMENTS:
pos2[3] (double)
Position vector, referred to origin at center of mass of the
Earth, corrected for aberration, components in AU
RETURNED
VALUE:
None.
GLOBALS
USED:
C_AUDAY novascon.c
FUNCTIONS
CALLED:
sqrt math.h
VER./DATE/
PROGRAMMER:
V1.0/01-93/TKB (USNO/NRL Optical Interfer.) Translate Fortran.
V1.1/08-93/WTH (USNO/AA) Update to C Standards.
V1.2/12-02/JAB (USNO/AA) Improve efficiency by removing calls to
'pow' and replacing final 'for' loop
with direct assignment.
V1.3/11-03/JAB (USNO/AA) Remove returned value.
V1.4/02-06/WKP (USNO/AA) Changed C to C_AUDAY.
NOTES:
1. This function is the C version of NOVAS Fortran routine
'aberat'.
2. If 'lighttime' = 0 on input, this function will compute it.
------------------------------------------------------------------------
*/
{
double p1mag, vemag, beta, dot, cosd, gammai, p, q, r;
if (lighttime == 0.0)
{
p1mag = sqrt (pos[0] * pos[0] + pos[1] * pos[1] + pos[2] *
pos[2]);
lighttime = p1mag / C_AUDAY;
}
else
p1mag = lighttime * C_AUDAY;
vemag = sqrt (ve[0] * ve[0] + ve[1] * ve[1] + ve[2] * ve[2]);
beta = vemag / C_AUDAY;
dot = pos[0] * ve[0] + pos[1] * ve[1] + pos[2] * ve[2];
cosd = dot / (p1mag * vemag);
gammai = sqrt (1.0 - beta * beta);
p = beta * cosd;
q = (1.0 + p / (1.0 + gammai)) * lighttime;
r = 1.0 + p;
pos2[0] = (gammai * pos[0] + q * ve[0]) / r;
pos2[1] = (gammai * pos[1] + q * ve[1]) / r;
pos2[2] = (gammai * pos[2] + q * ve[2]) / r;
return;
}
/********rad_vel */
void rad_vel (object *cel_object, double *pos, double *vel,
double *vel_obs, double d_obs_geo, double d_obs_sun,
double d_obj_sun,
double *rv)
/*
------------------------------------------------------------------------
PURPOSE:
Predicts the radial velocity of the observed object
as it would be measured by spectroscopic means. Radial
velocity is here defined as the radial velocity measure (z)
times the speed of light. For a solar system body, it applies
to a fictitious emitter at the center of the observed object,
assumed massless (no gravitational red shift), and does not
in general apply to reflected light. For stars, it includes
all effects, such as gravitational red shift, contained
in the catalog barycentric radial velocity measure, a scalar
derived from spectroscopy. Nearby stars with a known kinematic
velocity vector (obtained independently of spectroscopy) can be
treated like solar system objects.
REFERENCES:
Lindegren & Dravins (2003), Astronomy & Astrophysics 401,
1185-1201.
INPUT
ARGUMENTS:
*cel_object (struct object)
Specifies the celestial object of interest (structure defined
in novas.h).
*pos (double)
Geometric position vector of object with respect to observer,
corrected for light-time, in AU.
*vel (double)
Velocity vector of object with respect to solar system
barycenter, in AU/day.
*vel_obs (double)
Velocity vector of observer with respect to solar system
barycenter, in AU/day.
d_obs_geo (double)
Distance from observer to geocenter, in AU.
d_obs_sun (double)
Distance from observer to Sun, in AU.
d_obj_sun (double)
Distance from object to Sun, in AU.
OUTPUT
ARGUMENTS:
*rv (double)
The observed radial velocity measure times the speed of light,
in kilometers/second.
RETURNED
VALUE:
None.
GLOBALS
USED:
AU, C, GS, GE novascon.c
DEG2RAD novascon.c
FUNCTIONS
CALLED:
sqrt math.h
sin math.h
cos math.h
VER./DATE/
PROGRAMMER:
V1.0/09-06/JAB (USNO/AA)
V2.0/01-07/JAB (USNO/AA): Implement new algorithm
NOTES:
1. All the input arguments are BCRS quantities, expressed
with respect to the ICRS axes. 'vel' and 'vel_obs' are kinematic
velocities - derived from geometry or dynamics, not spectroscopy.
2. If the object is outside the solar system, the algorithm
used will be consistent with the IAU definition of stellar
radial velocity, specifically, the barycentric radial velocity
measure, which is derived from spectroscopy. In that case,
the vector 'vel' can be very approximate -- or, for distant stars
or galaxies, zero -- as it will be used only for a small geometric
correction that is proportional to proper motion.
3. Any of the distances (last three input arguments) can be set to
zero (0.0) if the corresponding general relativistic gravitational
potential term is not to be evaluated. These terms generally
are important only at the meter/second level. If 'd_obs_geo' and
'd_obs_sun' are both zero, an average value will be used for the
relativistic term for the observer, appropriate for an observer on
the surface of the Earth. 'd_obj_sun', if given, is used only for
solar system objects.
4. This function is the C version of NOVAS Fortran routine
'radvl'.
------------------------------------------------------------------------
*/
{
static short int first_call = 1;
short int i;
static double c2, toms, toms2;
double v[3], ra, dec, radvel, posmag, uk[3], v2, vo2, r, phigeo,
phisun, rel, rar, dcr, cosdec, du[3], zc, kv, zb1, kvobs, zobs1;
/*
Set up local constants the first time this function is called.
*/
if (first_call)
{
c2 = C * C;
toms = AU / 86400.0;
toms2 = toms * toms;
first_call = 0;
}
/*
Initialize variables needed for radial velocity calculation.
*/
for (i = 0; i < 3; i++)
{
v[i] = vel[i];
}
switch (cel_object->type)
{
case 2: /* Data for objects outside the solar system. */
ra = cel_object->star.ra;
dec = cel_object->star.dec;
radvel = cel_object->star.radialvelocity;
if (cel_object->star.parallax <= 0.0)
{
for (i = 0; i < 3; i++)
{
v[i] = 0.0;
}
}
break;
default: /* Data not used for solar system objects. */
ra = 0.0;
dec = 0.0;
radvel = 0.0;
}
/*
Compute length of position vector = distance to object in AU.
*/
posmag = sqrt (pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]);
/*
Compute unit vector toward object.
*/
for (i = 0; i < 3; i++)
{
uk[i] = pos[i] / posmag;
}
/*
Compute velocity-squared factors.
*/
v2 = (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) * toms2;
vo2 = (vel_obs[0] * vel_obs[0] + vel_obs[1] * vel_obs[1] +
vel_obs[2] * vel_obs[2]) * toms2;
/*
Compute geopotential at observer, unless observer is geocentric.
*/
r = d_obs_geo * AU;
if (r > 1.0e6)
{
phigeo = GE / r;
}
else
{
phigeo = 0.0;
}
/*
Compute solar potential at observer.
*/
r = d_obs_sun * AU;
if (r > 1.0e8)
{
phisun = GS / r;
}
else
{
phisun = 0.0;
}
/*
Compute relativistic potential and velocity factor for observer.
*/
if ((d_obs_geo != 0.0) || (d_obs_sun != 0.0))
{
/*
Lindegren & Dravins eq. (41), second factor in parentheses.
*/
rel = 1.0 - (phigeo + phisun) / c2 - 0.5 * vo2 / c2;
}
else
{
/*
Lindegren & Dravins eq. (42), inverse.
*/
rel = 1.0 - 1.550e-8;
}
/*
Complete radial velocity calculation.
*/
switch (cel_object->type)
{
case 2: /* Objects outside the solar system. */
/*
For stars, update barycentric radial velocity measure for change
in view angle.
*/
rar = ra * 15.0 * DEG2RAD;
dcr = dec * DEG2RAD;
cosdec = cos (dcr);
du[0] = uk[0] - (cosdec * cos (rar));
du[1] = uk[1] - (cosdec * sin (rar));
du[2] = uk[2] - sin (dcr);
zc = radvel * 1.0e3 +
(v[0] * du[0] + v[1] * du[1] + v[2] * du[2]) * toms;
/*
Compute observed radial velocity measure of a star (inverse of
Lindegren & Dravins eq. (41)).
*/
zb1 = 1.0 + zc / C;
kvobs = (uk[0] * vel_obs[0] + uk[1] * vel_obs[1] +
uk[2] * vel_obs[2]) * toms;
zobs1 = zb1 * rel / (1.0 + kvobs / C);
break;
case 0: /* Objects in the solar system */
case 1:
default:
/*
Compute solar potential at object, if within solar system.
*/
r = d_obj_sun * AU;
if ((r > 1.0e8) && (r < 1.0e16))
{
phisun = GS / r;
}
else
{
phisun = 0.0;
}
/*
Compute observed radial velocity measure of a planet or other
object -- including a nearby star -- where kinematic barycentric
velocity vector is known and gravitational red shift is negligible
(Lindegren & Dravins eq. (40), applied as per S. Klioner private
communication (2006)).
*/
kv = (uk[0] * vel[0] + uk[1] * vel[1] + uk[2] * vel[2]) * toms;
zb1 = (1.0 + kv / C ) / (1.0 - phisun / c2 - 0.5 * v2 / c2);
kvobs = (uk[0] * vel_obs[0] + uk[1] * vel_obs[1] +
uk[2] * vel_obs[2] ) * toms;
zobs1 = zb1 * rel / (1.0 + kvobs / C);
}
/*
Convert observed radial velocity measure to kilometers/second.
*/
*rv = (zobs1 - 1.0) * C / 1000.0;
return;
}
/********precession */
short int precession (double jd_tdb1, double *pos1, double jd_tdb2,
double *pos2)
/*
------------------------------------------------------------------------
PURPOSE:
Precesses equatorial rectangular coordinates from one epoch to
another. One of the two epochs must be J2000.0. The coordinates
are referred to the mean dynamical equator and equinox of the two
respective epochs.
REFERENCES:
Explanatory Supplement To The Astronomical Almanac, pp. 103-104.
Capitaine, N. et al. (2003), Astronomy And Astrophysics 412,
pp. 567-586.
Hilton, J. L. et al. (2006), IAU WG report, Celest. Mech., 94,
pp. 351-367.
INPUT
ARGUMENTS:
jd_tdb1 (double)
TDB Julian date of first epoch. See Note 1 below.
pos1[3] (double)
Position vector, geocentric equatorial rectangular coordinates,
referred to mean dynamical equator and equinox of first epoch.
jd_tdb2 (double)
TDB Julian date of second epoch. See Note 1 below.
OUTPUT
ARGUMENTS:
pos2[3] (double)
Position vector, geocentric equatorial rectangular coordinates,
referred to mean dynamical equator and equinox of second epoch.
RETURNED
VALUE:
(short int)
= 0 ... everything OK.
= 1 ... Precession not to or from J2000.0; 'jd_tdb1' or 'jd_tdb2'
not 2451545.0.
GLOBALS
USED:
T0, ASEC2RAD novascon.c
FUNCTIONS
CALLED:
fabs math.h
sin math.h
cos math.h
VER./DATE/
PROGRAMMER:
V1.0/01-93/TKB (USNO/NRL Optical Interfer.) Translate Fortran.
V1.1/08-93/WTH (USNO/AA) Update to C Standards.
V1.2/03-98/JAB (USNO/AA) Change function type from 'short int' to
'void'.
V1.3/12-99/JAB (USNO/AA) Precompute trig terms for greater
efficiency.
V2.0/10-03/JAB (USNO/AA) Update function for consistency with
IERS (2000) Conventions.
V2.1/01-05/JAB (USNO/AA) Update expressions for the precession
angles (extra significant digits).
V2.2/04-06/JAB (USNO/AA) Update model to 2006 IAU convention.
This is model "P03" of second reference.
V2.3/03-10/JAB (USNO/AA) Implement 'first-time' to fix bug when
'jd_tdb2' is 'T0' on first call to
function.
NOTES:
1. Either 'jd_tdb1' or 'jd_tdb2' must be 2451545.0 (J2000.0) TDB.
2. This function is the C version of NOVAS Fortran routine
'preces'.
------------------------------------------------------------------------
*/
{
static short int first_time = 1;
short int error = 0;
static double t_last = 0.0;
static double xx, yx, zx, xy, yy, zy, xz, yz, zz;
double eps0 = 84381.406;
double t, psia, omegaa, chia, sa, ca, sb, cb, sc, cc, sd, cd;
/*
Check to be sure that either 'jd_tdb1' or 'jd_tdb2' is equal to T0.
*/
if ((jd_tdb1 != T0) && (jd_tdb2 != T0))
return (error = 1);
/*
't' is time in TDB centuries between the two epochs.
*/
t = (jd_tdb2 - jd_tdb1) / 36525.0;
if (jd_tdb2 == T0)
t = -t;
if ((fabs (t - t_last) >= 1.0e-15) || (first_time == 1))
{
/*
Numerical coefficients of psi_a, omega_a, and chi_a, along with
epsilon_0, the obliquity at J2000.0, are 4-angle formulation from
Capitaine et al. (2003), eqs. (4), (37), & (39).
*/
psia = ((((- 0.0000000951 * t
+ 0.000132851 ) * t
- 0.00114045 ) * t
- 1.0790069 ) * t
+ 5038.481507 ) * t;
omegaa = ((((+ 0.0000003337 * t
- 0.000000467 ) * t
- 0.00772503 ) * t
+ 0.0512623 ) * t
- 0.025754 ) * t + eps0;
chia = ((((- 0.0000000560 * t
+ 0.000170663 ) * t
- 0.00121197 ) * t
- 2.3814292 ) * t
+ 10.556403 ) * t;
eps0 = eps0 * ASEC2RAD;
psia = psia * ASEC2RAD;
omegaa = omegaa * ASEC2RAD;
chia = chia * ASEC2RAD;
sa = sin (eps0);
ca = cos (eps0);
sb = sin (-psia);
cb = cos (-psia);
sc = sin (-omegaa);
cc = cos (-omegaa);
sd = sin (chia);
cd = cos (chia);
/*
Compute elements of precession rotation matrix equivalent to
R3(chi_a) R1(-omega_a) R3(-psi_a) R1(epsilon_0).
*/
xx = cd * cb - sb * sd * cc;
yx = cd * sb * ca + sd * cc * cb * ca - sa * sd * sc;
zx = cd * sb * sa + sd * cc * cb * sa + ca * sd * sc;
xy = -sd * cb - sb * cd * cc;
yy = -sd * sb * ca + cd * cc * cb * ca - sa * cd * sc;
zy = -sd * sb * sa + cd * cc * cb * sa + ca * cd * sc;
xz = sb * sc;
yz = -sc * cb * ca - sa * cc;
zz = -sc * cb * sa + cc * ca;
t_last = t;
first_time = 0;
}
if (jd_tdb2 == T0)
{
/*
Perform rotation from epoch to J2000.0.
*/
pos2[0] = xx * pos1[0] + xy * pos1[1] + xz * pos1[2];
pos2[1] = yx * pos1[0] + yy * pos1[1] + yz * pos1[2];
pos2[2] = zx * pos1[0] + zy * pos1[1] + zz * pos1[2];
}
else
{
/*
Perform rotation from J2000.0 to epoch.
*/
pos2[0] = xx * pos1[0] + yx * pos1[1] + zx * pos1[2];
pos2[1] = xy * pos1[0] + yy * pos1[1] + zy * pos1[2];
pos2[2] = xz * pos1[0] + yz * pos1[1] + zz * pos1[2];
}
return (error = 0);
}
/********nutation */
void nutation (double jd_tdb, short int direction, short int accuracy,
double *pos,
double *pos2)
/*
------------------------------------------------------------------------
PURPOSE:
Nutates equatorial rectangular coordinates from mean equator and
equinox of epoch to true equator and equinox of epoch. Inverse
transformation may be applied by setting flag 'direction'.
REFERENCES:
Explanatory Supplement To The Astronomical Almanac, pp. 114-115.
INPUT
ARGUMENTS:
jd_tdb (double)
TDB Julian date of epoch.
direction (short int)
Flag determining 'direction' of transformation;
direction = 0 transformation applied, mean to true.
direction != 0 inverse transformation applied, true to mean.
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
pos[3] (double)
Position vector, geocentric equatorial rectangular coordinates,
referred to mean equator and equinox of epoch.
OUTPUT
ARGUMENTS:
pos2[3] (double)
Position vector, geocentric equatorial rectangular coordinates,
referred to true equator and equinox of epoch.
RETURNED
VALUE:
None.
GLOBALS
USED:
DEG2RAD, ASEC2RAD novascon.c
FUNCTIONS
CALLED:
e_tilt novas.c
cos math.h
sin math.h
VER./DATE/
PROGRAMMER:
V1.0/01-93/TKB (USNO/NRL Optical Interfer.) Translate Fortran.
V1.1/08-93/WTH (USNO/AA) Update to C Standards.
V1.2/11-03/JAB (USNO/AA) Remove returned value.
V1.3/01-06/WKP (USNO/AA) Changed 'mode' to 'accuracy'.
NOTES:
1. This function is the C version of NOVAS Fortran routine
'nutate'.
------------------------------------------------------------------------
*/
{
double cobm, sobm, cobt, sobt, cpsi, spsi, xx, yx, zx, xy, yy, zy,
xz, yz, zz, oblm, oblt, eqeq, psi, eps;
/*
Call 'e_tilt' to get the obliquity and nutation angles.
*/
e_tilt (jd_tdb,accuracy, &oblm,&oblt,&eqeq,&psi,&eps);
cobm = cos (oblm * DEG2RAD);
sobm = sin (oblm * DEG2RAD);
cobt = cos (oblt * DEG2RAD);
sobt = sin (oblt * DEG2RAD);
cpsi = cos (psi * ASEC2RAD);
spsi = sin (psi * ASEC2RAD);
/*
Nutation rotation matrix follows.
*/
xx = cpsi;
yx = -spsi * cobm;
zx = -spsi * sobm;
xy = spsi * cobt;
yy = cpsi * cobm * cobt + sobm * sobt;
zy = cpsi * sobm * cobt - cobm * sobt;
xz = spsi * sobt;
yz = cpsi * cobm * sobt - sobm * cobt;
zz = cpsi * sobm * sobt + cobm * cobt;
if (!direction)
{
/*
Perform rotation.
*/
pos2[0] = xx * pos[0] + yx * pos[1] + zx * pos[2];
pos2[1] = xy * pos[0] + yy * pos[1] + zy * pos[2];
pos2[2] = xz * pos[0] + yz * pos[1] + zz * pos[2];
}
else
{
/*
Perform inverse rotation.
*/
pos2[0] = xx * pos[0] + xy * pos[1] + xz * pos[2];
pos2[1] = yx * pos[0] + yy * pos[1] + yz * pos[2];
pos2[2] = zx * pos[0] + zy * pos[1] + zz * pos[2];
}
return;
}
/********nutation_angles */
void nutation_angles (double t, short int accuracy,
double *dpsi, double *deps)
/*
------------------------------------------------------------------------
PURPOSE:
This function returns the values for nutation in longitude and
nutation in obliquity for a given TDB Julian date. The nutation
model selected depends upon the input value of 'accuracy'. See
notes below for important details.
REFERENCES:
Kaplan, G. (2005), US Naval Observatory Circular 179.
INPUT
ARGUMENTS:
t (double)
TDB time in Julian centuries since J2000.0
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
dpsi (double)
Nutation in longitude in arcseconds.
deps (double)
Nutation in obliquity in arcseconds.
RETURNED
VALUE:
None.
GLOBALS
USED:
T0, ASEC2RAD novascon.c
FUNCTIONS
CALLED:
iau2000a nutation.c
iau2000b nutation.c
nu2000k nutation.c
VER./DATE/
PROGRAMMER:
V1.0/12-04/JAB (USNO/AA)
V1.1/01-06/WKP (USNO/AA): Changed 'mode' to 'accuracy'.
V1.2/02-06/WKP (USNO/AA): Fixed units bug.
V1.3/01-07/JAB (USNO/AA): Implemented 'low_acc_choice' construct.
NOTES:
1. This function selects the nutation model depending first upon
the input value of 'accuracy'. If 'accuracy' = 0 (full accuracy),
the IAU 2000A nutation model is used. If 'accuracy' = 1 (reduced
accuracy, the model used depends upon the value of local
variable 'low_acc_choice', which is set below.
2. If local variable 'low_acc_choice' = 1 (the default), a
specially truncated version of IAU 2000A, called 'NU2000K' is
used. If 'low_acc_choice' = 2, the IAU 2000B nutation model is
used.
3. See the prologs of the nutation functions in file 'nutation.c'
for details concerning the models.
4. This function is the C version of NOVAS Fortran routine
'nod'.
------------------------------------------------------------------------
*/
{
/*
Set the value of 'low_acc_choice' according to the rules explained
under NOTES in the prolog.
*/
short int low_acc_choice = 1;
double t1;
t1 = t * 36525.0;
/*
High accuracy mode -- use IAU 2000A.
*/
if (accuracy == 0)
{
iau2000a (T0,t1, dpsi,deps);
}
/*
Low accuracy mode -- model depends upon value of 'low_acc_choice'.
*/
else
{
switch (low_acc_choice)
{
case 2: /* use IAU 2000B */
iau2000b (T0,t1, dpsi,deps);
break;
case 1: /* use NU2000K */
default:
nu2000k (T0,t1, dpsi,deps);
}
}
/*
Convert output to arcseconds.
*/
*dpsi /= ASEC2RAD;
*deps /= ASEC2RAD;
return;
}
/********fund_args */
void fund_args (double t,
double a[5])
/*
------------------------------------------------------------------------
PURPOSE:
To compute the fundamental arguments (mean elements) of the Sun
and Moon.
REFERENCES:
Simon et al. (1994) Astronomy and Astrophysics 282, 663-683,
esp. Sections 3.4-3.5.
INPUT
ARGUMENTS:
t (double)
TDB time in Julian centuries since J2000.0
OUTPUT
ARGUMENTS:
a[5] (double)
Fundamental arguments, in radians:
a[0] = l (mean anomaly of the Moon)
a[1] = l' (mean anomaly of the Sun)
a[2] = F (mean argument of the latitude of the Moon)
a[3] = D (mean elongation of the Moon from the Sun)
a[4] = Omega (mean longitude of the Moon's ascending node);
from Simon section 3.4(b.3),
precession = 5028.8200 arcsec/cy)
RETURNED
VALUE:
None.
GLOBALS
USED:
ASEC2RAD, ASEC360 novascon.c
FUNCTIONS
CALLED:
fmod math.h
VER./DATE/
PROGRAMMER:
V1.0/10-97/JAB (USNO/AA)
V1.1/07-98/JAB (USNO/AA): Place arguments in the range 0-TWOPI
radians.
V1.2/09-03/JAB (USNO/AA): Incorporate function 'norm_ang'.
V1.3/11-03/JAB (USNO/AA): Update with Simon et al. expressions.
V1.4/01-06/JAB (USNO/AA): Remove function 'norm_ang'; rewrite for
consistency with Fortran.
V1.5/02-11/WKP (USNO/AA): Clarified a[4] description in prolog.
NOTES:
1. This function is the C version of NOVAS Fortran routine
'funarg'.
------------------------------------------------------------------------
*/
{
a[0] = fmod (485868.249036 +
t * (1717915923.2178 +
t * ( 31.8792 +
t * ( 0.051635 +
t * ( - 0.00024470)))), ASEC360) * ASEC2RAD;
a[1] = fmod (1287104.79305 +
t * ( 129596581.0481 +
t * ( - 0.5532 +
t * ( 0.000136 +
t * ( - 0.00001149)))), ASEC360) * ASEC2RAD;
a[2] = fmod (335779.526232 +
t * (1739527262.8478 +
t * ( - 12.7512 +
t * ( - 0.001037 +
t * ( 0.00000417)))), ASEC360) * ASEC2RAD;
a[3] = fmod (1072260.70369 +
t * (1602961601.2090 +
t * ( - 6.3706 +
t * ( 0.006593 +
t * ( - 0.00003169)))), ASEC360) * ASEC2RAD;
a[4] = fmod (450160.398036 +
t * ( - 6962890.5431 +
t * ( 7.4722 +
t * ( 0.007702 +
t * ( - 0.00005939)))), ASEC360) * ASEC2RAD;
return;
}
/********mean_obliq */
double mean_obliq (double jd_tdb)
/*
------------------------------------------------------------------------
PURPOSE:
To compute the mean obliquity of the ecliptic.
REFERENCES:
Capitaine et al. (2003), Astronomy and Astrophysics 412, 567-586.
INPUT
ARGUMENTS:
jd_tdb (double)
TDB Julian Date.
OUTPUT
ARGUMENTS:
None.
RETURNED
VALUE:
(double)
Mean obliquity of the ecliptic in arcseconds.
GLOBALS
USED:
T0 novascon.c
FUNCTIONS
CALLED:
None.
VER./DATE/
PROGRAMMER:
V1.0/12-04/JAB (USNO/AA)
V2.0/04-06/JAB (USNO/AA) Update the expression for mean obliquity
using data from the reference.
NOTES:
None.
------------------------------------------------------------------------
*/
{
double t, epsilon;
/*
Compute time in Julian centuries from epoch J2000.0.
*/
t = (jd_tdb - T0) / 36525.0;
/*
Compute the mean obliquity in arcseconds. Use expression from the
reference's eq. (39) with obliquity at J2000.0 taken from eq. (37)
or Table 8.
*/
epsilon = (((( - 0.0000000434 * t
- 0.000000576 ) * t
+ 0.00200340 ) * t
- 0.0001831 ) * t
- 46.836769 ) * t + 84381.406;
return (epsilon);
}
/********vector2radec */
short int vector2radec (double *pos,
double *ra, double *dec)
/*
------------------------------------------------------------------------
PURPOSE:
Converts an vector in equatorial rectangular coordinates to
equatorial spherical coordinates.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
INPUT
ARGUMENTS:
pos[3] (double)
Position vector, equatorial rectangular coordinates.
OUTPUT
ARGUMENTS:
*rightascension (double)
Right ascension in hours.
*declination (double)
Declination in degrees.
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
= 1 ... All vector components are zero; 'ra' and 'dec' are
indeterminate.
= 2 ... Both pos[0] and pos[1] are zero, but pos[2] is nonzero;
'ra' is indeterminate.
GLOBALS
USED:
ASEC2RAD novascon.c
FUNCTIONS
CALLED:
sqrt math.h
atan2 math.h
VER./DATE/
PROGRAMMER:
V1.0/01-93/TKB (USNO/NRL Optical Interfer.) Translate Fortran.
V1.1/08-93/WTH (USNO/AA) Update to C Standards.
V1.2/01-04/JAB (USNO/AA) Remove calls to 'pow'.
NOTES:
1. This function is the C version of NOVAS Fortran routine
'angles'.
------------------------------------------------------------------------
*/
{
double xyproj;
xyproj = sqrt (pos[0] * pos[0] + pos[1] * pos[1]);
if ((xyproj == 0.0) && (pos[2] == 0))
{
*ra = 0.0;
*dec = 0.0;
return 1;
}
else if (xyproj == 0.0)
{
*ra = 0.0;
if (pos[2] < 0.0)
*dec = -90.0;
else
*dec = 90.0;
return 2;
}
else
{
*ra = atan2 (pos[1], pos[0]) / ASEC2RAD / 54000.0;
*dec = atan2 (pos[2], xyproj) / ASEC2RAD / 3600.0;
if (*ra < 0.0)
*ra += 24.0;
}
return 0;
}
/********radec2vector */
void radec2vector (double ra, double dec, double dist,
double *vector)
/*
------------------------------------------------------------------------
PURPOSE:
Converts equatorial spherical coordinates to a vector (equatorial
rectangular coordinates).
REFERENCES:
None.
INPUT
ARGUMENTS:
ra (double)
Right ascension (hours).
dec (double)
Declination (degrees).
dist (double)
Distance (AU)
OUTPUT
ARGUMENTS:
vector[3] (double)
Position vector, equatorial rectangular coordinates (AU).
RETURNED
VALUE:
(short int)
= 0 ... Everything OK.
GLOBALS
USED:
DEG2RAD novascon.c
FUNCTIONS
CALLED:
cos math.h
sin math.h
VER./DATE/
PROGRAMMER:
V1.0/05-92/TKB (USNO/NRL Optical Interfer.) Translate Fortran.
V1.1/08-93/WTH (USNO/AA) Update to C Standards.
V1.2/08-09/JLB (USNO/AA) Documented "dist" in prolog
NOTES:
None.
------------------------------------------------------------------------
*/
{
vector[0] = dist * cos (DEG2RAD * dec) * cos (DEG2RAD * 15.0 * ra);
vector[1] = dist * cos (DEG2RAD * dec) * sin (DEG2RAD * 15.0 * ra);
vector[2] = dist * sin (DEG2RAD * dec);
return;
}
/********starvectors */
void starvectors (cat_entry *star,
double *pos, double *vel)
/*
------------------------------------------------------------------------
PURPOSE:
Converts angular quantities for stars to vectors.
REFERENCES:
Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.
INPUT
ARGUMENTS:
*star (struct cat_entry)
Pointer to catalog entry structure containing ICRS catalog
data (defined in novas.h).
OUTPUT
ARGUMENTS:
pos[3] (double)
Position vector, equatorial rectangular coordinates,
components in AU.
vel[3] (double)
Velocity vector, equatorial rectangular coordinates,
components in AU/Day.
RETURNED
VALUE:
None.
GLOBALS
USED:
ASEC2RAD, DEG2RAD novascon.c
AU_KM, C novascon.c
FUNCTIONS
CALLED:
sin math.h
cos math.h
VER./DATE/
PROGRAMMER:
V1.0/01-93/TKB (USNO/NRL Optical Interfer.) Translate Fortran.
V1.1/08-93/WTH (USNO/AA) Updated to C programming standards.
V1.2/11-03/JAB (USNO/AA) Update to account for parallax in units
of milliarcseconds and proper motion
in milliarcseconds/year.
V1.3/01-06/WKP (USNO/AA) Changed 'a[3]' to 'd' and added Doppler
factor.
V1.4/02-06/WKP (USNO/AA) Corrected units of C to km/sec.
NOTES:
1. This function is the C version of NOVAS Fortran routine
'vectrs'.
------------------------------------------------------------------------
*/
{
double paralx, dist, r, d, cra, sra, cdc, sdc, pmr, pmd, rvl, k;
/*
If parallax is unknown, undetermined, or zero, set it to 1e-6
milliarcsecond, corresponding to a distance of 1 gigaparsec.
*/
paralx = star->parallax;
if (star->parallax <= 0.0)
paralx = 1.0e-6;
/*
Convert right ascension, declination, and parallax to position vector
in equatorial system with units of AU.
*/
dist = 1.0 / sin (paralx * 1.0e-3 * ASEC2RAD);
r = (star->ra) * 15.0 * DEG2RAD;
d = (star->dec) * DEG2RAD;
cra = cos (r);
sra = sin (r);
cdc = cos (d);
sdc = sin (d);
pos[0] = dist * cdc * cra;
pos[1] = dist * cdc * sra;
pos[2] = dist * sdc;
/*
Compute Doppler factor, which accounts for change in
light travel time to star.
*/
k = 1.0 / (1.0 - star->radialvelocity / C * 1000.0);
/*
Convert proper motion and radial velocity to orthogonal components of
motion with units of AU/Day.
*/
pmr = star->promora / (paralx * 365.25) * k;
pmd = star->promodec / (paralx * 365.25) * k;
rvl = star->radialvelocity * 86400.0 / AU_KM * k;
/*
Transform motion vector to equatorial system.
*/
vel[0] = - pmr * sra - pmd * sdc * cra + rvl * cdc * cra;
vel[1] = pmr * cra - pmd * sdc * sra + rvl * cdc * sra;
vel[2] = pmd * cdc + rvl * sdc;
return;
}
/********tdb2tt */
void tdb2tt (double tdb_jd,
double *tt_jd, double *secdiff)
/*
------------------------------------------------------------------------
PURPOSE:
Computes the Terrestrial Time (TT) or Terrestrial Dynamical Time
(TDT) Julian date corresponding to a Barycentric Dynamical Time
(TDB) Julian date.
REFERENCES:
Fairhead, L. & Bretagnon, P. (1990) Astron. & Astrophys. 229, 240.
Kaplan, G. (2005), US Naval Observatory Circular 179.
INPUT
ARGUMENTS:
tdb_jd (double)
TDB Julian date.
OUTPUT
ARGUMENTS:
*tt_jd (double)
TT Julian date.
*secdiff (double)
Difference 'tdb_jd'-'tt_jd', in seconds.
RETURNED
VALUE:
None.
GLOBALS
USED:
T0 novascon.c
FUNCTIONS
CALLED:
sin math.h
VER./DATE/
PROGRAMMER:
V1.0/07-92/TKB (USNO/NRL Optical Interfer.) Translate Fortran.
V1.1/08-93/WTH (USNO/AA) Update to C Standards.
V1.2/06-98/JAB (USNO/AA) Adopt new model (Explanatory Supplement
to the Astronomical Almanac, pp. 42-44
and p. 316.)
V1.3/11-03/JAB (USNO/AA) Changed variable names of the input
Julian dates to make more descriptive.
V1.4/01-07/JAB (USNO/AA) Adopt Fairhead & Bretagnon expression.
NOTES:
1. Expression used in this function is a truncated form of a
longer and more precise series given in the first reference. The
result is good to about 10 microseconds.
2. This function is the C version of NOVAS Fortran routine
'times'.
------------------------------------------------------------------------
*/
{
double t;
t = (tdb_jd - T0) / 36525.0;
/*
Expression given in USNO Circular 179, eq. 2.6.
*/
*secdiff = 0.001657 * sin ( 628.3076 * t + 6.2401)
+ 0.000022 * sin ( 575.3385 * t + 4.2970)
+ 0.000014 * sin (1256.6152 * t + 6.1969)
+ 0.000005 * sin ( 606.9777 * t + 4.0212)
+ 0.000005 * sin ( 52.9691 * t + 0.4444)
+ 0.000002 * sin ( 21.3299 * t + 5.5431)
+ 0.000010 * t * sin ( 628.3076 * t + 4.2490);
*tt_jd = tdb_jd - *secdiff / 86400.0;
return;
}
/********cio_ra */
short int cio_ra (double jd_tt, short int accuracy,
double *ra_cio)
/*
------------------------------------------------------------------------
PURPOSE:
This function computes the true right ascension of the celestial
intermediate origin (CIO) at a given TT Julian date. This is
-(equation of the origins).
REFERENCES:
Kaplan, G. (2005), US Naval Observatory Circular 179.
INPUT
ARGUMENTS:
jd_tt (double)
TT Julian date.
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*ra_cio (double)
Right ascension of the CIO, with respect to the true equinox
of date, in hours (+ or -).
RETURNED
VALUE:
(short int)
= 0 ... everything OK.
= 1 ... invalid value of 'accuracy'.
> 10 ... 10 + the error code from function 'cio_location'.
> 20 ... 20 + the error code from function 'cio_basis'.
GLOBALS
USED:
RAD2DEG, T0 novascon.c
FUNCTIONS
CALLED:
tdb2tt novas.c
cio_location novas.c
cio_basis novas.c
nutation novas.c
precession novas.c
frame_tie novas.c
atan2 math.h
VER./DATE/
PROGRAMMER:
V1.0/01-05/JAB (USNO/AA)
V1.1/12-05/WKP (USNO/AA) Corrected minor syntax problems.
V1.2/01-06/WKP (USNO/AA) Changed 'mode' to 'accuracy'.
V1.3/07-06/JAB (USNO/AA) Implement 'cio_location' construct.
NOTES:
1. This function is the C version of NOVAS Fortran routine
'ciora'.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int rs;
double unitx[3] = {1.0, 0.0, 0.0};
double jd_tdb, t, secdif, x[3], y[3], z[3], w1[3], w2[3],
eq[3], az, r_cio;
/*
Check for valid value of 'accuracy'.
*/
if ((accuracy < 0) || (accuracy > 1))
return (error = 1);
/*
'jd_tdb' is the TDB Julian date.
*/
jd_tdb = jd_tt;
tdb2tt (jd_tdb, &t,&secdif);
jd_tdb = jd_tt + secdif / 86400.0;
/*
Obtain the basis vectors, in the GCRS, for the celestial intermediate
system defined by the CIP (in the z direction) and the CIO (in the
x direction).
*/
if ((error = cio_location (jd_tdb,accuracy, &r_cio,&rs)) != 0)
{
*ra_cio = 0.0;
return (error += 10);
}
if ((error = cio_basis (jd_tdb,r_cio,rs,accuracy, x,y,z)) != 0)
{
return (error += 20);
}
/*
Compute the direction of the true equinox in the GCRS.
*/
nutation (jd_tdb,-1,accuracy,unitx, w1);
precession (jd_tdb,w1,T0, w2);
frame_tie (w2,-1, eq);
/*
Compute the RA-like coordinate of the true equinox in the celestial
intermediate system.
*/
az = atan2 (eq[0] * y[0] + eq[1] * y[1] + eq[2] * y[2],
eq[0] * x[0] + eq[1] * x[1] + eq[2] * x[2]) * RAD2DEG;
/*
The RA of the CIO is minus this coordinate.
*/
*ra_cio = -az / 15.0;
return (error);
}
/********cio_location */
short int cio_location (double jd_tdb, short int accuracy,
double *ra_cio, short int *ref_sys)
/*
------------------------------------------------------------------------
PURPOSE:
This function returns the location of the celestial
intermediate origin (CIO) for a given Julian date, as a
right ascension with respect to either the GCRS (geocentric ICRS)
origin or the true equinox of date. The CIO is always located on
the true equator (= intermediate equator) of date.
REFERENCES:
None.
INPUT
ARGUMENTS:
jd_tdb (double)
TDB Julian date.
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*ra_cio (double)
Right ascension of the CIO, in hours.
*ref_sys (short int)
Reference system in which right ascension is given
= 1 ... GCRS
= 2 ... True equator and equinox of date.
RETURNED
VALUE:
(short int)
= 0 ... everything OK.
= 1 ... unable to allocate memory for the 'cio' array.
> 10 ... 10 + the error code from function 'cio_array'.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
cio_array novas.c
ira_equinox novas.c
fopen stdio.h
fclose stdio.h
fabs math.h
calloc stdlib.h
VER./DATE/
PROGRAMMER:
V1.0/07-06/JAB (USNO/AA)
NOTES:
1. If an external file of CIO right ascensions is available,
it will be used and 'ref_sys' will be set to 1. Otherwise an
internal computation will be used and 'ref_sys' will be set to 2.
2. The external binary file of CIO right ascensions is assumed
to be named 'cio_ra.bin'. Utility program 'cio_file.c', provided
with the NOVAS-C package, creates this file from a text file also
provided with NOVAS-C.
3. This function is the C version of NOVAS Fortran routine
'cioloc'.
------------------------------------------------------------------------
*/
{
static short int first_call = 1;
static short int ref_sys_last = 0;
static short int use_file = 0;
short int error = 0;
long int n_pts = 6;
long int i, j;
static double t_last = 0.0;
static double ra_last;
double p, eq_origins;
size_t cio_size;
static ra_of_cio *cio;
static FILE *cio_file;
/*
Check if the input external binary file exists and can be read.
*/
if (first_call)
{
if ((cio_file = fopen ("cio_ra.bin", "rb")) == NULL)
{
use_file = 0;
}
else
{
use_file = 1;
fclose (cio_file);
}
}
/*
Check if previously computed RA value can be used.
*/
if ((fabs (jd_tdb - t_last) <= 1.0e-8))
{
*ra_cio = ra_last;
*ref_sys = ref_sys_last;
return (error = 0);
}
/*
Compute the RA of the CIO.
*/
switch (use_file)
{
/*
-----------------------------
Interpolate values from file.
-----------------------------
*/
case 1:
/*
Allocate memory for the array 'cio'. This array contains the values
to be interpolated, extracted from the CIO file.
*/
if (first_call)
{
cio_size = sizeof (ra_of_cio);
cio = (ra_of_cio *) calloc ((size_t) n_pts, cio_size);
if (cio == NULL)
return (error = 1);
else
first_call = 0;
}
/*
Get array of values to interpolate.
*/
if ((error = cio_array (jd_tdb,n_pts, cio)) != 0)
{
*ra_cio = 0.0;
return (error += 10);
}
/*
Perform Lagrangian interpolation for the RA at 'tdb_jd'.
*/
*ra_cio = 0.0;
for (j = 0L; j < n_pts; j++)
{
p = 1.0;
for (i = 0L; i < n_pts; i++)
{
if (i != j)
p *= ((jd_tdb - cio[i].jd_tdb) /
(cio[j].jd_tdb - cio[i].jd_tdb));
}
*ra_cio += (p * cio[j].ra_cio);
}
*ra_cio /= 54000.0;
*ref_sys = 1;
break;
/*
-------------------------
Use internal computation.
-------------------------
*/
case 0:
/*
Compute equation of the origins.
*/
if (first_call)
first_call = 0;
eq_origins = ira_equinox (jd_tdb,1,accuracy);
*ra_cio = -eq_origins;
*ref_sys = 2;
break;
}
t_last = jd_tdb;
ra_last = *ra_cio;
ref_sys_last = *ref_sys;
return (error);
}
/********cio_basis */
short int cio_basis (double jd_tdb, double ra_cio, short int ref_sys,
short int accuracy,
double *x, double *y, double *z)
/*
------------------------------------------------------------------------
PURPOSE:
To compute the orthonormal basis vectors, with
respect to the GCRS (geocentric ICRS), of the celestial
intermediate system defined by the celestial intermediate pole
(CIP) (in the z direction) and the celestial intermediate origin
(CIO) (in the x direction). A TDB Julian date and the right
ascension of the CIO at that date is required as input. The
right ascension of the CIO can be with respect to either the
GCRS origin or the true equinox of date -- different algorithms
are used in the two cases.
REFERENCES:
Kaplan, G. (2005), US Naval Observatory Circular 179.
INPUT
ARGUMENTS:
jd_tdb (double)
TDB Julian date of epoch.
ra_cio (double)
Right ascension of the CIO at epoch (hours).
ref_sys (short int)
Reference system in which right ascension is given (output
from function 'cio_location')
= 1 ... GCRS
= 2 ... True equator and equinox of date.
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
*x (double)
Unit vector toward the CIO, equatorial rectangular
coordinates, referred to the GCRS.
*y (double)
Unit vector toward the y-direction, equatorial rectangular
coordinates, referred to the GCRS.
*z (double)
Unit vector toward north celestial pole (CIP), equatorial
rectangular coordinates, referred to the GCRS.
RETURNED
VALUE:
(short int)
= 0 ... everything OK.
= 1 ... invalid value of input variable 'ref_sys'.
GLOBALS
USED:
T0, DEG2RAD novascon.c
FUNCTIONS
CALLED:
nutation novas.c
precession novas.c
frame_tie novas.c
fabs math.h
sin math.h
cos math.h
sqrt math.h
VER./DATE/
PROGRAMMER:
V1.0/08-04/JAB (USNO/AA)
V1.1/01-06/WKP (USNO/AA) Changed 'mode' to 'accuracy'.
V1.2/07-06/JAB (USNA/AA) Incorporate code to use 'ref_sys' input.
V1.3/06-08/WKP (USNO/AA) Changed value of direction argument in
calls to 'nutation' from 1 to -1 for
consistency.
NOTES:
1. This function effectively constructs the matrix C in eq. (3)
of the reference.
2. This function is the C version of NOVAS Fortran routine
'ciobas'.
------------------------------------------------------------------------
*/
{
static short int ref_sys_last = 0;
short int error = 0;
short int i;
static double t_last = 0.0;
static double xx[3], yy[3], zz[3];
double z0[3] = {0.0, 0.0, 1.0};
double w0[3], w1[3], w2[3], sinra, cosra, xmag;
/*
Compute unit vector z toward celestial pole.
*/
if (((fabs (jd_tdb - t_last) > 1.0e-8)) || (ref_sys != ref_sys_last))
{
nutation (jd_tdb,-1,accuracy,z0, w1);
precession (jd_tdb,w1,T0, w2);
frame_tie (w2,-1, zz);
t_last = jd_tdb;
ref_sys_last = ref_sys;
}
else
{
for (i = 0; i < 3; i++)
{
x[i] = xx[i];
y[i] = yy[i];
z[i] = zz[i];
}
return (error);
}
/*
Now compute unit vectors x and y. Method used depends on the
reference system in which right ascension of the CIO is given.
*/
switch (ref_sys)
{
/*
----------------------------
RA of CIO expressed in GCRS.
----------------------------
*/
case 1:
/*
Compute vector x toward CIO in GCRS.
*/
sinra = sin (ra_cio * 15.0 * DEG2RAD);
cosra = cos (ra_cio * 15.0 * DEG2RAD);
xx[0] = zz[2] * cosra;
xx[1] = zz[2] * sinra;
xx[2] = -zz[0] * cosra - zz[1] * sinra;
/*
Normalize vector x.
*/
xmag = sqrt (xx[0] * xx[0] + xx[1] * xx[1] + xx[2] * xx[2]);
xx[0] /= xmag;
xx[1] /= xmag;
xx[2] /= xmag;
/*
Compute unit vector y orthogonal to x and z (y = z cross x).
*/
yy[0] = zz[1] * xx[2] - zz[2] * xx[1];
yy[1] = zz[2] * xx[0] - zz[0] * xx[2];
yy[2] = zz[0] * xx[1] - zz[1] * xx[0];
break;
/*
----------------------------------------------------------
RA of CIO expressed in equator-and-equinox of date system.
----------------------------------------------------------
*/
case 2:
/*
Construct unit vector toward CIO in equator-and-equinox-of-date
system.
*/
w0[0] = cos (ra_cio * 15.0 * DEG2RAD);
w0[1] = sin (ra_cio * 15.0 * DEG2RAD);
w0[2] = 0.0;
/*
Rotate the vector into the GCRS to form unit vector x.
*/
nutation (jd_tdb,-1,accuracy,w0, w1);
precession (jd_tdb,w1,T0, w2);
frame_tie (w2,-1, xx);
/*
Compute unit vector y orthogonal to x and z (y = z cross x).
*/
yy[0] = zz[1] * xx[2] - zz[2] * xx[1];
yy[1] = zz[2] * xx[0] - zz[0] * xx[2];
yy[2] = zz[0] * xx[1] - zz[1] * xx[0];
break;
/*
---------------------------
Invalid value of 'ref_sys'.
---------------------------
*/
default:
for (i = 0; i < 3; i++)
{
xx[i] = 0.0;
yy[i] = 0.0;
zz[i] = 0.0;
}
error = 1;
break;
}
/*
Load the x, y, and z arrays.
*/
for (i = 0; i < 3; i++)
{
x[i] = xx[i];
y[i] = yy[i];
z[i] = zz[i];
}
return (error);
}
/********cio_array */
short int cio_array (double jd_tdb, long int n_pts,
ra_of_cio *cio)
/*
------------------------------------------------------------------------
PURPOSE:
Given an input TDB Julian date and the number of data points
desired, this function returns a set of Julian dates and
corresponding values of the GCRS right ascension of the celestial
intermediate origin (CIO). The range of dates is centered (at
least approximately) on the requested date. The function obtains
the data from an external data file.
REFERENCES:
None.
INPUT
ARGUMENTS:
jd_tdb (double)
TDB Julian date.
n_pts (long int)
Number of Julian dates and right ascension values requested
(not less than 2 or more than 20).
OUTPUT
ARGUMENTS:
*cio (struct ra_of_cio)
An time series (array) of the right ascension of the Celestial
Intermediate Origin (CIO) with respect to the GCRS (structure
defined in novas.h).
RETURNED
VALUE:
(short int)
= 0 ... everything OK
= 1 ... error opening or reading the 'cio_ra.bin' file.
= 2 ... 'jd_tdb' not in the range of the CIO file.
= 3 ... 'n_pts' out of range.
= 4 ... unable to allocate memory for the internal 't' array.
= 5 ... unable to allocate memory for the internal 'ra' array.
= 6 ... 'jd_tdb' is too close to either end of the CIO file;
unable to put 'n_pts' data points into the output
structure.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
fopen stdio.h
fread stdio.h
fclose stdio.h
abs math.h
free stdlib.h
calloc stdlib.h
fseek stdio.h
VER./DATE/
PROGRAMMER:
V1.0/09-04/JAB (USNO/AA)
V1.1/12-05/WKP (USNO/AA) Changed struct type from 'cio_ra' to
'ra_of_cio' to avoid conflicts.
V1.2/02-08/JAB (USNO/AA) Fix file-read strategy "Case 2" and
improve documentation.
NOTES:
1. This function assumes that binary, random-access file
'cio_ra.bin' has been created and is in the same directory as
your executable. This file is created by program 'cio_file.c',
included in the NOVAS-C package. On the first call to this
function, file 'cio_ra.bin' is opened in read mode.
------------------------------------------------------------------------
*/
{
static short int first_call = 1;
short int error = 0;
static long int last_index_rec = -50L;
static long int last_n_pts = 0L;
static long int header_size, record_size, n_recs;
long int min_pts = 2;
long int max_pts = 20;
long int del_n_pts, index_rec, half_int, lo_limit, hi_limit,
del_index, abs_del_index, bytes_to_lo, n_swap, n_read, i, j;
static double jd_beg, jd_end, t_int, *t, *ra;
double t_temp, ra_temp;
static size_t double_size, long_size;
static FILE *cio_file;
/*
Set the sizes of the file header and data records, open the CIO file,
and read the file header on the first call to this function.
*/
if (first_call)
{
double_size = sizeof (double);
long_size = sizeof (long int);
header_size = (long) ((size_t) 3 * double_size + long_size);
record_size = (long) ((size_t) 2 * double_size);
/*
Open the input (binary, random-access) file.
*/
if ((cio_file = fopen ("cio_ra.bin", "rb")) == NULL)
{
error = 1;
goto fail;
}
/*
Read the file header.
*/
if (1 != fread (&jd_beg, double_size, (size_t) 1, cio_file) ||
1 != fread (&jd_end, double_size, (size_t) 1, cio_file) ||
1 != fread (&t_int, double_size, (size_t) 1, cio_file) ||
1 != fread (&n_recs, long_size, (size_t) 1, cio_file))
{
error = 1;
goto fail;
}
}
/*
Check the input data against limits.
*/
if ((jd_tdb < jd_beg) || (jd_tdb > jd_end))
return (error = 2);
if ((n_pts < min_pts) || (n_pts > max_pts))
return (error = 3);
/*
Calculate the difference between the current value of 'n_pts' and
the last value of 'n_pts'.
*/
/* [Don Cross - 2021-12-29] Removed 'abs', because it was not necessary. */
del_n_pts = n_pts - last_n_pts;
/*
Allocate memory for the 't' and 'ra' arrays.
*/
if (del_n_pts != 0L)
{
if (!first_call)
{
free (t);
t = NULL;
free (ra);
ra = NULL;
}
t = (double *) calloc ((size_t) n_pts, double_size);
if (t == NULL)
{
error = 4;
goto fail;
}
ra = (double *) calloc ((size_t) n_pts, double_size);
if (ra == NULL)
{
error = 5;
goto fail;
}
first_call = 0;
}
/*
Calculate the record number of the record immediately preceding
the date of interest: the "index record".
*/
index_rec = (long int) ((jd_tdb - jd_beg) / t_int) + 1L;
/*
Test the range of 'n_pts' values centered on 'index_rec' to be sure
the range of values requested falls within the file limits.
*/
half_int = (n_pts / 2L) - 1L;
lo_limit = index_rec - half_int;
hi_limit = index_rec + (n_pts - half_int - 1L);
if ((lo_limit < 1L) || (hi_limit > n_recs))
return (error = 6);
/*
Compute the number of bytes from the beginning of the file to
the 'lo_limit'.
*/
bytes_to_lo = header_size + (lo_limit - 1L) * record_size;
/*
Compare the current index record with the previous index record.
*/
del_index = index_rec - last_index_rec;
/* [Don Cross - 2021-12-29] Reworked abs() with inline logic to eliminate warnings about long int. */
abs_del_index = (del_index < 0) ? -del_index : del_index;
/*
Determine the file read strategy.
*/
/*
Case 1: The input value of 'n_pts' changed since the last entry,
or there are no data in the current arrays that can be re-used in the
new arrays. In this case, read all new data points into the arrays.
*/
if ((abs_del_index > n_pts) || (del_n_pts != 0))
{
if (fseek (cio_file, bytes_to_lo, SEEK_SET))
{
error = 1;
goto fail;
}
for (i = 0L; i < n_pts; i++)
{
if (1 != fread (&t[i], double_size, (size_t) 1, cio_file) ||
1 != fread (&ra[i], double_size, (size_t) 1, cio_file))
{
error = 1;
goto fail;
}
}
}
/*
Case 2: The new index is close enough to the previous index that
there are some data in the current arrays that can be re-used
in the new arrays. The remaining data will be read from the CIO
file.
Note that if the new index is the same as the previous index (i.e.,
'del_index' == 0), neither Case 2a nor 2b is satisfied, and the
program flow goes directly to load the output arrays with the same
values as in the current arrays.
*/
else if ((abs_del_index <= n_pts) && (del_n_pts == 0))
{
/* [Don Cross - 2021-12-30] Reworked abs() with inline logic to eliminate warnings about long int. */
n_swap = n_pts - abs_del_index;
if (n_swap < 0)
n_swap = -n_swap;
n_read = abs_del_index;
/*
Case 2a: The new index is less than the previous one. Put the "old"
data at the end of the new arrays, and read "new" data into the
beginning of the new arrays.
*/
if (del_index < 0L)
{
for (i = 0L; i < n_swap; i++)
{
t_temp = t[i];
ra_temp = ra[i];
j = i + abs_del_index;
t[j] = t_temp;
ra[j] = ra_temp;
}
fseek (cio_file, bytes_to_lo, SEEK_SET);
for (i = 0L; i < n_read; i++)
{
if (1 != fread (&t[i], double_size, (size_t) 1, cio_file) ||
1 != fread (&ra[i], double_size, (size_t) 1, cio_file))
{
error = 1;
goto fail;
}
}
}
/*
Case 2b: The new index is greater than the previous one. Put the
"old" data at the beginning of the new arrays, and read "new" data
into the end of the new arrays.
*/
else if (del_index > 0L)
{
for (i = 0L; i < n_swap; i++)
{
j = i + abs_del_index;
t_temp = t[j];
ra_temp = ra[j];
t[i] = t_temp;
ra[i] = ra_temp;
}
if (fseek (cio_file, bytes_to_lo + (n_swap * record_size), SEEK_SET))
{
error = 1;
goto fail;
}
j = i++;
for (i = j; i < n_pts; i++)
{
if (1 != fread (&t[i], double_size, (size_t) 1, cio_file) ||
1 != fread (&ra[i], double_size, (size_t) 1, cio_file))
{
error = 1;
goto fail;
}
}
}
}
/*
Load the output 'cio' array with the values in the 't' and 'ra'
arrays.
Note that if the input value of 'n_pts' has not changed since the
last entry, all data in the current arrays can be re-used in
the new arrays. The if statements above are bypassed and the new
arrays are the same as the current arrays.
*/
for (i = 0L; i < n_pts; i++)
{
cio[i].jd_tdb = t[i];
cio[i].ra_cio = ra[i];
}
/*
Set values of 'last_index_rec' and 'last_n_pts'.
*/
last_index_rec = index_rec;
last_n_pts = n_pts;
fail:
if (error)
{
if (t != NULL)
{
free(t);
t = NULL;
}
if (ra != NULL)
{
free(ra);
ra = NULL;
}
if (cio_file != NULL)
{
fclose (cio_file);
cio_file = NULL;
}
}
return (error);
}
/********ira_equinox */
double ira_equinox (double jd_tdb, short int equinox,
short int accuracy)
/*
------------------------------------------------------------------------
PURPOSE:
To compute the intermediate right ascension of the equinox at
the input Julian date, using an analytical expression for the
accumulated precession in right ascension. For the true equinox,
the result is the equation of the origins.
REFERENCES:
Capitaine, N. et al. (2003), Astronomy and Astrophysics 412,
567-586, eq. (42).
INPUT
ARGUMENTS:
jd_tdb (double)
TDB Julian date.
equinox (short int)
Equinox selection flag:
= 0 ... mean equinox
= 1 ... true equinox.
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
None.
RETURNED
VALUE:
(double)
Intermediate right ascension of the equinox, in hours (+ or -).
If 'equinox' = 1 (i.e true equinox), then the returned value is
the equation of the origins.
GLOBALS
USED:
T0 novascon.c
FUNCTIONS
CALLED:
e_tilt novas.c
fabs math.h
VER./DATE/
PROGRAMMER:
V1.0/07-06/JAB (USNO/AA)
NOTES:
1. This function is the C version of NOVAS Fortran routine
'eqxra'.
------------------------------------------------------------------------
*/
{
static short int acc_last = 99;
static double t_last = 0.0;
static double eq_eq = 0.0;
double t, u, v, w, x, prec_ra, ra_eq;
/*
Compute time in Julian centuries.
*/
t = (jd_tdb - T0 ) / 36525.0;
/*
For the true equinox, obtain the equation of the equinoxes in time
seconds, which includes the 'complementary terms'.
*/
if (equinox == 1)
{
if (((fabs (jd_tdb - t_last)) > 1.0e-8) || (accuracy != acc_last))
{
e_tilt (jd_tdb,accuracy, &u, &v, &eq_eq, &w, &x);
t_last = jd_tdb;
acc_last = accuracy;
}
}
else
{
eq_eq = 0.0;
}
/*
Precession in RA in arcseconds taken from the reference.
*/
prec_ra = 0.014506 +
(((( - 0.0000000368 * t
- 0.000029956 ) * t
- 0.00000044 ) * t
+ 1.3915817 ) * t
+ 4612.156534 ) * t;
ra_eq = - (prec_ra / 15.0 + eq_eq) / 3600.0;
return (ra_eq);
}
/********ephemeris */
short int ephemeris (double jd[2], object *cel_obj, short int origin,
short int accuracy,
double *pos, double *vel)
/*
------------------------------------------------------------------------
PURPOSE:
Retrieves the position and velocity of a solar system body from
a fundamental ephemeris.
REFERENCES:
None.
INPUT
ARGUMENTS:
jd[2] (double)
TDB Julian date split into two parts, where the sum
jd[0] + jd[1] is the TDB Julian date.
*cel_obj (struct object)
Pointer to structure containing the designation of the body
of interest (defined in novas.h).
origin (short int)
Origin code
= 0 ... solar system barycenter
= 1 ... center of mass of the Sun
accuracy (short int)
Selection for accuracy
= 0 ... full accuracy
= 1 ... reduced accuracy
OUTPUT
ARGUMENTS:
pos[3] (double)
Position vector of the body at 'jd_tdb'; equatorial rectangular
coordinates in AU referred to the ICRS.
vel[3] (double)
Velocity vector of the body at 'jd_tdb'; equatorial rectangular
coordinates in AU/day referred to the ICRS.
RETURNED
VALUE:
(short int)
0 ... Everything OK.
1 ... Invalid value of 'origin'.
2 ... Invalid value of 'type' in 'cel_obj'.
3 ... Unable to allocate memory.
10+n ... where n is the error code from 'solarsystem'.
20+n ... where n is the error code from 'readeph'.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
solarsystem novas.c
solarsystem_hp novas.c
readeph readeph.c
strlen string.h
strcpy string.h
malloc stdlib.h
free stdlib.h
VER./DATE/
PROGRAMMER:
V1.0/06-97/JAB (USNO/AA)
V1.1/10-97/JAB (USNO/AA): Support error code from 'readeph'.
V1.2/08-98/JAB (USNO/AA): Add computation of barycentric
coordinates of a minor planet; support
new 'readeph' argument list.
V1.3/12-99/JAB (USNO/AA): Add error handling to call to
'solarsystem' (case 1).
V1.4/12-02/JAB (USNO/AA): Fix memory leaks with 'posvel' (case 1).
V1.5/10-08/JAB (USNO/AA): Incorporate higher-precision call to
the solar system ephemeris, primarily
to support light-time calculations.
V1.6/02-11/JLB (USNO/AA): Reformatted description of origin for
consistency with other documentation.
NOTES:
1. It is recommended that the input structure 'cel_obj' be
created using function 'make_object' in file novas.c.
------------------------------------------------------------------------
*/
{
char *mp_name;
int err = 0;
int mp_number;
short int error = 0;
short int ss_number, i;
double jd_tdb, *posvel, *sun_pos, *sun_vel;
size_t name_len;
/*
Check the value of 'origin'.
*/
if ((origin < 0) || (origin > 1))
return (error = 1);
/*
Invoke the appropriate ephemeris access software depending upon the
type of object.
*/
switch (cel_obj->type)
{
/*
Get the position and velocity of a major planet, Pluto, Sun, or Moon.
When high accuracy is specified, use function 'solarsystem_hp' rather
than 'solarsystem'.
*/
case 0:
ss_number = cel_obj->number;
if (accuracy == 0)
{
if ((error = solarsystem_hp (jd,ss_number,origin, pos,vel))
!= 0)
error += 10;
}
else
{
jd_tdb = jd[0] + jd[1];
if ((error = solarsystem (jd_tdb,ss_number,origin, pos,vel))
!= 0)
error += 10;
}
break;
/*
Get the position and velocity of a minor planet.
*/
case 1:
mp_number = (int) cel_obj->number;
name_len = (strlen (cel_obj->name) + 1L) * sizeof (char);
mp_name = (char *) malloc (name_len);
if (mp_name == NULL)
return (error = 3);
strcpy (mp_name, cel_obj->name);
/* The USNO minor planet software returns heliocentric positions and
velocities.
*/
jd_tdb = jd[0] + jd[1];
posvel = readeph (mp_number,mp_name,jd_tdb, &err);
if (posvel == NULL)
{
free (mp_name);
return (error = 3);
}
if (err != 0)
{
free (mp_name);
free (posvel);
return ((short int) (20 + err));
}
/* Barycentric coordinates of the minor planet, if desired, are
computed via the barycentric coordinates of the Sun, obtained
from the solar system ephemeris.
*/
if (origin == 0)
{
sun_pos = (double *) malloc (3L * sizeof (double));
if (sun_pos == NULL)
{
free (mp_name);
free (posvel);
return (error = 3);
}
sun_vel = (double *) malloc (3L * sizeof (double));
if (sun_vel == NULL)
{
free (mp_name);
free (posvel);
free (sun_pos);
return (error = 3);
}
if ((error = solarsystem (jd_tdb,10,0, sun_pos,sun_vel)) != 0)
{
free (mp_name);
free (posvel);
free (sun_pos);
free (sun_vel);
return (error += 10);
}
for (i = 0; i < 3; i++)
{
posvel[i] += sun_pos[i];
posvel[i+3] += sun_vel[i];
}
free (sun_pos);
free (sun_vel);
}
/*
Break up 'posvel' into separate position and velocity vectors.
*/
for (i = 0; i < 3; i++)
{
pos[i] = posvel[i];
vel[i] = posvel[i+3];
}
free (mp_name);
free (posvel);
break;
/*
Invalid type of object.
*/
default:
error = 2;
break;
}
return (error);
}
/********transform_hip */
void transform_hip (cat_entry *hipparcos,
cat_entry *hip_2000)
/*
------------------------------------------------------------------------
PURPOSE:
To convert Hipparcos catalog data at epoch J1991.25 to epoch
J2000.0, for use within NOVAS. To be used only for Hipparcos or
Tycho stars with linear space motion. Both input and output data
is in the ICRS.
REFERENCES:
None.
INPUT
ARGUMENTS:
*hipparcos (struct cat_entry)
An entry from the Hipparcos catalog, at epoch J1991.25, with
all members having Hipparcos catalog units. See Note 1
below (struct defined in novas.h).
OUTPUT
ARGUMENTS:
*hip_2000 (struct cat_entry)
The transformed input entry, at epoch J2000.0. See Note 2
below (struct defined in novas.h).
RETURNED
VALUE:
None.
GLOBALS
USED:
T0 novascon.c
FUNCTIONS
CALLED:
transform_cat novas.c
VER./DATE/
PROGRAMMER:
V1.0/03-98/JAB (USNO/AA)
V1.1/11-03/JAB (USNO/AA) Update for ICRS
NOTES:
1. Input (Hipparcos catalog) epoch and units:
Epoch: J1991.25
Right ascension (RA): degrees
Declination (Dec): degrees
Proper motion in RA: milliarcseconds per year
Proper motion in Dec: milliarcseconds per year
Parallax: milliarcseconds
Radial velocity: kilometers per second (not in catalog)
2. Output (modified Hipparcos) epoch and units:
Epoch: J2000.0
Right ascension: hours
Declination: degrees
Proper motion in RA: milliarcseconds per year
Proper motion in Dec: milliarcseconds per year
Parallax: milliarcseconds
Radial velocity: kilometers per second
3. This function is the C version of NOVAS Fortran routine
'gethip'.
------------------------------------------------------------------------
*/
{
const double epoch_hip = 2448349.0625;
cat_entry scratch;
/*
Set up a "scratch" catalog entry containing Hipparcos data in
"NOVAS units."
*/
strcpy (scratch.starname, hipparcos->starname);
scratch.starnumber = hipparcos->starnumber;
scratch.dec = hipparcos->dec;
scratch.promora = hipparcos->promora;
scratch.promodec = hipparcos->promodec;
scratch.parallax = hipparcos->parallax;
scratch.radialvelocity = hipparcos->radialvelocity;
strcpy (scratch.catalog, "SCR");
/*
Convert right ascension from degrees to hours.
*/
scratch.ra = hipparcos->ra / 15.0;
/*
Change the epoch of the Hipparcos data from J1991.25 to J2000.0.
*/
transform_cat (1,epoch_hip,&scratch,T0,"HP2", hip_2000);
return;
}
/********transform_cat */
short int transform_cat (short int option, double date_incat,
cat_entry *incat, double date_newcat,
char newcat_id[SIZE_OF_CAT_NAME],
cat_entry *newcat)
/*
------------------------------------------------------------------------
PURPOSE:
To transform a star's catalog quantities for a change of epoch
and/or equator and equinox. Also used to rotate catalog
quantities on the dynamical equator and equinox of J2000.0 to the
ICRS or vice versa.
REFERENCES:
None.
INPUT
ARGUMENTS:
option (short int)
Transformation option
= 1 ... change epoch; same equator and equinox
= 2 ... change equator and equinox; same epoch
= 3 ... change equator and equinox and epoch
= 4 ... change equator and equinox J2000.0 to ICRS
= 5 ... change ICRS to equator and equinox of J2000.0
date_incat (double)
TT Julian date, or year, of input catalog data.
*incat (struct cat_entry)
An entry from the input catalog, with units as given in
the struct definition (struct defined in novas.h).
date_newcat (double)
TT Julian date, or year, of transformed catalog data.
newcat_id[SIZE_OF_CAT_NAME] (char)
Catalog identifier ((SIZE_OF_CAT_NAME - 1) characters maximum)
e.g. HIP = Hipparcos, TY2 = Tycho-2.
OUTPUT
ARGUMENTS:
*newcat (struct cat_entry)
The transformed catalog entry, with units as given in
the struct definition (struct defined in novas.h).
RETURNED
VALUE:
= 0 ... Everything OK.
= 1 ... Invalid value of an input date for option 2 or 3 (see
Note 1 below).
= 2 ... length of 'newcat_id' out of bounds.
GLOBALS
USED:
SIZE_OF_CAT_NAME novas.h
T0, ASEC2RAD novascon.c
AU_KM, C novascon.c
FUNCTIONS
CALLED:
precession novas.c
frame_tie novas.c
sin math.h
cos math.h
sqrt math.h
atan2 math.h
asin math.h
strcpy string.h
VER./DATE/
PROGRAMMER:
V1.0/03-98/JAB (USNO/AA)
V1.1/11-03/JAB (USNO/AA) Update for ICRS; add options 4 and 5.
V1.2/06-08/WKP (USNO/AA) Changed 'a[3]' to 'd', added Doppler factor,
and corrected 'frame_tie' direction.
V1.3/06-08/JAB (USNO/AA) Error check on dates for options 2 and 3.
V1.4/02-11/WKP (USNO/AA) Implement SIZE_OF_CAT_NAME.
NOTES:
1. 'date_incat' and 'date_newcat' may be specified either as a
Julian date (e.g., 2433282.5) or a Julian year and fraction
(e.g., 1950.0). Values less than 10000 are assumed to be years.
For 'option' = 2 or 'option' = 3, either 'date_incat' or
'date_newcat' must be 2451545.0 or 2000.0 (J2000.0). For
'option' = 4 and 'option' = 5, 'date_incat' and 'date_newcat'
are ignored.
2. 'option' = 1 updates the star's data to account for the
star's space motion between the first and second dates, within a
fixed reference frame.
'option' = 2 applies a rotation of the reference frame
corresponding to precession between the first and second dates,
but leaves the star fixed in space.
'option' = 3 provides both transformations.
'option' = 4 and 'option' = 5 provide a a fixed rotation
about very small angles (<0.1 arcsecond) to take data from the
dynamical system of J2000.0 to the ICRS ('option' = 4) or vice
versa ('option' = 5).
3. For 'option' = 1, input data can be in any fixed reference
system. for 'option' = 2 or 'option' = 3, this function assumes
the input data is in the dynamical system and produces output in
the dynamical system. for 'option' = 4, the input data must be
on the dynamical equator and equinox of J2000.0. for
'option' = 5, the input data must be in the ICRS.
4. This function cannot be properly used to bring data from
old star catalogs into the modern system, because old catalogs
were compiled using a set of constants that are
incompatible with modern values. In particular, it should not
be used for catalogs whose positions and proper motions were
derived by assuming a precession constant significantly different
from the value implicit in function 'precession'.
5. This function is the C version of NOVAS Fortran routine
'catran'.
6. SIZE_OF_CAT_NAME is defined in novas.h.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int j;
double jd_incat, jd_newcat, paralx, dist, r, d, cra, sra, cdc,
sdc, k, pos1[3], term1, pmr, pmd, rvl, vel1[3], pos2[3],
vel2[3], xyproj;
/*
If necessary, compute Julian dates.
This function uses TDB Julian dates internally, but no
distinction between TDB and TT is necessary.
*/
if (date_incat < 10000.0)
jd_incat = T0 + (date_incat - 2000.0) * 365.25;
else
jd_incat = date_incat;
if (date_newcat < 10000.0)
jd_newcat = T0 + (date_newcat - 2000.0) * 365.25;
else
jd_newcat = date_newcat;
/*
Convert input angular components to vectors
If parallax is unknown, undetermined, or zero, set it to 1.0e-6
milliarcsecond, corresponding to a distance of 1 gigaparsec.
*/
paralx = incat->parallax;
if (paralx <= 0.0)
paralx = 1.0e-6;
/*
Convert right ascension, declination, and parallax to position
vector in equatorial system with units of AU.
*/
dist = 1.0 / sin (paralx * 1.0e-3 * ASEC2RAD);
r = incat->ra * 54000.0 * ASEC2RAD;
d = incat->dec * 3600.0 * ASEC2RAD;
cra = cos (r);
sra = sin (r);
cdc = cos (d);
sdc = sin (d);
pos1[0] = dist * cdc * cra;
pos1[1] = dist * cdc * sra;
pos1[2] = dist * sdc;
/*
Compute Doppler factor, which accounts for change in light
travel time to star.
*/
k = 1.0 / (1.0 - incat->radialvelocity / C * 1000.0);
/*
Convert proper motion and radial velocity to orthogonal components
of motion, in spherical polar system at star's original position,
with units of AU/day.
*/
term1 = paralx * 365.25;
pmr = incat->promora / term1 * k;
pmd = incat->promodec / term1 * k;
rvl = incat->radialvelocity * 86400.0 / AU_KM * k;
/*
Transform motion vector to equatorial system.
*/
vel1[0] = - pmr * sra - pmd * sdc * cra + rvl * cdc * cra;
vel1[1] = pmr * cra - pmd * sdc * sra + rvl * cdc * sra;
vel1[2] = pmd * cdc + rvl * sdc;
/*
Update star's position vector for space motion (only if 'option' = 1
or 'option' = 3).
*/
if ((option == 1) || (option == 3))
{
for (j = 0; j < 3; j++)
{
pos2[j] = pos1[j] + vel1[j] * (jd_newcat - jd_incat);
vel2[j] = vel1[j];
}
}
else
{
for (j = 0; j < 3; j++)
{
pos2[j] = pos1[j];
vel2[j] = vel1[j];
}
}
/*
Precess position and velocity vectors (only if 'option' = 2 or
'option' = 3).
*/
if ((option == 2) || (option == 3))
{
for (j = 0; j < 3; j++)
{
pos1[j] = pos2[j];
vel1[j] = vel2[j];
}
if ((error = precession (jd_incat,pos1,jd_newcat, pos2)) != 0)
return (error);
precession (jd_incat,vel1,jd_newcat, vel2);
}
/*
Rotate dynamical J2000.0 position and velocity vectors to ICRS
(only if 'option' = 4).
*/
if (option == 4)
{
frame_tie (pos1,-1, pos2);
frame_tie (vel1,-1, vel2);
}
/*
Rotate ICRS position and velocity vectors to dynamical J2000.0
(only if 'option' = 5).
*/
if (option == 5)
{
frame_tie (pos1,1, pos2);
frame_tie (vel1,1, vel2);
}
/*
Convert vectors back to angular components for output.
From updated position vector, obtain star's new position expressed
as angular quantities.
*/
xyproj = sqrt (pos2[0] * pos2[0] + pos2[1] * pos2[1]);
if (xyproj > 0.0)
r = atan2 (pos2[1], pos2[0]);
else
r = 0.0;
newcat->ra = r / ASEC2RAD / 54000.0;
if (newcat->ra < 0.0)
newcat->ra += 24.0;
if (newcat->ra >= 24.0)
newcat->ra -= 24.0;
d = atan2 (pos2[2], xyproj);
newcat->dec = d / ASEC2RAD / 3600.0;
dist = sqrt (pos2[0] * pos2[0] + pos2[1] * pos2[1] +
pos2[2] * pos2[2]);
paralx = asin (1.0 / dist) / ASEC2RAD * 1000.0;
newcat->parallax = paralx;
/*
Transform motion vector back to spherical polar system at star's
new position.
*/
cra = cos (r);
sra = sin (r);
cdc = cos (d);
sdc = sin (d);
pmr = - vel2[0] * sra + vel2[1] * cra;
pmd = - vel2[0] * cra * sdc - vel2[1] * sra * sdc + vel2[2] * cdc;
rvl = vel2[0] * cra * cdc + vel2[1] * sra * cdc + vel2[2] * sdc;
/*
Convert components of motion to from AU/day to normal catalog units.
*/
newcat->promora = pmr * paralx * 365.25 / k;
newcat->promodec = pmd * paralx * 365.25 / k;
newcat->radialvelocity = rvl * (AU_KM / 86400.0) / k;
/*
Take care of zero-parallax case.
*/
if (newcat->parallax <= 1.01e-6)
{
newcat->parallax = 0.0;
newcat->radialvelocity = incat->radialvelocity;
}
/*
Set the catalog identification code for the transformed catalog
entry.
*/
if ((short int) strlen (newcat_id) > SIZE_OF_CAT_NAME - 1)
return (2);
else
strcpy (newcat->catalog, newcat_id);
/*
Copy unchanged quantities from the input catalog entry to the
transformed catalog entry.
*/
strcpy (newcat->starname, incat->starname);
newcat->starnumber = incat->starnumber;
return (error);
}
/********limb_angle */
void limb_angle (double pos_obj[3], double pos_obs[3],
double *limb_ang, double *nadir_ang)
/*
------------------------------------------------------------------------
PURPOSE:
This function determines the angle of an object above or below
the Earth's limb (horizon). The geometric limb is computed,
assuming the Earth to be an airless sphere (no refraction or
oblateness is included). The observer can be on or above the
Earth. For an observer on the surface of the Earth, this
function returns the approximate unrefracted altitude.
REFERENCES:
None.
INPUT
ARGUMENTS:
pos_obj[3] (double)
Position vector of observed object, with respect to origin at
geocenter, components in AU.
pos_obs[3] (double)
Position vector of observer, with respect to origin at
geocenter, components in AU.
OUTPUT
ARGUMENTS:
*limb_ang (double)
Angle of observed object above (+) or below (-) limb
in degrees.
*nadir_ang (double)
Nadir angle of observed object as a fraction of apparent
radius of limb:
< 1.0 ... below the limb
= 1.0 ... on the limb
> 1.0 ... above the limb
RETURNED
VALUE:
None.
GLOBALS
USED:
TWOPI, ERAD novascon.c
AU, RAD2DEG novascon.c
FUNCTIONS
CALLED:
sqrt math.h
asin math.h
acos math.h
VER./DATE/
PROGRAMMER:
V1.0/09-04/JAB (USNO/AA)
NOTES:
1.This function is the C version of NOVAS Fortran routine
'limang'.
------------------------------------------------------------------------
*/
{
static short int first_entry = 1;
static double pi, halfpi, rade;
double disobj, disobs, aprad, zdlim, coszd, zdobj;
if (first_entry)
{
pi = TWOPI / 2.0;
halfpi = pi / 2.0;
rade = ERAD / AU;
first_entry = 0;
}
/*
Compute the distance to the object and the distance to the observer.
*/
disobj = sqrt (pos_obj[0] * pos_obj[0] +
pos_obj[1] * pos_obj[1] +
pos_obj[2] * pos_obj[2]);
disobs = sqrt (pos_obs[0] * pos_obs[0] +
pos_obs[1] * pos_obs[1] +
pos_obs[2] * pos_obs[2]);
/*
Compute apparent angular radius of Earth's limb.
*/
if (disobs >= rade)
{
aprad = asin (rade / disobs);
}
else
{
aprad = halfpi;
}
/*
Compute zenith distance of Earth's limb.
*/
zdlim = pi - aprad;
/*
Compute zenith distance of observed object.
*/
coszd = (pos_obj[0] * pos_obs[0] + pos_obj[1] * pos_obs[1] +
pos_obj[2] * pos_obs[2]) / (disobj * disobs);
if (coszd <= -1.0)
{
zdobj = pi;
}
else if (coszd >= 1.0)
{
zdobj = 0.0;
}
else
{
zdobj = acos (coszd);
}
/*
Angle of object wrt limb is difference in zenith distances.
*/
*limb_ang = (zdlim - zdobj) * RAD2DEG;
/*
Nadir angle of object as a fraction of angular radius of limb.
*/
*nadir_ang = (pi - zdobj) / aprad;
return;
}
/********refract */
double refract (on_surface *location, short int ref_option,
double zd_obs)
/*
------------------------------------------------------------------------
PURPOSE:
This function computes atmospheric refraction in zenith
distance. This version computes approximate refraction for
optical wavelengths.
REFERENCES:
Explanatory Supplement to the Astronomical Almanac, p. 144.
Bennett, G. (1982), Journal of Navigation (Royal Institute) 35,
pp. 255-259.
INPUT
ARGUMENTS:
*location (struct on_surface)
Pointer to structure containing observer's location. This
structure also contains weather data (optional) for the
observer's location (defined in novas.h).
ref_option (short int)
= 1 ... Use 'standard' atmospheric conditions.
= 2 ... Use atmospheric parameters input in the 'location'
structure.
zd_obs (double)
Observed zenith distance, in degrees.
OUTPUT
ARGUMENTS:
None.
RETURNED
VALUE:
(double)
Atmospheric refraction, in degrees.
GLOBALS
USED:
DEG2RAD novascon.c
FUNCTIONS
CALLED:
exp math.h
tan math.h
VER./DATE/
PROGRAMMER:
V1.0/06-98/JAB (USNO/AA)
NOTES:
1. This function can be used for planning observations or
telescope pointing, but should not be used for the reduction
of precise observations.
2. This function is the C version of NOVAS Fortran routine
'refrac'.
------------------------------------------------------------------------
*/
{
/*
's' is the approximate scale height of atmosphere in meters.
*/
const double s = 9.1e3;
double refr, p, t, h, r;
/*
Compute refraction only for zenith distances between 0.1 and
91 degrees.
*/
if ((zd_obs < 0.1) || (zd_obs > 91.0))
refr = 0.0;
else
{
/*
If observed weather data are available, use them. Otherwise, use
crude estimates of average conditions.
*/
if (ref_option == 2)
{
p = location->pressure;
t = location->temperature;
}
else
{
p = 1010.0 * exp (-location->height / s);
t = 10.0;
}
h = 90.0 - zd_obs;
r = 0.016667 / tan ((h + 7.31 / (h + 4.4)) * DEG2RAD);
refr = r * (0.28 * p / (t + 273.0));
}
return (refr);
}
/********julian_date */
double julian_date (short int year, short int month, short int day,
double hour)
/*
------------------------------------------------------------------------
PURPOSE:
This function will compute the Julian date for a given calendar
date (year, month, day, hour).
REFERENCES:
Fliegel, H. & Van Flandern, T. Comm. of the ACM, Vol. 11, No. 10,
October 1968, p. 657.
INPUT
ARGUMENTS:
year (short int)
Year.
month (short int)
Month number.
day (short int)
Day-of-month.
hour (double)
Hour-of-day.
OUTPUT
ARGUMENTS:
None.
RETURNED
VALUE:
(double)
Julian date.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
None.
VER./DATE/
PROGRAMMER:
V1.0/06-98/JAB (USNO/AA)
V1.1/03-08/WKP (USNO/AA) Updated prolog.
NOTES:
1. This function is the C version of NOVAS Fortran routine
'juldat'.
2. This function makes no checks for a valid input calendar
date.
3. Input calendar date must be Gregorian.
4. Input time value can be based on any UT-like time scale
(UTC, UT1, TT, etc.) - output Julian date will have the same basis.
------------------------------------------------------------------------
*/
{
long int jd12h;
double tjd;
jd12h = (long) day - 32075L + 1461L * ((long) year + 4800L
+ ((long) month - 14L) / 12L) / 4L
+ 367L * ((long) month - 2L - ((long) month - 14L) / 12L * 12L)
/ 12L - 3L * (((long) year + 4900L + ((long) month - 14L) / 12L)
/ 100L) / 4L;
tjd = (double) jd12h - 0.5 + hour / 24.0;
return (tjd);
}
/********cal_date */
void cal_date (double tjd,
short int *year, short int *month, short int *day,
double *hour)
/*
------------------------------------------------------------------------
PURPOSE:
This function will compute a date on the Gregorian calendar given
the Julian date.
REFERENCES:
Fliegel, H. & Van Flandern, T. Comm. of the ACM, Vol. 11, No. 10,
October 1968, p. 657.
INPUT
ARGUMENTS:
tjd (double)
Julian date.
OUTPUT
ARGUMENTS:
*year (short int)
Year.
*month (short int)
Month number.
*day (short int)
Day-of-month.
*hour (double)
Hour-of-day.
RETURNED
VALUE:
None.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
fmod math.h
VER./DATE/
PROGRAMMER:
V1.0/06-98/JAB (USNO/AA)
NOTES:
1. This routine valid for any 'jd' greater than zero.
2. Input Julian date can be based on any UT-like time scale
(UTC, UT1, TT, etc.) - output time value will have same basis.
3. This function is the C version of NOVAS Fortran routine
'caldat'.
------------------------------------------------------------------------
*/
{
long int jd, k, m, n;
double djd;
djd = tjd + 0.5;
jd = (long int) djd;
*hour = fmod (djd,1.0) * 24.0;
k = jd + 68569L;
n = 4L * k / 146097L;
k = k - (146097L * n + 3L) / 4L;
m = 4000L * (k + 1L) / 1461001L;
k = k - 1461L * m / 4L + 31L;
*month = (short int) (80L * k / 2447L);
*day = (short int) (k - 2447L * (long int) *month / 80L);
k = (long int) *month / 11L;
*month = (short int) ((long int) *month + 2L - 12L * k);
*year = (short int) (100L * (n - 49L) + m + k);
return;
}
/********norm_ang */
double norm_ang (double angle)
/*
------------------------------------------------------------------------
PURPOSE:
Normalize angle into the range 0 <= angle < (2 * pi).
REFERENCES:
None.
INPUT
ARGUMENTS:
angle (double)
Input angle (radians).
OUTPUT
ARGUMENTS:
None.
RETURNED
VALUE:
(double)
The input angle, normalized as described above (radians).
GLOBALS
USED:
TWOPI novascon.c
FUNCTIONS
CALLED:
fmod math.h
VER./DATE/
PROGRAMMER:
V1.0/09-03/JAB (USNO/AA)
NOTES:
None.
------------------------------------------------------------------------
*/
{
double a;
a = fmod (angle,TWOPI);
if (a < 0.0)
a += TWOPI;
return (a);
}
/********make_cat_entry */
short int make_cat_entry (const char *star_name,
const char *catalog,
long int star_num, double ra, double dec,
double pm_ra, double pm_dec, double parallax,
double rad_vel,
cat_entry *star)
/*
------------------------------------------------------------------------
PURPOSE:
To create a structure of type 'cat_entry' containing catalog
data for a star or "star-like" object.
REFERENCES:
None.
INPUT
ARGUMENTS:
star_name[SIZE_OF_OBJ_NAME] (char)
Object name ((SIZE_OF_OBJ_NAME - 1) characters maximum).
catalog[SIZE_OF_CAT_NAME] (char)
Catalog identifier ((SIZE_OF_CAT_NAME - 1) characters maximum)
e.g. HIP = Hipparcos, TY2 = Tycho-2.
star_num (long int)
Object number in the catalog.
ra (double)
Right ascension of the object (hours).
dec (double)
Declination of the object (degrees).
pm_ra (double)
Proper motion in right ascension (milliarcseconds/year).
pm_dec (double)
Proper motion in declination (milliarcseconds/year).
parallax (double)
Parallax (milliarcseconds).
rad_vel (double)
Radial velocity (kilometers/second).
OUTPUT
ARGUMENTS:
*star (struct cat_entry)
Structure containing the input data (structure defined in
novas.h).
RETURNED
VALUE:
(short int)
= 0 ... no problems.
= 1 ... length of 'star_name' out of bounds.
= 2 ... length of 'catalog' out of bounds.
GLOBALS
USED:
SIZE_OF_OBJ_NAME novas.h
SIZE_OF_CAT_NAME novas.h
FUNCTIONS
CALLED:
strlen string.h
strcpy string.h
VER./DATE/
PROGRAMMER:
V1.0/03-98/JAB (USNO/AA)
V1.1/08-04/JAB (USNO/AA)
V1.3/05-10/JAB (USNO/AA): Fix bugs in set-up of 'star->starname'
and 'star->catalog'.
NOTES:
1. This utility function creates a single data structure
encapsulating the input data.
2. SIZE_OF_OBJ_NAME and SIZE_OF_CAT_NAME are defined in novas.h.
------------------------------------------------------------------------
*/
{
/*
Set up the 'star' structure.
*/
if ((short int) strlen (star_name) > SIZE_OF_OBJ_NAME - 1)
return (1);
else
strcpy (star->starname, star_name);
if ((short int) strlen (catalog) > SIZE_OF_CAT_NAME - 1)
return (2);
else
strcpy (star->catalog, catalog);
star->starnumber = star_num;
star->ra = ra;
star->dec = dec;
star->promora = pm_ra;
star->promodec = pm_dec;
star->parallax = parallax;
star->radialvelocity = rad_vel;
return (0);
}
/********make_object */
short int make_object (short int type, short int number,
char name[/*SIZE_OF_OBJ_NAME*/],
cat_entry *star_data,
object *cel_obj)
/*
------------------------------------------------------------------------
PURPOSE:
Makes a structure of type 'object' - specifying a celestial
object - based on the input parameters.
REFERENCES:
None.
INPUT
ARGUMENTS:
type (short int)
Type of object
= 0 ... major planet, Pluto, Sun, or Moon
= 1 ... minor planet
= 2 ... object located outside the solar system
(e.g. star, galaxy, nebula, etc.)
number (short int)
Body number
For 'type' = 0: Mercury = 1,...,Pluto = 9, Sun = 10,
Moon = 11
For 'type' = 1: minor planet number
For 'type' = 2: set to 0 (zero)
name[SIZE_OF_OBJ_NAME] (char)
Name of the object ((SIZE_OF_OBJ_NAME - 1) characters maximum).
*star_data (struct cat_entry)
Structure containing basic astrometric data for any celestial
object located outside the solar system; the catalog
data for a star (defined in novas.h).
OUTPUT
ARGUMENTS:
struct object *cel_obj
Structure containing the object definition (defined in novas.h).
RETURNED
VALUE:
(short int)
= 0 ... everything OK
= 1 ... invalid value of 'type'
= 2 ... 'number' out of range
= 3 ... Initialization of 'cel_obj' failed (object name).
= 4 ... Initialization of 'cel_obj' failed (catalog name).
= 5 ... 'name' is out of string bounds.
GLOBALS
USED:
SIZE_OF_OBJ_NAME novas.h
SIZE_OF_CAT_NAME novas.h
FUNCTIONS
CALLED:
strlen string.h
strcpy string.h
toupper ctype.h
VER./DATE/
PROGRAMMER:
V1.0/06-97/JAB (USNO/AA)
V1.1/10-98/JAB (USNO/AA): Change object name to mixed case.
V1.2/12-02/JAB (USNO/AA): Fix bug by adding brackets in the "Set
the object number" 'if' block.
V1.3/07-04/JAB (USNO/AA): Change name from 'set_body'; expand to
include objects outside solar system.
V1.4/05-10/JAB (USNO/AA): Implement SIZE_OF_OBJ_NAME and
SIZE_OF_CAT_NAME; check length of
'name'.
V1.5/09-10/WKP (USNO/AA): Explicitly cast 'toupper' return value
to type 'char' to silence compiler.
NOTES:
1. SIZE_OF_OBJ_NAME and SIZE_OF_CAT_NAME are defined in novas.h.
------------------------------------------------------------------------
*/
{
short int error = 0;
short int i;
/*
Initialize the structure.
*/
cel_obj->type = 0;
cel_obj->number = 0;
strcpy (cel_obj->name, " ");
if ((short int) strlen (" ") > SIZE_OF_OBJ_NAME - 1)
return (error = 3);
else
{
strcpy (cel_obj->name, " ");
strcpy (cel_obj->star.starname, " ");
}
if ((short int) strlen (" ") > SIZE_OF_CAT_NAME - 1)
return (error = 4);
else
strcpy (cel_obj->star.catalog, " ");
cel_obj->star.starnumber = 0L;
cel_obj->star.ra = 0.0;
cel_obj->star.dec = 0.0;
cel_obj->star.promora = 0.0;
cel_obj->star.promodec = 0.0;
cel_obj->star.parallax = 0.0;
cel_obj->star.radialvelocity = 0.0;
/*
Set the object type.
*/
if ((type < 0) || (type > 2))
return (error = 1);
else
cel_obj->type = type;
/*
Set the object number.
*/
if (type == 0)
{
if ((number <= 0) || (number > 11))
return (error = 2);
}
else if (type == 1)
{
if (number <= 0)
return (error = 2);
}
else
{
number = 0;
}
cel_obj->number = number;
/*
Check length of 'name'. Set the object name in upper case.
*/
if ((short int) strlen (name) > SIZE_OF_OBJ_NAME - 1)
return (error = 5);
for (i = 0; i < SIZE_OF_OBJ_NAME - 1; i++)
{
cel_obj->name[i] = (char) toupper (name[i]);
if (name[i] == '\0')
break;
}
cel_obj->name[i] = '\0';
/*
Populate the astrometric-data structure if the object is 'type' = 2.
*/
if (type == 2)
{
strcpy (cel_obj->star.starname, star_data->starname);
strcpy (cel_obj->star.catalog, star_data->catalog);
cel_obj->star.starnumber = star_data->starnumber;
cel_obj->star.ra = star_data->ra;
cel_obj->star.dec = star_data->dec;
cel_obj->star.promora = star_data->promora;
cel_obj->star.promodec = star_data->promodec;
cel_obj->star.parallax = star_data->parallax;
cel_obj->star.radialvelocity = star_data->radialvelocity;
}
return (error);
}
/********make_observer */
short int make_observer (short int where, on_surface *obs_surface,
in_space *obs_space,
observer *obs)
/*
------------------------------------------------------------------------
PURPOSE:
Makes a structure of type 'observer' - specifying the location of
the observer.
REFERENCES:
None.
INPUT
ARGUMENTS:
where (short int)
Integer code specifying location of observer.
= 0: observer at geocenter
= 1: observer on surface of earth
= 2: observer on near-earth spacecraft
*obs_surface (struct on_surface)
Structure containing data for an observer's location on the
surface of the Earth; used when 'where' = 1 (defined in
novas.h).
*obs_space (struct in_space)
Data for an observer's location on a near-Earth spacecraft;
used when 'where' = 2 (defined in novas.h).
OUTPUT
ARGUMENTS:
struct observer *obs
Structure specifying the location of the observer (defined in
novas.h).
RETURNED
VALUE:
(short int)
= 0 ... everything OK.
= 1 ... input value of 'where' is out-of-range.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
None.
VER./DATE/
PROGRAMMER:
V1.0/08-04/JAB (USNO/AA)
NOTES:
None.
------------------------------------------------------------------------
*/
{
short int error = 0;
/*
Initialize the output structure.
*/
obs->where = where;
obs->on_surf.latitude = 0.0;
obs->on_surf.longitude = 0.0;
obs->on_surf.height = 0.0;
obs->on_surf.temperature = 0.0;
obs->on_surf.pressure = 0.0;
obs->near_earth.sc_pos[0] = 0.0;
obs->near_earth.sc_pos[1] = 0.0;
obs->near_earth.sc_pos[2] = 0.0;
obs->near_earth.sc_vel[0] = 0.0;
obs->near_earth.sc_vel[1] = 0.0;
obs->near_earth.sc_vel[2] = 0.0;
/*
Populate the output structure based on the value of 'where'.
*/
switch (where)
{
case (0): /* Geocentric */
break;
case (1): /* On surface of Earth */
obs->on_surf.latitude = obs_surface->latitude;
obs->on_surf.longitude = obs_surface->longitude;
obs->on_surf.height = obs_surface->height;
obs->on_surf.temperature = obs_surface->temperature;
obs->on_surf.pressure = obs_surface->pressure;
break;
case (2): /* In near-Earth spacecraft */
obs->near_earth.sc_pos[0] = obs_space->sc_pos[0];
obs->near_earth.sc_pos[1] = obs_space->sc_pos[1];
obs->near_earth.sc_pos[2] = obs_space->sc_pos[2];
obs->near_earth.sc_vel[0] = obs_space->sc_vel[0];
obs->near_earth.sc_vel[1] = obs_space->sc_vel[1];
obs->near_earth.sc_vel[2] = obs_space->sc_vel[2];
break;
default:
error = 1;
break;
}
return (error);
}
/********make_observer_at_geocenter */
void make_observer_at_geocenter (
observer *obs_at_geocenter)
/*
------------------------------------------------------------------------
PURPOSE:
Makes a structure of type 'observer' specifying an observer at the
geocenter.
REFERENCES:
None.
INPUT
ARGUMENTS:
None.
OUTPUT
ARGUMENTS:
struct observer *obs_at_geocenter
Structure specifying the location of the observer at the geocenter
(defined in novas.h).
RETURNED
VALUE:
None.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
make_in_space novas.c
make_on_surface novas.c
VER./DATE/
PROGRAMMER:
V1.0/03-07/JAB (USNO/AA)
NOTES:
1. The output data structure is one of the inputs to NOVAS-C
function 'place'.
------------------------------------------------------------------------
*/
{
double latitude = 0.0;
double longitude = 0.0;
double height = 0.0;
double temperature = 0.0;
double pressure = 0.0;
double sat_pos[3] = {0.0, 0.0, 0.0};
double sat_vel[3] = {0.0, 0.0, 0.0};
in_space sat_state;
on_surface surface_loc;
make_in_space (sat_pos,sat_vel, &sat_state);
make_on_surface (latitude,longitude,height,temperature,pressure,
&surface_loc);
obs_at_geocenter->where = 0;
obs_at_geocenter->on_surf = surface_loc;
obs_at_geocenter->near_earth = sat_state;
}
/********make_observer_on_surface */
void make_observer_on_surface (double latitude, double longitude,
double height, double temperature,
double pressure,
observer *obs_on_surface)
/*
------------------------------------------------------------------------
PURPOSE:
Makes a structure of type 'observer' specifying the location
of and weather for an observer on the surface of the Earth.
REFERENCES:
None.
INPUT
ARGUMENTS:
latitude (double)
Geodetic (ITRS) latitude in degrees; north positive.
longitude (double)
Geodetic (ITRS) longitude in degrees; east positive.
height (double)
Height of the observer (meters).
temperature (double)
Temperature (degrees Celsius).
pressure (double)
Atmospheric pressure (millibars).
OUTPUT
ARGUMENTS:
struct observer *obs_on_surface
Structure containing the location of and weather for an
observer on the surface of the Earth (defined in novas.h).
RETURNED
VALUE:
None.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
make_in_space novas.c
make_on_surface novas.c
VER./DATE/
PROGRAMMER:
V1.0/03-07/JAB (USNO/AA)
NOTES:
1. The output data structure is one of the inputs to NOVAS-C
function 'place'.
------------------------------------------------------------------------
*/
{
double sat_pos[3] = {0.0, 0.0, 0.0};
double sat_vel[3] = {0.0, 0.0, 0.0};
in_space sat_state;
on_surface surface_loc;
make_in_space (sat_pos,sat_vel, &sat_state);
make_on_surface (latitude,longitude,height,temperature,pressure,
&surface_loc);
obs_on_surface->where = 1;
obs_on_surface->on_surf = surface_loc;
obs_on_surface->near_earth = sat_state;
}
/********make_observer_in_space */
void make_observer_in_space (double sc_pos[3], double sc_vel[3],
observer *obs_in_space)
/*
------------------------------------------------------------------------
PURPOSE:
Makes a structure of type 'observer' specifying the position and
velocity of an observer situated on a near-Earth spacecraft.
REFERENCES:
None.
INPUT
ARGUMENTS:
sc_pos[3] (double)
Geocentric position vector (x, y, z) in km.
sc_vel[3] (double)
Geocentric velocity vector (x_dot, y_dot, z_dot) in km/s.
OUTPUT
ARGUMENTS:
struct observer *obs_in_space
Structure containing the position and velocity of an observer
situated on a near-Earth spacecraft (defined in novas.h).
RETURNED
VALUE:
None.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
make_in_space novas.c
make_on_surface novas.c
VER./DATE/
PROGRAMMER:
V1.0/03-07/JAB (USNO/AA)
NOTES:
1. Both input vectors are with respect to true equator and equinox
of date.
2. The output data structure is one of the inputs to NOVAS-C
function 'place'.
------------------------------------------------------------------------
*/
{
const double latitude = 0.0;
const double longitude = 0.0;
const double height = 0.0;
const double temperature = 0.0;
const double pressure = 0.0;
in_space sat_state;
on_surface surface_loc;
make_in_space (sc_pos,sc_vel, &sat_state);
make_on_surface (latitude,longitude,height,temperature,pressure,
&surface_loc);
obs_in_space->where = 2;
obs_in_space->on_surf = surface_loc;
obs_in_space->near_earth = sat_state;
}
/********make_on_surface */
void make_on_surface (double latitude, double longitude, double height,
double temperature, double pressure,
on_surface *obs_surface)
/*
------------------------------------------------------------------------
PURPOSE:
Makes a structure of type 'on_surface' - specifying the location
of and weather for an observer on the surface of the Earth.
REFERENCES:
None.
INPUT
ARGUMENTS:
latitude (double)
Geodetic (ITRS) latitude in degrees; north positive.
longitude (double)
Geodetic (ITRS) longitude in degrees; east positive.
height (double)
Height of the observer (meters).
temperature (double)
Temperature (degrees Celsius).
pressure (double)
Atmospheric pressure (millibars).
OUTPUT
ARGUMENTS:
struct on_surface *obs_surface
Structure containing the location of and weather for an
observer on the surface of the Earth (defined in novas.h).
RETURNED
VALUE:
None.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
None.
VER./DATE/
PROGRAMMER:
V1.0/08-04/JAB (USNO/AA)
NOTES:
None.
------------------------------------------------------------------------
*/
{
obs_surface->latitude = latitude;
obs_surface->longitude = longitude;
obs_surface->height = height;
obs_surface->temperature = temperature;
obs_surface->pressure = pressure;
}
/********make_in_space */
void make_in_space (double sc_pos[3], double sc_vel[3],
in_space *obs_space)
/*
------------------------------------------------------------------------
PURPOSE:
Makes a structure of type 'in_space' - specifying the position and
velocity of an observer situated on a near-Earth spacecraft.
REFERENCES:
None.
INPUT
ARGUMENTS:
sc_pos[3] (double)
Geocentric position vector (x, y, z) in km.
sc_vel[3] (double)
Geocentric velocity vector (x_dot, y_dot, z_dot) in km/s.
OUTPUT
ARGUMENTS:
struct in_space *obs_space
Structure containing the position and velocity of an observer
situated on a near-Earth spacecraft (defined in novas.h).
RETURNED
VALUE:
None.
GLOBALS
USED:
None.
FUNCTIONS
CALLED:
None.
VER./DATE/
PROGRAMMER:
V1.0/08-04/JAB (USNO/AA)
NOTES:
1. Both vectors with respect to true equator and equinox of date.
------------------------------------------------------------------------
*/
{
obs_space->sc_pos[0] = sc_pos[0];
obs_space->sc_pos[1] = sc_pos[1];
obs_space->sc_pos[2] = sc_pos[2];
obs_space->sc_vel[0] = sc_vel[0];
obs_space->sc_vel[1] = sc_vel[1];
obs_space->sc_vel[2] = sc_vel[2];
}
Reference in New Issue View Git Blame Copy Permalink