/*
Copyright (C) 2003 by Sean David Fleming

sean@ivec.org

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

The GNU GPL can also be found at http://www.gnu.org
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#ifndef __WIN32
#include <sys/times.h>
#endif

#include "gdis.h"
#include "coords.h"
#include "edit.h"
#include "matrix.h"
#include "opengl.h"
#include "render.h"
#include "select.h"
#include "shortcuts.h"
#include "interface.h"
#include "zone.h"

#define DEBUG 0

extern struct sysenv_pak sysenv;
extern struct elem_pak elements[];

gint current_colour[3];


/*********************/
/* copy an atom core */
/*********************/
/* return the number of items copied (ie 0-2) */
#define DEBUG_COPY_CORE 0
struct core_pak *copy_core(struct core_pak *core, struct model_pak *src,
                                                  struct model_pak *dest)
{
gint items=0;
gdouble vec[3];
struct core_pak *copyc;
struct shel_pak *copys;

/* checks */
g_assert(core != NULL);
g_assert(src != NULL);
g_assert(dest != NULL);

/* duplicate data structure */
copyc = dup_core(core);
items++;

/* setup status */
copyc->status = copyc->status & (~SELECT & ~SELECT);
copyc->orig = copyc->primary = TRUE;
copyc->primary_core = NULL;
VEC3SET(copyc->offset, 0.0, 0.0, 0.0);

/* coords, account for transformation matrices */
ARR3SET(vec, core->rx);
vecmat(dest->ilatmat, vec);
ARR3ADD(vec, dest->centroid);
ARR3SET(copyc->x, vec);

dest->cores = g_slist_prepend(dest->cores, copyc);

/* attached shell? */
if (copyc->shell)
  {
  copys = copyc->shell;
  items++;

/* main info */
  copys->status = copys->status & (~SELECT);
  copys->primary=copys->orig=TRUE; 
  copys->primary_shell = NULL;
  VEC3SET(copys->offset, 0.0, 0.0, 0.0);

/* coords, account for transformation matrices */
  ARR3SET(vec, copys->rx);
  vecmat(dest->ilatmat, vec);
  ARR3ADD(vec, dest->centroid);
  ARR3SET(copys->x, vec);

  dest->shels = g_slist_prepend(dest->shels, copys);
  }

return(copyc);
}

/********************/
/* duplicate a core */
/********************/
struct core_pak *dup_core(struct core_pak *orig)
{
struct core_pak *core;

/* checks */
g_assert(orig != NULL);

core = g_malloc(sizeof(struct core_pak));

memcpy(core, orig, sizeof(struct core_pak));

/* clear pointers */
core->bonds = NULL;
core->shell = NULL;
core->mol = NULL;
core->primary_core = NULL;

/* duplicate strings */
core->atom_label = g_strdup(orig->atom_label);
core->atom_type = g_strdup(orig->atom_type);
core->res_name = g_strdup(orig->res_name);
core->flags = g_strdup(orig->flags);

core->vibx_list = NULL;
core->viby_list = NULL;
core->vibz_list = NULL;

/* duplicate the shell */
if (orig->shell)
  {
  core->shell = dup_shel(orig->shell);
  (core->shell)->core = core;
  }

return(core);
}

/*********************/
/* shell duplication */
/*********************/
struct shel_pak *dup_shel(struct shel_pak *orig)
{
struct shel_pak *shel;

/* checks */
g_assert(orig != NULL);

shel = g_malloc(sizeof(struct shel_pak));

memcpy(shel, orig, sizeof(struct shel_pak));

/* clear pointers */
shel->core = NULL;
shel->primary_shell = NULL;
VEC3SET(shel->pic, 0, 0, 0);

/* duplicate strings */
shel->flags = g_strdup(orig->flags);

return(shel);
}

/****************************/
/* base atom initialization */
/****************************/
void core_init(gchar *elem, struct core_pak *core, struct model_pak *model)
{
gint code;
struct elem_pak elem_data;

/* attempt to match atom type with database */
if (model->protein)
  code = pdb_elem_type(elem);
else
  code = elem_test(elem);

/* general initialization */
core->atom_code = code;
core->atom_type = NULL;
core->atom_label = NULL;
core->res_name = g_strdup("");
core->res_no = 1;
core->chain = NULL;
core->atom_order = 0;

/* element related initialization */
get_elem_data(code, &elem_data, model);
core->bond_cutoff = elem_data.cova;
ARR3SET(core->colour, elem_data.colour);
core->colour[3] = 1.0;
}

