/****************************************************************************
* electrodynamics-scattering.c
* David Boozer
* 23 June 2008
* 30 July 2008 (modified)
*****************************************************************************
* Simulate scattering in 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 100
#define N (2*R + 1)
#define NUM_DOF (2*N+3)

#define SIGMA 5.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, t;
  } name;
  double dof[NUM_DOF];
};

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

double M, omega0, A, omega;

/****************************************************************************
 * 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 initialize_state (union s_type *s_ptr) {
  double temp;
  int i, j;

  s_ptr->Z = 0.0;
  s_ptr->P = 0.0;
  s_ptr->T = 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 harmonically bound particle, incoming wave */
void eqns_scattering (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 t = s_ptr->T;
  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 +
			    2*A*sin(omega*t));
    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;

  s_dot_ptr->T = 1.0;
}


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

void make_frames (void *f, union s_type *s_ptr, double tau) {
  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<=tau; 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
 ****************************************************************************/

double energy_harmonic (union s_type *s_ptr) {
  double p = s_ptr->P;
  double z = s_ptr->Z;

  return (p*p/M + M*omega0*omega0*z*z)/(M*A*A);
}

void figure_4a () {
  union s_type state;
  double tau, delta, detuning, t;
  int n;

  M = 5000.0;
  omega0 = 1/500.0;
  A = 0.1;
  detuning = 0.0;
  omega = omega0 + detuning/M;

  initialize_state (&state);

  tau = 50000.0;
  delta = 10.0;
  for (t=0.0, n=0; t<=tau; t+=delta, n++) {
    printf ("%f, %f\n", t, energy_harmonic (&state));
    evolve_state (eqns_scattering, &state, delta);
  }
}

void figure_4b () {
  union s_type state;
  double detuning, tau;

  M = 5000.0;
  omega0 = 1/500.0;
  A = 0.1;

  tau = 50000.0;
  for (detuning = -5.0; detuning <= 5.0; detuning += 0.2) {
    initialize_state (&state);
    omega = omega0 + detuning/M;
    evolve_state (eqns_scattering, &state, tau);
    printf ("%f, %f\n", detuning, energy_harmonic (&state));
    fprintf (stderr, "%f, %f\n", detuning, energy_harmonic (&state));
  }
}

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

int main () {
  union s_type state;
  double tau, delta, detuning;

  /*
  M = 5000.0;
  omega0 = 1/500.0;
  A = 0.1;
  detuning = 0.0;
  omega = omega0 + detuning/M;

  initialize_state (&state);

  make_frames (eqns_scattering, &state, 50000.0);
  */

  figure_4b ();
}