/****************************************************************************
* two-particle.c
* David Boozer
* 04 May 2008
*****************************************************************************
* Numerically integrate equation of motion for two particles.
****************************************************************************/

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

#define NUM_DOF 3

#define TAU 20000.0
#define DELTA 0.1

#define HISTORY_MAX 20000

#define M 5000.0
#define GAMMA (2/M)

/* note: unnamed struct in union may not work in some compilers */
union s_type {
  struct {
    double z, v, t;
  } name;
  double dof[NUM_DOF];
};

union s_type history[HISTORY_MAX];

#define Z name.z
#define V name.v
#define T name.t

/****************************************************************************
 * manipulate states
 ****************************************************************************/

void evolve_state (void (*f)(union s_type *s_ptr, union s_type *s_dot_ptr),
		   union s_type *s_ptr, double step) {
  
  union s_type s_dot_1, s_dot_2, s_dot_3, s_dot_4;
  union s_type s_temp;
  int i;

  (*f)(s_ptr, &s_dot_1);
  for (i=0; i<NUM_DOF; i++)
    s_temp.dof[i] = s_ptr->dof[i] + 0.5*step*s_dot_1.dof[i];
  (*f)(&s_temp, &s_dot_2);
  for (i=0; i<NUM_DOF; i++)
    s_temp.dof[i] = s_ptr->dof[i] + 0.5*step*s_dot_2.dof[i];
  (*f)(&s_temp, &s_dot_3);
  for (i=0; i<NUM_DOF; i++)
    s_temp.dof[i] = s_ptr->dof[i] + step*s_dot_3.dof[i];
  (*f)(&s_temp, &s_dot_4);
  for (i=0; i<NUM_DOF; i++)
    s_ptr->dof[i] += (step/6)*(s_dot_1.dof[i] + 2*s_dot_2.dof[i] +
			       s_dot_4.dof[i] + 2*s_dot_3.dof[i]);
}

void initialize_state (union s_type *s_ptr) {
  s_ptr->Z = 500.0;
  s_ptr->V = 0.0;
  s_ptr->T = 0.0;
}

/****************************************************************************
* equations of motion
****************************************************************************/

/* set *s_ptr to the state of the particle at time t */
void past (double t, union s_type *s_ptr) {
  double ym, yp, tm, tp;
  int n, k;

  n = (int) floor(t/DELTA);

  if (n < 0)
    initialize_state (s_ptr);
  else
    for (k=0; k<NUM_DOF; k++) {
      ym = history[n%HISTORY_MAX].dof[k];
      yp = history[(n+1)%HISTORY_MAX].dof[k];
      tm = n*DELTA;
      tp = (n+1)*DELTA;
      s_ptr->dof[k] = ((t - tm)*yp + (tp - t)*ym)/DELTA;
    }
}

/* return the retarded time corresponding to the event (t, x) */
double t_ret (double t, double x) {
  union s_type state_r;
  double tau, tau_old;
  int sign;
  
  tau = 0.0;
  do {
    past (t - tau, &state_r);
    sign = (x + state_r.Z > 0) ? 1 : -1;
    tau_old = tau;
    tau = tau - (tau - fabs(x + state_r.Z))/(1 + sign*state_r.V);
  } while (fabs(tau - tau_old) > DELTA);

  return (t - tau);
}

/* two-particle equation of motion */
void eqns_two_particles (union s_type *s_ptr, union s_type *s_dot_ptr) {
  union s_type state_r;
  double z = s_ptr->Z;
  double v = s_ptr->V;
  double t = s_ptr->T;
  double v_r;
  int sign;

  sign = (z > 0) ? 1 : -1;
  past (t_ret (t, z), &state_r);
  v_r = state_r.V;

  s_dot_ptr->V = -GAMMA*(v/(1 - v*v) + sign/(1 + sign*v_r));
  s_dot_ptr->Z = v;
  s_dot_ptr->T = 1.0;
}

int main () {
  union s_type state;
  double t;
  int n;

  initialize_state (&state);
  for (t=0.0, n=0; t<TAU; t+=DELTA, n++) {
    if (n%50 == 0) printf ("%f, %f, %f\n", state.T, state.V, state.Z);
    history[n%HISTORY_MAX] = state;
    history[(n+1)%HISTORY_MAX] = state;
    evolve_state (eqns_two_particles, &state, DELTA);
  }
}