/***************************/
/* atom creation primitive */
/***************************/
gpointer core_new(gchar *elem, gchar *label, struct model_pak *model)
{
struct core_pak *core;

g_assert(elem != NULL);
g_assert(model != NULL);

core = g_malloc(sizeof(struct core_pak));

/* element related initialization */
core_init(elem, core, model);

/* general initialization */
if (label)
  core->atom_label = g_strdup(label);
else
  core->atom_label = g_strdup(elem);

VEC4SET(core->x, 0.0, 0.0, 0.0, 1.0);
VEC4SET(core->rx, 0.0, 0.0, 0.0, 1.0);
VEC3SET(core->v, 0.0, 0.0, 0.0);
VEC3SET(core->offset, 0.0, 0.0, 0.0);
core->status = NORMAL;
core->primary = TRUE;
core->orig = TRUE;
core->ghost = FALSE;
core->breathe = FALSE;
core->growth = FALSE;
core->translate = FALSE;
core->render_mode = BALL_STICK;
core->render_wire = FALSE;
core->region = REGION1A;
core->radius = 0.0;
core->has_sof = FALSE;
core->sof = 1.0;
core->charge = 0.0;
core->lookup_charge = TRUE;
core->flags = NULL;
core->bonds = NULL;
core->mol = NULL;
core->shell = NULL;
core->primary_core = NULL;
core->vibx_list = NULL;
core->viby_list = NULL;
core->vibz_list = NULL;

/* NEW: hydrogen bond capabilities */
if (g_strrstr(elem, "H") != NULL)
  core->hydrogen_bond = TRUE;
else 
  core->hydrogen_bond = FALSE;

return(core);
}

/*******************************/
/* unified atom initialization */
/*******************************/
/* deprec */
#define DEBUG_NEW_CORE 0
struct core_pak *new_core(gchar *elem, struct model_pak *model)
{
struct core_pak *core;

core = core_new(elem, NULL, model);

return(core);
}

/****************************/
/* shell creation primitive */
/****************************/
gpointer shell_new(gchar *elem, gchar *label, struct model_pak *model)
{
gint code;
struct elem_pak elem_data;
struct shel_pak *shell;

g_assert(elem != NULL);
g_assert(model != NULL);

shell = g_malloc(sizeof(struct shel_pak));

/* attempt to match atom type with database */
code = elem_test(elem);
if (!code)
  printf("Warning: element [%s] not found.\n", elem);

/* init modifiable element data */
get_elem_data(code, &elem_data, model);
ARR3SET(shell->colour, elem_data.colour);

shell->atom_code = code;
if (label)
  shell->shell_label = g_strdup(label);
else
  shell->shell_label = g_strdup(elem);
VEC4SET(shell->x, 0.0, 0.0, 0.0, 1.0);
VEC4SET(shell->rx, 0.0, 0.0, 0.0, 1.0);
VEC3SET(shell->v, 0.0, 0.0, 0.0);
VEC3SET(shell->offset, 0.0, 0.0, 0.0);
shell->status = NORMAL;
shell->primary = TRUE;
shell->orig = TRUE;
shell->breathe = FALSE;
shell->translate = FALSE;
shell->region = REGION1A;
shell->core = NULL;
shell->primary_shell = NULL;
shell->radius = 0.0;
shell->charge = 0.0;
shell->has_sof = FALSE;
shell->sof = 1.0;
shell->lookup_charge = TRUE;
shell->flags = NULL;

return(shell);
}

/**********************************/
/* replacement shell init routine */
/**********************************/
/* deprec */
struct shel_pak *new_shell(gchar *elem, struct model_pak *model)
{
struct shel_pak *shell;

shell = shell_new(elem, NULL, model);

return(shell);
}

/***************************************/
/* fully duplicate a model's core list */
/***************************************/
GSList *dup_core_list(GSList *orig)
{
GSList *copy, *list;
struct core_pak *core;

copy = NULL;
for (list=orig ; list ; list=g_slist_next(list))
  {
  core = dup_core(list->data);
  copy = g_slist_prepend(copy, core);
  }
copy = g_slist_reverse(copy);

return(copy);
}

/****************************************/
/* fully duplicate a model's shell list */
/****************************************/
GSList *dup_shell_list(GSList *orig)
{
GSList *copy, *list;
struct shel_pak *shell;

copy = NULL;
for (list=orig ; list ; list=g_slist_next(list))
  {
  shell = dup_shel(list->data);
  copy = g_slist_prepend(copy, shell);
  }
copy = g_slist_reverse(copy);

return(copy);
}

