/*****************************************************************************
*        *****    jpl planetary and lunar ephemerides    *****     C ver.1.2 *
******************************************************************************
*                                                                            *
*  This program was written in standard fortran-77 and it was manually       *
*  translated to C language by Piotr A. Dybczynski (dybol@phys.amu.edu.pl),  *
*  subsequently revised heavily by Bill J Gray (pluto@gwi.net).              *
*                                                                            *
******************************************************************************
*                 Last modified: July 23, 1997 by PAD                        *
******************************************************************************
16 Mar 2001:  Revised by Bill J. Gray.  You can now use binary
ephemerides with either byte order ('big-endian' or 'small-endian');
the code checks to see if the data is in the "wrong" order for the
current platform,  and swaps bytes on-the-fly if needed.  (Yes,  this
can result in a slowdown... sometimes as much as 1%.  The function is
so mathematically intensive that the byte-swapping is the least of our
troubles.)  You can also use DE-200, 403, 404, 405,  or 406 without
recompiling (the constan( ) function now determines which ephemeris is
in use and its byte order).

Also,  I did some minor optimization of the interp( ) (Chebyshev
interpolation) function,  resulting in a bit of a speedup.

The code has been modified to be a separately linkable component,  with
details of the implementation encapsulated.
*****************************************************************************/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>

/**** include variable and type definitions, specific for this C version */

#include "jpleph.h"
#include "jpl_int.h"

#define TRUE 1
#define FALSE 0

double DLL_FUNC jpl_get_double( const void *ephem, const int value)
{
   return( *(double *)( (char *)ephem + value));
}

double DLL_FUNC jpl_get_long( const void *ephem, const int value)
{
   return( *(long *)( (char *)ephem + value));
}


