/*
 * Code to test atom.C
 */

#include <iostream.h>

#include "atom.h"
#include "bond.h"
#include "term.h"
#include "group.h"

#define X 0.076
#define DT 1.0E-6

group G;

int
main (int argc, char *argv[])
{
  double t;
  int i;
  atom *a1, *a2;
  term *t1, *t2, *t3;
  
  if (argc > 2)
    {
      a1 = new atom (C_SP3, 1.0, -X, 0.0, 0.0);
      a2 = new atom (C_SP3, 1.0, X, 0.0, 0.0);
    }
  else
    {
      a1 = new atom (C_SP3, 0.0, -X, 0.0, 0.0);
      a2 = new atom (C_SP3, 0.0, X, 0.0, 0.0);
    }
  t1 = length_term (a1, a2);
  t2 = spring_damper_term (a1, a2, 2.0e4, 0.3, 100.0);
  t3 = electrostatic_term (a1, a2);
  
  G.add_atom (a1);
  G.add_atom (a2);
  G.add_term (t1);
  G.add_term (t2);
  G.add_term (t3);
  
  for (t = 0.0, i = 0; t < 1.0 + 0.1 * DT; t += DT, i++)
    {
      G.physics_time_step (DT);
      /* if you want to apply external forces, do that now */
      
      if (argc > 1)
	if ((i % ((int) (0.01 / DT))) == 0)
	  cout << t << " "
	       << a1->x[0] << " "
	       << a2->x[0] << " "
	       << G.energy () << endl;
    }
}