/**************************/
/* remove duplicate cores */
/**************************/
/* NB: assumes zone_init() has been done */
#define DEBUG_REMOVE_DUPLICATES 0
#define DEBUG_REMOVE_DUPLICATES_MORE 0
void delete_duplicate_cores(struct model_pak *model)
{
gdouble vec[3];
GSList *list1, *list2;
gpointer zone;
struct core_pak *core1, *core2;
struct shel_pak *s1, *s2;

/* checks */
g_assert(model != NULL);
g_assert(model->zone_array != NULL);

#if DEBUG_REMOVE_DUPLICATES
printf("Initial cores: %d\n", g_slist_length(model->cores));
printf("Initial shels: %d\n", g_slist_length(model->shels));
#endif

/* enumerate all cores */
for (list1=model->cores ; list1 ; list1=g_slist_next(list1))
  {
  core1 = list1->data;
  if (core1->status & DELETED)
    continue;

/* enumerate cores in current locality */
  zone = zone_get(core1->x, model->zone_array);


#if DEBUG_REMOVE_DUPLICATES_MORE
printf(" + %s [%p] : [%p] :", core1->atom_label, core1, zone);
P3VEC(" ", core1->x);
#endif
/* should use zone_area_cores() as a very small coord difference */
/* can result in cores being put in different (neighbouring) zones */
/*
  for (list2=zone_cores(zone) ; list2 ; list2=g_slist_next(list2))
*/


  for (list2=zone_area_cores(1, zone, model->zone_array) ; list2 ; list2=g_slist_next(list2))

    {
    core2 = list2->data;
    if (core2->status & DELETED)
      continue;

/* avoid double counting */
    if (core1 >= core2)
      continue;

    if (core1->atom_code != core2->atom_code)
      continue;

#if DEBUG_REMOVE_DUPLICATES_MORE
printf(" - %s :", core2->atom_label);
P3VEC(" ", core2->x);
#endif

/* compute and test the minimum separation */
    ARR3SET(vec, core1->x);
    ARR3SUB(vec, core2->x);

    fractional_minsq(vec, model->periodic);

    if (VEC3MAGSQ(vec) < FRACTION_TOLERANCE)
      {
/* delete core2, unless primary AND core1 is non primary */
      if (core2->primary && !core1->primary)
        {
        core1->status |= DELETED;

/* FIXME - this may be a problem since we start the core2 loop */
/* requiring that core1 be undeleted */
#if DEBUG_REMOVE_DUPLICATES_MORE
printf(" * rm 1\n");
#endif
        }
      else
        {
        core2->status |= DELETED;

#if DEBUG_REMOVE_DUPLICATES_MORE
printf(" * rm 2\n");
#endif
        }
      }
    }
  }

/* commit before searching for duplicate shells, as a commit */
/* will delete some attached shells */
delete_commit(model);

for (list1=model->shels ; list1 ; list1=g_slist_next(list1))
  {
  s1 = list1->data;
  if (s1->status & DELETED)
    continue;

/* NEW - enumerate shells in the current locality */
  zone = zone_get(s1->x, model->zone_array); 
  for (list2=zone_shells(zone) ; list2 ; list2=g_slist_next(list2))
    {
    s2 = list2->data;

    if (s2->status & DELETED)
      continue;
    if (s1 == s2)
      continue;

    ARR3SET(vec, s1->x);
    ARR3SUB(vec, s2->x);

/* adjust for periodicity */
    fractional_minsq(vec, model->periodic);

    if (VEC3MAGSQ(vec) < FRACTION_TOLERANCE)
      {
/* delete shell2, unless primary AND shell1 is non primary */
      if (s2->primary && !s1->primary)
        s1->status |= DELETED;
      else
        s2->status |= DELETED;
      }
    }
  }
delete_commit(model);

#if DEBUG_REMOVE_DUPLICATES
printf("Final cores: %d\n", g_slist_length(model->cores));
printf("Final shels: %d\n", g_slist_length(model->shels));
#endif
}

/********************************************************/
/* set vector to minimum fractional magnitute (squared) */
/********************************************************/
#define DEBUG_FRAC_MINSQ 0
void fractional_minsq(gdouble *x, gint dim)
{
gint i, j;
gdouble whole;

g_assert(dim < 4);

#if DEBUG_FRAC_MINSQ
P3VEC("b4: ", x);
#endif

for (i=0 ; i<dim ; i++)
  {
/* clamp value -1 < x < 1 */ 
  x[i] = modf(x[i], &whole);

/* clamp to range (-0.5, 0.5) */
  j = (gint) (-2.0 * x[i]);
  x[i] += (gdouble) j;
  }

#if DEBUG_FRAC_MINSQ
P3VEC("af: ", x);
#endif
}

