text: /* crack natural satellite files from BDL */
  
  #include <math.h>
  
  #include "astro.h"
  #include "bdl.h"
  
  int read_bdl (FILE *fp, double jd, double *xp, double *yp, double *zp,
                char ynot[]) { return 0;}
  
  /* using a BDL planetary moon dataset defined in a struct in RAM and a
   * JD, find the x/y/z positions of each satellite. store in the given arrays,
   * assumed to have one entry per moon. values are planetocentric, +x east, +y
   * north, +z away from earth, all in au. corrected for light time.
   * files obtained from ftp://ftp.bdl.fr/pub/misc/satxyz.
   */
  void
  do_bdl (BDL_Dataset *dataset, double jd, double *xp, double *yp, double *zp)
  {
          int nsat = dataset->nsat;
          double djj = dataset->djj;
          unsigned *idn = dataset->idn;
          double *freq = dataset->freq;
          double *delt = dataset->delt;
  	double t0;
  	int i;
  
  	/* compute location of each satellite */
  	for (i = 0; i < nsat; i++) {
  	    int id = (int)floor((jd-djj)/delt[i]) + idn[i] - 2;
  	    double t1, anu, tau, tau2, at;
  	    double tbx, tby, tbz;
              double *cmx, *cfx, *cmy, *cfy, *cmz, *cfz;
              BDL_Record *r = & dataset->moonrecords[id];
  
              t0 = r->t0;
  	    t1 = floor(t0) + 0.5;
  	    anu = freq[i];
  	    tau = jd - t1;
  	    tau2 = tau * tau;
  	    at = tau*anu;
  
              cmx = & (r->cmx[0]); /* point at data in appropriate record */
              cfx = & (r->cfx[0]);
              cmy = & (r->cmy[0]);
              cfy = & (r->cfy[0]);
              cmz = & (r->cmz[0]);
              cfz = & (r->cfz[0]);
  
  	    tbx = cmx[0]+cmx[1]*tau+cmx[2]*sin(at+cfx[0])
  			    +cmx[3]*tau*sin(at+cfx[1])
  			    +cmx[4]*tau2*sin(at+cfx[2])
  			    +cmx[5]*sin(2*at+cfx[3]);
  	    tby = cmy[0]+cmy[1]*tau+cmy[2]*sin(at+cfy[0])
  			    +cmy[3]*tau*sin(at+cfy[1])
  			    +cmy[4]*tau2*sin(at+cfy[2])
  			    +cmy[5]*sin(2*at+cfy[3]);
  	    tbz = cmz[0]+cmz[1]*tau+cmz[2]*sin(at+cfz[0])
  			    +cmz[3]*tau*sin(at+cfz[1])
  			    +cmz[4]*tau2*sin(at+cfz[2])
  			    +cmz[5]*sin(2*at+cfz[3]);
  
  	    xp[i] = tbx*1000./149597870.;
  	    yp[i] = tby*1000./149597870.;
  	    zp[i] = tbz*1000./149597870.;
  	}
  }
  
  /* For RCS Only -- Do Not Edit */
  static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: bdl.c,v $ $Date: 2008/04/20 08:11:35 $ $Revision: 1.6 $ $Name:  $"};
  