text: /*
    Naval Observatory Vector Astrometry Software (NOVAS)
    C Edition, Version 3.1
  
    eph_manager.c: C version of JPL Ephemeris Manager for use with solsys1
  
    U. S. Naval Observatory
    Astronomical Applications Dept.
    Washington, DC
    http://www.usno.navy.mil/USNO/astronomical-applications
  */
  
  #ifndef _EPHMAN_
     #include "eph_manager.h"
  #endif
  
  #if 1
  
  // INDIGO hack, JPLEPH.421 is embeded into executable
  
  extern char binary_jpleph_start[];
  extern char binary_jpleph_end[];
  
  #if defined(INDIGO_LINUX)
  
  __asm__(".section \".rodata\", \"\"");
  __asm__("binary_jpleph_start:");
  __asm__(".incbin \"JPLEPH.421\"");
  __asm__("binary_jpleph_end:");
  __asm__(".previous");
  
  #define fopen(name, mode) fmemopen(binary_jpleph_start, binary_jpleph_end - binary_jpleph_start, mode)
  #endif
  
  #if defined(INDIGO_MACOS)
  
  __asm__(".section \".rodata\", \"\"");
  __asm__("_binary_jpleph_start:");
  __asm__(".incbin \"JPLEPH.421\"");
  __asm__("_binary_jpleph_end:");
  __asm__(".previous");
  
  #define fopen(name, mode) _fmemopen(binary_jpleph_start, binary_jpleph_end - binary_jpleph_start, mode)
  
  /*
   * Copyright (c) 2017  Joachim Nilsson <troglobit@gmail.com>
   *
   * Permission to use, copy, modify, and/or distribute this software for any
   * purpose with or without fee is hereby granted, provided that the above
   * copyright notice and this permission notice appear in all copies.
   *
   * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
   * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
   * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
   * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
   * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
   * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
   * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
   */
  
  #include <memory.h>
  
  struct ops {
  	char   *buf;
  	size_t  len, pos;
  };
  
  typedef struct ops ops_t;
  
  static int readfn(void *arg, char *buf, int len) {
  	int sz;
  	ops_t *ops = (ops_t *)arg;
  	sz = (int)(ops->len - ops->pos);
  	if (sz < 0)
  		sz = 0;
  	if (len > sz)
  		len = sz;
  	memcpy(buf, &ops->buf[ops->pos], len);
  	ops->pos += len;
  	return len;
  }
  
  static int writefn(void *arg, const char *buf, int len) {
  	int sz;
  	ops_t *ops = (ops_t *)arg;
  	sz = (int)(ops->len - ops->pos);
  	if (sz < 0)
  		sz = 0;
  	if (len > sz)
  		len = sz;
  	memcpy(&ops->buf[ops->pos], buf, len);
  	ops->pos += len;
  	return len;
  }
  
  static fpos_t seekfn(void *arg, fpos_t offset, int whence) {
  	fpos_t pos;
  	ops_t *ops = (ops_t *)arg;
  	switch (whence) {
  		case SEEK_SET:
  			pos = offset;
  			break;
  		case SEEK_END:
  			pos = ops->len + offset;
  			break;
  		case SEEK_CUR:
  			pos = ops->pos + offset;
  			break;
  		default:
  			return -1;
  	}
  	if (pos < 0 || (size_t)pos > ops->len) {
  		ops->pos = 0;
  		return -1;
  	}
  	ops->pos = pos;
  	return 0;
  }
  
  static int closefn(void *arg) {
  	free(arg);
  	return 0;
  }
  
  static FILE *_fmemopen(void *buf, size_t len, const char *type) {
  	ops_t *ops = malloc(sizeof(*ops));
  	if (!ops)
  		return NULL;
  	memset(ops, 0, sizeof(*ops));
  	ops->buf = buf;
  	ops->len = len;
  	ops->pos = 0;
  	return funopen(ops, readfn, writefn, seekfn, closefn);
  }
  
  #endif
  
  #endif
  
  
  /*
     Define global variables
  */
  
  short int KM;
  
  /*
     IPT and LPT defined as int to support 64 bit systems.
  */
  
  int IPT[3][12], LPT[3];
  
  long int NRL, NP, NV;
  long int RECORD_LENGTH;
  
  double SS[3], JPLAU, PC[18], VC[18], TWOT, EM_RATIO;
  double *BUFFER;
  
  FILE *EPHFILE = NULL;
  
  /********ephem_open */
  
  short int ephem_open (char *ephem_name,
  
                        double *jd_begin, double *jd_end,
                        short int *de_number)
  /*
  ------------------------------------------------------------------------
  
     PURPOSE:
        This function opens a JPL planetary ephemeris file and
        sets initial values.  This function must be called
        prior to calls to the other JPL ephemeris functions.
  
     REFERENCES:
        Standish, E.M. and Newhall, X X (1988). "The JPL Export
           Planetary Ephemeris"; JPL document dated 17 June 1988.
  
     INPUT
     ARGUMENTS:
        *ephem_name (char)
           Name of the direct-access ephemeris file.
  
     OUTPUT
     ARGUMENTS:
        *jd_begin (double)
           Beginning Julian date of the ephemeris file.
        *jd_end (double)
           Ending Julian date of the ephemeris file.
        *de_number (short int)
           DE number of the ephemeris file opened.
  
     RETURNED
     VALUE:
        (short int)
            0   ...file exists and is opened correctly.
            1   ...file does not exist/not found.
            2-10...error reading from file header.
            11  ...unable to set record length; ephemeris (DE number)
                   not in look-up table.
  
     GLOBALS
     USED:
        SS                eph_manager.h
        JPLAU             eph_manager.h
        PC                eph_manager.h
        VC                eph_manager.h
        TWOT              eph_manager.h
        EM_RATIO          eph_manager.h
        BUFFER            eph_manager.h
        IPT               eph_manager.h
        LPT               eph_manager.h
        NRL               eph_manager.h
        KM                eph_manager.h
        NP                eph_manager.h
        NV                eph_manager.h
        RECORD_LENGTH     eph_manager.h
        EPHFILE           eph_manager.h
  
     FUNCTIONS
     CALLED:
        fclose            stdio.h
        free              stdlib.h
        fopen             stdio.h
        fread             stdio.h
        calloc            stdlib.h
  
     VER./DATE/
     PROGRAMMER:
        V1.0/06-90/JAB (USNO/NA)
        V1.1/06-92/JAB (USNO/AA): Restructure and add initializations.
        V1.2/07-98/WTH (USNO/AA): Modified to open files for different
                                  ephemeris types. (200,403,404,405,406)
        V1.3/11-07/WKP (USNO/AA): Updated prolog.
        V1.4/09-10/WKP (USNO/AA): Changed ncon and denum variables and
                                  sizeof ipt array to type 'int' for
                                  64-bit system compatibility.
        V1.5/09-10/WTH (USNO/AA): Added support for DE421, default case
                                  for switch, close file on error.
        V1.6/10-10/WKP (USNO/AA): Renamed function to lowercase to
                                  comply with coding standards.
  
     NOTES:
        KM...flag defining physical units of the output states.
           = 1, km and km/sec
           = 0, AU and AU/day
        Default value is 0 (KM determines time unit for nutations.
                            Angle unit is always radians.)
  
  ------------------------------------------------------------------------
  */
  {
     char ttl[252], cnam[2400];
  
     short int i, j;
  
     int ncon, denum;
  
     if (EPHFILE)
     {
        fclose (EPHFILE);
        free (BUFFER);
     }
  
  /*
     Open file ephem_name.
  */
  
     if ((EPHFILE = fopen (ephem_name, "rb")) == NULL)
     {
        return 1;
     }
      else
     {
  
  /*
     File found. Set initializations and default values.
  */
  
        KM = 0;
  
        NRL = 0;
  
        NP = 2;
        NV = 3;
        TWOT = 0.0;
  
        for (i = 0; i < 18; i++)
        {
           PC[i] = 0.0;
           VC[i] = 0.0;
        }
  
        PC[0] = 1.0;
        VC[1] = 1.0;
  
  /*
     Read in values from the first record, aka the header.
  */
  
        if (fread (ttl, sizeof ttl, 1, EPHFILE) != 1)
        {
           fclose (EPHFILE);
           return 2;
        }
        if (fread (cnam, sizeof cnam, 1, EPHFILE) != 1)
        {
           fclose (EPHFILE);
           return 3;
        }
        if (fread (SS, sizeof SS, 1, EPHFILE) != 1)
        {
           fclose (EPHFILE);
           return 4;
        }
        if (fread (&ncon, sizeof ncon, 1, EPHFILE) != 1)
        {
           fclose (EPHFILE);
           return 5;
        }
        if (fread (&JPLAU, sizeof JPLAU, 1, EPHFILE) != 1)
        {
           fclose (EPHFILE);
           return 6;
        }
        if (fread (&EM_RATIO, sizeof EM_RATIO, 1, EPHFILE) != 1)
        {
           fclose (EPHFILE);
           return 7;
        }
        for (i = 0; i < 12; i++)
           for (j = 0; j < 3; j++)
              if (fread (&IPT[j][i], sizeof(int), 1, EPHFILE) != 1)
              {
                 fclose (EPHFILE);
                 return 8;
              }
        if (fread (&denum, sizeof denum, 1, EPHFILE) != 1)
        {
           fclose (EPHFILE);
           return 9;
        }
        if (fread (LPT, sizeof LPT, 1, EPHFILE) != 1)
        {
           fclose (EPHFILE);
           return 10;
        }
  
  /*
     Set the value of the record length according to what JPL ephemeris is
     being opened.
  */
  
        switch (denum)
        {
           case 200:
              RECORD_LENGTH = 6608;
              break;
           case 403: case 405:
           case 421:
              RECORD_LENGTH = 8144;
              break;
           case 404: case 406:
              RECORD_LENGTH = 5824;
              break;
  
  /*
     An unknown DE file was opened. Close the file and return an error
     code.
  */
  
           default:
              *jd_begin = 0.0;
              *jd_end = 0.0;
              *de_number = 0;
              fclose (EPHFILE);
              return 11;
              break;
        }
  
        BUFFER = (double *) calloc (RECORD_LENGTH / 8, sizeof(double));
  
        *de_number = (short int) denum;
        *jd_begin = SS[0];
        *jd_end = SS[1];
     }
  
     return 0;
  }
  
  /********ephem_close */
  
  short int ephem_close (void)
  /*
  ------------------------------------------------------------------------
  
     PURPOSE:
        This function closes a JPL planetary ephemeris file and frees the
        memory.
  
     REFERENCES:
        Standish, E.M. and Newhall, X X (1988). "The JPL Export
           Planetary Ephemeris"; JPL document dated 17 June 1988.
  
     INPUT
     ARGUMENTS:
        None.
  
     OUTPUT
     ARGUMENTS:
        None.
  
     RETURNED
     VALUE:
        (short int)
            0  ...file was already closed or closed correctly.
            EOF...error closing file.
  
     GLOBALS
     USED:
        BUFFER            eph_manager.h
        EPHFILE           eph_manager.h
  
     FUNCTIONS
     CALLED:
        fclose            stdio.h
        free              stdlib.h
  
     VER./DATE/
     PROGRAMMER:
        V1.0/11-07/WKP (USNO/AA)
        V1.1/09-10/WKP (USNO/AA): Explicitly cast fclose return value to
                                  type 'short int'.
        V1.2/10-10/WKP (USNO/AA): Renamed function to lowercase to
                                  comply with coding standards.
  
     NOTES:
        None.
  
  ------------------------------------------------------------------------
  */
  {
     short int error = 0;
  
     if (EPHFILE)
     {
        error =  (short int) fclose (EPHFILE);
        free (BUFFER);
     }
     return error;
  }
  
  /********planet_ephemeris */
  
  short int planet_ephemeris (double tjd[2], short int target,
                              short int center,
  
                              double *position, double *velocity)
  /*
  ------------------------------------------------------------------------
  
     PURPOSE:
        This function accesses the JPL planetary ephemeris to give the
        position and velocity of the target object with respect to the
        center object.
  
     REFERENCES:
        Standish, E.M. and Newhall, X X (1988). "The JPL Export
           Planetary Ephemeris"; JPL document dated 17 June 1988.
  
     INPUT
     ARGUMENTS:
        tjd[2] (double)
           Two-element array containing the Julian date, which may be
           split any way (although the first element is usually the
           "integer" part, and the second element is the "fractional"
           part).  Julian date is in the TDB or "T_eph" time scale.
        target (short int)
           Number of 'target' point.
        center (short int)
           Number of 'center' (origin) point.
           The numbering convention for 'target' and'center' is:
              0 = Mercury           7 = Neptune
              1 = Venus             8 = Pluto
              2 = Earth             9 = Moon
              3 = Mars             10 = Sun
              4 = Jupiter          11 = Solar system bary.
              5 = Saturn           12 = Earth-Moon bary.
              6 = Uranus           13 = Nutations (long int. and obliq.)
              (If nutations are desired, set 'target' = 13;
               'center' will be ignored on that call.)
  
     OUTPUT
     ARGUMENTS:
        *position (double)
           Position vector array of target relative to center, measured
           in AU.
        *velocity (double)
           Velocity vector array of target relative to center, measured
           in AU/day.
  
     RETURNED
     VALUE:
        (short int)
           0  ...everything OK.
           1,2...error returned from State.
  
     GLOBALS
     USED:
        EM_RATIO          eph_manager.h
  
     FUNCTIONS
     CALLED:
        state             eph_manager.h
  
     VER./DATE/
     PROGRAMMER:
        V1.0/03-93/WTH (USNO/AA): Convert FORTRAN to C.
        V1.1/07-93/WTH (USNO/AA): Update to C standards.
        V2.0/07-98/WTH (USNO/AA): Modified for ease of use and linearity.
        V3.0/11-06/JAB (USNO/AA): Allowed for use of input 'split' Julian
                                  date for higher precision.
        V3.1/11-07/WKP (USNO/AA): Updated prolog and error codes.
        V3.1/12-07/WKP (USNO/AA): Removed unreferenced variables.
        V3.2/10-10/WKP (USNO/AA): Renamed function to lowercase to
                                  comply with coding standards.
  
     NOTES:
        None.
  
  ------------------------------------------------------------------------
  */
  {
     short int i, error = 0, earth = 2, moon = 9;
     short int do_earth = 0, do_moon = 0;
  
     double jed[2];
     double pos_moon[3] = {0.0,0.0,0.0}, vel_moon[3] = {0.0,0.0,0.0},
            pos_earth[3] = {0.0,0.0,0.0}, vel_earth[3] = {0.0,0.0,0.0};
     double target_pos[3] = {0.0,0.0,0.0}, target_vel[3] = {0.0,0.0,0.0},
            center_pos[3] = {0.0,0.0,0.0}, center_vel[3] = {0.0,0.0,0.0};
  
  /*
     Initialize 'jed' for 'state' and set up component count.
  */
  
     jed[0] = tjd[0];
     jed[1] = tjd[1];
  
  /*
     Check for target point = center point.
  */
  
     if (target == center)
     {
        for (i = 0; i < 3; i++)
        {
           position[i] = 0.0;
           velocity[i] = 0.0;
        }
        return 0;
     }
  
  /*
     Check for instances of target or center being Earth or Moon,
     and for target or center being the Earth-Moon barycenter.
  */
  
     if ((target == earth) || (center == earth))
        do_moon = 1;
     if ((target == moon) || (center == moon))
        do_earth = 1;
     if ((target == 12) || (center == 12))
        do_earth = 1;
  
     if (do_earth)
     {
        error = state (jed,2, pos_earth,vel_earth);
        if (error)
           return error;
     }
  
     if (do_moon)
     {
        error = state (jed,9, pos_moon,vel_moon);
        if (error)
           return error;
     }
  
  /*
     Make call to State for target object.
  */
  
     if (target == 11)
     {
        for (i = 0; i < 3; i++)
        {
           target_pos[i] = 0.0;
           target_vel[i] = 0.0;
        }
     }
      else if (target == 12)
     {
        for (i = 0; i < 3; i++)
        {
           target_pos[i] = pos_earth[i];
           target_vel[i] = vel_earth[i];
        }
     }
      else
        error = state (jed,target, target_pos,target_vel);
  
     if (error)
        return error;
  
  /*
     Make call to State for center object.
  */
  
  /*
     If the requested center is the Solar System barycenter,
     then don't bother with a second call to State.
  */
  
     if (center == 11)
     {
        for (i = 0; i < 3; i++)
        {
           center_pos[i] = 0.0;
           center_vel[i] = 0.0;
        }
     }
  
  /*
     Center is Earth-Moon barycenter, which was already computed above.
  */
  
      else if (center == 12)
     {
        for (i = 0; i < 3; i++)
        {
           center_pos[i] = pos_earth[i];
           center_vel[i] = vel_earth[i];
        }
     }
      else
        error = state (jed,center, center_pos,center_vel);
  
     if (error)
        return error;
  
  /*
     Check for cases of Earth as target and Moon as center or vice versa.
  */
  
     if ((target == earth) && (center == moon))
     {
        for (i = 0; i < 3; i++)
        {
           position[i] = -center_pos[i];
           velocity[i] = -center_vel[i];
        }
        return 0;
     }
      else if ((target == moon) && (center == earth))
     {
        for (i = 0; i < 3; i++)
        {
           position[i] = target_pos[i];
           velocity[i] = target_vel[i];
        }
        return 0;
     }
  
  /*
     Check for Earth as target, or as center.
  */
  
      else if (target == earth)
     {
        for (i = 0; i < 3; i++)
        {
           target_pos[i] = target_pos[i] - (pos_moon[i] /
              (1.0 + EM_RATIO));
           target_vel[i] = target_vel[i] - (vel_moon[i] /
              (1.0 + EM_RATIO));
        }
     }
      else if (center == earth)
     {
        for (i = 0; i < 3; i++)
        {
           center_pos[i] = center_pos[i] - (pos_moon[i] /
              (1.0 + EM_RATIO));
           center_vel[i] = center_vel[i] - (vel_moon[i] /
              (1.0 + EM_RATIO));
        }
     }
  
  /*
     Check for Moon as target, or as center.
  */
  
      else if (target == moon)
     {
        for (i = 0; i < 3; i++)
        {
           target_pos[i] = (pos_earth[i] - (target_pos[i] /
              (1.0 + EM_RATIO))) + target_pos[i];
           target_vel[i] = (vel_earth[i] - (target_vel[i] /
              (1.0 + EM_RATIO))) + target_vel[i];
        }
     }
      else if (center == moon)
     {
        for (i = 0; i < 3; i++)
        {
           center_pos[i] = (pos_earth[i] - (center_pos[i] /
              (1.0 + EM_RATIO))) + center_pos[i];
           center_vel[i] = (vel_earth[i] - (center_vel[i] /
              (1.0 + EM_RATIO))) + center_vel[i];
        }
     }
  
  /*
     Compute position and velocity vectors.
  */
  
     for (i = 0; i < 3; i++)
     {
        position[i] = target_pos[i] - center_pos[i];
        velocity[i] = target_vel[i] - center_vel[i];
     }
  
     return 0;
  }
  
  
  /********state */
  
  short int state (double *jed, short int target,
  
                   double *target_pos, double *target_vel)
  /*
  ------------------------------------------------------------------------
  
     PURPOSE:
        This function reads and interpolates the JPL planetary
        ephemeris file.
  
     REFERENCES:
        Standish, E.M. and Newhall, X X (1988). "The JPL Export
           Planetary Ephemeris"; JPL document dated 17 June 1988.
  
     INPUT
     ARGUMENTS:
        *jed (double)
           2-element Julian date (TDB) at which interpolation is wanted.
           Any combination of jed[0]+jed[1] which falls within the time
           span on the file is a permissible epoch.  See Note 1 below.
        target (short int)
           The requested body to get data for from the ephemeris file.
           The designation of the astronomical bodies is:
                   0 = Mercury                    6 = Uranus
                   1 = Venus                      7 = Neptune
                   2 = Earth-Moon barycenter      8 = Pluto
                   3 = Mars                       9 = geocentric Moon
                   4 = Jupiter                   10 = Sun
                   5 = Saturn
  
     OUTPUT
     ARGUMENTS:
        *target_pos (double)
           The barycentric position vector array of the requested object,
           in AU.
           (If target object is the Moon, then the vector is geocentric.)
        *target_vel (double)
           The barycentric velocity vector array of the requested object,
           in AU/Day.
  
           Both vectors are referenced to the Earth mean equator and
           equinox of epoch.
  
     RETURNED
     VALUE:
        (short int)
           0...everything OK.
           1...error reading ephemeris file.
           2...epoch out of range.
  
     GLOBALS
     USED:
        KM                eph_manager.h
        EPHFILE           eph_manager.h
        IPT               eph_manager.h
        BUFFER            eph_manager.h
        NRL               eph_manager.h
        RECORD_LENGTH     eph_manager.h
        SS                eph_manager.h
        JPLAU             eph_manager.h
  
     FUNCTIONS
     CALLED:
        split             eph_manager.h
        fseek             stdio.h
        fread             stdio.h
        interpolate       eph_manager.h
        ephem_close       eph_manager.h
  
     VER./DATE/
     PROGRAMMER:
        V1.0/03-93/WTH (USNO/AA): Convert FORTRAN to C.
        V1.1/07-93/WTH (USNO/AA): Update to C standards.
        V2.0/07-98/WTH (USNO/AA): Modify to make position and velocity
                                  two distinct vector arrays.  Routine set
                                  to compute one state per call.
        V2.1/11-07/WKP (USNO/AA): Updated prolog.
        V2.2/10-10/WKP (USNO/AA): Renamed function to lowercase to
                                  comply with coding standards.
  
     NOTES:
        1. For ease in programming, the user may put the entire epoch in
           jed[0] and set jed[1] = 0. For maximum interpolation accuracy,
           set jed[0] = the most recent midnight at or before
           interpolation epoch, and set jed[1] = fractional part of a day
           elapsed between jed[0] and epoch. As an alternative, it may
           prove convenient to set jed[0] = some fixed epoch, such as
           start of the integration and jed[1] = elapsed interval between
           then and epoch.
  
  ------------------------------------------------------------------------
  */
  {
     short int i;
  
     long int nr, rec;
  
     double t[2], aufac = 1.0, jd[4], s;
  
  /*
     Set units based on value of the 'KM' flag.
  */
  
     if (KM)
        t[1] = SS[2] * 86400.0;
      else
     {
        t[1] = SS[2];
        aufac = 1.0 / JPLAU;
     }
  
  /*
     Check epoch.
  */
  
     s = jed[0] - 0.5;
     split (s, &jd[0]);
     split (jed[1], &jd[2]);
     jd[0] += jd[2] + 0.5;
     jd[1] += jd[3];
     split (jd[1], &jd[2]);
     jd[0] += jd[2];
  
  /*
     Return error code if date is out of range.
  */
  
     if ((jd[0] < SS[0]) || ((jd[0] + jd[3]) > SS[1]))
        return 2;
  
  /*
     Calculate record number and relative time interval.
  */
  
     nr = (long int) ((jd[0] - SS[0]) / SS[2]) + 3;
     if (jd[0] == SS[1])
        nr -= 2;
     t[0] = ((jd[0] - ((double) (nr-3) * SS[2] + SS[0])) + jd[3]) / SS[2];
  
  /*
     Read correct record if it is not already in memory.
  */
  
     if (nr != NRL)
     {
        NRL = nr;
        rec = (nr - 1) * RECORD_LENGTH;
        fseek (EPHFILE, rec, SEEK_SET);
        if (!fread (BUFFER, RECORD_LENGTH, 1, EPHFILE))
        {
           ephem_close ();
           return 1;
        }
     }
  
  /*
     Check and interpolate for requested body.
  */
  
     interpolate (&BUFFER[IPT[0][target]-1],t,IPT[1][target],
        IPT[2][target], target_pos,target_vel);
  
     for (i = 0; i < 3; i++)
     {
        target_pos[i] *= aufac;
        target_vel[i] *= aufac;
     }
  
     return 0;
  }
  
  /********interpolate */
  
  void interpolate (double *buf, double *t, long int ncf, long int na,
  
                    double *position, double *velocity)
  /*
  ------------------------------------------------------------------------
  
     PURPOSE:
        This function differentiates and interpolates a set of
        Chebyshev coefficients to give position and velocity.
  
     REFERENCES:
        Standish, E.M. and Newhall, X X (1988). "The JPL Export
           Planetary Ephemeris"; JPL document dated 17 June 1988.
  
     INPUT
     ARGUMENTS:
        *buf (double)
           Array of Chebyshev coefficients of position.
        *t (double)
           t[0] is fractional time interval covered by coefficients at
           which interpolation is desired (0 <= t[0] <= 1).
           t[1] is length of whole interval in input time units.
        ncf (long int)
           Number of coefficients per component.
        na (long int)
           Number of sets of coefficients in full array
           (i.e., number of sub-intervals in full interval).
  
     OUTPUT
     ARGUMENTS:
        *position (double)
           Position array of requested object.
        *velocity (double)
           Velocity array of requested object.
  
     RETURNED
     VALUE:
        None.
  
     GLOBALS
     USED:
        NP                eph_manager.h
        NV                eph_manager.h
        PC                eph_manager.h
        VC                eph_manager.h
        TWOT              eph_manager.h
  
     FUNCTIONS
     CALLED:
        fmod              math.h
  
     VER./DATE/
     PROGRAMMER:
        V1.0/03-93/WTH (USNO/AA): Convert FORTRAN to C.
        V1.1/07-93/WTH (USNO/AA): Update to C standards.
        V1.2/07-98/WTH (USNO/AA): Modify to make position and velocity
                                  two distinct vector arrays.
        V1.3/11-07/WKP (USNO/AA): Updated prolog.
        V1.4/12-07/WKP (USNO/AA): Changed ncf and na arguments from short
                                  int to long int.
        V1.5/10-10/WKP (USNO/AA): Renamed function to lowercase to
                                  comply with coding standards.
  
     NOTES:
        None.
  
  ------------------------------------------------------------------------
  */
  {
     long int i, j, k, l;
  
     double dna, dt1, temp, tc, vfac;
  
  /*
     Get correct sub-interval number for this set of coefficients and
     then get normalized Chebyshev time within that subinterval.
  */
  
     dna = (double) na;
     dt1 = (double) ((long int) t[0]);
     temp = dna * t[0];
     l = (long int) (temp - dt1);
  
  /*
     'tc' is the normalized Chebyshev time (-1 <= tc <= 1).
  */
  
     tc = 2.0 * (fmod (temp, 1.0) + dt1) - 1.0;
  
  /*
     Check to see whether Chebyshev time has changed, and compute new
     polynomial values if it has.  (The element PC[1] is the value of
     t1[tc] and hence contains the value of 'tc' on the previous call.)
  */
  
     if (tc != PC[1])
     {
        NP = 2;
        NV = 3;
        PC[1] = tc;
        TWOT = tc + tc;
     }
  
  /*
     Be sure that at least 'ncf' polynomials have been evaluated and
     are stored in the array 'PC'.
  */
  
     if (NP < ncf)
     {
        for (i = NP; i < ncf; i++)
           PC[i] = TWOT * PC[i-1] - PC[i-2];
        NP = ncf;
     }
  
  /*
     Interpolate to get position for each component.
  */
  
     for (i = 0; i < 3; i++)
     {
        position[i] = 0.0;
        for (j = ncf-1; j >= 0; j--)
        {
           k = j + (i * ncf) + (l * (3 * ncf));
           position[i] += PC[j] * buf[k];
        }
     }
  
  /*
     If velocity interpolation is desired, be sure enough derivative
     polynomials have been generated and stored.
  */
  
     vfac = (2.0 * dna) / t[1];
     VC[2] = 2.0 * TWOT;
     if (NV < ncf)
     {
        for (i = NV; i < ncf; i++)
           VC[i] = TWOT * VC[i-1] + PC[i-1] + PC[i-1] - VC[i-2];
        NV = ncf;
     }
  
  /*
     Interpolate to get velocity for each component.
  */
  
     for (i = 0; i < 3; i++)
     {
        velocity[i] = 0.0;
        for (j = ncf-1; j > 0; j--)
        {
           k = j + (i * ncf) + (l * (3 * ncf));
           velocity[i] += VC[j] * buf[k];
        }
        velocity[i] *= vfac;
     }
  
     return;
  }
  
  /********split */
  
  void split (double tt,
  
              double *fr)
  /*
  ------------------------------------------------------------------------
  
     PURPOSE:
        This function breaks up a double number into a double integer
        part and a fractional part.
  
     REFERENCES:
        Standish, E.M. and Newhall, X X (1988). "The JPL Export
           Planetary Ephemeris"; JPL document dated 17 June 1988.
  
     INPUT
     ARGUMENTS:
        tt (double)
           Input number.
  
     OUTPUT
     ARGUMENTS:
        *fr (double)
           2-element output array;
              fr[0] contains integer part,
              fr[1] contains fractional part.
           For negative input numbers,
              fr[0] contains the next more negative integer;
              fr[1] contains a positive fraction.
  
     RETURNED
     VALUE:
        None.
  
     GLOBALS
     USED:
        None.
  
     FUNCTIONS
     CALLED:
        None.
  
     VER./DATE/
     PROGRAMMER:
        V1.0/06-90/JAB (USNO/NA): CA coding standards
        V1.1/03-93/WTH (USNO/AA): Convert to C.
        V1.2/07-93/WTH (USNO/AA): Update to C standards.
        V1.3/10-10/WKP (USNO/AA): Renamed function to lowercase to
                                  comply with coding standards.
  
     NOTES:
        None.
  
  ------------------------------------------------------------------------
  */
  {
  
  /*
     Get integer and fractional parts.
  */
  
     fr[0] = (double)((long int) tt);
     fr[1] = tt - fr[0];
  
  /*
     Make adjustments for negative input number.
  */
  
     if ((tt >= 0.0) || (fr[1] == 0.0))
        return;
      else
     {
        fr[0] = fr[0] - 1.0;
        fr[1] = fr[1] + 1.0;
     }
  
     return;
  }
  