/*
 * General include file for MMTK C modules.
 *
 * Written by Konrad Hinsen
 * last revision: 2002-6-14
 */

#ifndef MMTK_H

#include "Python.h"
#include "Numeric/arrayobject.h"

/* MinGW doesn't have this */
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/* Configuration parameters */

/* max. number of threads that can try to access the same universe
   in parallel (energy evaluation threads don't count) */
#define MMTK_MAX_THREADS 10

/* Type definitions */

typedef double vector3[3];
typedef double tensor3[3][3];


/* Useful general macros */

#define sqr(x) ((x)*(x))
#define cube(x) ((x)*(x)*(x))
#define vector_add(v1, v2, f) \
   { int i; for (i=0; i<3; i++) (v1)[i] += f*(v2)[i]; }
#define vector_length_sq(r) (sqr(r[0])+sqr(r[1])+sqr(r[2]))
#define vector_length(r) (sqrt(vector_length_sq(r)))
#define dot(v1, v2) ((v1)[0]*(v2)[0] + (v1)[1]*(v2)[1] + (v1)[2]*(v2)[2])
#define cross(v, v1, v2) \
        { v[0] = (v1)[1]*(v2)[2]-(v1)[2]*(v2)[1]; \
          v[1] = (v1)[2]*(v2)[0]-(v1)[0]*(v2)[2]; \
          v[2] = (v1)[0]*(v2)[1]-(v1)[1]*(v2)[0];}

#define vector_scale(v, f) { v[0]*=f; v[1]*=f; v[2]*=f; }
#define vector_copy(v1, v2) { (v1)[0]=(v2)[0]; (v1)[1]=(v2)[1]; \
                              (v1)[2]=(v2)[2]; }
#define vector_changesign(v) { v[0]=-v[0]; v[1]=-v[1]; v[2]=-v[2]; }

#define tensor_product(t, v1, v2, f) \
   { int i,j; for (i=0;i<3;i++) for (j=0;j<3;j++) t[i][j]=f*(v1)[i]*(v2)[j]; }
#define symmetric_tensor_product(t, v1, v2, f) \
   { int i,j; for (i=0;i<3;i++) for (j=0;j<3;j++) \
		     t[i][j] = f*((v1)[i]*(v2)[j]+(v1)[j]*(v2)[i]); }
#define tensor_copy(t1, t2) \
   { int i,j; for (i=0;i<3;i++) for (j=0;j<3;j++) t1[i][j] = t2[i][j]; }
#define tensor_scale(t1, f) \
   { int i,j; for (i=0;i<3;i++) for (j=0;j<3;j++) t1[i][j] *= f; }
#define tensor_add(t1, t2, f) \
   { int i,j; for (i=0;i<3;i++) for (j=0;j<3;j++) t1[i][j] += f*t2[i][j]; }
#define tensor_changesign(t) { vector_changesign(t[0]); \
                               vector_changesign(t[1]); \
                               vector_changesign(t[2]); }
#define tensor_transpose(t) { double temp; \
                              temp=t[0][1]; t[0][1]=t[1][0]; t[1][0]=temp; \
			      temp=t[0][2]; t[0][2]=t[2][0]; t[2][0]=temp; \
			      temp=t[1][2]; t[1][2]=t[2][1]; t[2][1]=temp; }


/* Constants */

/* Numbers larger than 'undefined_limit' are considered to be undefined in
 * certain situations. This is used e.g. to represent unknown positions.
 * The value 'undefined' is used to mark numbers as undefined. The difference
 * between these two values should make sure that undefined values are
 * recognized across floating point representation boundaries.
 */
#define undefined_limit 1.e30;
#define undefined 1.e31;


#define MMTK_H
#endif