/*****************************************************************************
**           jpl_pleph( ephem,et,ntar,ncent,rrd,calc_velocity)              **
******************************************************************************
**                                                                          **
**    This subroutine reads the jpl planetary ephemeris                     **
**    and gives the position and velocity of the point 'ntarg'              **
**    with respect to 'ncent'.                                              **
**                                                                          **
**    Calling sequence parameters:                                          **
**                                                                          **
**      et = (double) julian ephemeris date at which interpolation          **
**           is wanted.                                                     **
**                                                                          **
**    ntarg = integer number of 'target' point.                             **
**                                                                          **
**    ncent = integer number of center point.                               **
**                                                                          **
**    The numbering convention for 'ntarg' and 'ncent' is:                  **
**                                                                          **
**            1 = mercury           8 = neptune                             **
**            2 = venus             9 = pluto                               **
**            3 = earth            10 = moon                                **
**            4 = mars             11 = sun                                 **
**            5 = jupiter          12 = solar-system barycenter             **
**            6 = saturn           13 = earth-moon barycenter               **
**            7 = uranus           14 = nutations (longitude and obliq)     **
**                                 15 = librations, if on eph. file         **
**                                                                          **
**            (If nutations are wanted, set ntarg = 14.                     **
**             For librations, set ntarg = 15. set ncent= 0)                **
**                                                                          **
**     rrd = output 6-element, double array of position and velocity        **
**           of point 'ntarg' relative to 'ncent'. The units are au and     **
**           au/day. For librations the units are radians and radians       **
**           per day. In the case of nutations the first four words of      **
**           rrd will be set to nutations and rates, having units of        **
**           radians and radians/day.                                       **
**                                                                          **
**           The option is available to have the units in km and km/sec.    **
**           for this, set km=TRUE at the begining of the program.          **
**                                                                          **
**     calc_velocity = integer flag;  if nonzero,  velocities will be       **
**           computed,  otherwise not.                                      **
**                                                                          **
*****************************************************************************/
int DLL_FUNC jpl_pleph( void *ephem, const double et, const int ntarg,
                      const int ncent, double rrd[], const int calc_velocity)
{
  struct jpl_eph_data *eph = (struct jpl_eph_data *)ephem;
  double pv[13][6];/* pv is the position/velocity array
                             NUMBERED FROM ZERO: 0=Mercury,1=Venus,...
                             8=Pluto,9=Moon,10=Sun,11=SSBary,12=EMBary
                             First 10 elements (0-9) are affected by
                             jpl_state(), all are adjusted here.         */


  int rval = 0, list_val = (calc_velocity ? 2 : 1);
  int i, k, list[12];    /* list is a vector denoting, for which "body"
                            ephemeris values should be calculated by
                            jpl_state():  0=Mercury,1=Venus,2=EMBary,...,
                            8=Pluto,  9=geocentric Moon, 10=nutations in
                            long. & obliq.  11= lunar librations  */

   for( i = 0; i < 6; ++i) rrd[i] = 0.0;

   if( ntarg == ncent) return( 0);

   for( i = 0; i < 12; i++) list[i] = 0;

/*   check for nutation call    */

   if( ntarg == 14)
      {
      if( eph->ipt[11][1] > 0) /* there is nutation on ephemeris */
         {
         list[10] = list_val;
         if( jpl_state( ephem, et, list, pv, rrd, 0))
            rval = -1;
         }
      else          /*  no nutations on the ephemeris file  */
         rval = -2;
      return( rval);
      }

/*  check for librations   */

   if( ntarg == 15)
      {
      if( eph->ipt[12][1] > 0) /* there are librations on ephemeris file */
         {
         list[11] = list_val;
         if( jpl_state( eph, et, list, pv, rrd, 0))
            rval = -3;
         for( i = 0; i < 6; ++i)
            rrd[i] = pv[10][i]; /* librations */
         }
      else          /*  no librations on the ephemeris file  */
         rval = -4;
      return( rval);
      }

/*  force barycentric output by 'state'     */

/*  set up proper entries in 'list' array for state call     */

   for( i = 0; i < 2; i++) /* list[] IS NUMBERED FROM ZERO ! */
      {
      k = ntarg-1;
      if( i == 1) k=ncent-1;            /* same for ntarg & ncent */
      if( k <= 9) list[k] = list_val;   /* Major planets */
      if( k == 9) list[2] = list_val;   /* for moon,  earth state is needed */
      if( k == 2) list[9] = list_val;   /* for earth,  moon state is needed */
      if( k == 12) list[2] = list_val;  /* EMBary state additionaly */
      }

/*   make call to state   */

   if( jpl_state( eph, et, list, pv, rrd, 1))
      rval = -5;
   /* Solar System barycentric Sun state goes to pv[10][] */
   if( ntarg == 11 || ncent == 11)
      for( i = 0; i < 6; i++)
         pv[10][i] = eph->pvsun[i];

   /* Solar System Barycenter coordinates & velocities equal to zero */
   if( ntarg == 12 || ncent == 12)
      for( i = 0; i < 6; i++)
         pv[11][i] = 0.0;

   /* Solar System barycentric EMBary state:  */
   if( ntarg == 13 || ncent == 13)
      for( i = 0; i < 6; i++)
         pv[12][i] = pv[2][i];

   /* if moon from earth or earth from moon ..... */
   if( (ntarg*ncent) == 30 && (ntarg+ncent) == 13)
      for( i = 0; i < 6; ++i) pv[2][i]=0.0;
   else
      {
      if( list[2])           /* calculate earth state from EMBary */
         for( i = 0; i < list[2] * 3; ++i)
            pv[2][i] -= pv[9][i]/(1.0+eph->emrat);

      if(list[9]) /* calculate Solar System barycentric moon state */
         for( i = 0; i < list[9] * 3; ++i)
            pv[9][i] += pv[2][i];
      }

   for( i = 0; i < list_val * 3; ++i)
      rrd[i] = pv[ntarg-1][i] - pv[ncent-1][i];

   return( rval);
}

