/*******************************************************************************
*                                                                              *
*                                   Viewmol                                    *
*                                                                              *
*                              M K U N I T C . C                               *
*                                                                              *
*                 Copyright (c) Joerg-R. Hill, December 2000                   *
*                                                                              *
********************************************************************************
*
* $Id: mkunitc.c,v 1.5 2000/12/10 15:11:57 jrh Exp $
* $Log: mkunitc.c,v $
* Revision 1.5  2000/12/10 15:11:57  jrh
* Release 2.3
*
* Revision 1.4  1999/05/24 01:26:34  jrh
* Release 2.2.1
*
* Revision 1.3  1999/02/07 21:52:54  jrh
* Release 2.2
*
* Revision 1.1  1996/12/10  18:42:14  jrh
* Initial revision
*
*/
#include<math.h>
#include<stdio.h>
#include "viewmol.h"

#define MIN(a, b) (a) < (b) ? (a) : (b)

int makeUnitCell(double, double, double, double, double, double, int);
void fractionalToCartesian(double, double, double, double, double, double, int,
                           int, int);
double *cartesianToFractional(double, double, double, double, double, double,
                              double *, double *, double *, int *);
void expandCell(void);

extern struct WINDOW windows[];
extern struct MOLECULE *molecules;
extern double forceScale;

extern int  makeConnectivity(struct MOLECULE *, int, int);
extern double bondLength(struct MOLECULE *, int, int);
extern double bondAngle(struct MOLECULE *, int, int, int);
extern double dist(double, double, double, double, double, double);
extern double angle(double, double, double, double, double, double,
                    double, double, double);
extern void millerPlane(void);
extern void *getmem(size_t, size_t);
extern void *expmem(void *, size_t, size_t);
extern void fremem(void **);

int makeUnitCell(double a, double b, double c, double alpha, double beta,
                 double gamma, int frac)
{
/* This function computes from the unit cell lengths
   and angles the cartesian coordinates of the eight
   corners and stores them in atoms[], if frac == TRUE
   the input coordinates are fractional coordinates
   and this function will convert them to cartesian
   coordinates */

  struct MOLECULE *mol;
  register double torad;
  register int n;

  if (windows[VIEWER].set >= 0)
    mol=&molecules[windows[VIEWER].set];
  else
    mol=&molecules[0];
  n=mol->na;
  mol->atoms=(struct ATOM *)expmem((void *)mol->atoms, n+9,
                                   sizeof(struct ATOM));
  
  torad=atan(1.0)/45.0;
  strcpy(mol->atoms[n].name, "Uc");
  mol->atoms[n].x=0.0;
  mol->atoms[n].y=0.0;
  mol->atoms[n].z=0.0;
  mol->atoms[n].ref=n;
  mol->atoms[n].mass=0.0;
  mol->atoms[n++].flags=ORIGINAL;
  strcpy(mol->atoms[n].name, "Uc");
  mol->atoms[n].x=a;
  mol->atoms[n].y=0.0;
  mol->atoms[n].z=0.0;
  mol->atoms[n].ref=n;
  mol->atoms[n].mass=0.0;
  mol->atoms[n++].flags=ORIGINAL;
  strcpy(mol->atoms[n].name, "Uc");
  mol->atoms[n].x=b*cos(torad*gamma);
  mol->atoms[n].y=sqrt(b*b-mol->atoms[n].x*mol->atoms[n].x);
  mol->atoms[n].z=0.0;
  mol->atoms[n].ref=n;
  mol->atoms[n].mass=0.0;
  mol->atoms[n++].flags=ORIGINAL;
  strcpy(mol->atoms[n].name, "Uc");
  mol->atoms[n].x=c*cos(torad*beta);
  mol->atoms[n].y=(b*c*cos(torad*alpha)-mol->atoms[n-1].x*mol->atoms[n].x)/mol->atoms[n-1].y;
  mol->atoms[n].z=sqrt(c*c-mol->atoms[n].x*mol->atoms[n].x-mol->atoms[n].y*mol->atoms[n].y);
  mol->atoms[n].ref=n;
  mol->atoms[n].mass=0.0;
  mol->atoms[n].flags=ORIGINAL;
  n++;
  strcpy(mol->atoms[n].name, "Uc");
  mol->atoms[n].x=mol->atoms[n-1].x+mol->atoms[n-3].x;
  mol->atoms[n].y=mol->atoms[n-1].y+mol->atoms[n-3].y;
  mol->atoms[n].z=mol->atoms[n-1].z+mol->atoms[n-3].z;
  mol->atoms[n].ref=n;
  mol->atoms[n].mass=0.0;
  mol->atoms[n].flags=ORIGINAL;
  n++;
  strcpy(mol->atoms[n].name, "Uc");
  mol->atoms[n].x=mol->atoms[n-1].x+mol->atoms[n-3].x;
  mol->atoms[n].y=mol->atoms[n-1].y+mol->atoms[n-3].y;
  mol->atoms[n].z=mol->atoms[n-1].z+mol->atoms[n-3].z;
  mol->atoms[n].ref=n;
  mol->atoms[n].mass=0.0;
  mol->atoms[n].flags=ORIGINAL;
  n++;
  strcpy(mol->atoms[n].name, "Uc");
  mol->atoms[n].x=mol->atoms[n-4].x+mol->atoms[n-3].x;
  mol->atoms[n].y=mol->atoms[n-4].y+mol->atoms[n-3].y;
  mol->atoms[n].z=mol->atoms[n-4].z+mol->atoms[n-3].z;
  mol->atoms[n].ref=n;
  mol->atoms[n].mass=0.0;
  mol->atoms[n].flags=ORIGINAL;
  n++;
  strcpy(mol->atoms[n].name, "Uc");
  mol->atoms[n].x=mol->atoms[n-6].x+mol->atoms[n-5].x;
  mol->atoms[n].y=mol->atoms[n-6].y+mol->atoms[n-5].y;
  mol->atoms[n].z=mol->atoms[n-6].z+mol->atoms[n-5].z;
  mol->atoms[n].ref=n;
  mol->atoms[n].mass=0.0;
  mol->atoms[n].flags=ORIGINAL;
  n++;

  if (frac) fractionalToCartesian(a, b, c, alpha, beta, gamma, TRUE, 0, mol->na);
  mol->na=n;
  return(TRUE);
}

