/*
 *  Groups of atoms, and interesting operations on them
 */

#include 
#include 

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

void
init_group (group * g)
{
  g->atom_count = g->bond_count = g->term_count = 0;
}

/* Everything that follows assumes that atoms, bonds, and terms are 
 * allocated in static arrays. Some day I'll do something more clever
 * than that, and all this code will need to be rewritten.
 */

void
add_atom (group * g, atom * a)
{
  if (g->atom_count < MAX_NUM_ATOMS)
    g->atoms[g->atom_count++] = a;
  else
    fprintf (stderr, "Too many atoms\n");
}

void
add_bond (group * g, bond * b)
{
  if (g->bond_count < MAX_NUM_BONDS)
    g->bonds[g->bond_count++] = b;
  else
    fprintf (stderr, "Too many bonds\n");
}

void
add_term (group * g, term * t)
{
  if (g->term_count < MAX_NUM_TERMS)
    g->terms[g->term_count++] = t;
  else
    fprintf (stderr, "Too many terms\n");
}

atom *
choose_atom (group * g, int index)
{
  return g->atoms[index];
}

bond *
choose_bond (group * g, int index)
{
  return g->bonds[index];
}

term *
choose_term (group * g, int index)
{
  return g->terms[index];
}

double
energy (group * g)
{
  double e;
  int i;

  for (i = 0, e = 0.0; i < g->atom_count; i++)
    e += kinetic_energy (g->atoms[i]);
  for (i = 0; i < g->term_count; i++)
    e += term_energy (g->terms[i]);
  return e;
}

void
physics_time_step (group * g, double dt)
{
  int i;

  for (i = 0; i < g->atom_count; i++)
    {
      iterate_motion (g->atoms[i], dt);
      g->atoms[i]->f[0] = g->atoms[i]->f[1] = g->atoms[i]->f[2] = 0.0;
    }
  for (i = 0; i < g->term_count; i++)
    {
      term_add_forces (g->terms[i]);
    }
}

void
select_subgroup (group * g, group * selected,
		 int (*selector) (atom * atm, int index))
{
  int i;

  for (i = 0; i < g->atom_count; i++)
    if ((*selector) (g->atoms[i], i))
      add_atom (selected, g->atoms[i]);
}

void
apply_to_atoms (group * g, int (*do_this) (atom * a))
{
  int i;

  for (i = 0; i < g->atom_count; i++)
    (*do_this) (g->atoms[i]);
}

void
apply_to_terms (group * g, int (*do_this) (term * t))
{
  int i;

  for (i = 0; i < g->term_count; i++)
    (*do_this) (g->terms[i]);
}

#define THRESH  1.8		/* maximum bond length in angstroms */
#define THRESH_SQ  (THRESH * THRESH)

void
infer_bonds_from_distances (group * g)
{
  int i, j;
  atom *p, *q;
  double diff, dsq;

  for (i = 0; i < g->bond_count; i++)
    free (g->bonds[i]);
  g->bond_count = 0;

  for (i = 0; i < g->bond_count; i++)
    {
      p = g->atoms[i];
      for (j = i + 1; j < g->bond_count; j++)
	{
	  q = g->atoms[j];
	  diff = p->x[0] - q->x[0];
	  dsq = diff * diff;
	  diff = p->x[1] - q->x[1];
	  dsq += diff * diff;
	  diff = p->x[2] - q->x[2];
	  dsq += diff * diff;
	  if (dsq <= THRESH_SQ)
	    {
	      bond *b;
	      b = malloc (sizeof (bond));
	      b->first = i;
	      b->second = j;
	      g->bonds[g->bond_count++] = b;
	    }
	}
    }
}