NOVAS_novas.c

/*
  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_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;

/*
   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;
         }

         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:             /* 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 / 3600.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];

/*
   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 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)
         return (error = 1);

/*
   Read the file header.
*/

      fread (&jd_beg, double_size, (size_t) 1, cio_file);
      fread (&jd_end, double_size, (size_t) 1, cio_file);
      fread (&t_int, double_size, (size_t) 1, cio_file);
      fread (&n_recs, long_size, (size_t) 1, cio_file);
   }

/*
   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'.
*/

   del_n_pts = abs (n_pts - last_n_pts);

/*
   Allocate memory for the 't' and 'ra' arrays.
*/

   if (del_n_pts != 0L)
   {
      if (!first_call)
      {
         free (t);
         free (ra);
      }

      t = (double *) calloc ((size_t) n_pts, double_size);
      if (t == NULL )
      {
         fclose (cio_file);
         return (error = 4);
      }

      ra = (double *) calloc ((size_t) n_pts, double_size);
      if (ra == NULL )
      {
         free (t);
         fclose (cio_file);
         return (error = 5);
      }

      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;
   abs_del_index = abs (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))
   {
      fseek (cio_file, bytes_to_lo, SEEK_SET);

      for (i = 0L; i < n_pts; i++)
      {
         fread (&t[i], double_size, (size_t) 1, cio_file);
         fread (&ra[i], double_size, (size_t) 1, cio_file);
      }
   }

/*
   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))
   {
      n_swap = abs (n_pts - abs_del_index);
      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++)
         {
            fread (&t[i], double_size, (size_t) 1, cio_file);
            fread (&ra[i], double_size, (size_t) 1, cio_file);
         }
      }

/*
   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;
         }

         fseek (cio_file, bytes_to_lo + (n_swap * record_size),
            SEEK_SET);

         j = i++;
         for (i = j; i < n_pts; i++)
         {
            fread (&t[i], double_size, (size_t) 1, cio_file);
            fread (&ra[i], double_size, (size_t) 1, cio_file);
         }
      }
   }

/*
   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;

   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 (char star_name[SIZE_OF_OBJ_NAME],
                          char catalog[SIZE_OF_CAT_NAME],
                          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];
}