/****************************************************************************
* finite-size.c
* David Boozer
* 15 April 2008
* 24 June 2008 (switch to new integrator)
*****************************************************************************
* Integrate equation of motion for a finite-size particle.
****************************************************************************/

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

#define ETA 1.0
#define T 20.0

#define DELTA 0.01
#define STEP_MAX 0.01

#define NUM_DOF 2

double history[10000];

union s_type {
  struct {
    double z, p;
  } name;
  double dof[NUM_DOF];
};

#define Z name.z
#define P name.p

/****************************************************************************
* evolve state in time from t to t+tau
****************************************************************************/

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

  num_steps = (int) ceil (tau/STEP_MAX);
  step = tau/num_steps;
  for (n=0; n<num_steps; n++) {
    (*f)(s_ptr, &s_dot_1, t);
    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, t);
    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, t);
    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, t);
    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]);
  }
}

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

double past_z (double t) {
  double u;
  int n;

  u = t/DELTA;
  n = (int) floor(u);
  if (n < 0)
    return 0.0;
  else
    return (u - n)*history[n+1] + (n + 1 - u)*history[n];
}

void f (union s_type *s_ptr, union s_type *s_dot_ptr, double t) {
  double sum, temp, s;

  s_dot_ptr->Z = s_ptr->P;

  for (s=0, sum=0.0; s<5.0; s+=DELTA) {
    temp = s_ptr->Z - past_z(t - s) - s;
    sum += exp(-temp*temp/4);
    temp = s_ptr->Z - past_z(t - s) + s;
    sum -= exp(-temp*temp/4);
  }
  s_dot_ptr->P = -sum*ETA*DELTA/(2*sqrt(M_PI));
}

/****************************************************************************
 * experiments
 ****************************************************************************/

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

  state.Z = 0.0;
  state.P = 0.3;

  for (t=0.0, n=0; t<T; t+=DELTA, n++) {
    history[n] = state.Z;
    evolve_state (f, &state, DELTA, t);
    printf ("%f, %f, %f\n", t, state.Z, state.P);
  }
  /*
  for (t=0; t<T; t+=0.01) printf ("%f, %f\n", t, past(t));
  */
}