#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include "include.h"
#include "improved.h"
#include "interact.h"

/*    
      Matrix elements for the longitudinal string tension using 
      an open string.
      See e-Print Archive: hep-th/9806231; revised interactions
      are discussed in e-Print Archive: hep-lat/9911035 

      The routine calculates v^+ P^-

      In addition, it is a function of the parameters:

        p_max = dimensionless link momentum cutoff:  
                P^+_link < v^+ sqrt(G^2 N) P_max

        kstep = p_max/K is the step size.

        ll = L sqrt{G^2 N}

      where K = kt/2 is the upper limit of the discretized momenta.
      The argument of long[] is the change in K of the state.
      
      memcmp(a,b,0) returns 0
      
      Test runs with n=2, 2-particle truncation suggest that
      it is better to not assume anything about large
      momentum of the wavefuntion for the purposes of calculating
      improved matrix elements.
      
      Future speed improvements may include recalculations
      of just the phase factors when possible.  Also, 
      calculations matrix elements for the 0->2 link field 
      interactions could be limited by ki+kj <= K. 

      The memory usage to store the matrix elements is pretty
      large.  It probably makes sense to calculate them dynamically,
      since they only occur on the ends of the strings.  This is
      what is done with the quark interactions.
      Also, I don't take any advantage of the symmetry of the matrices,
      and the code is rather ugly.
*/

size_t source_memory=0;    /* current memory use by source arrays, in bytes */
static element *ff[KKMAX];
static element *gg[KKMAX], *gg13[KKMAX];
static element *tau11[KKMAX], *tau02[KKMAX];
static element barepotential,kiki,b;

struct { 
  element p_max;      /* dimensionless cut-off on link P^+  */
  element beta;       /* beta for improved matrix elements */
  element ll;         /* sqrt(G^2 N) L longitudinal spacing */
  unsigned int kt;    /* twice the cut-off on discretized link momenta */
  element kstep;      /* 2.0*p_max/kt */
} longer;

#if IMPROVED >= 2
static element funpower(element x, element b){
  return (pow(x+0.5,b)-pow(x-0.5,b))/b;
}
#endif

