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

/*  This file calculates the overlap of the wavefunction
    for the longitudinal potential and glueball eigenstates.   */


int MAIN(){
  int i,vv,n;
  static struct sectors xy[2],lx,gx;
  clock_t t1;
  element beta;
  integer nval=5;
  doublereal *val,*vec=NULL,*lval=NULL,*lvec=NULL;
  partype part[NPMAX];
#define NOUT 100 /* Output file name length */
  char out[NOUT]="decay.out",*pout;
  FILE *fp;
#if HINDEX==1 /* have not implemented */

#elif HINDEX==2

  int     ht[HINDEX]={0,0},np=6,kt=26,cyclic=1,o=1,multi=8,multi2=9,nfig=8;
  element kmax=4.0,ll=0.0;
  int     npg=6,gkt=36,cyclicg=1,go=1,gmulti=3;

#if 0  /* test couplings */
  element  couplings[NPARAMS]={0.03, 1.0, 0.7,0.03,-0.005, 3.0,0.13,1.0,-2.0};
#elif 1 /* couplings used in first 3+1 paper hep-th/9810236 */
  element rescale=8.099820052, couplings[NPARAMS]={0.0517766953,1,0.7584261311,
	    -0.07623059914,-0.1458991891,8.889796764,0.1608104465,
            1.657995607,-3.268403913};
#endif
#endif

#ifdef BABBAGE  /* Routine to initialize Fortran on Babbage */
  if(hf_fint((char*)NULL)) exit(1);
#endif

#if 1  /* Ask for output file name*/
  printf("Output file name:  ");
  fgets(out,NOUT,stdin);
  if((pout=strchr(out,'\n'))!=NULL)*pout='\0';
  printf("%s\n",out);

#endif
  if((fp = fopen(out,"w")) == NULL){
     fprintf(stderr,__FILE__":  Can't open file %s in structure.c, exiting\n",out);
    goto cleanup;
  }

loop:

  printf("Input number of eigenvalues (0=exit):  ");
  scanf("%li",&nval);
  if(nval==0)goto cleanup;


  /******  Input parameters for the basis  *************/
#if 1   /* input parameters for basis (instead of defaults.) */
  printf("Input ht[%i], p truncation, kt=2*K, "
	       "cyclic, o, multiplet:\n",HINDEX);
  for(i=0;i<HINDEX;i++)scanf("%i",ht+i);
  scanf("%i %i %i %i %i",&np,&kt,&cyclic,&o,&multi);
  if(kt<0&&!cyclic){
    printf("Input P^+ cutoff kmax and long separation ll:\n");
    scanf("%lg %lg",&kmax,&ll);
  } else {
    printf("Input second multiplet:\n");
    scanf("%i",&multi2);
  }
#endif

#if 0 /******  Input coupling constants  *************/
  printf("Input %i couplings:  ",NPARAMS);
  for(i=0;i<NPARAMS;i++) scanf("%le",couplings+i);
  for(i=0;i<NPARAMS;i++) printf(" %.14g",couplings[i]);
  printf("\n");
#endif

/*  Print out parameters for basis and couplings */

#if HINDEX==1
  fprintf(fp,"****Wavefunction overlap for the 2+1 transverse lattice****\n");
  fprintf(fp,"Lowest %i Eigenvalue(s) for various momenta and n\n",nval);
  fprintf(fp,"Couplings:  m^2, G^2 N, lambda_1/a, lambda_2/a,"
	  " lambda_3/a tau\n");
#elif HINDEX==2
  fprintf(fp,"****Wavefunction overlap for the 3+1 transverse lattice****\n");
  fprintf(fp,"Couplings:  m^2, G^2 N, t/a^2, la_1/a^2, la_2/a^2\n");
  fprintf(fp,"           la_3/a^2, la_4/a^2, la_5/a^2, tau     \n");
#endif
  fprintf(fp,"For heavy sources (notation K<0), divide by sqrt(G^2 N).\n");
  for(i=0; i<NPARAMS; i++)fprintf(fp," " PDIGITS,couplings[i]);
  fprintf(fp,"\nOverall scale G^2 N/sigma = " PDIGITS "\n",rescale);
  fprintf(fp,"Winding (");
  for(i=0; i<HINDEX; i++)fprintf(fp," %i",ht[i]);
  fprintf(fp,") spectrum using (K,p) = (%i/2,%i)\n",kt,np);
#if MOM
  fprintf(fp," with multiplet=%i & %i ",multi,multi2);
  fprintf(fp,"and orientation symmetry indefinite.\n");
  if(!cyclic && type==TYPE_SOURCE)
    fprintf(fp,"\nHeavy sources, K_max = %f\n",kmax);
#else
  fprintf(fp," with multiplet=%i and orientation symmetry %i.\n",multi,o);
#endif

  /*  construct basis, any previously constructed basis in x is
      deallocated  */
  printf("Constructing basis and tables");
  t1=clock();
  makebasis(&gx,ht,gnp,gkt,cyclicg,go,gmulti);
  makebasis(xy,ht,np,kt,cyclic,o,multi);
  for(i=0, j=0; i<HINDEX; i++)j+=abs(ht[i]);
  if(j==0)j=2;
  for(i=0, k=j; k<=kt; i++, k+=2)
    makebasis(lx+i,ht,gnp,k,gcyclic,go,gmulti);
  nlx=i;
  if(o==0){
     fprintf(stderr,__FILE__":  Need nonzero o, exititing\n");
    exit(10);
  }
  if(type==TYPE_SOURCE&&!cyclic)
    makebasis(&y,ht,np,kt,cyclic,-o,multi);
  else{
     fprintf(stderr,__FILE__":  Need heavy source\n");
    goto cleanup;
    makebasis(xy+1,ht,np,kt,cyclic,-o,multi2);
  }
  printf(" in %f CPU seconds\n",
       (float)(clock()-t1)/(float) CLOCKS_PER_SEC);

  beta=mapcouplings(kt>gkt?kt:gkt,couplings);

  if(x.type==TYPE_SOURCE){
    maketenslist(ll,x.kt,kmax,beta); 
  } else {    
     fprintf(stderr,__FILE__":  This is supposed to be the heavy potential\n");
    goto cleanup;
  }
  printf(" in %f CPU seconds\n",
       (float)(clock()-t1)/(float) CLOCKS_PER_SEC);


  printf("Starting calculation of potential.\n");
  t1=clock();
  spectrum(xy,xy+1,&val,&vec,nval,couplings,nfig);
  printf("%f CPU seconds to calculate Hamiltonian and run Lanczos\n\n",
       (float)(clock()-t1)/(float) CLOCKS_PER_SEC);

  /*  Copy the resuling eigenvalues and eigenvectors into
      lval and lvec */

  for(vv=0; vv<nval; vv++)
    printf("Eigenvalue %21.14e\n",sqrt(rescale)*val[vv]);
  n=x.length+y.length;
  lvec=realloc(lvec,nval*n*sizeof(element));
  memcpy(lvec,vec,nval*n*sizeof(element));
  lval=realloc(lval,nval*sizeof(element));
  memcpy(lval,val,nval);

  printf("Starting calculation of spectrum.\n");
  t1=clock();
  spectrum(&gx,&gx,&val,&vec,nval,couplings,nfig);
  printf("%f CPU seconds to calculate Hamiltonian and run Lanczos\n\n",
       (float)(clock()-t1)/(float) CLOCKS_PER_SEC);

  for(vv=0; vv<nval; vv++)
    printf("Eigenvalue %21.14e\n",rescale*val[vv]);
  
  /* ground state overlap with various glueballs with various K.*/

     for(i=0, iist=0; i<x.nsectors; i++){
      temp=xy[0].basis[i]->locate;
      npi=x.np[i];
      for(ist=0; ist<x.basis[i]->length; ist++, temp+=npi)
	for(opi=0; opi<2; opi++)
	  if(x.inner[i][ist]!=0){
	    for(j=0, k=0; j<npi; j++)k+=((temp[j]&KMASK)<<1)+1;
	    for(jkt=0; jkt<nlx; jkt++)
	      if(lx[jkt].kt==k)break;
	    memcpy(part,temp,npi*sizeof(partype));
	    wellorder(part,npi,gcyclic,go,subgroup(gmulti));
	    for(j=0; j<lx[jkt].nsectors; j++)
	      if(lx[jkt].np[j]==npi)break;
	    for(jst=0; jst<lx[jkt].basis[j].length; jst++)
	    iist++;
	}
     }

  goto loop;    

cleanup:
  fclose(fp);
  return 0;
}