void fractionalToCartesian(double a, double b, double c, double alpha,
                           double beta, double gamma, int doHistory,
                           int from, int to)
{
  struct MOLECULE *mol;
  double torad, x2, x3, y2, y3, z3;
  register int i, j, k;

  if (windows[VIEWER].set >= 0)
    mol=&molecules[windows[VIEWER].set];
  else
    mol=&molecules[0];
  torad=atan(1.0)/45.0;
  alpha*=torad;
  beta*=torad;
  gamma*=torad;
  x2=b*cos(gamma);
  y2=sqrt(b*b-x2*x2);
  x3=c*cos(beta);
  y3=(b*c*cos(alpha)-x2*x3)/y2;
  z3=sqrt(c*c-x3*x3-y3*y3);
  for (i=from; i<to; i++)
  {
    mol->atoms[i].x = a*mol->atoms[i].x+x2*mol->atoms[i].y+x3*mol->atoms[i].z;
    mol->atoms[i].y =y2*mol->atoms[i].y+y3*mol->atoms[i].z;
    mol->atoms[i].z*=z3;
  }
  if (doHistory)
  {
    for (i=0; i<mol->nhist; i++)
    {
	k=i*(to+8);
      a=dist(mol->history[k+to].x, mol->history[k+to].y, mol->history[k+to].z,
             mol->history[k+to+1].x, mol->history[k+to+1].y, mol->history[k+to+1].z);
      b=dist(mol->history[k+to].x, mol->history[k+to].y, mol->history[k+to].z,
             mol->history[k+to+2].x, mol->history[k+to+2].y, mol->history[k+to+2].z);
      c=dist(mol->history[k+to].x, mol->history[k+to].y, mol->history[k+to].z,
             mol->history[k+to+3].x, mol->history[k+to+3].y, mol->history[k+to+3].z);
	alpha=torad*angle(mol->history[k+to+2].x, mol->history[k+to+2].y, mol->history[k+to+2].z,
                        mol->history[k+to].x, mol->history[k+to].y, mol->history[k+to].z,
                        mol->history[k+to+3].x, mol->history[k+to+3].y, mol->history[k+to+3].z);
      beta=torad*angle(mol->history[k+to+1].x, mol->history[k+to+1].y, mol->history[k+to+1].z,
                       mol->history[k+to].x, mol->history[k+to].y, mol->history[k+to].z,
                       mol->history[k+to+3].x, mol->history[k+to+3].y, mol->history[k+to+3].z);
      gamma=torad*angle(mol->history[k+to+1].x, mol->history[k+to+1].y, mol->history[k+to+1].z,
                        mol->history[k+to].x, mol->history[k+to].y, mol->history[k+to].z,
                        mol->history[k+to+2].x, mol->history[k+to+2].y, mol->history[k+to+2].z);
      x2=b*cos(gamma);
      y2=sqrt(b*b-x2*x2);
      x3=c*cos(beta);
      y3=(b*c*cos(alpha)-x2*x3)/y2;
      z3=sqrt(c*c-x3*x3-y3*y3);
      for (j=from; j<to; j++)
      {
        mol->history[k+j].x = a*mol->history[k+j].x
                            +x2*mol->history[k+j].y
                            +x3*mol->history[k+j].z;
        mol->history[k+j].y =y2*mol->history[k+j].y
                            +y3*mol->history[k+j].z;
        mol->history[k+j].z*=z3;
      }
    }
  }
}

