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

/*    Matrix elements for the longitudinal string tension using an open string

    This is the interaction for the open string, it returns a 
    complex number. x[] is of length NPARAMS.

    In addition, it is a function of the parameters:

      kmax = K_{max}/(v^+ sqrt{G^2 N})

      kstep = 2*kmax/abs(kt)

      ll = L sqrt{G^2 N}

    where abs(kt), kt<0, is the upper limit of the discretized momenta
    and K = kt/2.  The argument of long[] is the change in kt/2 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<=abskt/2. 
   
*/

static element ff[KKMAX][KKMAX];
static complex_e gg[KKMAX][KKMAX],gg13[KKMAX][KKMAX],
		      tt[KKMAX][KKMAX],tt13[KKMAX][KKMAX];
static element barepotential,kstep,kmax,beta,kiki,b;
static int abskt;

#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 llin, int ktin, element kmaxin, element betain){
  int i,j,l,ipmj;
  complex_e phase;
  element ll;
#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;
  const element pi=3.1415926535897932384626433832795;
  ll=llin; kmax=kmaxin; abskt=abs(ktin);
  if(abskt>2*KKMAX){
    fprintf(stderr,"ktin in maketenslist too large, exiting\n");
    abskt=0;
    return;
  }
  kstep=2.0*kmax/((element) abskt);
  beta=betain;

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

  barepotential=0.25*fabs(ll)*kstep/kmax;

  /* All interactions (except for self-inertias)
     The sqrt(funpower(ki... in the demoninator
     is the normalization of the x^beta wavefunction */


  for(i=0; 2*i+1<=abskt; i++){
    ki=(element) i+0.5;

  /* first set imaginary part equal to zero */

    for(j=0; 2*j+1<=abskt; j++){
      gg[i][j].im=0.0; gg13[i][j].im=0.0;
      tt[i][j].im=0.0; tt13[i][j].im=0.0;
    }

    /* Do one triangle of matrix */

#if IMPROVED<2 /* DLCQ matrix elements  */

    for(j=0; j<=i; j++){
      kj=(element) j+0.5;
#if NEWTAU
      tt[i][j].re = 1.0/(8.0*kmax*pi*sqrt(ki*kj));
      if(i!=j) tt[i][j].re /= ki-kj;
      tt13[i][j].re= 1.0/(8.0*kmax*pi*sqrt(ki*kj)*(ki+kj));
#else /* not NEWTAU */
      tt[i][j].re = kstep/(8.0*pi*kmax*sqrt(ki*kj));
      tt13[i][j].re=tt[i][j].re;
#endif /* NEWTAU */
      if(i!=j){
	gg[i][j].re = -(ki+kj)/(8.0*kmax*pi*sqrt(ki*kj)*(ki-kj)*(ki-kj));
	gg13[i][j].re = (ki-kj)/(8.0*kmax*pi*sqrt(ki*kj)*(ki+kj)*(ki+kj));
      }
    }

#else /* improved matrix elements */

    for(j=0; j<=i; j++){
      kj=(element) j+0.5;
      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
#if NEWTAU
      tt[i][j].re = funpower(ki,beta+0.5)*funpower(kj,beta+0.5)/
	(8.0*kmax*pi*sqrt(funpower(ki,2.0*beta+1.0)*
			  funpower(kj,2.0*beta+1.0)));
      if(i!=j)tt[i][j].re /= ki-kj;
      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));
      /* The i is the second monentum and j is the first, thus a sign: */
      tt13[i][j].re = ansvtau/(8.0*kmax*pi*sqrt(ans));
#else /* not NEWTAU */
      tt[i][j].re = kstep*funpower(ki,beta+0.5)*funpower(kj,beta+0.5)/
	(8.0*pi*kmax*sqrt(funpower(ki,2.0*beta+1.0)*
			  funpower(kj,2.0*beta+1.0)));
      DGAUS8(VTAU_TO_FORTRAN,&x,&m,&err,&ansvtau,&ierr);
      if(i==0)ansvtau += 1.0/((beta+0.5)*pow(kj+0.5,beta+0.5));
      if(j==0)ansvtau += 1.0/((beta+0.5)*pow(ki+0.5,beta+0.5));
      tt13[i][j].re = kstep*ansvtau/(8.0*pi*kmax*sqrt(ans));