void maketenslist(element ll, int kt, element p_max, element beta){
  int i,j,l;
#if IMPROVED>=2
  integer ierr;
  /* err is the error of the integrator */
  element err=1.0e-10,zero=0.0,ans;
  element ansinst,ansvtau;
#endif
  element x,ki,kj,m;
  size_t temp;
  const element pi=3.1415926535897932384626433832795;

  /* initialize global parameters */
  longer.ll=ll; longer.p_max=p_max; longer.kt=kt; longer.beta=beta;

  if(kt>2*KKMAX){
    fprintf(stderr,__FILE__":  kt in maketenslist too large, exiting\n");
    return;
  }
  longer.kstep=2.0*p_max/((element) kt);

  /*  This is for the "bare" potential. */

  barepotential=0.25*fabs(ll);

  /* All interactions (except for self-inertias)
     The sqrt(funpower(ki... in the demoninator
     is the normalization of the x^beta wavefunction */
  
  free_source();  /* we malloc all the memory each time */

  for(i=0; 2*i+(GLUE_PERIODIC?0:1)<=kt; i++){
    ki=(element) i+(GLUE_PERIODIC?0:0.5);
    
    /* allocate memory */
      temp=(kt/2+1)*sizeof(element);
      source_memory+=4*temp;
      gg[i]=malloc(temp);
      if(gg[i]==NULL){
	fprintf(stderr,__FILE__":%i  out of memory\n",__LINE__);
	return;
      }
      gg13[i]=malloc(temp);
      if(gg13[i]==NULL){
	fprintf(stderr,__FILE__":%i  out of memory\n",__LINE__);
	return;
      }
      tau11[i]=malloc(temp);
      if(tau11[i]==NULL){
	fprintf(stderr,__FILE__":%i  out of memory\n",__LINE__);
	return;
      }
      tau02[i]=malloc(temp);
      if(tau02[i]==NULL){
	fprintf(stderr,__FILE__":%i  out of memory\n",__LINE__);
	return;
      }

    /* Do one triangle of matrix */

    for(j=0; 2*j+(GLUE_PERIODIC?0:1)<=kt; j++){
      kj=(element) j+(GLUE_PERIODIC?0:0.5);
#if IMPROVED<2 /* DLCQ matrix elements  */
      tau11[i][j] = 1.0/(8.0*longer.kstep*pi*sqrt(ki*kj));
#else /* improved matrix elements */
      kiki=ki; b=kj;
      /* Here we compute ans, the norm^2 of our basis states. */ 
      x=(ki-0.5)/(ki+kj); m=(ki+0.5)/(ki+kj);
      DGAUS8(SEENORM_TO_FORTRAN,&x,&m,&err,&ans,&ierr);
#if 0 /* print norm^2 of basis states */
      printf("ki=%g kj=%g beta=%.16g norm=%.16g\n",ki,kj,beta,ans);
#endif
      tau11[i][j] = funpower(ki,beta+0.5)*funpower(kj,beta+0.5)/
	(8.0*longer.kstep*pi*sqrt(funpower(ki,2.0*beta+1.0)*
			  funpower(kj,2.0*beta+1.0)));
#endif  /* improved matrix elements */
      if(i==j){
	tau11[i][j] *= 0.5*ll*longer.kstep;
      } else 
	tau11[i][j] *= sin(0.5*ll*longer.kstep*(ki-kj))/(ki-kj);
 
#if IMPROVED<2 /* DLCQ matrix elements  */
      tau02[i][j] = sin(0.5*ll*longer.kstep*(ki+kj))/
	(8.0*pi*longer.kstep*sqrt(ki*kj)*(ki+kj));
#else /* improved matrix elements */
      DGAUS8(VNEWTAU_TO_FORTRAN,&x,&m,&err,&ansvtau,&ierr);
      if(i==j && i==0)
        ansvtau += 2.0/(beta+0.5);
      else if(i==0)
	ansvtau += log((kj+0.5)/(kj-0.5))/((beta+0.5)*pow(kj+0.5,beta+0.5));
      else if(j==0)
	ansvtau += log((ki+0.5)/(ki-0.5))/((beta+0.5)*pow(ki+0.5,beta+0.5));
      tau02[i][j] = ansvtau/(8.0*pi*longer.kstep*sqrt(ans));
      /* At the boundary, Taylor expand the sine... */
      if(i==0 && j==0)
	tau02[i][j] *= 0.5*ll*longer.kstep;
      else
	tau02[i][j] *= sin(0.5*ll*longer.kstep*(ki+kj));
#endif  /* improved matrix elements */
      if(i!=j){
#if IMPROVED<2 /* DLCQ matrix elements  */
	gg[i][j] = -(ki+kj)/(8.0*longer.kstep*pi*sqrt(ki*kj)*(ki-kj)*(ki-kj));
#else /* improved matrix elements */
	gg[i][j] = -(funpower(ki,beta+1.5)*funpower(kj,beta+0.5)+
			funpower(ki,beta+0.5)*funpower(kj,beta+1.5))/
	  (8.0*longer.kstep*pi*sqrt(funpower(ki,2.0*beta+1.0)*
			    funpower(kj,2.0*beta+1.0))*(ki-kj)*(ki-kj));
#endif  /* improved matrix elements */
#if IMPROVED<2 /* DLCQ matrix elements  */
	gg13[i][j] = (ki-kj)/(8.0*longer.kstep*pi*sqrt(ki*kj)*(ki+kj)*(ki+kj));
#else /* improved matrix elements */
	DGAUS8(VINST_TO_FORTRAN,&x,&m,&err,&ansinst,&ierr);
	if(i==0)
	  ansinst += log((kj+0.5)/(kj-0.5))/((beta+0.5)*pow(kj+0.5,beta+0.5));
	else if(j==0)
	  ansinst -= log((ki+0.5)/(ki-0.5))/((beta+0.5)*pow(ki+0.5,beta+0.5));
	/* The i is the second monentum and j is the first, thus a sign: */
	gg13[i][j] = -ansinst/(8.0*longer.kstep*pi*sqrt(ans));
#endif  /* improved matrix elements */
      }
    }

    /* Do Diagonal parts */
    gg[i][i]=0.0;
    gg13[i][i]=0.0;
  }

  /* Calculate the self-intertias  */

  for(l=0; l<=kt/2; l++){
    temp=((kt-2*l)/2+1)*sizeof(element);
    source_memory+=temp;
    ff[l]=malloc(temp);
    if(ff[l]==NULL){
      fprintf(stderr,__FILE__":%i  out of memory\n",__LINE__);
	return;
    }

    for(j=0; 2*j+(GLUE_PERIODIC?0:1)<=kt-2*l; j++){
      kj=(element)j+(GLUE_PERIODIC?0:0.5);
      m=((element) l+kj)/kj;
#if IMPROVED < 2 /* DLCQ matrix elements */
      /* This is currently done with Burkardt's technique
	 (That is, the ll exponential and the x^beta are not included.)  */
      /* note that gg[i][j] is negative. */
      for(i=0, x=0.0; i<=l+j; i++){
	if(GLUE_PERIODIC && (i==0 || j==0))continue;
	x -= gg[i][j];
      }
      if(l>0)
	ff[l][j] = x+sqrt(m)/(4.0*pi*longer.kstep*(element) l);
      else
	ff[l][j] = x+1.0/(8.0*pi*longer.kstep*kj);
#if 0 /*debugging print */
      printf("ff calculation l=%i j=%i x = %e\n",l,j,x);
#endif 

#else            /* improved matrix elements */

      for(i=0, x=0.0; i<=l+j; i++){
	ki = (element)i+0.5;
	/* note that gg[i][j] is negative.    */
	x -= pow(ki/kj,beta)*gg[i][j];
      }
      b=kj*funpower(kj,beta+0.5)/funpower(kj,beta+1.5);
      DGAUS8(SELFINERTIA_TO_FORTRAN,&zero,&m,&err,&ans,&ierr);
#if 0 /*debugging print */
      printf("j=%i integration answer = %e\n",j,ans);
      printf("     sum answer* 8.0 pi longer.kstep = %e\n",
	     8.0*x*pi*longer.kstep);
#endif 
      if(m>1.0+err)
	ans += (-1.*(1 + b)*m)/(-1. + m) + 
	  (1.*pow(m,0.5 + beta))/(0.5 + 1.*beta) + 
	  (-0.25 + b*(0.25 + 0.5*beta) + 0.5*beta)*
	  log(pow(-1. + m,2));
      else
	ans += -1. - 1.*b + 2./(1. + 2.*beta);
      ff[l][j] = (x-ans*funpower(kj,beta+1.5)/
	    (8.0*pi*longer.kstep*pow(kj,1.5)*sqrt(funpower(kj,2.0*beta+1.0))));
#endif /* end improved matrix elements */
#if 0 /*debugging print */
      printf("b=%21.14e m=%21.14e beta=%21.14e\n",b,m,beta);
#endif
#if 0 /* DEBUGGING print, print out self-inertias */
      printf(" self inertia:  kstep*ff[%i][%i]=%21.14e\n",
	     l,j,longer.kstep*ff[l][j]);
#endif
    }

#if 0  /* DEBUGGING print: plot integrand in region near z=1 */
    printf(" Integrand near z=1:\n");
    for(ki=1.0-1.5e-4; ki<1.0+1.5e-4; ki+=2.0*1.5e-4/200.)
      printf(" %21.14e %21.14e\n",ki,selfinertia(&ki));
#endif
  }

#if 0  /* debugging print */
  for(i=0; 2*i+(GLUE_PERIODIC?0:1)<=kt && i<5; i++)
    for(j=0; 2*j+(GLUE_PERIODIC?0:1)<=kt && j<5; j++)
      printf("   tau11[%i][%i] = %g and %g\n",i,j,tau11[i][j],tau11[i][j]);
#endif
}

