#include <math.h>
#include <stdio.h>
#include <stdarg.h>
#include <stdlib.h>
#include <strings.h>

extern const char *__progname;

static FILE *df;

static int hdrksize;
static int ksize;
static int ncoeffs;

static char ttl[3][84];
static char cnam[400][6];
static double ss[3];
static int ncon;
static double au;
static double emrat;
static int ipt[12][3];
static int numde;
static int lpt[3];
static double cval[400];
static double et;
static int ntarg;
static int nctr;
static int ncoord;
static int * const tvec[3] = { &ntarg, &nctr, &ncoord };
static double xi;
static int bary = 0;
static int km = 0;
static double pvsun[6];

static char *line = 0;
static int linelen = 0;
static int linealloc = 0;
static int lineno = 0;

static void plaint(const char *, ...)
	__attribute__((__format__(__printf__,1,2)));
static void plaint(const char *fmt, ...)
{
 va_list ap;

 va_start(ap,fmt);
 fprintf(stderr,"%s: %d: ",__progname,lineno);
 vfprintf(stderr,fmt,ap);
 fprintf(stderr,"\n");
 va_end(ap);
}

static char *getline(int mustget)
{
 int c;

 lineno ++;
 linelen = 0;
 while (1)
  { c = getchar();
    if (c == EOF)
     { if (mustget)
	{ plaint("unexpected EOF");
	  exit(1);
	}
       if (linelen > 0) plaint("partial last line ignored");
       return(0);
     }
    if (linealloc >= linelen) line = realloc(line,linealloc=linelen+16);
    if (c == '\n')
     { line[linelen] = '\0';
       return(line);
     }
    line[linelen++] = c;
  }
}

static void doheader(void)
{
 hdrksize = -1;
 while (1)
  { getline(1);
    if (!strcmp(line,"EOT")) return;
    if (!strncmp(line,"KSIZE=",6)) hdrksize = atoi(line+6);
  }
}

static void compute_ksize(void)
{
 int mxx;
 int i;
 int nd;

 mxx = 0;
 for (i=1;i<12;i++) if (ipt[i][0] > ipt[mxx][0]) mxx = i;
 nd = (mxx == 11) ? 2 : 3;
 ncoeffs = ipt[mxx][0] + (nd*ipt[mxx][1]*ipt[mxx][2]) - 1;
 ksize = 2 * ncoeffs;
 printf("KHI = %6d, KMX = %6d, IPT(*,KHI) = %8d%8d%8d\n",
	mxx+1,ipt[mxx][0],ipt[mxx][0],ipt[mxx][1],ipt[mxx][2]);
 printf("ND = %6d, KSIZE = %6d\n",nd,ksize);
}

static void split(double in, double *out)
{
 double ipart;
 double fpart;

 fpart = modf(in,&ipart);
 if (fpart < 0)
  { fpart += 1;
    ipart -= 1;
  }
 out[0] = ipart;
 out[1] = fpart;
}

static void Xprintf(const char *, ...)
	__attribute__((__format__(__printf__,1,2)));
static void Xprintf(const char *fmt, ...)
{
#if 0
 va_list ap;
 va_start(ap,fmt);
 vprintf(fmt,ap);
 va_end(ap);
#else
 fmt = fmt;
#endif
}

