/*******************************************************************************
*                                                                              *
*                                   Viewmol                                    *
*                                                                              *
*                            M O D I F Y G E O . C                             *
*                                                                              *
*                 Copyright (c) Joerg-R. Hill, October 2003                    *
*                                                                              *
********************************************************************************
*
* $Id: modifygeo.c,v 1.4 2003/11/07 11:07:25 jrh Exp $
* $Log: modifygeo.c,v $
* Revision 1.4  2003/11/07 11:07:25  jrh
* Release 2.4
*
* Revision 1.3  2000/12/10 15:12:10  jrh
* Release 2.3
*
* Revision 1.2  1999/05/24 01:28:58  jrh
* Release 2.2.1
*
* Revision 1.1  1999/02/07 21:53:09  jrh
* Initial revision
*
*
*/
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<Xm/Xm.h>
#include<GL/gl.h>
#include "viewmol.h"
#include "dialog.h"

extern double angle(double, double, double, double, double, double,double,double,double);
extern double bondLength(struct ATOM *, int, int);
extern double bondAngle(struct ATOM *, int, int, int);
extern double torsionAngle(struct ATOM *, int, int, int, int);
extern int messgb(Widget, int, char *, struct PushButtonRow *, int);
extern void GetMessageBoxButton(Widget, XtPointer, caddr_t);
extern struct SAVE *saveGeometry(struct MOLECULE *);
extern void restoreGeometry(struct SAVE *, int);
extern int  calcInternal(struct MOLECULE *, int);
extern void makeConnectivity(struct MOLECULE *, int, int);
extern char *getStringResource(Widget, char *);
extern void vectorProduct(double *, double *, double *);                               
extern double scalarProduct(double *, double *); 
extern void *getmem(size_t, size_t);
extern void *expmem(void *, size_t, size_t);
extern void fremem(void **);
extern void redraw(int);
extern void setAtom(struct MOLECULE *);
extern void deleteMolecule(Widget, caddr_t, XmAnyCallbackStruct *);
extern void clearGeometry(Widget, caddr_t, XmAnyCallbackStruct *);
extern void setMenuItem(int, int);

int findBondPartners(int *, int, int);
int findAttachedAtoms(int *, int, int);
void changeBondLength(int, int, int *, double);
int  changeAngle(int, int, int *, double);
void makeRotationMatrix(double, double, double, double, double matrix[3][3]);
void addUndo(struct MOLECULE *);

extern struct MOLECULE *molecules;
extern struct WINDOW windows[];
extern struct UNDO *undo;
extern Widget topShell, undoButton;
extern double bndfac;
extern int debug, nUndo;
extern int showWarning;

