/**=============================================================**/
/**  function jpleph interfaces to David Hoffman's C-version of **/
/**  processing JPL EPHEMERIS data ephem_read.c                 **/
/**                                                             **/
/**  See also the Fortran procedure PLEPH from NOVAS            **/
/**                                                             **/
/**  Author: Jaap Spies                                         **/
/**=============================================================**/


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

#include "ephem_read.h"
#ifndef TYPES_DEFINED
#include "ephem_types.h"
#endif

#define JPLFILE "JPLEPH"
static char *ephemerisFile = NULL;
                
/*
The IAU recommandation
#define EMRAT (1.0 / 0.0123000345)
but JPL uses:
*/
#define EMRAT	81.3005600
#define AU	149597870.691
#define SEC_PER_DAY	(24.0*60.0*60.0)

#define DEBUG	0


/* Set the ephemeris file to read */
void setEphemeris(const char *file)
{
  static char fileName[100];
  if(DEBUG)
    printf("jpleph.c: the pointer *file is %s\n",file);
  strcpy(fileName, file);
  ephemerisFile = fileName;
  if(DEBUG)
       printf("jpleph.c: setEphemeris has ephemerisFile = %s \n",ephemerisFile);
}

int jpleph ( double jdt, int targ, int cent, double *pos, double *vel)
{
   int i;
   stateType  State[13];
   static int firstTime = 1;

   if(DEBUG)
	printf("jpleph called: %7.3f  targ: %d  cent %d \n", jdt, targ, cent);

   /* Initialize the ephemeris.................................................*/

   if (firstTime)
     {
       firstTime = 0;
       if (ephemerisFile == NULL)
	 Initialize_Ephemeris(JPLFILE);
       else
	 Initialize_Ephemeris(ephemerisFile);
     }

   /* Compute the desired ephemeris data.......................................*/

   if (cent == SS_BARY)
   {
       for (i=0; i<3;i++)
       {
           State[SS_BARY].Position[i] = 0.0;
           State[SS_BARY].Velocity[i] = 0.0;
       }
   }
   else
        Interpolate_State (jdt, SUN, &State[SUN]);

   if (targ == MOON)
       Interpolate_State (jdt, EARTH, &State[EARTH] );
   else if ( targ == EARTH)
       Interpolate_State (jdt, MOON, &State[MOON] );

   Interpolate_State( jdt , targ , &State[targ] );


   if (targ == MOON || targ == EARTH)
   {
       for(i=0; i<3; i++)
       {
           State[EARTH].Position[i] -= State[MOON].Position[i] / (1.0 + EMRAT);
           State[EARTH].Velocity[i] -= State[MOON].Velocity[i] / (1.0 + EMRAT);
       }
   }
   if (targ == MOON)
       for (i=0; i<3; i++)
       {
           State[targ].Position[i] += State[EARTH].Position[i];
           State[targ].Velocity[i] += State[EARTH].Velocity[i];
       }

/*
solsys wants position in AU and velocity in AU/DAY
*/

   for (i=0; i<3; i++)
   {
       pos[i] = State[targ].Position[i] - State[cent].Position[i];
       pos[i] /= AU;
       vel[i] = State[targ].Velocity[i] - State[cent].Velocity[i];
       vel[i] = vel[i] * SEC_PER_DAY / AU;
   }
   if(DEBUG)
       for (i=0; i<3; i++)
       {
           printf(" %7.3f \n",pos[i] );
           printf(" %7.3f \n", vel[i]);
       }
   return 0;
}