#endif /* NEWTAU */
      if(i!=j){
	gg[i][j].re = -(funpower(ki,beta+1.5)*funpower(kj,beta+0.5)+
			funpower(ki,beta+0.5)*funpower(kj,beta+1.5))/
	  (8.0*kmax*pi*sqrt(funpower(ki,2.0*beta+1.0)*
			    funpower(kj,2.0*beta+1.0))*(ki-kj)*(ki-kj));
	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].re = -ansinst/(8.0*kmax*pi*sqrt(ans));
      }
    }

#endif  /* improved matrix elements */

    /* Do Diagonal parts */

    gg[i][i].re=0.0;
    gg13[i][i].re=0.0;

#if 0  /* debugging print */
    for(j=0; j<=i; j++){
      printf("new improved ki=%g kj=%g",ki,kj);
      printf(" norm=%g vinst=%g vtau=%g\n",ans,ansinst,ansvtau);
      printf("                 gg13=%g tt13=%g\n",
	     gg13[i][j].re,tt13[i][j].re);
    }
#endif
  }

  /* Calculate the self-intertias  */

  for(l=0; l<=abskt>>1; l++){
    for(j=0; 2*j+1<=abskt-2*l; j++){
      kj=(element)j+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.  Only one triangle of the
       gg matrix is now filled. */
      for(i=0, x=0.0; i<=l+j; i++)
	x -= (i>j?gg[i][j].re:gg[j][i].re);
      if(l>0)
	ff[l][j] = x+sqrt(m)/(4.0*pi*kmax*(element) l);
      else
	ff[l][j] = x+1.0/(8.0*pi*kmax*kj);
#if 0 /*debugging print */
      printf("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. Also, only one
	 triangle of the gg matrix is now filled.     */
	x -= pow(ki/kj,beta)*(i>j?gg[i][j].re:gg[j][i].re);
      }
      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 kmax = %e\n",8.0*x*pi*kmax);
#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*kmax*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:  kmax*ff[%i][%i]=%21.14e\n",l,j,kmax*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
  }

  /*  Add phases for nonzero ll and 
      Fill up other triangle of array 

      In the following "ipmj" means i-j or i+j. 

      We have to handle the endpoints (when there
      is a possible divergence) of the new 
      tau interactions in a special way */

#if NEWTAU
  for(j=0; 2*j+1<=abskt; j++)
    tt[j][j].im=0.5*ll*kstep*tt[j][j].re;
#endif
  for(ipmj=1; 2*ipmj+1<=abskt; ipmj++){
    phasetocomplex(&phase,0.5*ll*kstep*(element) ipmj);
    for(j=0; 2*(i=ipmj+j)+1<=abskt; j++){
      multiply(gg[i]+j,&phase);
      gg[j][i].re=gg[i][j].re;
      gg[j][i].im=-gg[i][j].im;
      multiply(tt[i]+j,&phase);
#if NEWTAU
      tt[j][i].re=tt[i][j].re;
      tt[j][i].im=tt[i][j].im;
#else /* not NEWTAU */
      tt[j][i].re=tt[i][j].re;
      tt[j][i].im=-tt[i][j].im;
#endif /* NEWTAU */
    }
  }

  /*  for the 0->2 link interactions, we only calculate 
      phases when the total link
      momentum is less than the cutoff.   */

  for(ipmj=0; 2*ipmj+2<=abskt; ipmj++){
    phasetocomplex(&phase,0.5*ll*kstep*(element) (ipmj+1));
    for(j=0, i=ipmj; j<=i; j++,i--){
      multiply(gg13[i]+j,&phase);
      gg13[j][i].re = -gg13[i][j].re;
      gg13[j][i].im = -gg13[i][j].im;
      /* We have to handle the endpoint of the new 
         tau interactions in a special way for 
         improved matrix elements */
      if(ipmj==0 && NEWTAU && IMPROVED>=2)
	tt13[i][i].im=0.5*ll*kstep*tt13[i][i].re;
      else
	multiply(tt13[i]+j,&phase);
      tt13[j][i].re = tt13[i][j].re;
      tt13[j][i].im = tt13[i][j].im;
    }
  }

#if 0  /* debugging print */
  for(i=0; 2*i+1<=abskt && i<5; i++)
    for(j=0; 2*j+1<=abskt && j<5; j++)
      printf("   tt[%i][%i] = %g %g\n",i,j,tt[i][j].re,tt[i][j].im);
#endif
}

/*  function to mulitply two complex numbers 
    the result is returned in x. */

