#include "ephlib.h"

#include <math.h>
#include <stdio.h>
#include <errno.h>
#include <stdarg.h>
#include <stdlib.h>
/*
#include <strings.h>
*/

typedef struct eph EPH;

struct eph {
  FILE *f;
  int ncoeffs;
  double ss[3];
  int ncon;
  double au;
  double emrat;
  int ipt[12][3];
  int numde;
  int lpt[3];
  double pvsun[6];
  double buf[1500];
  double nrl;
  double pc[18];
  int np;
  int nv;
  double twot;
  double vc[18];
  } ;

static void compute_ncoeffs(EPH *e)
{
 int mxx;
 int i;
 int nd;

 mxx = 0;
 for (i=1;i<12;i++) if (e->ipt[i][0] > e->ipt[mxx][0]) mxx = i;
 nd = (mxx == 11) ? 2 : 3;
 e->ncoeffs = e->ipt[mxx][0] + (nd*e->ipt[mxx][1]*e->ipt[mxx][2]) - 1;
}

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 interp(EPH *e, double *buf, double *t, int ncf, int ncm, int na, int ifl, double *pv)
{
 double dt1;
 double temp;
 double tc;
 double djunk;
 int l;
 int i;
 int j;
 double vfac;

 modf(t[0],&dt1);
 temp = na * t[0];
 l = ((int)(temp - dt1)) + 1;
 tc = (2 * (modf(temp,&djunk) + dt1)) - 1;
 if (tc != e->pc[1])
  { e->np = 2;
    e->nv = 3;
    e->pc[1] = tc;
    e->twot = tc + tc;
  }
 if (e->np < ncf)
  { for (i=e->np+1;i<=ncf;i++)
     { e->pc[i-1] = (e->twot * e->pc[i-2]) - e->pc[i-3];
     }
    e->np = ncf;
  }
 for (i=0;i<ncm;i++)
  { pv[i] = 0;
    for (j=ncf-1;j>=0;j--)
     { pv[i] += e->pc[j] * buf[((((l-1)*ncm)+i)*ncf)+j];
     }
  }
 if (ifl <= 1) return;
 vfac = (2 * na) / t[1];
 e->vc[2] = e->twot + e->twot;
 if (e->nv < ncf)
  { for (i=e->nv+1;i<=ncf;i++)
     { e->vc[i-1] = (e->twot * e->vc[i-2]) + e->pc[i-2] + e->pc[i-2] - e->vc[i-3];
     }
    e->nv = ncf;
  }
 for (i=0;i<ncm;i++)
  { pv[ncm+i] = 0;
    for (j=ncf-1;j>0;j--)
     { pv[ncm+i] += e->vc[j] * buf[((((l-1)*ncm)+i)*ncf)+j];
     }
    pv[ncm+i] *= vfac;
  }
}

static int state(EPH *e, double date, int *list, double (*pv)[6])
{
 double pjd[2];
 double t[2];
 double aufac;
 int nr;
 int i;
 int j;

 split(date-.5,&pjd[0]);
 pjd[0] += .5;
 if ((pjd[0]+pjd[1] < e->ss[0]) || (pjd[0]+pjd[1] > e->ss[1]))
  { return(EPHL_ERR_RANGE);
  }
 nr = floor((pjd[0] - e->ss[0]) / e->ss[2]) + 3;
 if (pjd[0] == e->ss[1]) nr --;
 t[0] = ((pjd[0] - (((nr-3) * e->ss[2]) + e->ss[0])) + pjd[1]) / e->ss[2];
 if (nr != e->nrl)
  { e->nrl = nr;
    fseek(e->f,(nr-1)*e->ncoeffs*8,SEEK_SET);
    fread(&e->buf[0],8,e->ncoeffs,e->f);
  }
 t[1] = e->ss[2];
 aufac = 1 / e->au;
 interp(e,&e->buf[e->ipt[10][0]-1],&t[0],e->ipt[10][1],3,e->ipt[10][2],2,&e->pvsun[0]);
 for (i=0;i<6;i++) e->pvsun[i] *= aufac;
 for (i=0;i<10;i++)
  { if (list[i] == 0) continue;
    interp(e,&e->buf[e->ipt[i][0]-1],&t[0],e->ipt[i][1],3,e->ipt[i][2],list[i],&pv[i][0]);
    for (j=0;j<6;j++) pv[i][j] *= aufac;
  }
 return(0);
}

void *eph_open(const char *fn)
{
 EPH *e;
 int err;
 int i;

 if (sizeof(int) != 4)
  { fprintf(stderr,"sizeof(int) != 4; ephlib needs work\n");
    exit(1);
  }
 if (sizeof(double) != 8)
  { fprintf(stderr,"sizeof(double) != 8; ephlib needs work\n");
    exit(1);
  }
 e = malloc(sizeof(EPH));
 if (e == 0) return(0);
 e->f = fopen(fn,"r");
 if (e->f == 0)
  { err = errno;
    free(e);
    errno = err;
    return(0);
  }
 fseek(e->f,(3*84)/*ttl*/+(400*6)/*cnam*/,SEEK_SET);
 fread(&e->ss[0],8,3,e->f);
 fread(&e->ncon,4,1,e->f);
 fread(&e->au,8,1,e->f);
 fread(&e->emrat,8,1,e->f);
 fread(&e->ipt[0][0],4,12*3,e->f);
 fread(&e->numde,4,1,e->f);
 fread(&e->lpt[0],4,3,e->f);
 compute_ncoeffs(e);
 e->nrl = -1;
 for (i=0;i<18;i++)
  { e->pc[i] = 0;
    e->vc[i] = 0;
  }
 e->pc[0] = 1;
 e->np = 2;
 e->nv = 3;
 e->twot = 0;
 e->vc[1] = 1;
 return(e);
}

void eph_close(void *eph)
{
 EPH *e;

 e = eph;
 fclose(e->f);
 free(e);
}

int eph_lookup(void *eph, double date, int targ, int cent, double *out)
{
 EPH *e;
 int i;
 int list[12];
 double pv[13][6];

 e = eph;
 if (targ == cent)
  { for (i=0;i<6;i++) out[i] = 0;
    return(0);
  }
 for (i=0;i<12;i++) list[i] = 0;
 for (i=0;i<2;i++)
  { int k;
    k = i ? cent : targ;
    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;
  }
 i = state(e,date,&list[0],&pv[0]);
 if (i) return(i);
 if ((targ == 11) || (cent == 11)) for (i=0;i<6;i++) pv[10][i] = e->pvsun[i];
 if ((targ == 12) || (cent == 12)) for (i=0;i<6;i++) pv[11][i] = 0;
 if ((targ == 13) || (cent == 13)) for (i=0;i<6;i++) pv[12][i] = pv[2][i];
 if ( ((targ == 10) && (cent == 3)) ||
      ((targ == 3) && (cent == 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+e->emrat);
    if (list[9] == 2) for (i=0;i<6;i++) pv[9][i] += pv[2][i];
  }
 for (i=0;i<6;i++) out[i] = pv[targ-1][i] - pv[cent-1][i];
 return(0);
}