static void interp(double *buf, double *t, int ncf, int ncm, int na, int ifl, double *pv)
{
 double dt1;
 double temp;
 double tc;
 double djunk;
 static double pc[18] = { 1, 0 };
 static int np = 2;
 static int nv = 3;
 static double twot = 0;
 static double vc[18] = { [1] = 1 };
 int l;
 int i;
 int j;
 double vfac;

 Xprintf("INTERP:\nT = %26.12f%26.12f\nNCF = %2d, NCM = %2d, NA = %2d, IFL = %2d\n",t[0],t[1],ncf,ncm,na,ifl);
 modf(t[0],&dt1);
 temp = na * t[0];
 l = ((int)(temp - dt1)) + 1;
 Xprintf("DNA = %26.12f, DT1 = %26.12f\nTEMP = %26.12f, L = %6d\n",(double)na,dt1,temp,l);
 tc = (2 * (modf(temp,&djunk) + dt1)) - 1;
 Xprintf("TC = %26.12f\n",tc);
 if (tc != pc[1])
  { Xprintf(" UPDATING NP/NV/PC(2)/TWOT\n");
    np = 2;
    nv = 3;
    pc[1] = tc;
    twot = tc + tc;
  }
 if (np < ncf)
  { Xprintf(" UPDATING PC\n");
    for (i=np+1;i<=ncf;i++)
     { pc[i-1] = (twot * pc[i-2]) - pc[i-3];
       Xprintf("PC(%2d) = %26.12f\n",i,pc[i-1]);
     }
    np = ncf;
  }
 for (i=0;i<ncm;i++)
  { pv[i] = 0;
    for (j=ncf-1;j>=0;j--)
     { pv[i] += pc[j] * buf[((((l-1)*ncm)+i)*ncf)+j];
       Xprintf("INTERP PV(%2d,1) = %26.12f\n       BUF(%2d,%2d,%2d) = %26.12f\n",
		i+1,pv[i],j+1,i+1,l,buf[((((l-1)*ncm)+i)*ncf)+j]);
     }
  }
 if (ifl <= 1) return;
 vfac = (2 * na) / t[1];
 vc[2] = twot + twot;
 Xprintf("VFAC = %26.12f, VC(3) = %26.12f\n",vfac,vc[2]);
 if (nv < ncf)
  { for (i=nv+1;i<=ncf;i++)
     { vc[i-1] = (twot * vc[i-2]) + pc[i-2] + pc[i-2] - vc[i-3];
       Xprintf("INTERP VC(%2d) = %26.12f\n",i,vc[i-1]);
       Xprintf("       TWOT = %26.12f\n",twot);
       Xprintf("       VC(%2d) = %26.12f\n",i-1,vc[i-2]);
       Xprintf("       PC(%2d) = %26.12f\n",i-1,pc[i-2]);
       Xprintf("       VC(%2d) = %26.12f\n",i-2,vc[i-3]);
     }
    nv = ncf;
  }
 for (i=0;i<ncm;i++)
  { pv[ncm+i] = 0;
    for (j=ncf-1;j>0;j--)
     { pv[ncm+i] += vc[j] * buf[((((l-1)*ncm)+i)*ncf)+j];
       Xprintf("INTERP PV(%2d,2) = %26.12f\n       BUF(%2d,%2d,%2d) = %26.12f\n",
		i+1,pv[ncm+i],j+1,i+1,l,buf[((((l-1)*ncm)+i)*ncf)+j]);
     }
    pv[ncm+i] *= vfac;
    Xprintf("IFINAL PV(%2d,2) = %26.12f\n",i+1,pv[ncm+i]);
  }
 Xprintf(" INTERP EXIT, PV =\n");
 l = 0;
 for (i=0;i<ncm;i++)
  { for (j=0;j<ifl;j++)
     { if (l > 2)
	{ Xprintf("\n");
	  l = 0;
	}
       Xprintf("%26.12f",pv[(j*ncm)+i]);
       l ++;
     }
  }
 Xprintf("\n");
}

static void state(double *et2, int *list, double (*pv)[6], double *pnut)
{
 double s;
 double pjd[4];
 double t[2];
 double aufac;
 static double buf[1500];
 static int nrl = -1;
 int nr;
 int i;
 int j;

 Xprintf(" STATE ENTERED\nET2 = %26.12f%26.12f\n",et2[0],et2[1]);
 s = et2[0] - .5;
 split(s,&pjd[0]);
 split(et2[1],&pjd[2]);
 pjd[0] += pjd[2] + .5;
 pjd[1] += pjd[3];
 split(pjd[1],&pjd[2]);
 pjd[0] += pjd[2];
 Xprintf("PJD = %26.12f%26.12f\n      %26.12f%26.12f\n",pjd[0],pjd[1],pjd[2],pjd[3]);
 if ((pjd[0]+pjd[3] < ss[0]) || (pjd[0]+pjd[3] > ss[1]))
  { plaint(" ***  Requested JED, %12.2f, not within ephemeris limits %12.2f - %12.2f",et2[0]+et2[1],ss[0],ss[1]);
    exit(1);
  }
 nr = floor((pjd[0] - ss[0]) / ss[2]) + 3;
 if (pjd[0] == ss[1]) nr --;
 t[0] = ((pjd[0] - (((nr-3) * ss[2]) + ss[0])) + pjd[3]) / ss[2];
 Xprintf("NR = %6d, T(1) = %26.12f\n",nr,t[0]);
 if (nr != nrl)
  { nrl = nr;
    Xprintf(" READING\n");
    fseek(df,(nr-1)*ksize*4,SEEK_SET);
    fread(&buf[0],8,ncoeffs,df);
  }
 if (km)
  { t[1] = ss[2] * 86400;
    aufac = 1;
  }
 else
  { t[1] = ss[2];
    aufac = 1 / au;
  }
 Xprintf("T(2) = %26.12f, AUFAC = %26.12f\n",t[1],aufac);
 Xprintf(" INTERP FOR SSBARY SUN (11)\n");
 interp(&buf[ipt[10][0]-1],&t[0],ipt[10][1],3,ipt[10][2],2,&pvsun[0]);
 for (i=0;i<6;i++) pvsun[i] *= aufac;
 Xprintf("PVSUN = \n%26.12f%26.12f%26.12f\n%26.12f%26.12f%26.12f\n",
	pvsun[0],pvsun[1],pvsun[2],pvsun[3],pvsun[4],pvsun[5]);
 for (i=0;i<10;i++)
  { if (list[i] == 0) continue;
    Xprintf("INTERP FOR %2d\n",i+1);
    interp(&buf[ipt[i][0]-1],&t[0],ipt[i][1],3,ipt[i][2],list[i],&pv[i][0]);
    for (j=0;j<6;j++)
     { if ((i <= 8) && !bary)
	{ pv[i][j] = (pv[i][j] * aufac) - pvsun[j];
	}
       else
	{ pv[i][j] *= aufac;
	}
     }
    Xprintf("PV(%2d) = \n%26.12f%26.12f%26.12f\n%26.12f%26.12f%26.12f\n",
	i+1,pv[i][0],pv[i][1],pv[i][2],pv[i][3],pv[i][4],pv[i][5]);
  }
 if ((list[10] > 0) && (ipt[11][1] > 0))
  { interp(&buf[ipt[11][0]-1],&t[0],ipt[11][1],3,ipt[11][2],list[10],pnut);
  }
 if ((list[11] > 0) && (ipt[12][1] > 0))
  { interp(&buf[ipt[12][0]-1],&t[0],ipt[12][1],3,ipt[12][2],list[11],&pv[10][0]);
  }
 Xprintf(" STATE EXIT, PV =\n");
 for (i=0;i<12;i++)
  { Xprintf("%26.12f%26.12f%26.12f\n%26.12f%26.12f%26.12f\n",
	pv[i][0],pv[i][1],pv[i][2],pv[i][3],pv[i][4],pv[i][5]);
  }
}