double *cartesianToFractional(double a, double b, double c, double alpha,
                              double beta, double gamma, double *xmin,
                              double *ymin, double *zmin, int *natoms)
{
  struct MOLECULE *mol;
  double torad, x2, x3, y2, y3, z3;
  double *frac;
  register int i, j, n=0;

/* Transform cartesian coordinates to fractional coordinates */

  if (windows[VIEWER].set >= 0)
    mol=&molecules[windows[VIEWER].set];
  else
    mol=&molecules[0];
  torad=atan(1.0)/45.0;
  alpha*=torad;
  beta*=torad;
  gamma*=torad;
  x2=b*cos(gamma);
  y2=sqrt(b*b-x2*x2);
  x3=c*cos(beta);
  y3=(b*c*cos(alpha)-x2*x3)/y2;
  z3=1.0/sqrt(c*c-x3*x3-y3*y3);
  y2=1.0/y2;
  frac=(double *)getmem((size_t)(3*mol->na), sizeof(double));
  j=0;
  *xmin=0.0;
  *ymin=0.0;
  *zmin=0.0;
  for (i=0; i<mol->na-8; i++)
  {
    if (mol->atoms[i].flags & ORIGINAL)
    {
      frac[j+2]=mol->atoms[i].z*z3;
      frac[j+1]=(mol->atoms[i].y-frac[j+2]*y3)*y2;
      frac[j]  =(mol->atoms[i].x-frac[j+1]*x2-frac[j+2]*x3)/a;
      *xmin=MIN(*xmin, frac[j]);
      *ymin=MIN(*ymin, frac[j+1]);
      *zmin=MIN(*zmin, frac[j+2]);
      j+=3;
      n++;
    }
  }
  j=0;
  for (i=0; i<n; i++)
  {
    if (*xmin < 0.0) frac[j]  -=(*xmin);
    if (*ymin < 0.0) frac[j+1]-=(*ymin);
    if (*zmin < 0.0) frac[j+2]-=(*zmin);
    j+=3;
  }
  *natoms=n;
  return(frac);
}