void multiply(complex_e *x,complex_e *y){
  element temp;
  temp = x->re * y->im + x->im * y->re;
  x->re = x->re * y->re - x->im * y->im;
  x->im=temp;
}

void phasetocomplex(complex_e *p, element x){
  p->re = cos(x); p->im = sin(x);
}

/*  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.*beta)*
      (-5. + 2.*beta + b*(-1. + 2.*beta)) + 
     (-1. + 4.*pow(beta,2))*(-9. + 2.*beta + b*(-3. + 2.*beta))*
      (-1. + *z));
  else
    return (0.5*(-3. - 2.*beta*(-1. + *z) + *z + 4.*pow(*z,0.5 + beta) - 
       2.*pow(*z,1.5 + beta) - 
       1.*b*(1. + 2.*beta*(-1. + *z) + *z - 2.*pow(*z,0.5 + 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*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),beta-0.5)*(1-2.0* *x)*log(MINK/MAXK);
  if(ki<0.75)z -= pow(*x,beta-0.5)*log((kj+0.5)/(kj-0.5));
  if(kj<0.75)z += pow(1.0- *x,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),beta-0.5)*(MINK-MAXK)
      -pow(*x,beta-0.5) - pow(1.0- *x,beta-0.5);
  } else {
    z = pow(*x *(1.0- *x),beta-0.5)*log(MINK/MAXK);
    if(ki<0.75)z -= pow(*x,beta-0.5)*log((kj+0.5)/(kj-0.5));
    if(kj<0.75)z -= pow(1.0- *x,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),beta-0.5)*(MINK-MAXK);
  if(ki<0.75)z -= pow(*x,beta-0.5);
  if(kj<0.75)z -= pow(1.0- *x,beta-0.5);
  return z;
}

/*  This calculates matrix elements of the Hamiltonian for 
    an open string, using variables defined above. 
    Both HINDEX==1 and HINDEX==2 are done together.
*/