void modifyGeometry(int which, char *text)
{
  static struct PushButtonRow buttons[] = {{"ok", GetMessageBoxButton,
                                            (XtPointer)0, NULL}};
  struct MOLECULE *mol;
  double factor, value, diff;
  int *move, *tmove, nmove, m, save;
  char *word;
  register int i;

  if (windows[VIEWER].set >= 0)
    mol=&molecules[windows[VIEWER].set];
  else
    mol=&molecules[0];
  move=(int *)getmem(mol->na, sizeof(int));
  move[0]=mol->internals[which].atoms[0];
  if (mol->atoms[move[0]].nbonds > 1)
    nmove=findAttachedAtoms(move, 1, mol->internals[which].atoms[1]);
  else
    nmove=1;
  if (nmove < 0 && (mol->internals[which].type == ANGLE || mol->internals[which].type == TORSION))
  {
    move[0]=mol->internals[which].atoms[1];
    nmove=findAttachedAtoms(move, 1, mol->internals[which].atoms[2]);
  }
  if (nmove < 0)
  {
    if (showWarning)
    {
      word=getStringResource(topShell, "notChangeable");           
      messgb(topShell, 1, word, buttons, 1);
    }
    fremem((void *)&move);
    return;
  }

  addUndo(mol);

  value=atof(text);
  switch (mol->internals[which].type)
  {
    case BONDLENGTH:  factor=1.0;
                      if (strstr(text, "pm")) factor=0.01;
                      if (strstr(text, "bohr") || strstr(text, "au")) factor=0.52917706;
                      changeBondLength(which, nmove, move, factor*value);
                      break;
    case ANGLE:       (void)changeAngle(which, nmove, move, value);
                      break;
    case TORSION:     diff=value-torsionAngle(mol->atoms,
                                              mol->internals[which].atoms[0],
                                              mol->internals[which].atoms[1],
                                              mol->internals[which].atoms[2],
                                              mol->internals[which].atoms[3]);
                      (void)changeAngle(which, nmove, move, value);
                      move[0]=mol->internals[which].atoms[1];
                      if ((nmove=findBondPartners(move, 1, mol->internals[which].atoms[0])) > 1)
                      {
                        tmove=(int *)getmem(mol->na, sizeof(int));
                        save=mol->internals[which].atoms[0];
                        for (i=1; i<nmove; i++)
                        {
                          if (move[i] != mol->internals[which].atoms[2])
                          {
                            tmove[0]=move[i];
                            if ((m=findAttachedAtoms(tmove, 1, move[0])) > 0)
                            {
                              value=torsionAngle(mol->atoms, move[i], move[0],
                                                 mol->internals[which].atoms[2],
                                                 mol->internals[which].atoms[3])+diff;
                              mol->internals[which].atoms[0]=move[i];
                              (void)changeAngle(which, m, tmove, value);
                            }
                            }
                        }
                        mol->internals[which].atoms[0]=save;
                        fremem((void *)&tmove);
                      }
                      break;
  }
  fremem((void *)&move);
  makeConnectivity(mol, FALSE, TRUE);
  for (i=0; i<mol->ninternal; i++)
    (void)calcInternal(mol, i);
}

int findBondPartners(int *move, int n, int except)
{
  struct MOLECULE *mol;
  int found;
  register int i, j, k, l;

  /* Walk through list of bonds and add all atoms connected to
     the atoms in move[] to move[], ignore the bond between
     move[0] and except and all hydrogen bonds */

  if (windows[VIEWER].set >= 0)
    mol=&molecules[windows[VIEWER].set];
  else
    mol=&molecules[0];
  l=n;
  for (i=0; i<mol->nb; i++)
  {
    for (j=0; j<n; j++)
    {
      if (mol->bonds[i].order != (-1)) /* ignore hydrogen bonds */
      {
        if ((mol->bonds[i].first == move[0] &&
             mol->bonds[i].second == except) ||
            (mol->bonds[i].first == except &&
             mol->bonds[i].second == move[0])) continue;
        if (mol->bonds[i].first == move[j])
        {
          found=FALSE;
          for (k=0; k<l; k++)
            if (mol->bonds[i].second == move[k]) found=TRUE;
          if (!found) move[l++]=mol->bonds[i].second;
        }
        if (mol->bonds[i].second == move[j])
        {
          found=FALSE;
          for (k=0; k<l; k++)
            if (mol->bonds[i].first == move[k]) found=TRUE;
          if (!found) move[l++]=mol->bonds[i].first;
         }
      }
    }
  }
  return(l);
}

int findAttachedAtoms(int *move, int n, int except)
{
  int nmove=1, nmoveold=0, save=except;
  register int i;

  while (nmove != nmoveold)
  {
    nmoveold=nmove;
    nmove=findBondPartners(move, nmoveold, except);
  }
  for (i=1; i<nmove; i++)
  {
    if (move[i] == save)
    {
      nmove=(-nmove);
      break;
    }
  }
  if (debug)
  {
    printf("Atoms to move: ");
    if (nmove < 0)
      printf("ring detected, no atoms to move.");
    else
    {
      for (i=0; i<nmove; i++)
        printf("%d ", move[i]);
    }
    printf("\n");
  }
  return(nmove);
}

