#include <stdio.h>
#include <math.h>

struct POLYCO {
  char psrname[80];
  char date[80];
  double utc;
  double tmid;
  double dm;
  double doppler;
  double lfrms;
  double rphase;
  double f0;
  double fobs;
  int span;
  int obsno;
  int nc;
  double *coeff;
} ;

int read_polycoset(FILE *polycofile, struct POLYCO *polyco) /*includefile*/
{
  char coefficient[25];
  int i,j;

  /* read in the first two lines of the polyco set (header info) */
  if((fscanf(polycofile,"%10s%9s%12lf%20lf%21lf%7lf%7lf",
      &polyco->psrname,&polyco->date,&polyco->utc,&polyco->tmid,&polyco->dm,
      &polyco->doppler,&polyco->lfrms))== -1) return(NULL);
  if ((fscanf(polycofile,"%20lf%18lf%5i%6i%5i%9lf\n",&polyco->rphase,
        &polyco->f0,&polyco->obsno,&polyco->span,&polyco->nc,
	&polyco->fobs)) == -1) return(NULL);

  free(polyco->coeff);
  polyco->coeff = (double *) malloc((polyco->nc)*sizeof(double));

  i=0;
  while (i<polyco->nc) {
    /* look for coefficients of the form .....D... (TEMPO - f77 style) */
    if ((fscanf(polycofile,"%25s",coefficient))==-1) return(NULL);
    if (strstr(coefficient,"D") != NULL) {
      /* replace the D with an e and copy string to coefficient array */
      for (j=0; j<strlen(coefficient); j++) {
	if(coefficient[j]=='D') coefficient[j]='e';
      }
      polyco->coeff[i]=atof(coefficient);
      i++;
    }
  }
  return(polyco->nc);
}

void get_nearest_polyco(char *filename, double mjd, struct POLYCO *polyco) /* includefile */
{
  double diff,best;
  int i,set;

  FILE *polycofile;
  polycofile=fopen(filename,"r");

  best=1.0e32;
  i=0;
  while (read_polycoset(polycofile, polyco) != NULL) 
    {
      diff=abs((polyco->tmid)-mjd);
      if (diff<best) {
	set=i;
	best=diff;
      }
      i++;
    }
  fclose(polycofile);

  polycofile=fopen(filename,"r");
  for(i=1; i<=set; i++) read_polycoset(polycofile,polyco);
  fclose(polycofile);
}

double polyco_period(double mjd, struct POLYCO polyco) /*includefile*/
{
  int i;
  double dt,freq,factor;
  dt=(mjd-polyco.tmid)*1440.0;
  factor=polyco.coeff[1];
  for (i=2; i<polyco.nc; i++) factor+=dt*i*polyco.coeff[i]; 
  freq=polyco.f0+(factor/60.0);
  return(1.0/freq);
}

double polyco_phase(double mjd, struct POLYCO polyco) /*includefile*/
{
  int i;
  double dt,phase;
  dt=(mjd-polyco.tmid)*1440.0;
  phase=polyco.rphase +(dt*60.0*polyco.f0)+polyco.coeff[0];
  for (i=1; i<polyco.nc; i++) phase+=pow(dt,i)*polyco.coeff[i]; 
  return(phase-floor(phase));
}

main (int argc, char *argv[]) 
{
  struct POLYCO p;
  double mjd;
  if (argc < 3) {
    fprintf(stderr,"usage: readpolyco filename mjd\n");
    exit(1);
  }
  mjd=atof(argv[2]);
  get_nearest_polyco(argv[1],atof(argv[2]),&p);
  printf("%s TMID %f MJD P0 %.10f SEC\n",argv[1],p.tmid,1.0/p.f0);
  printf("T %s MJD P(t) %.10f SEC PHI %.10f DM %f cm-3 pc\n",
	 argv[2],polyco_period(mjd,p),polyco_phase(mjd,p),p.dm);
}