void interactiontens(complex_e *x, partype *left, partype *right, 
		 int npl, int npr){
  int i,j,ktl=0,ktr=0,l[NPMAX],r[NPMAX],jj;
  element *row,a;
#if HINDEX==1
  int hl[NPMAX],hr[NPMAX];
  for(i=0; i<npl; i++){l[i] = hitol(left[i],hl+i); ktl += 2*l[i]+1;}
  for(i=0; i<npr; i++){r[i] = hitol(right[i],hr+i); ktr += 2*r[i]+1;}
#elif HINDEX==2
  int hl[NPMAX][HINDEX],hr[NPMAX][HINDEX];
  for(i=0; i<npl; i++){l[i] = hitol(left[i],hl[i]); ktl += 2*l[i]+1;}
  for(i=0; i<npr; i++){r[i] = hitol(right[i],hr[i]); ktr += 2*r[i]+1;}
#endif
  /*printf("ktl=%i ktr=%i\n",ktl,ktr);*/
  for(i=0; i<NPARAMS; i++){x[i].re=0.0; x[i].im=0.0;}

  /*  The first part calculates the interior of the string and
      follows the closed loop case */
#if 1      /* for debugging only -- turn off light-light interactions */
  if(ktl==ktr){

#if HINDEX==1  /* 2+1 dimensional case */

    if(npl==npr && npl>1){
      for(j=0; j<npl-1; j++){
	/* If there are any mismatches on the left 
           hand side, then there is no hope for further matches */
	if(j>0 && left[j-1]!=right[j-1])break;
	for(jj=j+2; jj<npl; jj++)
	  if(left[jj]!=right[jj])break;
	if(jj==npl){
	  row=twoelements[l[j]+l[j+1]]+four*((l[j]+l[j+1]+1)*l[j]+r[j]);
	  if(hl[j] == hr[j]){
            x[0].re+=row[0]; 
	    x[1].re+=row[1];
	    if(hl[j]==hr[j+1])
	      x[3].re+=row[3];
	    else
	      x[2].re+=row[3];
	  } 
	  if(hl[j] != hl[j+1]){
	    /* Note that the instantaneous-annihilation graph
	       can also act in the two particle sector, unlike
               in the closed loop case. */
	      x[1].re-=row[2]; 
	    if(hl[j]==hr[j])
	      x[2].re+=row[3]; 
	    else
	      x[3].re+=row[3]; 
	  } 
	}
      }
      /* The link field "pinching" terms are suppressed in the open 
	 string.                                          */
    } else if(npl+2 == npr){
      for(j=0; j<npl; j++){
	if(j>0 && left[j-1]!=right[j-1])break;
	for(jj=j+1; jj<npl; jj++)
	  if(left[jj]!=right[jj+2])break;
	if(jj==npl){
	  if(hl[j] == hr[j+2]){
	    row=threeelements[l[j]-1]+
	      ((((2*l[j]-r[j]+1)*r[j])>>1)+r[j+1])*two;
            x[1].re+=row[0];
	    if(hl[j]==hr[j+1])
	      x[3].re+=row[1];
	    else
	      x[2].re+=row[1];
	  }
	  if(hl[j] == hr[j]){
	    row=threeelements[l[j]-1]+
	      ((((2*l[j]-r[j+2]+1)*r[j+2])>>1)+r[j+1])*two;
            x[1].re+=row[0];
	    if(hl[j]==hr[j+1])
	      x[3].re+=row[1];
	    else
	      x[2].re+=row[1];
	  }
	} 
      } 
    }
    else if(npr+2 == npl){
      for(j=0; j<npr; j++){
	if(j>0 && left[j-1]!=right[j-1])break;
	for(jj=j+1; jj<npr; jj++)
	  if(left[jj+2]!=right[jj])break;
	if(jj==npr){
	  if(hr[j] == hl[j+2]){
	    row=threeelements[r[j]-1]+
	      ((((2*r[j]-l[j]+1)*l[j])>>1)+l[j+1])*two;
            x[1].re+=row[0];
	    if(hr[j]==hl[j+1])
	      x[3].re+=row[1];
	    else
	      x[2].re+=row[1];
	  }
	  if(hr[j] == hl[j]){
	    row=threeelements[r[j]-1]+
	      ((((2*r[j]-l[j+2]+1)*l[j+2])>>1)+l[j+1])*two;
            x[1].re+=row[0];
	    if(hr[j]==hl[j+1])
	      x[3].re+=row[1];
	    else
	      x[2].re+=row[1];
	  }
	} 
    }
    }

#elif HINDEX==2 /* 3+1 dimensional case */

    if(npl==npr && npl>1){
      for(j=0; j<npl-1; j++){
	/* If there are any mismatches on the left 
           hand side, then there is no hope for further matches */
	if(j>0 && left[j-1]!=right[j-1])break;
	for(jj=j+2; jj<npl; jj++)
	  if(left[jj]!=right[jj])break;
	if(jj==npl){
	  row=twoelements[l[j]+l[j+1]]+four*((l[j]+l[j+1]+1)*l[j]+r[j]);
          if(EQUAL(hl[j],hr[j])){
	    x[0].re+=row[0]; x[1].re+=row[1];
            if(EQUAL(hl[j],hl[j+1]))
	      x[4].re+=row[3];
	    else if(MINUS(hl[j],hl[j+1]))
	      x[3].re+=row[3];
	    else
	      x[6].re+=row[3];
	  } 
          if(MINUS(hl[j],hl[j+1])){
	    /* Note that the instantaneous-annihilation graph
	       can also act in the two particle sector, unlike
               in the closed loop case. */
	    x[1].re-=row[2];
            if(EQUAL(hl[j],hr[j]))
	      x[3].re+=row[3];
	    else if(MINUS(hl[j],hr[j]))
	      x[4].re+=row[3];
	    else
	      x[6].re+=row[3];
          } 
          if(EQUAL(hl[j],hr[j+1])&&ORTHO(hl[j],hl[j+1]))
	    x[2].re-=row[3];
	}
      /* The "pinching" terms are suppressed in the open string although
         one could introduce terms like (\phi^\d M \phi)^2  */
      } 
    }
    else if(npl+2 == npr){
      for(j=0; j<npl; j++){
	if(j>0 && left[j-1]!=right[j-1])break;
	for(jj=j+1; jj<npl; jj++)
	  if(left[jj]!=right[jj+2])break;
	if(jj==npl){
	  if(EQUAL(hl[j],hr[j+2])){
	    row=threeelements[l[j]-1]+
	      ((((2*l[j]-r[j]+1)*r[j])>>1)+r[j+1])*two;
	    x[1].re+=row[0];
	    if(EQUAL(hl[j],hr[j+1]))
	      x[4].re+=row[1];
	    else if(MINUS(hl[j],hr[j+1]))
	      x[3].re+=row[1];
	    else
	      x[6].re+=row[1];
	  }
	  if(EQUAL(hl[j],hr[j])){
	    row=threeelements[l[j]-1]+
	      ((((2*l[j]-r[j+2]+1)*r[j+2])>>1)+r[j+1])*two;
	    x[1].re+=row[0];
	    if(EQUAL(hl[j],hr[j+1]))
	      x[4].re+=row[1];
	  else if(MINUS(hl[j],hr[j+1]))
	    x[3].re+=row[1];
	    else
	    x[6].re+=row[1];
	  }
	  if(EQUAL(hl[j],hr[j+1])&&ORTHO(hl[j],hr[j])){
	    row=threeelements[l[j]-1]+((((2*l[j]-r[j]+1)*r[j])>>1)+r[j+1])*two;
	    x[2].re-=row[1];
	  }
	}
      }  
    }
    else if(npr+2 == npl){
      for(j=0; j<npr; j++){
	if(j>0 && left[j-1]!=right[j-1])break;
	for(jj=j+1; jj<npr; jj++)
	  if(left[jj+2]!=right[jj])break;
	if(jj==npr){
	  if(EQUAL(hr[j],hl[j+2])){
	    row=threeelements[r[j]-1]+
	      ((((2*r[j]-l[j]+1)*l[j])>>1)+l[j+1])*two;
            x[1].re+=row[0];
            if(EQUAL(hr[j],hl[j+1]))
	      x[4].re+=row[1];
	    else if(MINUS(hr[j],hl[j+1]))
	      x[3].re+=row[1];
	    else
	      x[6].re+=row[1];
	  }
	  if(EQUAL(hr[j],hl[j])){
	    row=threeelements[r[j]-1]+
	      ((((2*r[j]-l[j+2]+1)*l[j+2])>>1)+l[j+1])*two;
            x[1].re+=row[0];
            if(EQUAL(hr[j],hl[j+1]))
	      x[4].re+=row[1];
	    else if(MINUS(hr[j],hl[j+1]))
	      x[3].re+=row[1];
	    else
	      x[6].re+=row[1];
	  }
	  if(EQUAL(hr[j],hl[j+1])&&ORTHO(hl[j],hr[j])){
	    row=threeelements[r[j]-1]+((((2*r[j]-l[j]+1)*l[j])>>1)+l[j+1])*two;
	    x[2].re-=row[1];
	  } 
	}
      } 
    }

#endif  /* end 3+1 dimensional case */

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

        For the interior of the string, we divide by an extra factor 
	of 2 since we want P^- and not 2 P^+ P^- */
    a=0.5/kmax;
    for(i=0; i<NPARAMS; i++)x[i].re *= a;
  }
#else
  fprintf(stderr,"link-link interactions turned off.\n");
#endif  /* end of debugging flag to turn off interior interactions */

  /*  For the entire string, we divide by kstep; this
      should be done by the calling Hamiltonian routine.  */

  /*  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 HINDEX==1
#define TAU 5
#elif HINDEX==2
#define TAU 8
#endif

  if(npl==0 && npr==0)
    x[1].re += barepotential;

  else if(npl==npr){
    if(!memcmp(left,right,npl*sizeof(partype))){
      x[1].re += 0.5/kmax*kstep*kstep*(element) ktl/2.0;
      /* This is half of the mass term on each end
	 since the light-light interactions include the other
	 half.                         */
      if(beta>1.0e-12)
#if IMPROVED<2 /* DLCQ matrix elements  */
      x[0].re += 0.25*(1.0/(l[0]+0.5)+1.0/(l[npl-1]+0.5))/kmax;
#else           /* improved matrix elements */
      x[0].re += 0.25*(funpower(l[0]+0.5,2*beta)/funpower(l[0]+0.5,2*beta+1.0)+
	      funpower(l[npl-1]+0.5,2*beta)/funpower(l[npl-1]+0.5,2*beta+1.0))
	       /kmax;
#endif 
#if 0  /* debugging print statement */
      printf("npl=%i ff[%i][%i]=%e ff[%i][%i]=%e\n",npl,
	     ktl>>1,l[0],ff[ktl>>1][l[0]],
	     ktl>>1,l[npl-1],ff[ktl>>1][l[npl-1]]);
#endif
      x[1].re += ff[(abskt-ktl)>>1][l[0]]+ff[(abskt-ktl)>>1][l[npl-1]];
    }
    if(!memcmp(left+1,right+1,(npl-1)*sizeof(partype))){
      x[1].re += gg[l[0]][r[0]].re;
      x[1].im += gg[l[0]][r[0]].im;
#if !NEWTAU
      x[TAU].re += tt[l[0]][r[0]].re;
      x[TAU].im += tt[l[0]][r[0]].im;
#endif  /* !NEWTAU */
    }
    if(!memcmp(left,right,(npl-1)*sizeof(partype))){
      x[1].re += gg[l[npl-1]][r[npr-1]].re;
      x[1].im -= gg[l[npl-1]][r[npr-1]].im;
#if !NEWTAU
      x[TAU].re += tt[l[npl-1]][r[npr-1]].re;
      x[TAU].im -= tt[l[npl-1]][r[npr-1]].im;      
#endif  /* !NEWTAU */
    }
    if(NEWTAU&&npl==1){
      x[TAU+1].re += tt[l[0]][r[0]].im;
    }

  } else if(npl==npr+2){
    if(!memcmp(left+2,right,npr*sizeof(partype))){
      x[1].re += gg13[l[1]][l[0]].re;
      x[1].im += gg13[l[1]][l[0]].im;
#if !NEWTAU
      x[TAU].re += tt13[l[0]][l[1]].re;
      x[TAU].im += tt13[l[0]][l[1]].im;      
#endif  /* !NEWTAU */
    }
    if(!memcmp(left,right,npr*sizeof(partype))){
      x[1].re += gg13[l[npl-2]][l[npl-1]].re;
      x[1].im -= gg13[l[npl-2]][l[npl-1]].im;
#if !NEWTAU
      x[TAU].re += tt13[l[npl-1]][l[npl-2]].re;
      x[TAU].im -= tt13[l[npl-1]][l[npl-2]].im;      
#endif  /* !NEWTAU */
    }
    if(NEWTAU&&npl==2)
      x[TAU].re += tt13[l[0]][l[1]].im;

  } else if(npr==npl+2){
    if(!memcmp(left,right+2,npl*sizeof(partype))){
      x[1].re += gg13[r[1]][r[0]].re;
      x[1].im -= gg13[r[1]][r[0]].im;
#if !NEWTAU
      x[TAU].re += tt13[r[0]][r[1]].re;
      x[TAU].im -= tt13[r[0]][r[1]].im;            
#endif  /* !NEWTAU */
    }
    if(!memcmp(left,right,npl*sizeof(partype))){
      x[1].re += gg13[r[npr-2]][r[npr-1]].re;
      x[1].im += gg13[r[npr-2]][r[npr-1]].im;
#if !NEWTAU
      x[TAU].re += tt13[r[npr-1]][r[npr-2]].re;
      x[TAU].im += tt13[r[npr-1]][r[npr-2]].im;            
#endif  /* !NEWTAU */
    } 
    if(NEWTAU&&npr==2)
      x[TAU].re += tt13[r[0]][r[1]].im;
  }