/*****************************************************************************
**                     interp(buf,t,ncf,ncm,na,ifl,pv)                      **
******************************************************************************
**                                                                          **
**    this subroutine differentiates and interpolates a                     **
**    set of chebyshev coefficients to give position and velocity           **
**                                                                          **
**    calling sequence parameters:                                          **
**                                                                          **
**      input:                                                              **
**                                                                          **
**      iinfo   stores certain chunks of interpolation info,  in hopes      **
**              that if you call again with similar parameters,  the        **
**              function won't have to re-compute all coefficients/data.    **
**                                                                          **
**       coef   1st location of array of d.p. chebyshev coefficients        **
**              of position                                                 **
**                                                                          **
**          t   t[0] is double fractional time in interval covered by       **
**              coefficients at which interpolation is wanted               **
**              (0 <= t[0] <= 1).  t[1] is dp length of whole               **
**              interval in input time units.                               **
**                                                                          **
**        ncf   # of coefficients per component                             **
**                                                                          **
**        ncm   # of components per set of coefficients                     **
**                                                                          **
**         na   # of sets of coefficients in full array                     **
**              (i.e., # of sub-intervals in full interval)                 **
**                                                                          **
**         ifl  integer flag: =1 for positions only                         **
**                            =2 for pos and vel                            **
**                                                                          **
**                                                                          **
**      output:                                                             **
**                                                                          **
**    posvel   interpolated quantities requested.  dimension                **
**              expected is posvel[ncm*ifl], double precision.              **
**                                                                          **
*****************************************************************************/
static void interp( struct interpolation_info *iinfo,
                 const double coef[], const double t[2], const int ncf,
                 const int ncm, const int na, const int ifl, double posvel[])
{
  double dna, dt1, temp, tc, vfac, temp1;
  int l, i, j;

/*  entry point. get correct sub-interval number for this set
    of coefficients and then get normalized chebyshev time
    within that subinterval.                                             */

  dna = (double)na;
  modf( t[0], &dt1);
  temp = dna * t[0];
  l = (int)(temp - dt1);

/*  tc is the normalized chebyshev time (-1 <= tc <= 1)    */

  tc = 2.0 * (modf( temp, &temp1) + dt1) - 1.0;

/*  check to see whether chebyshev time has changed,
    and compute new polynomial values if it has.
    (the element iinfo->pc[1] is the value of t1[tc] and hence
    contains the value of tc on the previous call.)     */

  if(tc != iinfo->pc[1])
    {
      iinfo->np = 2;
      iinfo->nv = 3;
      iinfo->pc[1] = tc;
      iinfo->twot = tc+tc;
    }

/*  be sure that at least 'ncf' polynomials have been evaluated
    and are stored in the array 'iinfo->pc'.    */

  if(iinfo->np < ncf)
    {
    double *pc_ptr = iinfo->pc + iinfo->np;

    for(i=ncf - iinfo->np; i; i--, pc_ptr++)
       *pc_ptr = iinfo->twot * pc_ptr[-1] - pc_ptr[-2];
    iinfo->np=ncf;
    }

/*  interpolate to get position for each component  */

  for( i = 0; i < ncm; ++i)        /* ncm is a number of coordinates */
     {
     const double *coeff_ptr = coef + ncf * (i + l * ncm + 1);
     const double *pc_ptr = iinfo->pc + ncf;

     posvel[i]=0.0;
     for( j = ncf; j; j--)
        posvel[i] += (*--pc_ptr) * (*--coeff_ptr);
     }

   if(ifl <= 1) return;

/*  if velocity interpolation is wanted, be sure enough
    derivative polynomials have been generated and stored.    */

  vfac=(dna+dna)/t[1];
  iinfo->vc[2] = iinfo->twot + iinfo->twot;
  if( iinfo->nv < ncf)
     {
     double *vc_ptr = iinfo->vc + iinfo->nv;
     const double *pc_ptr = iinfo->pc + iinfo->nv - 1;

     for( i = ncf - iinfo->nv; i; i--, vc_ptr++, pc_ptr++)
        *vc_ptr = iinfo->twot * vc_ptr[-1] + *pc_ptr + *pc_ptr - vc_ptr[-2];
     iinfo->nv = ncf;
     }

/*  interpolate to get velocity for each component    */

   for( i = 0; i < ncm; ++i)
      {
      double tval = 0.;
      const double *coeff_ptr = coef + ncf * (i + l * ncm + 1);
      const double *vc_ptr = iinfo->vc + ncf;

      for( j = ncf; j; j--)
         tval += (*--vc_ptr) * (*--coeff_ptr);
      posvel[ i + ncm] = tval * vfac;
      }
   return;
}

/* swap_long_integer() and swap_double() are used when reading a binary
ephemeris that was created on a machine with 'opposite' byte order to
the currently-used machine (signalled by the 'swap_bytes' flag in the
jpl_eph_data structure).  In such cases,  every double and integer
value read from the ephemeris must be byte-swapped by these two functions. */

#define SWAP_MACRO( A, B, TEMP)   { TEMP = A;  A = B;  B = TEMP; }

