JPL Ephemeris function declarations

Collaboration diagram for JPL Ephemeris function declarations:

Functions
void *	jpl_init_ephemeris (const char *ephemeris_filename, char nam[][6], double *val)
void	jpl_close_ephemeris (void *ephem)
int	jpl_state (void *ephem, const double et, const int list[12], double pv[][6], double nut[4], const int bary)
int	jpl_pleph (void *ephem, const double et, const int ntarg, const int ncent, double rrd[], const int calc_velocity)
double	jpl_get_double (const void *ephem, const int value)
double	jpl_get_long (const void *ephem, const int value)
int	make_sub_ephem (const void *ephem, const char *sub_filename, const double start_jd, const double end_jd)

Detailed Description

Function Documentation

void* 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;
    size_t count;

    if( !ifile)
    {
        errno = -JPLEPHEM_ERROR_NOTFOUND;
        return nullptr;
    }
    /* 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);
        errno = -JPLEPHEM_ERROR_INSUFFICIENT_MEMORY;
        return nullptr;
    }
    rval->ifile = ifile;
    count = fread( title, 84, 1, ifile);
    if (count)
    {
        fseek( ifile, 2652L, SEEK_SET);      /* skip title & constant name data */
        count = 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);
        count = fread( val, (size_t)rval->ncon, sizeof( double), ifile);
        if( count && 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++)
            count = fread( nam[i], 6, 1, ifile);
    }
    return( rval);
}

void jpl_close_ephemeris

(

void *

ephem

)

{
    struct jpl_eph_data *eph = (struct jpl_eph_data *)ephem;

    fclose( eph->ifile);
    free( ephem);
}

int 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 JPLEPHEM_ERROR_OUTOFRANGE;

    /*   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);
        size_t count = fread( buf, (size_t)eph->ncoeff, sizeof( double), eph->ifile);
        if(count && eph->swap_bytes)
            swap_double( buf, eph->ncoeff);
    }
    // Choose between KM and AU
    if (KM)
    {
        aufac = 1.;
        t[1] = eph->ephem_step * 86400.;
    }
    else
    {
        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);
}

int 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;
            rval = jpl_state( ephem, et, list, pv, rrd, 0);
        }
        else          /*  no nutations on the ephemeris file  */
            rval = JPLEPHEM_ERROR_NUTATIONS;
        return( rval);
    }

    /*  check for librations   */

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

    rval = jpl_state( eph, et, list, pv, rrd, 1);
    /* 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;
}

double jpl_get_double	(	const void *	ephem,
		const int	value
	)

{
    return( *(double *)( (char *)ephem + value));
}

double jpl_get_long	(	const void *	ephem,
		const int	value
	)

{
    return( *(long *)( (char *)ephem + value));
}

int make_sub_ephem	(	const void *	ephem,
		const char *	sub_filename,
		const double	start_jd,
		const double	end_jd
	)