/*****************************************/
/*  Micro-structure V1.0                 */
/*****************************************/
/*  1. Nano-materials by Voronoi cells   */
/*  2. Associated Grain Boundaries       */
/*  3. Spherical Voids                   */
/*  4. Elliptical (un)through cracks     */
/*                                       */
/*                         Apr. 9, 1999  */
/*                          Ju Li (MIT)  */
/*****************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include <sys/time.h>
#include <sys/types.h>

#ifndef PI
#define PI 3.14159265358979323846
#endif
#define MAX_STRLEN 128
#define TP_STRLEN 2
#define TINY (1.E-9)
#define MAX_COMPRESSION 1.1
#define TIME_TO_NOTICE (0.)
#define MAX_GRAINS 1024
#define MAX_VOIDS  128
#define MAX_CRACKS 2
#define min(x,y) ((x)<(y)?(x):(y))
#define max(x,y) ((x)>(y)?(x):(y))
#define square(x) ((x)*(x))

static int iseed, npa, np, npm, npr, npgb, nvoids, ncracks,
    use_config_H, *tag, *gb, nc[3], nv[5], ex[3];
static double CL, min_d, min_r, *mass, *s, *massm, *sm, r[4],
    sv[3*MAX_GRAINS], ov[4*MAX_GRAINS], dv[4*MAX_GRAINS],
    cuto[3], cuts[3][3], H[3][3], HI[3][3],
    voids[MAX_VOIDS][4], cracks[MAX_CRACKS][6],
    R[MAX_GRAINS][3][3];
static char buf[MAX_STRLEN], fn_config[MAX_STRLEN],
    voronoi_cell[MAX_STRLEN], fn_config_gb[MAX_STRLEN],
    **tp, **tpm;

static void write_config (char filename[], int np, int GB_only);
static void ds (double s1[], double s2[], double d[]);
static double r2 (double si[], double sj[]);
static void waitfor (double seconds);
static double frandom();
static void veceqv (double a[3], double b[3]);
static void vecadd (double a[3], double b[3], double c[3]);
static void vecminus (double a[3], double b[3], double c[3]);
static double dot (double a[3], double b[3]);
static double normalize (double a[3]);
static void cross (double a[3], double b[3], double c[3]);
static void random_unit_vec (double x[3]);
static void rotate (double a[4], double b[3], double c[3]);
static void mateqv (double A[3][3], double B[3][3]);
static void matran (double A[3][3], double B[3][3]);
static void matadd (double A[3][3], double B[3][3], double C[3][3]);
static void matsub (double A[3][3], double B[3][3], double C[3][3]);
static void matmul (double A[3][3], double B[3][3], double C[3][3]);
static void vec_mul_mat (double a[3], double B[3][3], double c[3]);
static double matinv (double A[3][3], double B[3][3]);


int main (int argc, char *argv[])
{
    int i, j, k, l, m, n;
    double cc, a[3], b[4], c[3], d[4];
    FILE *config;
    
    for (buf[0]=(char)0,i=0; i<argc; i++)
	strcat(strcat(buf,argv[i]),(i==argc-1)?"":" ");
    printf ("\nProgram \"%s\" starts (pid = %d).\n\n",
	    buf, getpid());
    
    /* read instructions from "mic" file */
    fscanf (stdin, "Load configuration from file ('NULL'->start "
	    "anew); use its H?:\n");
    fscanf (stdin, "%s %s\n\n", fn_config, buf);
    use_config_H = (buf[0]=='1') || (buf[0]=='y') || (buf[0]=='Y');
    
    fscanf (stdin, "Default unit cell (number of atoms, "
	    "characteristic length in A):\n");
    fscanf (stdin, "%d %lf\n", &npa, &CL);
    fscanf (stdin, "Bravais lattice [Cl]:\n");
    fscanf (stdin, "%lf %lf %lf\n", &H[0][0], &H[0][1], &H[0][2]);
    fscanf (stdin, "%lf %lf %lf\n", &H[1][0], &H[1][1], &H[1][2]);
    fscanf (stdin, "%lf %lf %lf\n", &H[2][0], &H[2][1], &H[2][2]);
    fscanf (stdin, "Default multiplications of the unit cell:\n");
    fscanf (stdin, "%d %d %d\n", &nc[0], &nc[1], &nc[2]);

    /* determine the number of particles of the reference system */
    if ( strcasecmp(fn_config, "NULL") )
    {
	if ((config=fopen(fn_config,"r")) == NULL)
	{
	    printf("Error: cannot open config file \"%s\".\n",
		   fn_config);
	    exit(1);
	}
	printf ("The reference configuration is loaded from \"%s\".\n",
		fn_config);
	fscanf (config, "Number of particles = %d\n\n", &np);
    }
    else np = npa*nc[0]*nc[1]*nc[2];
    
    printf ("Number of particles prior to micro-structural transformation "
	    "= %d.\n\n", np);
    npm = ceil(MAX_COMPRESSION*np);
    
    /* allocate reference system properties */
    s     = (double *) malloc (3*np*sizeof(double));
    mass  = (double *) malloc (np*sizeof(double));
    tp    = (char **)  malloc (np*sizeof(char *));
    tp[0] = (char *)   malloc (np*(TP_STRLEN+1)*sizeof(char));
    for (i=1; i<np; i++) tp[i] = tp[i-1] + TP_STRLEN + 1;
    /* allocate transformed system properties */
    sm     = (double *) malloc (3*npm*sizeof(double));
    massm  = (double *) malloc (npm*sizeof(double));
    tpm    = (char **)  malloc (npm*sizeof(char *));
    tpm[0] = (char *)   malloc (npm*(TP_STRLEN+1)*sizeof(char));
    for (i=1; i<npm; i++) tpm[i] = tpm[i-1] + TP_STRLEN + 1;
    tag = (int *) malloc (npm*sizeof(int)); /* species tag */
    gb  = (int *) malloc (npm*sizeof(int)); /* gb atoms index */

    /* load in default system lattice */
    matinv (H, HI);
    fscanf (stdin, "TP Mass(amu)     x [Cl]       y [Cl]       z [Cl]\n");
    for (i=0; i<3*npa; i+=3)
    {  /* what is read is not s yet, but x */
	fscanf (stdin, "%2s %lf %lf %lf %lf\n", tp[i/3],
		&mass[i/3], &s[i], &s[i+1], &s[i+2]);
	vec_mul_mat (&s[i], HI, &s[i]);
	s[i]   /= nc[0];
	s[i+1] /= nc[1];
	s[i+2] /= nc[2];
    }
    fscanf (stdin, "\n");

    /* convert H[][] in stdin from CL unit to A */
    for (i=0; i<3; i++)
	for (j=0; j<3; j++)
	    H[i][j] *= nc[i] * CL;
    
    /* actual lattice and particle properties */
    if ( strcasecmp(fn_config, "NULL") )
    {  /* if loading the configuration from a file */
	fscanf (config, "H(1,1) = %lf A\n", &HI[0][0]);
	fscanf (config, "H(1,2) = %lf A\n", &HI[0][1]);
	fscanf (config, "H(1,3) = %lf A\n", &HI[0][2]);
	fscanf (config, "nc(1)  = %d\n\n",  &nv[0]);
	fscanf (config, "H(2,1) = %lf A\n", &HI[1][0]);
	fscanf (config, "H(2,2) = %lf A\n", &HI[1][1]);
	fscanf (config, "H(2,3) = %lf A\n", &HI[1][2]);
	fscanf (config, "nc(2)  = %d\n\n",  &nv[1]);
	fscanf (config, "H(3,1) = %lf A\n", &HI[2][0]);
	fscanf (config, "H(3,2) = %lf A\n", &HI[2][1]);
	fscanf (config, "H(3,3) = %lf A\n", &HI[2][2]);
	fscanf (config, "nc(3)  = %d\n\n",  &nv[2]);
	if ( use_config_H )
	{ /* H and nc from config file will be used */
	    mateqv (HI, H);
	    nc[0] = nv[0];
	    nc[1] = nv[1];
	    nc[2] = nv[2];
	}
	fscanf (config, "TP Mass(amu)     sx        sy        sz"
		"     sx1(1/ns)  sy1(1/ns)  sz1(1/ns)\n");
	for (i=0; i<3*np; i+=3)
	    fscanf (config, "%2s %lf %lf %lf %lf %lf %lf %lf\n",
		    tp[i/3], &mass[i/3], &s[i], &s[i+1], &s[i+2],
		    &cc,  &cc, &cc);
	fclose (config);
    }
    else /* build our own crystalline reference system based on stdin */
	for (n=i=0; i<nc[0]; i++)
	    for (j=0; j<nc[1]; j++)
		for (k=0; k<nc[2]; k++)
		    for (l=0; l<npa; l++,n++)
		    {
			s[3*n]   = s[3*l]   + (double)i/nc[0];
			s[3*n+1] = s[3*l+1] + (double)j/nc[1];
			s[3*n+2] = s[3*l+2] + (double)k/nc[2];
			mass[n]  = mass[l];
			for (m=0; m<TP_STRLEN+1; m++) tp[n][m] = tp[l][m];
		    }
    
    /* H is ready now */
    printf ("              | %f %f %f |\n", H[0][0], H[0][1], H[0][2]);
    printf ("Reference H = | %f %f %f | A, nc = (%d %d %d)\n",
	    H[1][0], H[1][1], H[1][2], nc[0], nc[1], nc[2]);
    printf ("              | %f %f %f |\n\n", H[2][0], H[2][1], H[2][2]);
    matinv (H, HI);

    /* need randomness from now on */
    fscanf (stdin, "Random number generator seed (0->time seed):\n");
    fscanf (stdin, "%d\n\n", &iseed); 
    if (iseed == 0) iseed = time(NULL);
    srandom ((unsigned)iseed);
    
    /* Voronoi sites can be arranged as a crystal too */
    fscanf (stdin, "Supercell partition into Voronoi site unit cells:\n");
    fscanf (stdin, "%d %d %d\n\n", &nv[0], &nv[1], &nv[2]);

    /* its unit cell */
    fscanf (stdin, "Voronoi sites unit cell (sc,bcc,fcc,random*):\n");
    fscanf (stdin, "%s\n\n", voronoi_cell);
    if ( !strcasecmp(voronoi_cell, "sc") )
    { /* simple cubic */
	nv[3] = 1;  /* number of sites in a unit cell */
	sv[0] = 0.5; sv[1] = 0.5; sv[2] = 0.5;
    }
    else if ( !strcasecmp(voronoi_cell, "bcc") )
    {
	nv[3] = 2;
	sv[0] = 0.25; sv[1] = 0.25; sv[2] = 0.25;
	sv[3] = 0.75; sv[4] = 0.75; sv[5] = 0.75;
    }
    else if ( !strcasecmp(voronoi_cell, "fcc") )
    {
	nv[3] = 4;
	sv[0] = 0.25;  sv[1]  = 0.25;  sv[2]  = 0.25;
	sv[3] = 0.75;  sv[4]  = 0.75;  sv[5]  = 0.25;
	sv[6] = 0.25;  sv[7]  = 0.75;  sv[8]  = 0.75;
	sv[9] = 0.75;  sv[10] = 0.25;  sv[11] = 0.75;
    }
    else if ( !strncasecmp(voronoi_cell, "random", 6) )
    {   /* e.g. "random16" means a unit cell with 16 random sites */
	nv[3] = (strlen(voronoi_cell)==6)? 1 : atoi(voronoi_cell+6);
	/* assign the sites random positions */
	for (i=0; i<3*nv[3]; i++) sv[i] = frandom();
    }
    else
    {
	printf ("Error: invalid Voronoi cell type \"%s\".\n",
		voronoi_cell);
	exit(1);
    }
    
    /* build Voronoi site super-lattice */
    for (n=i=0; i<nv[0]; i++)
	for (j=0; j<nv[1]; j++)
	    for (k=0; k<nv[2]; k++)
		for (l=0; l<nv[3]; l++,n++)
		{   /* not yet "normalized" by nv[0-2] */
		    sv[3*n]   = sv[3*l]   + i;
		    sv[3*n+1] = sv[3*l+1] + j;
		    sv[3*n+2] = sv[3*l+2] + k;
		    /* assign them random orientations */
		    random_unit_vec (&ov[4*n]);
		    ov[4*n+3] = frandom()*360.;
		}

    fscanf (stdin, "Number of non-random Voronoi site rotations:\n");
    fscanf (stdin, "%d\n", &n);
    fscanf (stdin, "Site   Axis x,y,z       Degrees\n");
    for (i=0; i<n; i++)
    { /* cater to special requests */
	fscanf (stdin, "%d ", &j);
	fscanf (stdin, "%lf %lf %lf %lf\n",
		&ov[4*j], &ov[4*j+1], &ov[4*j+2], &ov[4*j+3]);
	normalize (&ov[4*j]);
    }
    fscanf (stdin, "\n");
    
    /* total number of Voronoi sites */
    nv[4] = nv[0]*nv[1]*nv[2]*nv[3];
    for (n=0; n<nv[4]; n++)
    {   /* "normalize" origin of Voronoi sites */
	sv[3*n]   /= nv[0];
	sv[3*n+1] /= nv[1];
	sv[3*n+2] /= nv[2];
	/* calculate the rotation matrix for each grain */
	for (i=0; i<3; i++)
	{  /* s[i]=1 denotes vector H[i][] */
	    rotate (&ov[4*n], &H[i][0], a);
	    vec_mul_mat (a, HI, &R[n][i][0]);
	}
	printf ("Grain %d, center at (%f %f %f):\n", n, sv[3*n],
		sv[3*n+1], sv[3*n+2]);
	printf ("    | %f %f %f |\n", R[n][0][0], R[n][0][1], R[n][0][2]);
	printf ("R = | %f %f %f |\n", R[n][1][0], R[n][1][1], R[n][1][2]);
	printf ("    | %f %f %f |\n\n", R[n][2][0], R[n][2][1], R[n][2][2]);
    }
    
    fscanf (stdin, "Minimum atomic distance [A], GB thickness/2 [A]:\n");
    fscanf (stdin, "%lf %lf\n\n", &min_r, &min_d);
    if (min_d < min_r)
    {
	printf ("GB thickness/2 (%.4f) must >= minimum atomic distance"
		" (%.4f)\nfor GB particle list construction -> min_r "
		"modified.\n\n", min_d, min_r);
	min_d = min_r;
	waitfor (TIME_TO_NOTICE);
    }

    /* fill the grains with atoms after on-site rotation */
    for (npr=npgb=n=0; n<nv[4]; n++)
    {  /* first calculate the distances between sites */ 
	for (m=0; m<nv[4]; m++)
	{
	    if (m==n) continue;
	    ds (&sv[3*n], &sv[3*m], c);
	    vec_mul_mat (c, H, &dv[4*m]);
	    dv[4*m+3] = normalize (&dv[4*m]);
	}
	for (i=0; i<np; i++)
	{
	    ds (&sv[3*n], &s[3*i], a);
	    vec_mul_mat (a, R[n], a);
	    vecadd (&sv[3*n], a, &sm[3*npr]);
	    vec_mul_mat (a, H, b);
	    b[3] = dot(b, b);
	    for (tag[npr]=m=0; m<nv[4]; m++)
	    { /* tag = 0 - normal, >0 - GB, -1 - get rid of. */
		if (m==n) continue;
		if (b[3] >= r2(&sv[3*m], &sm[3*npr])) break;
		/* those tagged are close to GB */
		if (dot(&dv[4*m], b) > 0.5*dv[4*m+3]-min_d) tag[npr]++;
	    }
	    if (m == nv[4])
	    { /* this particle truly belongs to this Voronoi site */
		while (sm[3*npr  ] > 1) sm[3*npr  ]--;
		while (sm[3*npr  ] < 0) sm[3*npr  ]++;
		while (sm[3*npr+1] > 1) sm[3*npr+1]--;
		while (sm[3*npr+1] < 0) sm[3*npr+1]++;
		while (sm[3*npr+2] > 1) sm[3*npr+2]--;
		while (sm[3*npr+2] < 0) sm[3*npr+2]++;
		massm[npr] = mass[i];	
		for (j=0; j<TP_STRLEN+1; j++) tpm[npr][j] = tp[i][j];
		if (tag[npr] > 0) gb[npgb++] = npr;
		if (++npr >= npm)
		{
		    printf ("Error: npm = %d reached.\n", npm);
		    exit(1);
		}
	    } 
	} /* go over all atoms */
    } /* go over all Voronoi sites */

    for (i=0; i<npgb; i++)
    { /* randomly delete GB atoms that are too close */
	if (tag[gb[i]] > 0)
	    for (j=i+1; j<npgb; j++)
		if (tag[gb[j]] > 0)
		    if (r2(sm+3*gb[i], sm+3*gb[j]) < min_d*min_d)
			if (frandom() < 0.5) goto delete_i;
			else tag[gb[j]] = -1;
	continue;
    delete_i: tag[gb[i]] = -1;
    }
    
    /* create voids in the configuration */
    fscanf (stdin, "Create voids (sx sy sz radius [A]):\n");
    for (nvoids=0; ; nvoids++)
    {
	while ( ((buf[0]=getc(stdin)) != '\n') &&
		(buf[0] != (char)EOF) &&
		(buf[0] != 'n') &&
		(buf[0] != 'N') &&
		(buf[0] != '(') );
	if (buf[0] == '(')
	{
	    if (nvoids >= MAX_VOIDS)
	    {
		printf ("Error: MAX_VOIDS = %d exceeded.\n", MAX_VOIDS);
		exit(1);
	    }
	    fscanf (stdin, "%lf %lf %lf %lf)",
		    &voids[nvoids][0], &voids[nvoids][1],
		    &voids[nvoids][2], &voids[nvoids][3]);
	    for (i=0; i<npr; i++)
		if ( r2(&voids[nvoids][0], &sm[3*i]) <
		     square(voids[nvoids][3]) ) tag[i] = -1;
	}
	else
	{
	    if ( (buf[0] == 'n') || (buf[0] == 'N') )
		while ( ((buf[0]=getc(stdin) != '\n')) &&
			(buf[0] != (char)EOF) );
	    break;
	}
    }
    fscanf (stdin, "\n");
    
    /* create cracks in the configuration */
    fscanf (stdin, "Create elliptical cracks (sx sy sz d_sx d_sy d_sz):\n");
    for (ncracks=0; ; ncracks++)
    {
	while ( ((buf[0]=getc(stdin)) != '\n') &&
		(buf[0] != (char)EOF) &&
		(buf[0] != 'n') &&
		(buf[0] != 'N') &&
		(buf[0] != '(') );
	if (buf[0] == '(')
	{
	    if (ncracks >= MAX_CRACKS)
	    {
		printf ("Error: MAX_CRACKS = %d exceeded.\n", MAX_CRACKS);
		exit(1);
	    }
	    fscanf (stdin, "%lf %lf %lf ", &cracks[ncracks][0],
		    &cracks[ncracks][1], &cracks[ncracks][2]);
	    for (i=0; i<3; i++)
	    { /* read in the minor axis */
		fscanf (stdin, "%s", buf);
		cracks[ncracks][3+i] = 
		    (!strncasecmp(buf,"inf",3))? -1. : fabs(atof(buf));
	    }
	    fscanf (stdin, ")");
	    printf ("%f %f %f %f %f %f\n", cracks[ncracks][0],
		    cracks[ncracks][1], cracks[ncracks][2],
		    cracks[ncracks][3], cracks[ncracks][4],
		    cracks[ncracks][5]);
	    for (i=0; i<npr; i++)
	    {
		ds (&cracks[ncracks][0], &sm[3*i], a);
		for (cc=j=0; j<3; j++)
		    if (cracks[ncracks][3+j] > 0)
			cc += square(a[j]/cracks[ncracks][3+j]);
		if (cc <= 1) tag[i] = -1;
	    }
	}
	else
	{
	    if ( (buf[0] == 'n') || (buf[0] == 'N') )
		while ( ((buf[0]=getc(stdin) != '\n')) &&
			(buf[0] != (char)EOF) );
	    break;
	}
    }
    fscanf (stdin, "\n");
    
    fscanf (stdin, "Cut out s-chunk and replicate periodically:\n");
    fscanf (stdin, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf "
	    "%d %d %d\n", &cuto[0], &cuto[1], &cuto[2],
	    &cuts[0][0], &cuts[0][1], &cuts[0][2],
	    &cuts[1][0], &cuts[1][1], &cuts[1][2],
	    &cuts[2][0], &cuts[2][1], &cuts[2][2],
	    &ex[0], &ex[1], &ex[2]);
    matinv (cuts, HI);
    for (npgb=j=i=0; i<npr; i++)
    {
	vecminus (&sm[3*i], cuto, a);
	vec_mul_mat (a, HI, &sm[3*i]);
	if ( (sm[3*i  ]<0) || (sm[3*i  ]>=1) ||
	     (sm[3*i+1]<0) || (sm[3*i+1]>=1) ||
	     (sm[3*i+2]<0) || (sm[3*i+2]>=1) ) tag[i]=-1;
	if (tag[i] >= 0)
	{  
	    j++;  /* atom count */
	    if (tag[i] > 0) npgb++; /* GB atom count */
	}
    }
    matmul (cuts, H, HI);
    mateqv (HI, H);
    nc[0] = ceil(cuts[0][0]*nc[0]);
    nc[1] = ceil(cuts[1][1]*nc[1]);
    nc[2] = ceil(cuts[2][2]*nc[2]);
    
    fscanf (stdin, "Save entire (->*) and GB (default=*.gb) "
	    "configurations on files:\n");
    fscanf (stdin, "%s %s\n", fn_config, fn_config_gb);
    if ( !strcasecmp(fn_config, "default") )
	sprintf (fn_config, "config");
    if ( strcasecmp(fn_config, "NULL") )
    {
	printf ("Total number of atoms after transformation = %d "
		"(%.2f%%).\n", j, 100.*(j-np)/np); 
	write_config (fn_config, j, 0);
	printf ("Entire configuration saved on \"%s\".\n\n",
		fn_config);
    }
    
    if (npgb > 0)
    {
	printf ("Number of GB atoms (|d| < %.4f A) = %d (%.2f%%).\n",
		min_d, npgb, 100.*npgb/j ); 
	if ( !strcasecmp(fn_config_gb, "default") )
	    sprintf (fn_config_gb, "%s.gb", fn_config);
	if ( strncasecmp(fn_config, "NULL", 4) )
	{
	    write_config (fn_config_gb, npgb, 1);
	    printf ("GB configuration saved on \"%s\".\n\n",
		    fn_config_gb);
	}
    }
    return(0);
} /* end main() */


