/*---------------------------------------------------*/
/* getarrtime.c                                      */
/* a c translation of the Jodrell Bank Observatory   */
/* pulse arrival time program from Crab Notes        */
/* compile using eg gcc -o getarrtime getarrtime.c   */
/*---------------------------------------------------*/ 

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

int main ()
{

  double  epoch, secs, v, vdot,freq;
  double treq, ti, p, pdot, vddot;
  double turns, preq, tarr;
  char q='y',eph[20];
 
  printf("Enter values from pulsar ephemeris...... \n");
  
  printf("Epoch of ephemeris (MJD)       :");
  scanf("%lf",&epoch);
  printf("Epoch of ephemeris (secs)      :");
  scanf("%lf",&secs);
  epoch = epoch + secs/86400.0;
  
  printf("Frequency (nu) in Hz           :");
  scanf("%lf",&v);
  printf("Freqdot   (nudot) in e-15Hz/sec:");
  scanf("%lf",&vdot);
  vdot = vdot* 1e-15;
  
  while (q=='y') {
    printf("\nRequired arrival time (mjd):");
    scanf("%lf",&treq);
    
    /* Calculate arrival time */
    ti = (treq-epoch)*86400.0;
    p = 1.0/v;
    pdot = -vdot/(v*v);
    vddot = 2.0*pdot*pdot/(p*p*p);
    freq = v + ti*vdot + ti*ti/2.0*vddot;
    preq = 1.0/freq;
    
    turns = v*ti + vdot*ti*ti/2.0 + vddot*ti*ti*ti/6.0;

    /*    printf("ti = %lf, p = %lf, pdot = %f, vddot = %lf, freq = %lf\n",
          ti,p,pdot,vddot,freq);
    printf("turns %lf\n",turns);
    */

    turns = turns - (int)turns;
    tarr = treq + (1.0 - turns) *preq/86400.0;
      
    /*    printf("tarr = %lf\n",tarr); */
    printf(    
"\nFirst pulse arrival time after required epoch (mjd & secs):%d %lf\n",
        (int)tarr, (tarr - (int)tarr)*86400.0);
    printf("Preq(secs) and nu(Hz): %.9lf %.9lf\n\n",preq,freq);
  
    /* Clear input buffer! */
    q=getchar();
    
    printf("Another epoch (n/y):");
    q=getchar();
    
    q=tolower(q);
    }
  exit(0);
    
  }