#endif  /* end of debugging flag to turn off end of string interactions */
  return;
}



#if 1==0  /* older versions of maketenslist and selfinertia */

void maketenslist(element llin, int ktin, element kmaxin, element betain){
  int i,j,l;
  integer ierr;
  element c=1.0;
  /* err is the error of the integrator */
  element x,ki,kj,m,zero=0.0,err=1.0e-10,ans;
  const element pi=3.1415926535897932384626433832795;
  ll=llin; kmax=kmaxin; abskt=abs(ktin);
  kstep=2.0*kmax/((element) abskt);
  tailpow=2.5;
  beta=betain;
  for(i=0; 2*i<=abskt; i++){
    ki=0.5*ll*kstep*(element) i;
    longlist[i].re=cos(ki);
    longlist[i].im=sin(ki);
  }
  /* self-intertias; This is currently done with Burkardt's technique
     (That is, the ll exponential and the x^beta are not included.)  */
  for(l=0; l<=abskt>>1; l++)
    for(j=0; 2*j+1<=abskt-2*l; j++){
      x=0.0; kj=(element)j+0.5;
#if 1==1 /* switch between new and old self inertias */
      m=((element) l+kj)/kj;
      w=c*pow(kj*kstep,tailpow);
      for(i=0; i<=l+j; i++)
	if(i!=j){
	  ki = (element)i+0.5;
	  x += pow(ki/kj,beta-0.5)*(ki+kj)/
	    ((ki-kj)*(ki-kj)*(1.0+w*pow(ki/kj,tailpow)));
	}
      DGAUS8(SELFINERTIA_TO_FORTRAN,&zero,&m,&err,&ans,&ierr);
#if 1==0 /*debugging print */
      printf("  integration answer, m=%e w=%e integral= %e\n",m,w,ans);
      printf("  sum answer = %e\n",x);
#endif 
      if(m>1.0)
	ans += (pow(m,0.5 + beta))/(0.5 + beta) - 
	  (2.*m)/((-1. + m)*(1. + w)) + 
	  ((beta + beta*w - tailpow*w)*2.0*log(fabs(-1. + m)))/
	  ((1. + w)*(1. + w));
      else
	ans += (2.*(-2.*beta + w))/((1. + 2.*beta)*(1. + w));
      ff[l][j] = (x-ans)*(1.0+w)/(8.0*pi*kmax*kj);
#if 1==0 /*debugging print */
      printf("w=%e m=%e beta=%e tailpow=%21.14e\n",w,m,beta,tailpow);
#endif
#else
      for(i=0; i<=l+j; i++)
	if(i!=j){
	  ki = (element)i+0.5;
	  x += (ki+kj)/(sqrt(ki*kj)*(ki-kj)*(ki-kj));
	}
      ff[l][j] = x/(8.0*pi*kmax);
#endif
#if 1==0 /* DEBUGGING print, print out self-inertias */
      printf(" self inertia:  kmax*ff[%i][%i]=%21.14e\n",l,j,kmax*ff[l][j]);
#endif
    }
  /* Other interactions */
  for(i=0; 2*i+1<=abskt; i++){
    ki=(element) i+0.5;  
    for(j=0; 2*j+1<=abskt; j++){
      kj=(element) j+0.5;    
      if(i!=j) 
	gg[i][j].re = -(ki+kj)/(8.0*kmax*pi*sqrt(ki*kj)*(ki-kj)*(ki-kj));
      else
	gg[i][j].re = 0.0;
      gg13[i][j].re = (ki-kj)/(8.0*kmax*pi*sqrt(ki*kj)*(ki+kj)*(ki+kj));
      tt[i][j].re = kstep/(8.0*pi*kmax*sqrt(ki*kj));
    }
  }
#if 1==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

}