void changeBondLength(int which, int nmove, int *move, double newLength)
{
  struct MOLECULE *mol;
  double oldLength, factor;
  double vector[3];
  int i, j;

  if (windows[VIEWER].set >= 0)
    mol=&molecules[windows[VIEWER].set];
  else
    mol=&molecules[0];
  i=mol->internals[which].atoms[0];
  j=mol->internals[which].atoms[1];
  oldLength=bondLength(mol->atoms, i, j);
  factor=newLength/oldLength-1.0;
  vector[0]=factor*(mol->atoms[i].x-mol->atoms[j].x);
  vector[1]=factor*(mol->atoms[i].y-mol->atoms[j].y);
  vector[2]=factor*(mol->atoms[i].z-mol->atoms[j].z);
  mol->internals[which].value=bndfac*newLength;
  mol->internals[which].x=(mol->atoms[i].x+vector[0]+mol->atoms[j].x)*0.5+0.05;
  mol->internals[which].y=(mol->atoms[i].y+vector[1]+mol->atoms[j].y)*0.5+0.05;
  mol->internals[which].z=(mol->atoms[i].z+vector[2]+mol->atoms[j].z)*0.5+0.05;
  for (i=0; i<nmove; i++)
  {
    mol->atoms[move[i]].x+=vector[0];
    mol->atoms[move[i]].y+=vector[1];
    mol->atoms[move[i]].z+=vector[2];
  }
}

int changeAngle(int which, int nmove, int *move, double newAngle)
{
  struct MOLECULE *mol;
  double oldAngle, matrix[3][3];
  double v1[3], v2[3], v3[3];
  register double x, y, z;
  register int i, j, k;

  if (windows[VIEWER].set >= 0)
    mol=&molecules[windows[VIEWER].set];
  else
    mol=&molecules[0];
  i=mol->internals[which].atoms[0];
  j=mol->internals[which].atoms[1];
  k=mol->internals[which].atoms[2];
  if (mol->internals[which].type == ANGLE)
  {
    oldAngle=bondAngle(mol->atoms, i, j, k);
    v1[0]=mol->atoms[i].x-mol->atoms[j].x;
    v1[1]=mol->atoms[i].y-mol->atoms[j].y;
    v1[2]=mol->atoms[i].z-mol->atoms[j].z;
    v2[0]=mol->atoms[k].x-mol->atoms[j].x;
    v2[1]=mol->atoms[k].y-mol->atoms[j].y;
    v2[2]=mol->atoms[k].z-mol->atoms[j].z;
    vectorProduct(v1, v2, v3);

/* The two bonds are co-linear, find an arbitrary vector perpendicular to the first bond */

    if (fabs(v3[0]) < 1.0e-6 && fabs(v3[1]) < 1.0e-6 && fabs(v3[2]) < 1.0e-6)
    {
      v3[0]=-v1[1]/v1[0];
      v3[1]=1.0;
      v3[2]=0.0;
    }
  }
  else  /* torsion angle */
  {
    oldAngle=torsionAngle(mol->atoms, i, j, k,
                          mol->internals[which].atoms[3]);
/* Torsion angle is undefined (one of the bond angles in 180 deg) */
    if (oldAngle > 360.) return(FALSE);
    v3[0]=mol->atoms[k].x-mol->atoms[j].x;
    v3[1]=mol->atoms[k].y-mol->atoms[j].y;
    v3[2]=mol->atoms[k].z-mol->atoms[j].z;
  }
  makeRotationMatrix(newAngle-oldAngle, v3[0], v3[1], v3[2], matrix);
  for (i=0; i<nmove; i++)
  {
    x=mol->atoms[move[i]].x-mol->atoms[j].x;
    y=mol->atoms[move[i]].y-mol->atoms[j].y;
    z=mol->atoms[move[i]].z-mol->atoms[j].z;
    mol->atoms[move[i]].x=x*matrix[0][0]+y*matrix[0][1]+z*matrix[0][2]+mol->atoms[j].x;
    mol->atoms[move[i]].y=x*matrix[1][0]+y*matrix[1][1]+z*matrix[1][2]+mol->atoms[j].y;
    mol->atoms[move[i]].z=x*matrix[2][0]+y*matrix[2][1]+z*matrix[2][2]+mol->atoms[j].z;
  }
  mol->internals[which].value=newAngle;
  return(TRUE);
}

