/* Read and interpolate de118i ephemeris.  */#include "kep.h"/* Define this or not in the makefile.  Use kepjpl.c when not defined.  */#if SSYSTEM/* Interpolator degree */#define NINTERP 12/* Include constants for appropriate ephemeris */#include "de118i.h"/* long rec0 = 0L; *//* long recsiz = 584L; */double JDi, JDb;extern double au;  /* Kilometers per au */extern double emrat; /* Earth/Moon mass ratio *//* System header files required for I/O.   Define IBMPC, etc. 1 or 0 in kep.h.  */#ifdef _MSC_VER#if _MSC_VER >= 1000#include <stdlib.h>#endif#include <sys\types.h>#include <sys\stat.h>#include <fcntl.h>#include <io.h>#endif#if UNIX#include <sys/types.h>#include <sys/stat.h>#include <fcntl.h>#include <unistd.h>#endif#ifndef O_BINARY#define O_BINARY 0#endifextern double pearthb[];static long bcnt = 0;extern struct orbit earth;extern double coseps, sineps, J2000;int defd = 0;extern int defd;extern char defile[];double zatan2(), asin(), sqrt();/* heliocentric velocity of given object (geocentric, if moon) */double pobjb[3];double pobjh[3];double vobjb[3];double vobjh[3];double psunb[3];double vsunb[3];int kepjpl(J, e, rect, polar)double J, rect[], polar[];struct orbit *e;{double JD, ratio, x, y, z;double pm[3], vm[3];int nobj, i;floop:/* Open up the ephemeris file.  */if( defd <= 0 )	{	if( defd < 0 )		{		printf( "Enter DE file name ? " );		gets( &defile[0] );		}	defd = open( &defile[0], O_BINARY | O_RDONLY, S_IREAD );	if( defd <= 0 )		{		printf( "Can't find DE file <%s>\n", &defile[0] );		defd = -1;		goto floop;		}	printf( "Opened %s\n", &defile[0] );/* rec0 is assumed initialized properly in de245.h or de200.h or de102.h */	bcnt = rec0;	lseek( defd, bcnt, SEEK_SET );	read( defd, &JDb, 8 );/* The number of days covered by one record is the difference between   the dates on two adjacent records.  */	lseek( defd, bcnt+recsiz, SEEK_SET );	read( defd, &JDi, 8 );	JDi -= JDb;/*#if DEBUG*/#if 1	printf( "First date in file = %.8E\n", JDb );#endif	}nobj = objnum;if( nobj == 3 )	nobj = 10;	/* moon */if( nobj == 0 )	nobj = 3;	/* sun */JD = tdb(J);/* Sun */#if DEBUGprintf( "psunb, vsunb\n" );#endifif( jpl( JD, 11, psunb, vsunb ) )	{	printf( "Error, DE file boundary violation\n" );	exit(0);	}if( (e == &earth) || (nobj == 3) || (nobj == 10) )	{/* Moon, geocentric */#if DEBUG	printf( "pgmoon, vgmoon\n" );#endif	jpl( JD, 10, pm, vm );	if( e != &earth )		{		for( i=0; i<3; i++ )			{			vobjb[i] = vm[i];			pobjb[i] = pm[i];			vobjh[i] = vm[i];			pobjh[i] = pm[i];			}		goto output;		}/* Earth-Moon barycenter */#if DEBUG	printf( "pemb, vemb\n" );#endif	jpl( JD, 3, pobjb, vobjb );/* Heliocentric EMB */	for( i=0; i<3; i++ )		{		pobjh[i] = pobjb[i] - psunb[i];		vobjh[i] = vobjb[i] - vsunb[i];		}/* Heliocentric Earth */	ratio = 1.0/(emrat+1.0);	for( i=0; i<3; i++ )		{		pobjh[i] = pobjh[i] - ratio * pm[i];		vobjh[i] = vobjh[i] - ratio * vm[i];		}/* Earth position and velocity re solar system barycenter */#if DEBUG	printf( "barycentric earth = pemb - pgmoon/(1+emratio)\n" );#endif	for( i=0; i<3; i++ )		{		pearthb[i] = pobjh[i] + psunb[i];		vearth[i] = vobjh[i] + vsunb[i];		}	jvearth = JD;	}else	{/* Selected object, not Earth or Moon */#if DEBUG	printf( "pobjb, vobjb\n" );#endif	jpl( JD, nobj, pobjb, vobjb );	for( i=0; i<3; i++ )		{		pobjh[i] = pobjb[i] - psunb[i];		vobjh[i] = vobjb[i] - vsunb[i];		}#if DEBUG	printf( "pobjb-psunb\n" );	prvec(pobjh);	prvec(vobjh); /* heliocentric velocity */#endif	}output:/* Write out the heliocentric object coordinates */for( i=0; i<3; i++ )	rect[i] = pobjh[i];epsiln(J2000);x = pobjh[0];y = pobjh[1];z = pobjh[2];/* radius distance */ratio = x*x + y*y + z*z;polar[2] = sqrt(ratio);/* Convert from equatorial to ecliptic coordinates */ratio  =  coseps * y  +  sineps * z;z  = -sineps * y  +  coseps * z;y = ratio;/* convert from rectangular to polar */polar[0] = zatan2( x, y );polar[1] = asin( z/polar[2] );return(0);}#if DEBUGvoid prvec( vec )double vec[];{int k;for( k=0; k<3; k++ )	printf( "%17.9E ", vec[k] );printf( "\n" );}#endif#define DOUBLE doublevoid divdif();DOUBLE difpol();int jpl( JD, iobj, p, v )double JD;int iobj;double p[], v[];{int i, k;long recno;double JD0, x, t;double rec[73];double pcofs[NINTERP+2], vcofs[NINTERP+2];double pdiffs[NINTERP+2], vdiffs[NINTERP+2];/*iobj -= 1;*//* First record number later than requested date.  */x = ((JD - JDb) / JDi) + 1.0;recno = (long) x;JD0 = JDb + (double )recno * JDi;bcnt = (recno - NINTERP) * recsiz + rec0;#if DEBUGprintf( "recno = %ld, bcnt = %ld\n", recno, bcnt );#endift = (JD - JD0)/JDi; /* time scaled to fraction of subrecord */for( k=0; k<3; k++ )	{	  lseek( defd, bcnt, SEEK_SET );	  for( i=0; i<=NINTERP; i++ )	    {	      if( read( defd, rec, (int) recsiz ) != (int) recsiz )		return(-1);	      vcofs[i] = rec[6*iobj + 2*k + 1];	      pcofs[i] = rec[6*iobj + 2*k + 2];	    }	  divdif (vcofs, NINTERP, vdiffs);	  divdif (pcofs, NINTERP, pdiffs);	  p[k] = difpol (pdiffs, NINTERP, t);	  v[k] = difpol (vdiffs, NINTERP, t);	}/* The last record read should have the date JD0.  */if( rec[0] != JD0 )  printf("Record has wrong date %.16e.\n", rec[0]);#if DEBUGprvec(p);prvec(v);#endifreturn(0);}/* Interpolation routines taken from de118i.  */#define DOUBLE double#define One 1.0/* Compute zeroth through kth backward differences * of the data in the input array */void divdif(vec , k, diffn)DOUBLE vec[]; /* input array of k+1 data items */DOUBLE *diffn; /* output array of ith differences */int k;{DOUBLE diftbl[NINTERP+1];DOUBLE *p, *q;DOUBLE y;int i, o;/* Copy the given data (zeroth difference) into temp array */p = diftbl;q = vec;for( i=0; i<=k; i++ )	*p++ = *q++;/* On the first outer loop, k-1 first differences are calculated. * These overwrite the original data in the temp array. */o = k;for( o=k; o>0; o-- )	{	p = diftbl;	q = p;	for( i=0; i<o; i++ )		{		y = *p++;		*q++ = *p - y;		}	*diffn++ = *p; /* copy out the last (undifferenced) item */#if DEBUG	printf( "%.5e ", *p );#endif	}#if DEBUG	printf( "%.5e\n", *(q-1) );#endif*diffn++ = *(q-1);}/* Update array of differences, given new data value. * diffn is an array of k+1 differences, starting with the * zeroth difference (the previous original data value). */void dupdate( diffn, k, f )register DOUBLE *diffn;  /* input and output array of differences */int k; /* max order of differences */DOUBLE f; /* new data point (zeroth difference) */{DOUBLE new, old;int i;new = f;for( i=0; i<k; i++ )	{	old = *diffn;	*diffn++ = new;#if DEBUG	printf( "%.5e ", new );#endif	new = new - old;	}#if DEBUG	printf( "%.5e\n", new );#endif*diffn++ = new;}/* Evaluate the interpolating polynomial * *              (x - x ) *                    n    1 * P(x) = f  +  --------  D f   +  ... *         n       h         n * *     (x - x )(x - x   )...(x - x     ) *           n       n-1          n+2-k    k-1 *  +  ---------------------------------  D    f *                   k-1                        n *                  h     (k-1)! * * *         j *  where D denotes the jth backward difference, see dupdate(), and * *  f   =   f( x , y(x ) )  is the interpolated derivative y'(x ) . *   n          n     n                                        n * * The subroutine argument t is linearly scaled so that t = 1.0 * will evaluate the polynomial at x = x_n + h, * t = 0.0 corresponds to x = x_n, etc. */DOUBLE difpol( diffn, k, t )DOUBLE *diffn;int k; /* differences go up to order k-1 */DOUBLE t; /* scaled argument */{DOUBLE f, fac, s, u;int i;f = *diffn++; /* the zeroth difference = nth data point */u = One;/*s = x/h - n;*//*s = One; */	/* to evaluate the polynomial at x = xn + h */s = t;fac = One;for( i=1; i<k; i++ )	{	if( s == 0 )		break;	u *= s / fac;	f += u  *  (*diffn++);	fac += One;	s += One;	}return( f );}#endif  /* SSYSTEM */