/****************************************************************************
* electrodynamics.c
* David Boozer
* 23 June 2008
* 30 July 2008 (modified)
*****************************************************************************
* Simulate toy model of electrodynamics.
*****************************************************************************
* How to convert sequence of frames into a movie:
*
* $ makejpgs
* $ mencoder \*.jpg -mf on:fps=25 -ovc lavc
****************************************************************************/

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

#define R 2000
#define N (2*R + 1)
#define NUM_DOF (2*N+2)

#define SIGMA 250.0
#define STEP_MAX 0.1
#define TOTAL_CHARGE 1.0

/* charge density profile of particle */
#define F(x) (exp(-(x)*(x)/(2*SIGMA*SIGMA))/(SIGMA*sqrt(2*M_PI)))

/* first derivative of charge density profile */
#define F1(x) (-(x)*F(x)/(SIGMA*SIGMA))

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

#define PHI name.phi
#define PI name.pi
#define Z name.z
#define P name.p

double M, omega0;

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

void evolve_state (void (*f)(union s_type *s_ptr, union s_type *s_dot_ptr),
		   union s_type *s_ptr, double tau) {
  
  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);
    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 reverse_state (union s_type *s_ptr) {
  int n;

  for (n=0; n<N; n++) s_ptr->PI[n] = -s_ptr->PI[n];
  s_ptr->P = -s_ptr->P;
}

void initialize_state (union s_type *s_ptr) {
  double temp;
  int i, j;

  s_ptr->Z = 0.0;
  s_ptr->P = 0.0;
  for (i=0; i<N; i++) {
    for (j=0, temp=0.0; j<N; j++) temp += abs(i-j)*F(j-R);
    s_ptr->PHI[i] = temp;
    s_ptr->PI[i] = 0.0;
  }
}

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

/* equations of motion for a simple harmonic oscillator */
void eqns_sho (union s_type *s_ptr, union s_type *s_dot_ptr) {
  s_dot_ptr->Z = s_ptr->P;
  s_dot_ptr->P = -s_ptr->Z - s_ptr->P;
}

/* equations of motion for a harmonically bound particle */
void eqns_harmonic (union s_type *s_ptr, union s_type *s_dot_ptr) {
  double *phi = s_ptr->PHI;
  double *pi = s_ptr->PI;
  double z = s_ptr->Z;
  double p = s_ptr->P;
  double temp;
  int n;

  for (n=0; n<N; n++) {
    if (n == 0)
      s_dot_ptr->PI[n] = 2*(phi[1] - phi[0] - pi[0] + TOTAL_CHARGE);
    else if (n == N-1)
      s_dot_ptr->PI[n] = 2*(phi[N-2] - phi[N-1] - pi[N-1] + TOTAL_CHARGE);
    else
      s_dot_ptr->PI[n] = phi[n+1] - 2*phi[n] + phi[n-1] - 2*F(z-(n-R));

    s_dot_ptr->PHI[n] = pi[n];
  }

  s_dot_ptr->Z = p/M;

  for (n=0, temp=0.0; n<N; n++) temp -= phi[n]*F1(z-(n-R));
  s_dot_ptr->P = 2*temp - M*omega0*omega0*z;
}

/****************************************************************************
 * output files
 ****************************************************************************/

void make_frames (void *f, union s_type *s_ptr, double T) {
  FILE *file_field, *file_particle;
  char filename[32];
  double *pi = s_ptr->PI;
  double *phi = s_ptr->PHI;
  double *p = &s_ptr->P;
  double *z = &s_ptr->Z;
  double E, B, Erad, Es[N], temp;
  double t, delta;
  int i, n;

  delta = 4.0;
  file_particle = fopen ("particle.txt", "w");
  for (t=0.0, n=0; t<=T; t+=delta, n++) {
    /* calculate static electric field */
    Es[0] = -TOTAL_CHARGE;
    for (i=1; i<N; i++) Es[i] = Es[i-1] + F(*z-(i-1-R)) + F(*z-(i-R));

    /* print out field data */
    sprintf (filename, "field-%004d.txt\0", (int) floor(delta*n));
    file_field = fopen (filename, "w");
    for (i=1; i<N-1; i++) {
      E = 0.5*(phi[i+1] - phi[i-1]);
      B = pi[i];
      Erad = E - Es[i];
      fprintf (file_field, "%d, %f, %f, %f, %f, %f, %f\n",
	       i-R, E, B, phi[i], 0.5*(Erad*Erad + B*B), -Erad*B, t);
    }
    fclose (file_field);
    
    /* print out particle data */
    if (n%5 == 0) fprintf (file_particle, "%f, %f, %f\n", t, *p/M, *z);

    evolve_state (f, s_ptr, delta);
  }
  fclose (file_particle);
}

/****************************************************************************
 * paper figures
 ****************************************************************************/

void figure_2 () {
  union s_type state;

  M = 2000.0;
  omega0 = 1/50.0;
  initialize_state (&state);
  state.P = 0.3*M;

  make_frames (eqns_harmonic, &state, 5000.0);
}

void figure_3 () {
  union s_type state;

  M = 500.0;
  omega0 = 0.0;
  initialize_state (&state);
  state.P = 0.3*M;

  make_frames (eqns_harmonic, &state, 1000.0);
}

void figure_6 () {
  union s_type state;

  M = 500.0;
  omega0 = 0.0;
  initialize_state (&state);
  state.P = 0.3*M;

  make_frames (eqns_harmonic, &state, 5000.0);
}

/****************************************************************************
 * time reversal experiments
 ****************************************************************************/

void time_reverse_sho () {
  union s_type state;
  int n;

  state.Z = 0.0;
  state.P = 1.0;

  for (n=0; n<20; n++) {
    printf ("%f, %f, %f\n", 0.1*n, state.Z, state.P);
    evolve_state (eqns_sho, &state, 0.1);
  }
  reverse_state (&state);
  for (; n<=40; n++) {
    printf ("%f, %f, %f\n", 0.1*n, state.Z, state.P);
    evolve_state (eqns_sho, &state, 0.1);
  }
}

void time_reverse_toymodel () {
  union s_type state;
  int n;

  for (n=0; n<200; n++) {
    printf ("%f, %f, %f\n", 1.0*n, state.Z, state.P/M);
    evolve_state (eqns_harmonic, &state, 1.0);
  }
  reverse_state (&state);
  for (; n<400.0; n++) {
    printf ("%f, %f, %f\n", 1.0*n, state.Z, state.P/M);
    evolve_state (eqns_harmonic, &state, 1.0);
  }
}

/****************************************************************************
 * main
 ****************************************************************************/

int main () {
  figure_6 ();
}