#ifndef INTERACT_INCLUDED
#define INTERACT_INCLUDED

#include "basis.h"

/*
   Define number and order of coupling constants
   
   In principle, the assignment of numbers to these variables
   should be arbitrary:  work like an enumeration
*/

#if HINDEX==1
#define NPARAMS 7 /* number of coupling constants */
#define MM 0
#define GGN 1
#define LAMBDA_1 2
#define LAMBDA_2 3
#define LAMBDA_3 4
#define TAU_1 5
#define TAU_2 6
#elif HINDEX==2
#define NPARAMS 16 /* number of coupling constants */
#define MM 0
#define GGN 1
#define BETA 2
#define LAMBDA_1 3
#define LAMBDA_2 4
#define LAMBDA_3 5
#define LAMBDA_4 6
#define LAMBDA_5 7
#define TAU_1 8
#define TAU_2 9
#define KAPPA_S 10
#define KAPPA_A 11
#define MU_F 12
#define MMF_0 13
#define MMF_1 14
#define MMF_2 15
#endif

/* Flag for choices of phase factor convention for
   nonzero P_1.  The 1 and 3 choices are probably faster,
   and 2 or 3 are necessary for nonzero winding number.
   Setting it to be zero will cause the automatic 
   choice of 1 for n=0 and 3 for n>0.
   These are the fastest two choices.  */

#if HINDEX==1
#define PHASE_CONVENTION 0
#elif HINDEX==2
#define PHASE_CONVENTION 0
#endif

/* 
   Outside of the routines interact_*,
   one should not assume that there are 4-point interactions.
   This is not yet fully implimented
*/

#define GLUE_POINTS 4  /* Number of particles in glue interactions */
#define QUARK_POINTS 4 /* Number of particles in quark interactions */

struct {
  int nptrunc;         /* number of particles in truncation */
  unsigned int kt;     /* The total momentum.
		          For antiperiodic links, this is twice the momentum.
                          Otherwise, it is the momentum rounded down. */
  int  momflag;        /* flag for nonzero momentum/angles */ 
  element apperp[HINDEX+2];   /* angles and momenta, length HINDEX+2 */
} working;

/* --------------- Macro Definitions --------------------------------*/

#if HINDEX==2
/*
  Macros for comparing two vectors.
*/

#define EQUAL(A,B) (A[0]==B[0] && A[1]==B[1])
#define MINUS(A,B) (A[0]==-B[0] && A[1]==-B[1])
#define ORTHO(A,B) (A[0]*B[0]== -A[1]*B[1])
#endif

/*
  macros for complex times real multiply    

        ZTIMES(A,B,C):  A += B C     
       ZTIMESC(A,B,C):  A += B* C (complex conjugate of B)
          ZZPROD(A,B):  Z = A B

  Z, A and B are pointers while C is a number
*/

#define ZTIMES(A,B,C) do{A->re+=(B)->re*(C); A->im+=(B)->im*(C);}while(0)
#define ZTIMESC(A,B,C) do{A->re+=(B)->re*(C); A->im-=(B)->im*(C);}while(0)
#define ZZPROD(Z,A,B) do{(Z)->re=(A)->re*(B)->re-(A)->im*(B)->im; \
                         (Z)->im=(A)->im*(B)->re+(A)->re*(B)->im; \
                      }while(0)

/* convert quark and gluon momenta to real */

#define QE(A) ((element) ((A)&KMASK)+(QUARK_PERIODIC?0.0:0.5))
#define GE(A) ((element) ((A)&KMASK)+(GLUE_PERIODIC?0.0:0.5))

/* 
   order the four momenta for an interaction 
   The two symmetries are:  L1<->R1, L2<->R2 and L1<->R2, R1<->L2.
*/

#define sort_k(KA,KB,L1,L2,R1,R2) do{           \
  partype a=(L1<R2?L1:R2), b=(R1<L2?R1:L2);     \
  if(a>=b){*KA=a; *KB=b; }else{ *KA=b; *KB=a; } \
}while(0)

/******************** from files: interact_*.c  **************************/

void interact_glueball(element *x, partype *left, int npl, partype *right, 
		   int npr, element *couplings);
void interact_glueballc(complex_element *x, partype *left, int npl, 
		    partype *right, int npr, element *couplings);
void interact_glueballcu(complex_element *x, partype *left, int npl, 
		    partype *right, int npr, element *couplings);
void interact_glueballcv(complex_element *x, partype *left, int npl, 
		    partype *right, int npr, element *couplings);
void interact_source(complex_element *x, partype *left, int npl, 
		     partype *right, int npr, element *couplings);
void print_source_interaction_memory(FILE *fp);
void free_source(void);  /* free memory used by source interactions */
char *coupling_string(int t); /* string for coupling */
void interact_heavy_light(complex_element *x, partype *left, int npl, 
		    partype *right, int npr, element *couplings);
void interact_meson(complex_element *x, partype *left, int npl, 
		     partype *right, int npr, element *couplings);

void interact_glue(element *x, partype *left, int npl, partype *right, 
		   int npr, element *couplings);
void interact_gluec(complex_element *x, partype *left, int npl, 
		    partype *right, int npr, element *couplings);
void interact_quark(complex_element *x, partype *left, int npl, 
		     partype *right, int npr, element *couplings, int ngl);
void interact_quarkc(complex_element *x, partype *left, int npl, 
		     partype *right, int npr, element *couplings, int ngl);
void interact_heavy(complex_element *x, partype *left, int npl, 
		     partype *right, int npr, element *couplings, 
		     partype added_k);

                 /* initialize parameters needed for the 
                    interaction calculations */
void init_params(int nptrunc, unsigned int kt, int momflag, element *apperp);

        /* find matrix element of hamiltonian, assuming cyclic permutations */
void makephaselist(element *,int); /*make a list of angles for interactionc */
void shuffle(full_basis *,element *, element *,int);       
                              /*order 0++* and 2++ by particle content */
void maketenslist(element llin, int ktin, element p_maxin, element betain);
                                   /* initialize longitudinal separations */
element selfinertia(element *); /* integrand for self-inertia calculations */
element seenorm(element *);  /* norm for improved 0->2 interactions */
element vinst(element *);  /* instantaneous improved 0->2 interaction */
element vtau(element *);   /* tau improved 0->2 interaction */
element vnewtau(element *);   /*  new tau improved 0->2 interaction */

#endif