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

#include <stdio.h>
#include <malloc.h>
#include <string.h>

#include "debug.h"
#include "atom.h"
#include "bond.h"
#include "term.h"
#include "group.h"
#include "getword.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;

  printf ("%d\n", g->atom_count);
  for (i = 0; i < g->atom_count; i++)
    {
      a = choose_atom (g, i);
      printf ("%s ", element_name (a->species));
      printf ("%f %f %f\n", a->x[0], a->x[1], a->x[2]);
    }
}

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

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

int
main (void)
{
  g = malloc (sizeof (group));
  init_getword ();
  while (!getword ())
    {
      if (strcmp (word, "atom") == 0)
	{
	  atom *a;
	  enum SPECIES s;
	  double ch, x, y, z;

	  getword ();
	  s = lookup_species (word);
	  ch = getdouble ();
	  x = getdouble ();
	  y = getdouble ();
	  z = getdouble ();
	  DBGPRINTF1("New atom at %0X\n", (unsigned int) a);
	  a = malloc (sizeof (atom));
	  init_atom (a, s, ch, x, y, z);
	  add_atom (g, a);
	}
      else if (strcmp (word, "bond") == 0)
	{
	  bond *b;
	  int f, s;

	  f = getint ();
	  s = getint ();
	  b = malloc (sizeof (bond));
	  DBGPRINTF1("New bond at %0X\n", (unsigned int) b);
	  init_bond (b, f, s);
	  add_bond (g, b);
	}
      else if (strcmp (word, "length") == 0)
	{
	  term *t;
	  int a1, a2;

	  a1 = getint ();
	  a2 = getint ();
	  t = malloc (sizeof (term));
	  DBGPRINTF1("New length term at %0X\n", (unsigned int) t);
	  length_term (t, choose_atom (g, a1), choose_atom (g, a2));
	  add_term (g, t);
	}
      else if (strcmp (word, "angle") == 0)
	{
	  term *t;
	  int a1, a2, a3;

	  t = malloc (sizeof (term));
	  DBGPRINTF1("New angle term at %0X\n", (unsigned int) t);
	  a1 = getint ();
	  a2 = getint ();
	  a3 = getint ();
	  angle_term (t, choose_atom (g, a1), choose_atom (g, a2),
		      choose_atom (g, a3));
	  add_term (g, t);
	}
      else if (strcmp (word, "torsion") == 0)
	{
	  term *t;
	  int a1, a2, a3, a4;

	  t = malloc (sizeof (term));
	  DBGPRINTF1("New torsion term at %0X\n", (unsigned int) t);
	  a1 = getint ();
	  a2 = getint ();
	  a3 = getint ();
	  a4 = getint ();
	  torsion_term (t, choose_atom (g, a1), choose_atom (g, a2),
			choose_atom (g, a3), choose_atom (g, a4));
	  add_term (g, t);
	}
      else if (strcmp (word, "vdw") == 0)
	{
	  term *t;
	  int a1, a2;

	  t = malloc (sizeof (term));
	  DBGPRINTF1("New vdw term at %0X\n", (unsigned int) t);
	  a1 = getint ();
	  a2 = getint ();
	  vdw_term (t, choose_atom (g, a1), choose_atom (g, a2));
	  add_term (g, t);
	}
      else if (strcmp (word, "electrostatic") == 0)
	{
	  term *t;
	  int a1, a2;

	  t = malloc (sizeof (term));
	  DBGPRINTF1("New electrostatic term at %0X\n", (unsigned int) t);
	  a1 = getint ();
	  a2 = getint ();
	  electrostatic_term (t, choose_atom (g, a1), choose_atom (g, a2));
	  add_term (g, t);
	}
      else if (strcmp (word, "spring-damper") == 0)
	{
	  term *t;
	  int a1, a2;
	  double k, r0, d;

	  t = malloc (sizeof (term));
	  DBGPRINTF1("New spring-damper term at %0X\n", (unsigned int) t);
	  a1 = getint ();
	  a2 = getint ();
	  k = getdouble ();
	  r0 = getdouble ();
	  d = getdouble ();
	  spring_damper_term (t, choose_atom (g, a1), choose_atom (g, a2),
			      k, r0, d);
	  add_term (g, t);
	}
      else if (strcmp (word, "time-step") == 0)
	{
	  dt = getdouble ();
	}
      else if (strcmp (word, "every") == 0)
	{
	  every_modulo = getint ();
	}
      else if (strcmp (word, "run") == 0)
	{
	  int n;
	  double t0, t1;

	  t0 = getdouble ();
	  t1 = getdouble ();
	  t1 += 0.5 * dt;
	  for (n = 0, t = t0; t < t1; t += dt, n++)
	    {
	      if ((n % every_modulo) == 0)
		printout ();
	      physics_time_step (g, 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++];
	  getword ();
	  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;
	      p->which_atom = getint ();
	    }
	  else if (strcmp (word, "position") == 0)
	    {
	      p->token = POSITION;
	      p->which_atom = getint ();
	    }
	  else if (strcmp (word, "velocity") == 0)
	    {
	      p->token = VELOCITY;
	      p->which_atom = getint ();
	    }
	  else if (strcmp (word, "force") == 0)
	    {
	      p->token = FORCE;
	      p->which_atom = getint ();
	    }
	  else if (strcmp (word, "xyz") == 0)
	    {
	      p->token = XYZ;
	    }
	  else if (strcmp (word, "newline") == 0)
	    {
	      p->token = NEWLINE;
	    }
	  else
	    {
	      fprintf (stderr, "Bad print statement: %s\n", word);
	      num_printables--;
	    }
	}
      else
	fprintf (stderr, "Bad command: %s\n", word);
    }

  /* At this point, it'd be nice to deallocate all those
   * atoms, bonds, and terms.
   */

  return 0;
}