/*
  find phase factor for given heavy interaction
*/

void source_phase(complex_element *z, partype *left, int npl,
		partype *right, int npr){
  int i;
  element phase=0;
  for(i=0; i<npl; i++)phase += GE(left[i]);
  for(i=0; i<npr; i++)phase -= GE(right[i]);
  /* working.kt is twice the momentum */
  phase *= 0.5*longer.ll*longer.kstep; 
  z->re=cos(phase);
  z->im=sin(phase);
}


/*  This is the pole-free function to be integrated over
    numerically in order to obtain the self-inertias.

    There can be roundoff error problems near the z=1 pole.
    The two choices of integrand are based on the roundoff
    errors in the mathematica program.  The first choice is
    a linear Taylor expansion about the pole and the 
    second is the full expression.

    It uses the external variables b and beta. */

element selfinertia(element *z){
  if(fabs((*z)-1.0)<1.0e-4)
    return 0.02083333333333333*(6.*(1. + 2.*longer.beta)*
      (-5. + 2.*longer.beta + b*(-1. + 2.*longer.beta)) + 
     (-1. + 4.*pow(longer.beta,2))*(-9. + 2.*longer.beta + b*(-3. + 2.*longer.beta))*
      (-1. + *z));
  else
    return (0.5*(-3. - 2.*longer.beta*(-1. + *z) + *z + 4.*pow(*z,0.5 + longer.beta) - 
       2.*pow(*z,1.5 + longer.beta) - 
       1.*b*(1. + 2.*longer.beta*(-1. + *z) + *z - 2.*pow(*z,0.5 + longer.beta))))/
   pow(-1. + *z,2);
}

