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

#include <iostream.h>

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

group::group ()
{
  atom_count = bond_count = 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 group::
add_atom (atom * a)
{
  if (atom_count < MAX_NUM_ATOMS)
    atoms[atom_count++] = a;
  else
    cerr << "Too many atoms" << endl;
}

void group::
add_bond (bond * b)
{
  if (bond_count < MAX_NUM_BONDS)
    bonds[bond_count++] = b;
  else
    cerr << "Too many bonds" << endl;
}

void group::
add_term (term * t)
{
  if (term_count < MAX_NUM_TERMS)
    terms[term_count++] = t;
  else
    cerr << "Too many terms" << endl;
}

atom *group::
choose_atom (int index)
{
  return atoms[index];
}

bond *group::
choose_bond (int index)
{
  return bonds[index];
}

term *group::
choose_term (int index)
{
  return terms[index];
}

double group::
energy ()
{
  double e;
  int i;
  
  for (i = 0, e = 0.0; i < atom_count; i++)
    e += atoms[i]->kinetic_energy ();
  for (i = 0; i < term_count; i++)
    e += terms[i]->energy_function ();
  return e;
}

void group::
physics_time_step (double dt)
{
  int i;
  
  for (i = 0; i < atom_count; i++)
    {
      atoms[i]->iterate_motion (dt);
      atoms[i]->f[0] = atoms[i]->f[1] = atoms[i]->f[2] = 0.0;
    }
  for (i = 0; i < term_count; i++)
    {
      terms[i]->add_forces ();
    }
}

void group::
select_subgroup (group * selected,
		 int (*selector) (atom * atm, int index))
{
  int i;
  
  for (i = 0; i < atom_count; i++)
    if ((*selector) (atoms[i], i))
      selected->add_atom (atoms[i]);
}

void group::
apply_to_atoms (int (*do_this) (atom * a))
{
  int i;
  
  for (i = 0; i < atom_count; i++)
    (*do_this) (atoms[i]);
}

void group::
apply_to_terms (int (*do_this) (term * t))
{
  int i;
  
  for (i = 0; i < term_count; i++)
    (*do_this) (terms[i]);
}

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

void group::
infer_bonds_from_distances ()
{
  int i, j;
  atom *p, *q;
  double diff, dsq;
  
  /* get rid of any existing bonds */
  for (i = 0; i < bond_count; i++)
    delete bonds[i];
  bond_count = 0;
  
  for (i = 0; i < bond_count; i++)
    {
      p = atoms[i];
      for (j = i + 1; j < bond_count; j++)
	{
	  q = 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)
	    bonds[bond_count++] = new bond (i, j);
	}
    }
}