static void swap_long_integer( void *ptr)
{
   char *tptr = (char *)ptr, tchar;

   SWAP_MACRO( tptr[0], tptr[3], tchar);
   SWAP_MACRO( tptr[1], tptr[2], tchar);
}

static void swap_double( void *ptr, long count)
{
   char *tptr = (char *)ptr, tchar;

   while( count--)
      {
      SWAP_MACRO( tptr[0], tptr[7], tchar);
      SWAP_MACRO( tptr[1], tptr[6], tchar);
      SWAP_MACRO( tptr[2], tptr[5], tchar);
      SWAP_MACRO( tptr[3], tptr[4], tchar);
      tptr += 8;
      }
}

/*****************************************************************************
**                        jpl_state(ephem,et2,list,pv,nut,bary)             **
******************************************************************************
** This subroutine reads and interpolates the jpl planetary ephemeris file  **
**                                                                          **
**    Calling sequence parameters:                                          **
**                                                                          **
**    Input:                                                                **
**                                                                          **
**        et2[] double, 2-element JED epoch at which interpolation          **
**              is wanted.  Any combination of et2[0]+et2[1] which falls    **
**              within the time span on the file is a permissible epoch.    **
**                                                                          **
**               a. for ease in programming, the user may put the           **
**                  entire epoch in et2[0] and set et2[1]=0.0               **
**                                                                          **
**               b. for maximum interpolation accuracy, set et2[0] =        **
**                  the most recent midnight at or before interpolation     **
**                  epoch and set et2[1] = fractional part of a day         **
**                  elapsed between et2[0] and epoch.                       **
**                                                                          **
**               c. as an alternative, it may prove convenient to set       **
**                  et2[0] = some fixed epoch, such as start of integration,**
**                  and et2[1] = elapsed interval between then and epoch.   **
**                                                                          **
**       list   12-element integer array specifying what interpolation      **
**              is wanted for each of the "bodies" on the file.             **
**                                                                          **
**                        list[i]=0, no interpolation for body i            **
**                               =1, position only                          **
**                               =2, position and velocity                  **
**                                                                          **
**              the designation of the astronomical bodies by i is:         **
**                                                                          **
**                        i = 0: mercury                                    **
**                          = 1: venus                                      **
**                          = 2: earth-moon barycenter                      **
**                          = 3: mars                                       **
**                          = 4: jupiter                                    **
**                          = 5: saturn                                     **
**                          = 6: uranus                                     **
**                          = 7: neptune                                    **
**                          = 8: pluto                                      **
**                          = 9: geocentric moon                            **
**                          =10: nutations in longitude and obliquity       **
**                          =11: lunar librations (if on file)              **
**                                                                          **
**    output:                                                               **
**                                                                          **
**    pv[][6]   double array that will contain requested interpolated       **
**              quantities.  The body specified by list[i] will have its    **
**              state in the array starting at pv[i][0]  (on any given      **
**              call, only those words in 'pv' which are affected by the    **
**              first 10 'list' entries (and by list(11) if librations are  **
**              on the file) are set.  The rest of the 'pv' array           **
**              is untouched.)  The order of components in pv[][] is:       **
**              pv[][0]=x,....pv[][5]=dz.                                   **
**                                                                          **
**              All output vectors are referenced to the earth mean         **
**              equator and equinox of epoch. The moon state is always      **
**              geocentric; the other nine states are either heliocentric   **
**              or solar-system barycentric, depending on the setting of    **
**              global variables (see below).                               **
**                                                                          **
**              Lunar librations, if on file, are put into pv[10][k] if     **
**              list[11] is 1 or 2.                                         **
**                                                                          **
**        nut   dp 4-word array that will contain nutations and rates,      **
**              depending on the setting of list[10].  the order of         **
**              quantities in nut is:                                       **
**                                                                          **
**                       d psi  (nutation in longitude)                     **
**                       d epsilon (nutation in obliquity)                  **
**                       d psi dot                                          **
**                       d epsilon dot                                      **
**                                                                          **
*****************************************************************************/
int DLL_FUNC jpl_state( void *ephem, const double et, const int list[12],
                          double pv[][6], double nut[4], const int bary)
{
  struct jpl_eph_data *eph = (struct jpl_eph_data *)ephem;
  int i,j, n_intervals;
  long int nr;
  double prev_midnight, time_of_day;
  double *buf = eph->cache;
  double s,t[2],aufac;
  struct interpolation_info *iinfo = (struct interpolation_info *)eph->iinfo;


/*  ********** main entry point **********  */

  s=et - 0.5;
  prev_midnight = floor( s);
  time_of_day = s - prev_midnight;
  prev_midnight += 0.5;

/* here prev_midnight contains last midnight before epoch desired (in JED: *.5)
   and time_of_day contains the remaining, fractional part of the epoch    */

/*   error return for epoch out of range  */

  if( et < eph->ephem_start || et > eph->ephem_end)
     return( -1);

/*   calculate record # and relative time in interval   */

   nr = (long)((prev_midnight-eph->ephem_start)/eph->ephem_step)+2;
        /* add 2 to adjus for the first two records containing header data */
   if( prev_midnight == eph->ephem_end)
      nr--;
   t[0]=( prev_midnight-( (1.0*nr-2.0)*eph->ephem_step+eph->ephem_start) +
           time_of_day )/eph->ephem_step;

/*   read correct record if not in core (static vector buf[])   */

   if( nr != eph->curr_cache_loc)
      {
      eph->curr_cache_loc = nr;
      fseek( eph->ifile, nr * eph->recsize, SEEK_SET);
      fread( buf, (size_t)eph->ncoeff, sizeof( double), eph->ifile);
      if( eph->swap_bytes)
         swap_double( buf, eph->ncoeff);
      }
   t[1] = eph->ephem_step;
   aufac = 1.0 / eph->au;

   for( n_intervals = 1; n_intervals <= 8; n_intervals *= 2)
      for( i = 0; i < 11; i++)
         if( n_intervals == eph->ipt[i][2] && (list[i] || i == 10))
            {
            int flag = ((i == 10) ? 2 : list[i]);
            double *dest = ((i == 10) ? eph->pvsun : pv[i]);

            interp( iinfo, &buf[eph->ipt[i][0]-1], t, (int)eph->ipt[i][1], 3,
                                    n_intervals, flag, dest);
                               /* gotta convert units */
            for( j = 0; j < flag * 3; j++)
               dest[j] *= aufac;
            }

    if( !bary)                             /* gotta correct everybody for */
       for( i = 0; i < 9; i++)            /* the solar system barycenter */
          for( j = 0; j < list[i] * 3; j++)
              pv[i][j] -= eph->pvsun[j];

/*  do nutations if requested (and if on file)    */

      if(list[10] > 0 && eph->ipt[11][1] > 0)
         interp( iinfo, &buf[eph->ipt[11][0]-1], t, (int)eph->ipt[11][1], 2,
                              (int)eph->ipt[11][2], list[10], nut);

/*  get librations if requested (and if on file)    */

      if(list[11] > 0 && eph->ipt[12][1] > 0)
         {
         double pefau[6];

         interp( iinfo, &buf[eph->ipt[12][0]-1], t, (int)eph->ipt[12][1], 3,
                              (int)eph->ipt[12][2], list[11], pefau);
         for( j = 0; j < 6; ++j) pv[10][j]=pefau[j];
         }
   return( 0);
}