/* It is useful to define the following quantities when 
   computing improved matrix elements for the 0->2 link
   interactions.  If we make the substitution 
       ki = x * k  and kj = (1-x) * k 
   Then the boundaries in the k integration are given
   by the following:     

   The OR makes it handle the *x=0 correctly. */

#define MAXK ((*x*(ki+kj-1.0)<(ki-0.5)||*x==1.0)?\
              (ki-0.5)/(*x):(kj-0.5)/(1.0-*x))

#define MINK ((*x*(ki+kj+1.0)<(ki+0.5))?\
              (kj+0.5)/(1.0-*x):(ki+0.5)/(*x))

/* this is a function to be integrated over for the
   norm of a state                              */

element seenorm(element *x){
  element ki=kiki, kj=b;
  return 0.5*pow(*x *(1.0- *x),2*longer.beta)*(pow(MINK,2)-pow(MAXK,2));
}

/* this is a pole-free function to be integrated over for the
   heavy-light instantaneous interaction               */

element vinst(element *x){
element ki=kiki, kj=b, z;
  if(ki<0.75 && kj<0.75)return 0.0;
  z = pow(*x *(1.0- *x),longer.beta-0.5)*(1-2.0* *x)*log(MINK/MAXK);
  if(ki<0.75)z -= pow(*x,longer.beta-0.5)*log((kj+0.5)/(kj-0.5));
  if(kj<0.75)z += pow(1.0- *x,longer.beta-0.5)*log((ki+0.5)/(ki-0.5));
  return z;
}

/* this is a pole-free function to be integrated over for the
   new 0->2 links tau interaction               */

element vnewtau(element *x){
element ki=kiki, kj=b, z;
  if(ki<0.75 && kj<0.75){
   z = pow(*x *(1.0- *x),longer.beta-0.5)*(MINK-MAXK)
      -pow(*x,longer.beta-0.5) - pow(1.0- *x,longer.beta-0.5);
  } else {
    z = pow(*x *(1.0- *x),longer.beta-0.5)*log(MINK/MAXK);
    if(ki<0.75)z -= pow(*x,longer.beta-0.5)*log((kj+0.5)/(kj-0.5));
    if(kj<0.75)z -= pow(1.0- *x,longer.beta-0.5)*log((ki+0.5)/(ki-0.5));
  }
  return z;
}

/* this is a pole-free function to be integrated over for the
   heavy-light tau interaction               */

element vtau(element *x){
  element ki=kiki, kj=b, z;
  z = pow(*x *(1.0- *x),longer.beta-0.5)*(MINK-MAXK);
  if(ki<0.75)z -= pow(*x,longer.beta-0.5);
  if(kj<0.75)z -= pow(1.0- *x,longer.beta-0.5);
  return z;
}

/*
  This calculates the matrix elements of the Hamiltonian v^+ P^- for 
  an open string, using variables defined above. 

  At the end, the answer is divided by sqrt(G^2 N).
*/

