/* Here's an example where we link the libmol.a stuff to a little interpreter
 * of sorts. This understands the 'script' file.
 */

#include <iostream.h>

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

enum WHAT_TO_PRINT
{
  /* things that can be printed */
  TIME, TOTAL_ENERGY, ENERGY, POSITION, VELOCITY, FORCE, XYZ, NEWLINE
};

struct to_print
{
  enum WHAT_TO_PRINT token;
  int which_atom;
};

#define MAX_PRINTABLES 200
int num_printables = 0;
struct to_print printables[MAX_PRINTABLES];

group g;

int every_modulo = 1;
double dt, t;

enum SPECIES lookup_species (char *w)
{
  if (strcmp (w, "c-sp3") == 0) return C_SP3;
  if (strcmp (w, "c-sp2") == 0) return C_SP2;
  if (strcmp (w, "c-sp" ) == 0) return C_SP;
  if (strcmp (w, "h"    ) == 0) return H;
  if (strcmp (w, "n-sp3") == 0) return N_SP3;
  if (strcmp (w, "n-sp2") == 0) return N_SP2;
  if (strcmp (w, "n-sp" ) == 0) return N_SP;
  if (strcmp (w, "o-sp3") == 0) return O_SP3;
  if (strcmp (w, "o-sp2") == 0) return O_SP2;
  return C_SP3;
}

char *
element_name (enum SPECIES s)
{
  switch (s)
    {
    case C_SP3:
    case C_SP2:
    case C_SP:
      return "C";
    case N_SP3:
    case N_SP2:
    case N_SP:
      return "N";
    case O_SP3:
    case O_SP2:
      return "O";
    case H:
      return "H";
    }
  return "?";
}

/* The idea here is to print the kinds of XYZ frames that
 * Al Globus uses to make animations.
 */
void
print_xyz (void)
{
  int i;
  atom *a;

  cout << g.atom_count << endl;
  for (i = 0; i < g.atom_count; i++)
    {
      a = g.choose_atom (i);
      cout << element_name (a->species) << " "
	   << a->x[0] << " "
	   << a->x[1] << " "
	   << a->x[2] << endl;
    }
}

void printout (void)
{
  int i;
  atom *a;

  for (i = 0; i < num_printables; i++)
    {
      a = g.choose_atom (printables[i].which_atom);
      switch (printables[i].token)
	{
	case TIME:
	  cout << t << " ";
	  break;
	case TOTAL_ENERGY:
	  cout << g.energy() << " ";
	  break;
	case ENERGY:
	  cout << a->kinetic_energy() << " ";
	  break;
	case POSITION:
	  cout << "["
	       << a->x[0] << " "
	       << a->x[1] << " "
	       << a->x[2] << "] ";
	  break;
	case VELOCITY:
	  cout << "["
	       << a->v[0] << " "
	       << a->v[1] << " "
	       << a->v[2] << "] ";
	  break;
	case FORCE:
	  cout << "["
	       << a->f[0] << " "
	       << a->f[1] << " "
	       << a->f[2] << "] ";
	  break;
	case XYZ:
	  print_xyz ();
	  break;
	case NEWLINE:
	  cout << endl;
	  break;
	}
    }
}

int main (void)
{
  char word[200];

  while (cin >> word)
    {
      if (word[0] == '#')
	{
	  cin.getline (word, 200, '\n');
	}
      else if (strcmp (word, "atom") == 0)
	{
	  enum SPECIES s;
	  double ch, x, y, z;

	  cin >> word;
	  s = lookup_species (word);
	  cin >> ch >> x >> y >> z;
	  g.add_atom (new atom (s, ch, x, y, z));
	}
      else if (strcmp (word, "bond") == 0)
	{
	  int f, s;
	  cin >> f >> s;
	  g.add_bond (new bond(f, s));
	}
      else if (strcmp (word, "length") == 0)
	{
	  term *t;
	  int a1, a2;
	  cin >> a1 >> a2;
	  t = length_term (g.choose_atom (a1), g.choose_atom (a2));
	  g.add_term (t);
	}
      else if (strcmp (word, "angle") == 0)
	{
	  term *t;
	  int a1, a2, a3;
	  cin >> a1 >> a2 >> a3;
	  t = angle_term (g.choose_atom (a1), g.choose_atom (a2),
			  g.choose_atom (a3));
	  g.add_term (t);
	}
      else if (strcmp (word, "torsion") == 0)
	{
	  term *t;
	  int a1, a2, a3, a4;
	  cin >> a1 >> a2 >> a3 >> a4;
	  t = torsion_term (g.choose_atom (a1), g.choose_atom (a2),
			    g.choose_atom (a3), g.choose_atom (a4));
	  g.add_term (t);
	}
      else if (strcmp (word, "vdw") == 0)
	{
	  term *t;
	  int a1, a2;
	  cin >> a1 >> a2;
	  t = vdw_term (g.choose_atom (a1), g.choose_atom (a2));
	  g.add_term (t);
	}
      else if (strcmp (word, "electrostatic") == 0)
	{
	  term *t;
	  int a1, a2;
	  cin >> a1 >> a2;
	  t = electrostatic_term (g.choose_atom (a1), g.choose_atom (a2));
	  g.add_term (t);
	}
      else if (strcmp (word, "spring-damper") == 0)
	{
	  term *t;
	  int a1, a2;
	  double k, r0, d;
	  cin >> a1 >> a2 >> k >> r0 >> d;
	  t = spring_damper_term (g.choose_atom (a1), g.choose_atom (a2),
				  k, r0, d);
	  g.add_term (t);
	}
      else if (strcmp (word, "time-step") == 0)
	{
	  cin >> dt;
	}
      else if (strcmp (word, "every") == 0)
	{
	  cin >> every_modulo;
	}
      else if (strcmp (word, "run") == 0)
	{
	  int n;
	  double t0, t1;
	  cin >> t0 >> t1;
	  t1 += 0.5 * dt;
	  for (n = 0, t = t0; t < t1; t += dt, n++)
	    {
	      if ((n % every_modulo) == 0)
		printout ();
	      g.physics_time_step (dt);
	    }
	}
      else if (strcmp (word, "newline") == 0)
	{
	  printables[num_printables++].token = NEWLINE;
	}
      else if (strcmp (word, "print") == 0)
	{
	  struct to_print *p = &printables[num_printables++];
	  cin >> word;
	  if      (strcmp (word, "time") == 0)
	    {
	      p->token = TIME;
	    }
	  else if (strcmp (word, "total-energy") == 0)
	    {
	      p->token = TOTAL_ENERGY;
	    }
	  else if (strcmp (word, "energy") == 0)
	    {
	      p->token = ENERGY;
	      cin >> p->which_atom;
	    }
	  else if (strcmp (word, "position") == 0)
	    {
	      p->token = POSITION;
	      cin >> p->which_atom;
	    }
	  else if (strcmp (word, "velocity") == 0)
	    {
	      p->token = VELOCITY;
	      cin >> p->which_atom;
	    }
	  else if (strcmp (word, "force") == 0)
	    {
	      p->token = FORCE;
	      cin >> p->which_atom;
	    }
	  else if (strcmp (word, "xyz") == 0)
	    {
	      p->token = XYZ;
	    }
	  else if (strcmp (word, "newline") == 0)
	    {
	      p->token = NEWLINE;
	    }
	  else
	    {
	      cerr << "Bad print statement: " << word << endl;
	      num_printables--;
	    }
	}
      else cerr << "Bad command: " << word << endl;
    }
}