void expandCell(void)
{
  struct NEW
  {
    double x;
    double y;
    double z;
    int    ref;
  } *new;
  struct MOLECULE *mol;
  struct ELEMENT *uc;
  double a, b, c, alpha, beta, gamma;
  double xmin, ymin, zmin;
  double *frac, cell[3][8], rad;
  int n, nnew;
  register int i, j, k, l, m;

  if (windows[VIEWER].set >= 0)
    mol=&molecules[windows[VIEWER].set];
  else
    mol=&molecules[0];

  if (mol->nmiller != 0)
  {
    mol->na-=mol->nmiller;
    mol->atoms=(struct ATOM *)expmem((void *)mol->atoms, (size_t)mol->na, sizeof(struct ATOM));
  }
  n=mol->na-8;
  a=bondLength(mol, n, n+1);
  b=bondLength(mol, n, n+2);
  c=bondLength(mol, n, n+3);
  alpha=bondAngle(mol, n+2, n, n+3);
  beta=bondAngle(mol, n+1, n, n+3);
  gamma=bondAngle(mol, n+1, n, n+2);
  n=0;
  for (i=0; i<mol->na; i++)
    if (mol->atoms[i].flags & ORIGINAL) n++;
  n-=8;
  frac=cartesianToFractional(a, b, c, alpha, beta, gamma, &xmin, &ymin, &zmin,
                             &n);
  new=(struct NEW *)getmem((size_t)(n*ceil(mol->cellFactor[0])*ceil(mol->cellFactor[1])*ceil(mol->cellFactor[2])),
                           sizeof(struct NEW));
  nnew=0;
  for (i=0; i<(int)ceil(mol->cellFactor[0]); i++)
  {
    for (j=0; j<(int)ceil(mol->cellFactor[1]); j++)
    {
      for (k=0; k<(int)ceil(mol->cellFactor[2]); k++)
      {
        if (i == 0 && j == 0 && k == 0) continue;
        l=0;
        for (m=0; m<n; m++)
        {
          new[nnew].x=frac[l++]+(double)i;
          new[nnew].y=frac[l++]+(double)j;
          new[nnew].z=frac[l++]+(double)k;
          if (new[nnew].x <= mol->cellFactor[0] && new[nnew].y <= mol->cellFactor[1] &&
              new[nnew].z <= mol->cellFactor[2]) new[nnew++].ref=m;
        }
      }
    }
  }
  fremem((void *)&frac);

  /* Save unit cell corners since code assumes that the last eight atoms
     are the unit cell corners if a unit cell is present */
  uc=mol->atoms[mol->na-8].element;
  rad=mol->atoms[mol->na-8].rad;
  for (i=mol->na-8; i<mol->na; i++)
  {
    cell[0][i-mol->na+8]=mol->atoms[i].x;
    cell[1][i-mol->na+8]=mol->atoms[i].y;
    cell[2][i-mol->na+8]=mol->atoms[i].z;
  }

  /* Extend/shrink atoms structure */
  mol->atoms=(struct ATOM *)expmem((void *)mol->atoms, (size_t)(n+nnew+8),
                                   sizeof(struct ATOM));
  j=0;
  for (i=n; i<n+nnew; i++)
  {
    mol->atoms[i].x=new[j].x+xmin;
    mol->atoms[i].y=new[j].y+ymin;
    mol->atoms[i].z=new[j].z+zmin;
    k=new[j++].ref;
    mol->atoms[i].rad=mol->atoms[k].rad;
    mol->atoms[i].radScale=mol->atoms[k].radScale;
    mol->atoms[i].mass=mol->atoms[k].mass;
    mol->atoms[i].neutronScatterfac=mol->atoms[k].neutronScatterfac;
    mol->atoms[i].nbonds=mol->atoms[k].nbonds;
    mol->atoms[i].nelectrons=mol->atoms[k].nelectrons;
    mol->atoms[i].ref=k;
    mol->atoms[i].flags=mol->atoms[k].flags & ~ORIGINAL;
    mol->atoms[i].element=mol->atoms[k].element;
    mol->atoms[i].basis=mol->atoms[k].basis;
    strcpy(mol->atoms[i].name, mol->atoms[k].name);
    strcpy(mol->atoms[i].basisname, mol->atoms[k].basisname);
  }
  fremem((void *)&new);
  fractionalToCartesian(a, b, c, alpha, beta, gamma, FALSE, n, n+nnew);
  mol->na=n+nnew;
  j=0;
  for (i=mol->na; i<mol->na+8; i++)
  {
    mol->atoms[i].x=cell[0][i-mol->na];
    mol->atoms[i].y=cell[1][i-mol->na];
    mol->atoms[i].z=cell[2][i-mol->na];
    mol->atoms[i].element=uc;
    mol->atoms[i].rad=rad;
    mol->atoms[i].radScale=1.0;
    mol->atoms[i].ref=n+j;
    strcpy(mol->atoms[i].name, "Uc");
    j++;
  }
  mol->na+=8;
  a=forceScale;
  makeConnectivity(mol, TRUE, TRUE);
  forceScale=a;
  if (mol->nmiller != 0)
  {
    mol->na+=mol->nmiller;
    mol->atoms=(struct ATOM *)expmem((void *)mol->atoms, (size_t)mol->na, sizeof(struct ATOM));
    millerPlane();
  }
}