static void dpleph(double *et2, int ntarg, int ncent, double *r)
{
 int i;
 int list[12];
 double pv[13][6];
 int bsave;

 for (i=0;i<6;i++) r[i] = 0;
 Xprintf("DPLEPH ENTRY\nET2 = %26.12f%26.12f\nNTARG = %2d, NCENT = %2d\n",
	et2[0],et2[1],ntarg,ncent);
 if (ntarg == ncent) return;
 for (i=0;i<12;i++) list[i] = 0;
 if (ntarg == 14)
  { if (ipt[11][1] > 0)
     { list[10] = 2;
       state(et2,&list[0],&pv[0],r);
       return;
     }
    else
     { plaint(" *****  NO NUTATIONS ON THE EPHEMERIS FILE  *****");
       exit(1);
     }
  }
 if (ntarg == 15)
  { if (ipt[12][2] > 0)
     { list[11] = 2;
       state(et2,&list[0],&pv[0],r);
       for (i=0;i<6;i++) r[i] = pv[10][i];
       return;
     }
    else
     { plaint(" *****  NO LIBRATIONS ON THE EPHEMERIS FILE  *****");
       exit(1);
     }
  }
 bsave = bary;
 bary = 1;
 for (i=0;i<2;i++)
  { int k;
    k = i ? ncent : ntarg;
    if (k <= 10) list[k-1] = 2;
    if (k == 10) list[2] = 2;
    if (k == 3) list[9] = 2;
    if (k == 13) list[2] = 2;
  }
 for (i=0;i<13;i++)
  { int j;
    for (j=0;j<6;j++) pv[i][j] = 0;
  }
 Xprintf("CALLING STATE\nET2 = %26.12f%26.12f\nLIST = %2d%2d%2d%2d%2d%2d%2d%2d%2d%2d%2d%2d\n",
	et2[0],et2[1],list[0],list[1],list[2],list[3],list[4],list[5],
	list[6],list[7],list[8],list[9],list[10],list[11]);
 state(et2,&list[0],&pv[0],r);
 if ((ntarg == 11) || (ncent == 11))
  { for (i=0;i<6;i++) pv[10][i] = pvsun[i];
  }
 if ((ntarg == 12) || (ncent == 12))
  { for (i=0;i<6;i++) pv[11][i] = 0;
  }
 if ((ntarg == 13) || (ncent == 13))
  { for (i=0;i<6;i++) pv[12][i] = pv[2][i];
  }
 if ( ((ntarg == 10) && (ncent == 3)) ||
      ((ntarg == 3) && (ncent == 10)) )
  { for (i=0;i<6;i++) pv[2][i] = 0;
  }
 else
  { if (list[2] == 2)
     { for (i=0;i<6;i++) pv[2][i] -= pv[9][i] / (1+emrat);
     }
    if (list[9] == 2)
     { for (i=0;i<6;i++) pv[9][i] += pv[2][i];
     }
  }
 for (i=0;i<6;i++) r[i] = pv[ntarg-1][i] - pv[ncent-1][i];
 bary = bsave;
 Xprintf("DPLEPH EXIT, RRD =\n%26.12f%26.12f%26.12f\n%26.12f%26.12f%26.12f\n",
	r[0],r[1],r[2],r[3],r[4],r[5]);
}