/* save configuration in extended replicas */
void write_config (char filename[], int np, int GB_only)
{
    int i, j, k, l;
    FILE *config;
    config = fopen (filename, "w");
    fprintf (config, "Number of particles = %d\n\n",
	     np*ex[0]*ex[1]*ex[2]);
    fprintf (config, "H(1,1) = %f A\n", ex[0]*H[0][0]);
    fprintf (config, "H(1,2) = %f A\n", ex[0]*H[0][1]);
    fprintf (config, "H(1,3) = %f A\n", ex[0]*H[0][2]);
    fprintf (config, "nc(1)  = %d\n\n", ex[0]*nc[0]);
    fprintf (config, "H(2,1) = %f A\n", ex[1]*H[1][0]);
    fprintf (config, "H(2,2) = %f A\n", ex[1]*H[1][1]);
    fprintf (config, "H(2,3) = %f A\n", ex[1]*H[1][2]);
    fprintf (config, "nc(2)  = %d\n\n", ex[1]*nc[1]);
    fprintf (config, "H(3,1) = %f A\n", ex[2]*H[2][0]);
    fprintf (config, "H(3,2) = %f A\n", ex[2]*H[2][1]);
    fprintf (config, "H(3,3) = %f A\n", ex[2]*H[2][2]);
    fprintf (config, "nc(3)  = %d\n\n", ex[2]*nc[2]);
    fprintf (config, "TP Mass(amu)     sx        sy        sz ");
    fprintf (config, "    sx1(1/ns)  sy1(1/ns)  sz1(1/ns)\n");
    for (i=0; i<npr; i++)
	if ( (tag[i]>0) || ((tag[i]==0) && (!GB_only)) )
	    for (j=0; j<ex[0]; j++)
		for (k=0; k<ex[1]; k++)
		    for (l=0; l<ex[2]; l++)
			fprintf (config, "%s %f  %f  %f  %f  %.6f  %.6f  "
				 "%.6f\n", tpm[i], massm[i],
				 (sm[3*i]   + j) / ex[0],
				 (sm[3*i+1] + k) / ex[1],
				 (sm[3*i+2] + l) / ex[2], 0., 0., 0.);
    fclose (config);
    return;
} /* end write_config() */


