/*
 * Hacking atoms
 * Units: mass in 10^-27 kg, distance in nanometers, time in picoseconds,
 * force in piconewtons, energy in milli-attojoules (10^-21 J)
 */

#include <stdio.h>

#include "debug.h"
#include "atom.h"

void
iterate_motion (atom * a, double time_step)
{
  double m, a0, a1, a2;

  m = 0.5 * time_step / a->mass;

  a0 = m * a->f[0];
  a1 = m * a->f[1];
  a2 = m * a->f[2];

  a->v[0] += a0;
  a->v[1] += a1;
  a->v[2] += a2;

  a->x[0] += a->v[0] * time_step;
  a->x[1] += a->v[1] * time_step;
  a->x[2] += a->v[2] * time_step;

  a->v[0] += a0;
  a->v[1] += a1;
  a->v[2] += a2;
}

/* it might seem silly, but it helps keep main() a little simpler */
void
zero_force (atom * a)
{
  a->f[0] = 0.0;
  a->f[1] = 0.0;
  a->f[2] = 0.0;
}

double
kinetic_energy (atom * a)
{
  double vsq, e;

  vsq = a->v[0] * a->v[0] + a->v[1] * a->v[1] +
    a->v[2] * a->v[2];
  e = 0.5 * a->mass * vsq;
  DBGPRINTF1 ("Kinetic energy: %f\n", e);
  return e;
}

static struct SPECIES_INFO
{
  enum SPECIES species;
  double mass, evdw, rvdw;	/* 10^-27 kg, maJ, 10^-10 m */
}
species_info[] =
{
  /*       10^-27 kg  10^-21 J  nm */
  {
    C_SP3, 19.925, 0.357, 0.190,
  }
  ,
  {
    C_SP2, 19.925, 0.357, 0.194,
  }
  ,
  {
    C_SP, 19.925, 0.357, 0.194,
  }
  ,
  {
    H, 1.674, 0.382, 0.150,
  }
  ,
  {
    N_SP3, 23.251, 0.447, 0.182,
  }
  ,
  {
    O_SP3, 26.565, 0.406, 0.174,
  }
  ,
  {
    O_SP2, 26.565, 0.536, 0.174
  }
};

static int num_species = sizeof (species_info) / sizeof (species_info[0]);

void
init_atom (atom * a, enum SPECIES species, double charge,
	   double x0, double x1, double x2)
{
  int i, found;

  a->species = species;
  a->charge = charge;
  for (i = found = 0; i < num_species && !found; i++)
    if (species_info[i].species == species)
      {
	a->mass = species_info[i].mass;
	a->evdw = species_info[i].evdw;
	a->rvdw = species_info[i].rvdw;
	found = 1;
      }

  a->x[0] = x0;
  a->x[1] = x1;
  a->x[2] = x2;
  a->v[0] = a->v[1] = a->v[2] = 0.0;
  a->f[0] = a->f[1] = a->f[2] = 0.0;
}