static void pleph(double et, int ntarg, int ncent, double *r)
{
 double et2[2];

 et2[0] = et;
 et2[1] = 0;
 dpleph(&et2[0],ntarg,ncent,r);
}

#if 0

static int txt_inc(char *digs, int n)
{
 for (n--;(n>=0)&&(digs[n]=='9');n--) digs[n] = '0';
 if (n < 0) return(1);
 digs[n] ++;
 return(0);
}

static int txt_dec(char *digs, int n)
{
 for (n--;(n>=0)&&(digs[n]=='0');n--) digs[n] = '9';
 if (n < 0) return(1);
 digs[n] --;
 return(0);
}

static void txt_exp_inc(char *exp)
{
 if ((exp[1] == '-') && (exp[2] == '0') && (exp[3] == '1'))
  { exp[1] = '+';
    exp[3] = '0';
  }
 else if (exp[1] == '-')
  { txt_dec(exp+2,2);
  }
 else if (exp[1] == '+')
  { txt_inc(exp+2,2);
  }
}

static void shift_exp_txt_round(char *t)
{
 char *dot;
 char *exp;

 dot = index(t,'.');
 exp = index(t,'E');
 if (!dot || !exp) return;
 bcopy(dot+1,dot+2,(exp-dot)-2);
 dot[1] = dot[-1];
 dot[-1] = '0';
 txt_exp_inc(exp);
 if (txt_inc(dot+1,(exp-dot)-1))
  { dot[1] = 1;
    txt_exp_inc(exp);
  }
}

#endif

int main(void);
int main(void)
{
 int c;
 int i;
 double r[6];
 double err;

 if (sizeof(int) != 4)
  { plaint("sizeof(int) != 4");
    exit(1);
  }
 if (sizeof(double) != 8)
  { plaint("sizeof(double) != 8");
    exit(1);
  }
 doheader();
 df = fopen("JPLEPH","r");
 fread(&ttl[0][0],1,3*84,df);
 fread(&cnam[0][0],1,400*6,df);
 fread(&ss[0],8,3,df);
 fread(&ncon,4,1,df);
 fread(&au,8,1,df);
 fread(&emrat,8,1,df);
 fread(&ipt[0][0],4,12*3,df);
 fread(&numde,4,1,df);
 fread(&lpt[0],4,3,df);
 printf("  JPL TEST-EPHEMERIS program.  Last modified July 1997.\n");
 compute_ksize();
 if ((hdrksize >= 0) && (hdrksize != ksize))
  { plaint("warning: ksize in header (%d) != ksize from file (%d)",hdrksize,ksize);
  }
 printf("KSIZE = %6d, NCOEFFS = %6d, IRECSZ = %6d\n",ksize,ncoeffs,ksize*4);
 printf("\n%14.2f%14.2f%14.2f\n",ss[0],ss[1],ss[2]);
 fseek(df,ksize*4,SEEK_SET);
 fread(&cval[0],8,400,df);
#if 0
 for (i=0;i<ncon;i++)
  { char b[64];
    sprintf(&b[0],"%24.16E",cval[i]);
    if (strcmp(&b[0],"  0.0000000000000000E+00"))
     { shift_exp_txt_round(&b[0]);
     }
    printf("%s  %6.6s%s",(i&1)?"":"\n",&cnam[i][0],&b[0]);
  }
 printf("\n");
#endif
 printf("   line -- jed --   t#   c#   x#   --- jpl value"
	" ---   --- user value --    -- difference --\n\n");
 while (getline(0))
  { if (strlen(line) < 15+10+(3*3)+1)
     { plaint("input line too short: %s",line);
       continue;
     }
    c = line[15+10];
    line[15+10] = '\0';
    et = atof(line+15);
    line[15+10] = c;
    for (i=0;i<3;i++)
     { c = line[15+10+(3*i)+3];
       line[15+10+(3*i)+3] = '\0';
       *tvec[i] = atoi(line+15+10+(3*i));
       line[15+10+(3*i)+3] = c;
     }
    xi = atof(line+15+10+(3*3));
    if ((et < ss[0]) || (et > ss[1])) continue;
    pleph(et,ntarg,nctr,&r[0]);
    err = fabs(r[ncoord-1]-xi);
    if ((ntarg == 15) && (ncoord == 3)) err /= (.23 * (et - 2451545));
    if (((lineno % 100) == 0) || (err > 1e-13))
     { printf("%6d%10.1f%5d%5d%5d%20.13f%20.13f%20.13f%s\n",
		lineno,et,ntarg,nctr,ncoord,xi,r[ncoord-1],err,
		(err>1e-13)?"  **** ERROR":"");
     }
  }
 exit(0);
}