void interact_source(complex_element *x, partype *left, int npl, 
		     partype *right, int npr, element *couplings){
  int i,ktl=0,ktr=0;
  partype left_flip[NPMAX],right_flip[NPMAX];

  x->re=0.0; x->im=0.0;

#if 1      /* for debugging only -- turn off glue-glue interactions */

  for(i=0; i<npl; i++)ktl += left[i]&KMASK;
  for(i=0; i<npr; i++)ktr += right[i]&KMASK;
  if(!GLUE_PERIODIC){
    ktl=2*ktl+npl;
    ktr=2*ktr+npr;
  }
  if(ktl==ktr){
    interact_glue(&(x->re),left,npl,right,npr,couplings);

/*     In the one gluon case, the mass term is given among
       the heavy-light interactions, below.  */

    x->re /= longer.kstep; /* kt is twice total K */
  }
#else
   fprintf(stderr,__FILE__":  link-link interactions turned off.\n");
#endif  /* end of debugging flag to turn off interior interactions */

  /*  In the following, we calculate the contibutions from the ends of
      the strings.   */

#if 1   /* debugging flag to turn off ends of the strings */

	/* The only modification for different numbers of
           dimensions is the position of the tau interaction.
           Matching "spectators" is used to determine if the
           interaction is nonzero.  */

   if(npl==0 && npr==0){
     x->re += barepotential*couplings[GGN];
   } else {
     interact_heavy(x,left,npl,right,npr,couplings,0);
     orientation(left_flip,left,TYPE_SOURCE,npl);
     orientation(right_flip,right,TYPE_SOURCE,npr);
     /* left and right are switched to get the complex conjugate */ 
     interact_heavy(x,right_flip,npr,left_flip,npl,couplings,0);
     /* "pinching" terms tau_1 and tau_2, these are particular
        to the heavy source system. */

     if(npl==1 && npr==1)
       x->re+=couplings[TAU_2]*tau11[left[0]&KMASK][right[0]&KMASK];
     else if(npl==2 && npr==0)
       x->re+=couplings[TAU_1]*tau02[left[0]&KMASK][left[1]&KMASK];
     else if(npl==0 && npr==2)
       x->re+=couplings[TAU_1]*tau02[right[0]&KMASK][right[1]&KMASK];
   }
      
#endif  /* end of debugging flag to turn off end of string interactions */

  x->re /= sqrt(couplings[GGN]);
  x->im /= sqrt(couplings[GGN]);

  return;
}
/*
  This calculates matrix elements of the Hamiltonian v^+ P^- for 
  a heavy-light system, using variables defined above. 

  At the end, the answer is divided by sqrt(G^2 N).
*/

void interact_heavy_light(complex_element *x, partype *left, int npl, 
			  partype *right, int npr, element *couplings){
  int i,ktl=0,ktr=0;
  partype left_flip[NPMAX],right_flip[NPMAX];
  
  x->re=0.0; x->im=0.0;
  
  for(i=0; i<npl; i++)ktl += left[i]&KMASK;
  for(i=0; i<npr; i++)ktr += right[i]&KMASK;
  if(!GLUE_PERIODIC){
    ktl=2*ktl+npl;
    ktr=2*ktr+npr;
  }
  if(ktl==ktr){
    
#if 1      /* for debugging only -- turn off glue-glue interactions */
    if(left[npl-1]==right[npr-1]) /* check quarks */
      interact_glue(&(x->re),left,npl-1,right,npr-1,couplings);
#else
    fprintf(stderr,__FILE__":  link-link interactions turned off.\n");
#endif  /* end of debugging flag to turn off interior interactions */
    
#if 1   /* debugging flag to turn off light end of the string */
    /* Do the anti-quark interactions by transforming the state */
    
    charge_conjugation(left_flip,left,TYPE_MESON,npl);
    charge_conjugation(right_flip,right,TYPE_MESON,npr);
    interact_quark(x,left_flip,npl,right_flip,npr,couplings,npl-1);
    
#endif  /* end of debugging flag to turn off end of string interactions */
    
    x->re /= longer.p_max;
  }
  
#if 1   /* debugging flag to turn off heavy end of the string */
  
  /* The only modification for different numbers of
     dimensions is the position of the tau interaction.
     Matching "spectators" is used to determine if the
     interaction is nonzero.  */
  
  if(npl==1 && npr==1){
    printf(__FILE__":  have to figure this out\n");
    exit(-1);
  } else 
    if(left[npl-1]==right[npr-1])  /* make sure quarks match */
      interact_heavy(x,left,npl-1,right,npr-1,couplings,
		     2*(left[npl-1]&KMASK)+(QUARK_PERIODIC?0:1));
  
  /*  printf("x=%f%+f*i\n",x->re,x->im);*/
  
#endif  /* end of debugging flag to turn off end of string interactions */
  
  x->re /= sqrt(couplings[GGN]);
  x->im /= sqrt(couplings[GGN]);
  
  return;
}