/****************************************************************************
**    jpl_init_ephemeris( ephemeris_filename, nam, val, n_constants)       **
*****************************************************************************
**                                                                         **
**    this function does the initial prep work for use of binary JPL       **
**    ephemerides.                                                         **
**      const char *ephemeris_filename = full path/filename of the binary  **
**          ephemeris (on the Willmann-Bell CDs,  this is UNIX.200, 405,   **
**          or 406)
**      char nam[][6] = array of constant names (max 6 characters each)    **
**          You can pass nam=NULL if you don't care about the names        **
**      double *val = array of values of constants                         **
**          You can pass val=NULL if you don't care about the constants    **
**      Return value is a pointer to the jpl_eph_data structure            **
**      NULL is returned if the file isn't opened or memory isn't alloced  **
****************************************************************************/
void * DLL_FUNC jpl_init_ephemeris( const char *ephemeris_filename,
                          char nam[][6], double *val)
{
   int i, j;
   long de_version;
   char title[84];
   FILE *ifile = fopen( ephemeris_filename, "rb");
   struct jpl_eph_data *rval;
   struct interpolation_info *iinfo;

   if( !ifile)
      return( NULL);
               /* Rather than do three separate allocations,  everything   */
               /* we need is allocated in _one_ chunk,  then parceled out. */
               /* This looks a little weird,  but it does simplify error   */
               /* handling and cleanup.                                    */
   rval = (struct jpl_eph_data *)calloc( sizeof( struct jpl_eph_data)
                        + sizeof( struct interpolation_info)
                        + MAX_KERNEL_SIZE * 4, 1);
   if( !rval)
      {
      fclose( ifile);
      return( NULL);
      }
   rval->ifile = ifile;
   fread( title, 84, 1, ifile);
   fseek( ifile, 2652L, SEEK_SET);      /* skip title & constant name data */
   fread( rval, JPL_HEADER_SIZE, 1, ifile);

   de_version = atoi( title + 26);
         /* Once upon a time,  the kernel size was determined from the */
         /* DE version.  This was not a terrible idea,  except that it */
         /* meant that when the code faced a new version,  it broke.   */
         /* Thus,  the 'OLD_CODE' section was replaced with some logic */
         /* that figures out the kernel size.  The 'OLD_CODE' section  */
         /* should therefore be of only historical interest.           */
#ifdef OLD_CODE
   switch( de_version)
      {
      case 200:
      case 202:
         rval->kernel_size = 1652;
         break;
      case 403:
      case 405:
         rval->kernel_size = 2036;
         break;
      case 404:
      case 406:
         rval->kernel_size = 1456;
         break;
      }
#endif
          /* The 'iinfo' struct is right after the 'jpl_eph_data' struct: */
   iinfo = (struct interpolation_info *)(rval + 1);
   rval->iinfo = (void *)iinfo;
   iinfo->np = 2;
   iinfo->nv = 3;
   iinfo->pc[0] = 1.0;
   iinfo->pc[1] = 0.0;
   iinfo->vc[1] = 1.0;
   rval->curr_cache_loc = -1L;
          /* The 'cache' data is right after the 'iinfo' struct: */
   rval->cache = (double *)( iinfo + 1);
          /* A small piece of trickery:  in the binary file,  data is stored */
          /* for ipt[0...11],  then the ephemeris version,  then the         */
          /* remaining ipt[12] data.  A little switching is required to get  */
          /* the correct order. */
   for( i = 0; i < 3; i++)
      rval->ipt[12][i] = rval->ipt[12][i + 1];
   rval->ephemeris_version = de_version;

   rval->swap_bytes = ( rval->ncon < 0 || rval->ncon > 65536L);
   if( rval->swap_bytes)     /* byte order is wrong for current platform */
      {
      swap_double( &rval->ephem_start, 1);
      swap_double( &rval->ephem_end, 1);
      swap_double( &rval->ephem_step, 1);
      swap_long_integer( &rval->ncon);
      swap_double( &rval->au, 1);
      swap_double( &rval->emrat, 1);
      for( j = 0; j < 3; j++)
         for( i = 0; i < 13; i++)
            swap_long_integer( &rval->ipt[i][j]);
      }
   rval->kernel_size = 4;
   for( i = 0; i < 13; i++)
      rval->kernel_size +=
                       rval->ipt[i][1] * rval->ipt[i][2] * ((i == 11) ? 4 : 6);
// printf( "Kernel size = %d\n", rval->kernel_size);
   rval->recsize = rval->kernel_size * 4L;
   rval->ncoeff = rval->kernel_size / 2L;

   if( val)
      {
      fseek( ifile, rval->recsize, SEEK_SET);
      fread( val, (size_t)rval->ncon, sizeof( double), ifile);
      if( rval->swap_bytes)     /* gotta swap the constants,  too */
         swap_double( val, rval->ncon);
      }

   if( nam)
      {
      fseek( ifile, 84L * 3L, SEEK_SET);   /* just after the 3 'title' lines */
      for( i = 0; i < (int)rval->ncon; i++)
         fread( nam[i], 6, 1, ifile);
      }
   return( rval);
}

/****************************************************************************
**    jpl_close_ephemeris( ephem)                                          **
*****************************************************************************
**                                                                         **
**    this function closes files and frees up memory allocated by the      **
**    jpl_init_ephemeris( ) function.                                      **
****************************************************************************/
void DLL_FUNC jpl_close_ephemeris( void *ephem)
{
   struct jpl_eph_data *eph = (struct jpl_eph_data *)ephem;

   fclose( eph->ifile);
   free( ephem);
}
/*************************** THE END ***************************************/