void makeRotationMatrix(double angle, double x, double y, double z, double matrix[3][3])
{
  double r, ca, sa, xx, xy, xz, yy, yz, zz;
  register double ca1;

  r=1.0/sqrt(x*x+y*y+z*z);
  x*=r;
  y*=r;
  z*=r;
  angle*=atan(1.0)/45.0;
  ca=cos(angle);
  sa=sin(angle);
  xx=x*x;
  xy=x*y;
  xz=x*z;
  yy=y*y;
  yz=y*z;
  zz=z*z;
  ca1=1.0-ca;
  matrix[0][0]=ca1*xx+ca;
  matrix[1][0]=ca1*xy-sa*z;
  matrix[2][0]=ca1*xz+sa*y;
  matrix[0][1]=ca1*xy+sa*z;
  matrix[1][1]=ca1*yy+ca;
  matrix[2][1]=ca1*yz-sa*x;
  matrix[0][2]=ca1*xz-sa*y;
  matrix[1][2]=ca1*yz+sa*x;
  matrix[2][2]=ca1*zz+ca;
}

void addUndo(struct MOLECULE *mol)
{
  if (nUndo++ == 0)
    undo=(struct UNDO *)getmem((size_t)nUndo, sizeof(struct UNDO));
  else
    undo=(struct UNDO *)expmem((void *)undo, (size_t)nUndo, sizeof(struct UNDO));
  undo[nUndo-1].coord=saveGeometry(mol);
  undo[nUndo-1].mol=windows[VIEWER].set;
  undo[nUndo-1].na=mol->na;
  setMenuItem(VIEWER2_UNDO, True);
  if (undoButton != NULL) XtVaSetValues(undoButton, XmNsensitive, True, NULL);
}

void undoGeometry(Widget w, caddr_t dummy, XmAnyCallbackStruct *data)
{
  int imolSave;
  register int i, imol;

  if (undo != NULL)
  {
    imolSave=windows[VIEWER].set;
    imol=windows[VIEWER].set=undo[--nUndo].mol;
    if (imol == (-1))
      deleteMolecule((Widget)0, (caddr_t)TRUE, (XmAnyCallbackStruct *)0);
    else
    {
      if (undo[nUndo].na >= molecules[imol].na)
      {
          molecules[imol].atoms=(struct ATOM *)expmem((void *)molecules[imol].atoms,
                                  (size_t)undo[nUndo].na, sizeof(struct ATOM));
      }
      molecules[imol].na=undo[nUndo].na;
      restoreGeometry(undo[nUndo].coord, imol);
      for (i=0; i<molecules[imol].ninternal; i++)
      {
        if (molecules[imol].internals[i].atoms[0] >= molecules[imol].na ||
            molecules[imol].internals[i].atoms[1] >= molecules[imol].na ||
            molecules[imol].internals[i].atoms[2] >= molecules[imol].na ||
            molecules[imol].internals[i].atoms[3] >= molecules[imol].na)
          clearGeometry((Widget)0, (caddr_t)i, (XmAnyCallbackStruct *)0);
        else
          {
          (void)calcInternal(&molecules[imol], i);
            if (molecules[imol].internals[i].type == TORSION &&
                molecules[imol].internals[i].value > 360.)
            clearGeometry((Widget)0, (caddr_t)i, (XmAnyCallbackStruct *)0);
        }
      }
      makeConnectivity(&molecules[imol], FALSE, TRUE);
      setAtom(&molecules[imol]);
    }
    if (nUndo)
      undo=(struct UNDO *)expmem((void *)undo, (size_t)nUndo, sizeof(struct UNDO));
    else
    {
      fremem((void *)&undo);
      setMenuItem(VIEWER2_UNDO, False);
      if (undoButton != NULL) XtVaSetValues(undoButton, XmNsensitive, False, NULL);
    }
    windows[VIEWER].set=imolSave;
    redraw(VIEWER);
  }
}