/*************************************/
/* clamp coords to fractional limits */
/*************************************/
#define DEBUG_CLAMP 0
void fractional_clamp(gdouble *vec, gint *mov, gint dim)
{
gint i;
gdouble ip;

#if DEBUG_CLAMP
P3VEC("inp: ", vec);
#endif

/* NB: init ALL 3 coords as isolated molecules call this */
/* routine & it is expected that mov to be set to 0,0,0 */
mov[0] = mov[1] = mov[2] = 0;

/*
for (i=0 ; i<dim ; i++)
*/
for (i=dim ; i-- ; )
  {
  mov[i] = -vec[i];
  if (vec[i] < 0.0)
    {
/* increment the move for -ve's with exception -1.0, -2.0, -3.0,... */
    if (modf(vec[i], &ip) != 0.0)
      mov[i]++;
    }
  else
    mov[i] = -vec[i];

  vec[i] += (gdouble) mov[i];
  }

#if DEBUG_CLAMP
P3VEC("out: ", vec);
#endif
}

/********************************************/
/* confine a list of atoms to the unit cell */
/********************************************/
#define DEBUG_PBC_ATOMS 0
void coords_confine_cores(GSList *cores, struct model_pak *model)
{
gint xlat[3];
gdouble mov[3];
GSList *list;
struct core_pak *core;
struct shel_pak *shell;

g_assert(model != NULL);

if (!model->periodic)
  return;

/* translate cores to within the cell */
for (list=cores ; list ; list=g_slist_next(list))
  {
  core = list->data;
  fractional_clamp(core->x, xlat, model->periodic);

/* move shell */
  if (core->shell)
    {
    shell = core->shell;
    ARR3SET(mov, xlat);
    ARR3ADD(shell->x, mov);
    }
  }
}

/**********************************************/
/* confine a molecule's centroid to unit cell */
/**********************************************/
void coords_confine_centroid(struct mol_pak *mol, struct model_pak *model)
{
gint xlat[3];
gdouble mov[3];
GSList *list;
struct core_pak *core;

/* calc moves required to bring centroid within pbc */
fractional_clamp(mol->centroid, xlat, model->periodic);
ARR3SET(mov, xlat);

/* apply to all atoms/shells in this molecule */
for (list=mol->cores ; list ; list=g_slist_next(list))
  {
  core = list->data;
  ARR3ADD(core->x, mov);
  if (core->shell)
    {
    ARR3ADD((core->shell)->x, mov);
    }
  }
}

/*******************************************/
/* draw a box - contents are the selection */
/*******************************************/
void update_box(gint x, gint y, struct model_pak *data, gint call_type)
{

switch (call_type)
  {
  case START:
/* don't clean the selection (allows multiple box selections) */
/* setup the box selection object */
    data->box_on = TRUE;
    data->select_box[0] = x;
    data->select_box[1] = y;
  case UPDATE:
    data->select_box[2] = x;
    data->select_box[3] = y;
    break;
  }
}

/**********************************************************/
/* check if a given label matches with element/atom label */
/**********************************************************/
#define DEBUG_SHELL_MATCH 0
gint shel_match(const gchar *label, struct shel_pak *shell)
{
gint code;

#if DEBUG_SHELL_MATCH
printf("[%s] : [%s]\n", label, shell->shell_label);
#endif

code = elem_symbol_test(label);

/* if input label doesnt match the element symbol length - it means the */
/* user has put in something like H1 - compare this with the atom label */
if (code)
  {
  if (g_ascii_strcasecmp(label, elements[shell->atom_code].symbol) != 0)
    {
    if (g_ascii_strcasecmp(shell->shell_label, label) == 0)
      return(1);
    }
  else
    return(1);
  }
else
  return(1);

#if DEBUG_SHELL_MATCH
printf("rejected.\n");
#endif

return(0);
}

/**********************************************************/
/* check if a given label matches with element/atom label */
/**********************************************************/
#define DEBUG_ATOM_MATCH 0
gint core_match(const gchar *label, struct core_pak *core)
{
gint code;

#if DEBUG_ATOM_MATCH
printf("[%s] : [%s]\n", label, core->atom_label);
#endif

code = elem_symbol_test(label);

/* if input label doesnt match the element symbol length - it means the */
/* user has put in something like H1 - compare this with the atom label */
if (code)
  {
  if (g_ascii_strcasecmp(label, elements[core->atom_code].symbol) != 0)
    {
    if (g_ascii_strcasecmp(core->atom_label, label) == 0)
      return(1);
    }
  else
    return(1);
  }
else
  return(1);

#if DEBUG_ATOM_MATCH
printf("rejected.\n");
#endif

return(0);
}