/*
  Interactions with heavy source at one end of the string 

   x is not initialized

   On input, kt_light is twice the quark K;
   it becomes twice the total light field K.
*/

void interact_heavy(complex_element *x, partype *left, int npl, 
		     partype *right, int npr, element *couplings, 
		     partype kt_light){
  int i,j,l[NPMAX],r[NPMAX];
  complex_element z;

  for(i=0; i<npl; i++){
    l[i] = left[i]&KMASK; 
    kt_light += 2*l[i]+(GLUE_PERIODIC?0:1);
  }
  for(i=0; i<npr; i++)r[i] = right[i]&KMASK;

  if(npl==npr){

    /* diagonal interactions */
    for(j=0; j<npl; j++)if(left[j]!=right[j])break;
    if(j==npl){
      /* this is half the  "kinetic" term for this end of the string */
      x->re += 0.25*couplings[GGN]*longer.kstep
	*0.5*(element) kt_light; /* kt_light is twice K */
      /* This is half of the mass term on each end
	 since the light-light interactions include the other half. */
      if(longer.beta>ELEMENT_EPSILON)
#if IMPROVED<2 /* DLCQ matrix elements  */
      x->re += 0.25*couplings[MM]/(GE(left[0])*longer.kstep);
#else           /* improved matrix elements */
      x->re += 0.25*couplings[MM]*funpower(GE(left[0]),2*longer.beta)/
				  (funpower(GE(left[0]),2*longer.beta+1.0)*
				   longer.kstep);
#endif 
#if 0  /* debugging print statement */
      printf("npl=%i longer.kt=%i kt_light=%i ff[%i][%i]=%e\n",npl,
	     longer.kt,kt_light,
	     (longer.kt-kt_light)>>1,l[0],ff[(longer.kt-kt_light)/2][l[0]]);
#endif
      x->re += couplings[GGN]*ff[(longer.kt-kt_light)/2][l[0]];

    }

    /*   non-diagonal interactions */
    for(j=1; j<npl; j++)if(left[j]!=right[j])break;
    if(j==npl){
      source_phase(&z,left,1,right,1);
      ZTIMES(x,&z,gg[l[0]][r[0]]*couplings[GGN]);
    }

  } else if(npl==npr+2){

    for(j=0; j<npr; j++)if(left[j+2]!=right[j])break;   
    if(j==npr){
      source_phase(&z,left,2,right,0);
      ZTIMES(x,&z,gg13[l[1]][l[0]]*couplings[GGN]);
   }
    
  } else if(npr==npl+2){

    for(j=0; j<npl; j++)if(left[j]!=right[j+2])break;   
    if(j==npl){
      source_phase(&z,left,0,right,2);
      ZTIMES(x,&z,gg13[r[1]][r[0]]*couplings[GGN]);
     }

  }

#if 0 /* debug print */
  printf("npl=%i npr=%i x=%f%+f*i\n",npl,npr,x->re,x->im);
#endif

  return;
}

/* This prints out the memory used by the heavy source interactions */

void print_source_interaction_memory(FILE *fp){
  fprintf(fp,"MM heavy source interactions use " SIZE_T_PRINT " bytes\n",
	      source_memory);
}

/* Free memory used by the heavy source interactions */

void free_source(void){
  int j;
  for(j=0; j<KKMAX; j++){
    free(ff[j]); ff[j]=NULL;
    free(gg[j]); gg[j]=NULL;
    free(gg13[j]); gg13[j]=NULL;
    free(tau11[j]); tau11[j]=NULL;
    free(tau02[j]); tau02[j]=NULL;
  }
  source_memory=0;
}