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

// d^2 x/dtheta^2 + Q_O^2 x = eps * cos(m*theta) * x^p

double eps; 
int Nturns;  // number of turns for each particle
int Nsteps;  // number of steps per turn
int p; 	     // power of x in perturbation driving term
int m;       // azimuthal harmonic driving perturbation 
double xiO, xpiO; 
double QO;   // tune of unperturbed oscillator 
int npart;   // number of particles to track 
double rat; 

FILE *f; 
//======= 
int process_cmd(int argc, char **argv)

{
  int narg, i;

  narg = 10+1;
  if(argc!=narg)
    {
      fprintf(stderr,
      "Usage: %s Nturns, Nsteps, eps, m, pow, xO, x'O, QO,",
      argv[0]);
      fprintf(stderr, " npart, rat\n");
      fprintf(stderr, " Nturns:\t num turns per particle\n");
      fprintf(stderr, " Nsteps:\t num integ steps/turn\n");
      fprintf(stderr, " eps:\t epsilon\n");
      fprintf(stderr, " m:\t azimuthal harmonic\n");
      fprintf(stderr, " pow:\t power of x in drive term\n");
      fprintf(stderr, " xiO:\t init x for 1st particle\n");
      fprintf(stderr, " xpiO:\t init x' for 1st particle\n");
      fprintf(stderr, " Q_O for unperturbed oscillator\n");
      fprintf(stderr, " npart:\t num of particles\n");
      fprintf(stderr, " ratio:\t multiplier of [xiO, xpiO]");
      fprintf(stderr, "for succeeding particles\n");
      return -1;
    }
  i = 1;
  Nturns = atoi(argv[i++]);
  Nsteps = atoi(argv[i++]);
  eps    = atof(argv[i++]);
  m      = atoi(argv[i++]);
  p      = atoi(argv[i++]);
  xiO    = atof(argv[i++]);
  xpiO   = atof(argv[i++]);
  QO     = atof(argv[i++]);
  npart  = atoi(argv[i++]);
  rat    = atof(argv[i++]);

  return 0;
}

//-------
void dumppar(void)

   {
    fprintf(f, "# %s:\t %d \n", "Nturns", Nturns);
    fprintf(f, "# %s:\t %d \n", "Nsteps", Nsteps);
    fprintf(f, "# %s:\t %f \n", "eps", eps);
    fprintf(f, "# %s:\t %d \n", "m", m);
    fprintf(f, "# %s:\t %d \n", "p", p);
    fprintf(f, "# %s:\t %f \n", "xiO", xiO);
    fprintf(f, "# %s:\t %f \n", "xpiO", xpiO);
    fprintf(f, "# %s:\t %f \n", "QO", QO);
    fprintf(f, "# %s:\t %d \n", "npart", npart);
    fprintf(f, "# %s:\t %f \n", "rat", rat);
   }
//-------
void trackit(void)
   {
    double x, xp, xi, xpi, xt, xtp;
    double fdrive;
    double mustep;
    double ma, mb, mc, md; // linear matrix
    int nturn, ipart, jstep, kpow;
    double pi;
    double mtheta, dmtheta;
    double a, b, c, d;

    pi = acos(-1.0);
    mustep = 2.0*pi*QO/((double) Nsteps);
    a = cos(mustep);
    b = sin(mustep);
    c = -b;
    d = a;
    dmtheta = 2.*pi*m/((double) Nsteps);

    xi  = xiO;
    xpi = xpiO;
    for(ipart=0;ipart<npart;ipart++)
      {
       x = xi;
       xp = xpi;
       fprintf(f, "%6d %6d %e %e\n", ipart, 0, x, xp);

    for(nturn=1;nturn<Nturns;nturn++)
{
 mtheta = 0;
 for(jstep=0;jstep<Nsteps;jstep++)
   {
    xt = a*x+b*xp;
    xtp = c*x+d*xp;
    fdrive = eps*cos(mtheta)*dmtheta;
    for(kpow=0;kpow<p;kpow++)
{
 fdrive = fdrive*xt;
}
    x = xt;
    xp = xtp+fdrive;
    mtheta += dmtheta;
   }
fprintf(f, "%6d %6d %e %e\n", ipart, nturn, x, xp);
}
    // xi  = xi*rat;
    // xpi = xpi*rat;
    xi  = xiO*(1.0+(rat-1.0)*ipart);
    xpi = xpiO*(1.0+(rat-1.0)*ipart);
   }
}
//-------
int main(int argc, char** argv)
{
   if(process_cmd(argc, argv)) return -1;

   f = fopen("track.dat", "w");
   dumppar();
   trackit();
   fclose (f);

   return 0;
}