/*  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 tailpow, w, and beta. */

element selfinertia(element *z){
  if(fabs((*z)-1.0)<1.0e-4)
    return (0.0833333333333333*(3.*(1. + w)*
        (-3. - 2.*(5. - 2.*tailpow + 2.*pow(tailpow,2))*w + 
          (-11. + 4.*tailpow + 4.*pow(tailpow,2))*pow(w,2) - 
          4.*pow(w,3) + 4.*pow(beta,2)*pow(1. + w,2) - 
          4.*beta*(1. + w)*(1. + w + 2.*tailpow*w)) + 
       (3. + (15. - 11.*tailpow + 12.0*pow(tailpow,2) - 
             4.*pow(tailpow,3))*w + 
          (27. - 22.*tailpow + 16.*pow(tailpow,3))*pow(w,2) - 
          1.*(-21. + 11.*tailpow + 12.0*pow(tailpow,2) + 
             4.*pow(tailpow,3))*pow(w,3) + 6.*pow(w,4) + 
          4.*pow(beta,3)*pow(1. + w,3) - 
          12.0*pow(beta,2)*pow(1. + w,2)*
           (1. + w + tailpow*w) - 
          1.*beta*(1. + w)*(1. + 
             2.*(7. - 12.0*tailpow + 6.*pow(tailpow,2))*
              w + (25. - 24.0*tailpow - 
                12.0*pow(tailpow,2))*pow(w,2) + 
             12.0*pow(w,3)))*(-1. + (*z))))/pow(1. + w,4);
  else
    return -2./((1. + w)*(-1. + (*z))*(-1. + (*z))) - 
      (2.*(beta + beta*w - tailpow*w))/((1. + w)*(1. + w)*(-1. + (*z))) - 
      1.*pow((*z),-0.5 + beta) + (pow((*z),-0.5 + beta)*(1. + (*z)))/
      ((-1. + (*z))*(-1. + (*z))*(1. + w*pow((*z),tailpow)));
}


#endif