#ifndef mymath_h
#define mymath_h
/** \file mymath.h
 * Header file for \ref mymath.c
 * 
 * Contains declarations of the functions implemented in \ref mymath.c, shorter
 * definitions of mathematical constants such as PI, square root of 2 SQRT2
 * and hard-coded constructs for calculating maximum MAX(), row majors CalcRowMajor2D(),
 * CalcRowMajor3D(), determinants RDET2(), RDET3(),  scalar products real RSP3()
 * and complex CSP3re(), CSP3im(), euclidian norm real RNORMSQ3() and complex CNORMSQ3(),
 * complex multiplication CCMULTre(), CCMULTim(), complex multiplication with scalar
 * RCMULTre(), RCMULTim().
 * 
  Project: ParallelCarParrinello
  Jan Hamaekers
  2000

  File: mymath.h
  $Id: mymath.h,v 1.15 2007-03-29 13:35:51 foo Exp $
*/

// use double precision fft when we have it
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#if defined _BSD_SOURCE || defined _XOPEN_SOURCE
//! short form for pi from math.h
# define PI M_PI
//! short form for square root of 2 from math.h
# define SQRT2 M_SQRT2
#else /* generische Form */
//! short form for pi
# define PI (acos(-1.0))
//! short form for square root of 2
# define SQRT2 (sqrt(2.))
#endif

#include "defs.h"

#ifdef HAVE_DFFTW_H
#include "dfftw.h"
#else
#include "fftw.h"
#endif

#define MAX(a,b) ((a) > (b) ? (a) : (b))															//!< returns maximum of a or b
#define CalcRowMajor3D(R0,R1,R2,N0,N1,N2) ((R2)+(N2)*((R1)+(N1)*(R0)))//!< calculates row major of 3x3 matrix
#define CalcRowMajor2D(R0,R1,N0,N1) ((R1)+(N1)*(R0))									//!< calculates row major of 2x2 matrix
#define RSP3(a,b) ((a)[0]*(b)[0] + (a)[1]*(b)[1] + (a)[2]*(b)[2])			//!< scalar product of two 3-dim vectors
#define RNORMSQ3(a) ((a)[0]*(a)[0] + (a)[1]*(a)[1] + (a)[2]*(a)[2])		//!< squared euclidian norm
#define RDET3(a) ((a)[0]*(a)[4]*(a)[8] + (a)[3]*(a)[7]*(a)[2] + (a)[6]*(a)[1]*(a)[5] - (a)[2]*(a)[4]*(a)[6] - (a)[5]*(a)[7]*(a)[0] - (a)[8]*(a)[1]*(a)[3])	//!< hard-coded determinant of a 3x3 matrix
#define RDET2(a0,a1,a2,a3) ((a0)*(a3)-(a1)*(a2))											//!< hard-coded determinant of a 2x2 matrix
#define CCMULTre(a,b) ((a).re*(b).re - (a).im*(b).im)									//!< real part of a complex multiplication
#define CCMULTim(a,b) ((a).re*(b).im + (a).im*(b).re)									//!< imaginary part of a complex multiplication
#define RCMULTre(a,b) ((a).re*(b))																		//!< real part of a complex number scaled by a real number
#define RCMULTim(a,b) ((a).im*(b))																		//!< imaginary part of a complex number scaled by a real number
#define CSP3re(a,b) (CCMULTre((a)[0],(b)[0]) + CCMULTre((a)[1],(b)[1]) + CCMULTre((a)[2],(b)[2])) //!< real part of a scalar product of two 3x3 complex vectors
#define CSP3im(a,b) (CCMULTim((a)[0],(b)[0]) + CCMULTim((a)[1],(b)[1]) + CCMULTim((a)[2],(b)[2])) //!< imaginary part of a scalar product of two 3x3 complex vectors
#define CNORMSQ3(a) ((a).re[0]*(a).re[0] + (a).re[1]*(a).re[1] + (a).re[2]*(a).re[2] + (a).im[0]*(a).im[0] + (a).im[1]*(a).im[1] + (a).im[2]*(a).im[2])	//!< square of complex euclidian norm

inline double tpow(double, int);
inline int Rest(int n, int m);
inline void RTranspose3(double *A);
inline void RMatMat33(double C[NDIM*NDIM], const double A[NDIM*NDIM], const double B[NDIM*NDIM]);
inline void RMat33Vec3(double C[NDIM], const double M[NDIM*NDIM], const double V[NDIM]);
inline int RMatReci3(double B[NDIM_NDIM], const double A[NDIM_NDIM]);
inline void RVec3Mat33(double C[NDIM], const double V[NDIM], const double M[NDIM*NDIM]);
inline void VP3(double V[NDIM], double A[NDIM], double B[NDIM]);
/* Skalarprodukt */
inline double SP(const double *a, const double *b, const int n);
/* Multiplikation mit Skalar */
inline void SM(double *a, const double c, const int n);
/* Nullvektor erzeugen */
inline void NV(double *a, int n);
inline double dSum(int n, double *dx, int incx);
inline double Simps(int n, double *f, double h);
inline double derf(double x);
/* Initialisiere a array[3] mit b - c Orte mit periodisch */
double Dist(const double *a, const double *b, const int n);
inline void SetArrayToDouble0(double *a, int n);
void PrintCMat330(fftw_complex M[NDIM_NDIM]);
void PrintRMat330(fftw_real M[NDIM_NDIM]);
void PrintCVec30(fftw_complex M[NDIM]);
void PrintRVec30(fftw_real M[NDIM]);
void RotateToAlign(fftw_real Q[NDIM_NDIM], fftw_real matrix[NDIM_NDIM], fftw_real vector[NDIM]);
#endif