/* ds() and r2() both take first argument as origin */
void ds (double s1[], double s2[], double d[])
{
    d[0] = s2[0] - s1[0];
    d[1] = s2[1] - s1[1];
    d[2] = s2[2] - s1[2];
    /* unique j image in H-parallelepiped */
    while (d[0] > 0.5)  d[0]--;
    while (d[0] < -0.5) d[0]++; 
    while (d[1] > 0.5)  d[1]--;
    while (d[1] < -0.5) d[1]++;
    while (d[2] > 0.5)  d[2]--;
    while (d[2] < -0.5) d[2]++; 
    return;
} /* end ds() */


/* (r[0],r[1],r[2]) is x_{ji}, r[3] is |x_{ji}|^2 */
/* r[] should be a cached global for efficiency.  */
static double r2 (double si[], double sj[])
{
    double d[3];
    ds (si, sj, d);
    r[0] = d[0]*H[0][0] + d[1]*H[1][0] + d[2]*H[2][0];
    r[1] = d[0]*H[0][1] + d[1]*H[1][1] + d[2]*H[2][1];
    r[2] = d[0]*H[0][2] + d[1]*H[1][2] + d[2]*H[2][2];
    return (r[3] = r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
} /* end r2() */


void waitfor (double seconds)
{
    struct timeval time;
    time.tv_sec = floor(seconds);
    time.tv_usec = rint((seconds-time.tv_sec)*1E6);
    select (0, NULL, NULL, NULL, &time);
    return;
} /* end waitfor() */


/* return pseudo random number on (0,1) */
double frandom()
{
    return ((random()+0.5)/pow(2.,31.));
} /* end frandom() */


void veceqv (double a[3], double b[3])
{
    b[0] = a[0];
    b[1] = a[1];
    b[2] = a[2];
    return;
} /* end veceqv() */


void vecadd (double a[3], double b[3], double c[3])
{
    c[0] = a[0] + b[0];
    c[1] = a[1] + b[1];
    c[2] = a[2] + b[2];
    return;
} /* end vecadd() */


void vecminus (double a[3], double b[3], double c[3])
{
    c[0] = a[0] - b[0];
    c[1] = a[1] - b[1];
    c[2] = a[2] - b[2];
    return;
} /* end vecminus() */


double dot (double a[3], double b[3])
{
    return (a[0]*b[0]+a[1]*b[1]+a[2]*b[2]);
} /* end dot() */


double normalize (double a[3])
{
    double length = sqrt(dot(a,a));
    if ( length < TINY )
    {
	printf ("normalize: length = %e < TINY = %e\n",
		length, TINY);
	exit(1);
    }
    a[0] /= length;
    a[1] /= length;
    a[2] /= length;
    return (length);
} /* end normalize() */


void cross (double a[3], double b[3], double c[3])
{
    c[0] = a[1]*b[2] - a[2]*b[1];
    c[1] = a[2]*b[0] - a[0]*b[2];
    c[2] = a[0]*b[1] - a[1]*b[0];
    return;
} /* end cross() */


void random_unit_vec (double x[3])
{
    int i;
    double norm;
    while(1)
    {
	x[0] = frandom() - 0.5;
	x[1] = frandom() - 0.5;
	x[2] = frandom() - 0.5;
	if ( (norm=dot(x,x)) < 0.25 ) break;
    }
    norm = sqrt(norm);
    x[0] /= norm;
    x[1] /= norm;
    x[2] /= norm;
    return;
} /* end random_unit_vec() */


/* rotate b in a[0-2] by a[3] degrees; a[0-2] should be normalized */
void rotate (double a[4], double b[3], double c[3])
{
    double d[3], e[3], dr[3], er[3], alpha, beta, ba, bd, be;
    random_unit_vec (d);
    cross (a, d, e);
    normalize (e);
    cross (e, a, d);
    alpha = cos(a[3]/180.*PI);
    beta  = sin(a[3]/180.*PI);
    dr[0] = alpha*d[0] + beta*e[0];
    dr[1] = alpha*d[1] + beta*e[1];
    dr[2] = alpha*d[2] + beta*e[2];
    er[0] = -beta*d[0] + alpha*e[0];
    er[1] = -beta*d[1] + alpha*e[1];
    er[2] = -beta*d[2] + alpha*e[2];
    ba = dot(b,a); bd = dot(b,d); be = dot(b,e); 
    c[0] = ba*a[0] + bd*dr[0] + be*er[0];
    c[1] = ba*a[1] + bd*dr[1] + be*er[1];
    c[2] = ba*a[2] + bd*dr[2] + be*er[2];
    return;
} /* end rotate() */


void mateqv(double A[3][3], double B[3][3])
{
    B[0][0] = A[0][0];
    B[0][1] = A[0][1];
    B[0][2] = A[0][2];
    B[1][0] = A[1][0];
    B[1][1] = A[1][1];
    B[1][2] = A[1][2];
    B[2][0] = A[2][0];
    B[2][1] = A[2][1];
    B[2][2] = A[2][2];
    return;
} /* end mateqv() */


void matran(double A[3][3], double B[3][3])
{
    B[0][0] = A[0][0];
    B[0][1] = A[1][0];
    B[0][2] = A[2][0];
    B[1][0] = A[0][1];
    B[1][1] = A[1][1];
    B[1][2] = A[2][1];
    B[2][0] = A[0][2];
    B[2][1] = A[1][2];
    B[2][2] = A[2][2];
    return;
} /* end matran() */


void matadd(double A[3][3], double B[3][3], double C[3][3])
{
    C[0][0] = A[0][0] + B[0][0];
    C[0][1] = A[0][1] + B[0][1];
    C[0][2] = A[0][2] + B[0][2];
    C[1][0] = A[1][0] + B[1][0];
    C[1][1] = A[1][1] + B[1][1];       
    C[1][2] = A[1][2] + B[1][2];       
    C[2][0] = A[2][0] + B[2][0];           
    C[2][1] = A[2][1] + B[2][1];      
    C[2][2] = A[2][2] + B[2][2];
    return;
} /* end matadd() */


void matsub(double A[3][3], double B[3][3], double C[3][3])
{
    C[0][0] = A[0][0] - B[0][0];
    C[0][1] = A[0][1] - B[0][1];
    C[0][2] = A[0][2] - B[0][2]; 
    C[1][0] = A[1][0] - B[1][0];         
    C[1][1] = A[1][1] - B[1][1];       
    C[1][2] = A[1][2] - B[1][2];       
    C[2][0] = A[2][0] - B[2][0];           
    C[2][1] = A[2][1] - B[2][1];      
    C[2][2] = A[2][2] - B[2][2];
    return;
} /* end matsub() */


void matmul(double A[3][3], double B[3][3], double C[3][3])
{
    C[0][0] = A[0][0]*B[0][0] + A[0][1]*B[1][0] + A[0][2]*B[2][0];
    C[0][1] = A[0][0]*B[0][1] + A[0][1]*B[1][1] + A[0][2]*B[2][1];
    C[0][2] = A[0][0]*B[0][2] + A[0][1]*B[1][2] + A[0][2]*B[2][2];
    C[1][0] = A[1][0]*B[0][0] + A[1][1]*B[1][0] + A[1][2]*B[2][0];
    C[1][1] = A[1][0]*B[0][1] + A[1][1]*B[1][1] + A[1][2]*B[2][1];
    C[1][2] = A[1][0]*B[0][2] + A[1][1]*B[1][2] + A[1][2]*B[2][2];
    C[2][0] = A[2][0]*B[0][0] + A[2][1]*B[1][0] + A[2][2]*B[2][0];
    C[2][1] = A[2][0]*B[0][1] + A[2][1]*B[1][1] + A[2][2]*B[2][1];
    C[2][2] = A[2][0]*B[0][2] + A[2][1]*B[1][2] + A[2][2]*B[2][2];
    return;
} /* end matmul() */


/* c = aB: storage of c could be a */
void vec_mul_mat (double a[3], double B[3][3], double c[3])
{
    double d[3];
    d[0] = a[0]*B[0][0] + a[1]*B[1][0] + a[2]*B[2][0];
    d[1] = a[0]*B[0][1] + a[1]*B[1][1] + a[2]*B[2][1];
    d[2] = a[0]*B[0][2] + a[1]*B[1][2] + a[2]*B[2][2];
    veceqv (d, c);
    return;
} /* end vec_mul_mat() */

    
double matinv(double A[3][3], double B[3][3])
{
    double D11,D22,D33,D12,D21,D13,D31,D32,D23,determinant;
    D11=A[1][1]*A[2][2]-A[1][2]*A[2][1];
    D22=A[2][2]*A[0][0]-A[2][0]*A[0][2];
    D33=A[0][0]*A[1][1]-A[0][1]*A[1][0];
    D12=A[1][2]*A[2][0]-A[1][0]*A[2][2];
    D23=A[2][0]*A[0][1]-A[2][1]*A[0][0];
    D31=A[0][1]*A[1][2]-A[0][2]*A[1][1];
    D13=A[1][0]*A[2][1]-A[2][0]*A[1][1];
    D21=A[2][1]*A[0][2]-A[0][1]*A[2][2];
    D32=A[0][2]*A[1][0]-A[1][2]*A[0][0];
    determinant=A[0][0]*D11+A[0][1]*D12+A[0][2]*D13;
    B[0][0]=D11/determinant;
    B[1][1]=D22/determinant;
    B[2][2]=D33/determinant;
    B[0][1]=D21/determinant;
    B[1][2]=D32/determinant;
    B[2][0]=D13/determinant;
    B[1][0]=D12/determinant;
    B[2][1]=D23/determinant;
    B[0][2]=D31/determinant;
    return (determinant);
} /* end matinv() */