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

/*  This file calculates various properties of the 
    wavefunction associated with a given eigenstate
    or states.   */

#define MOM 0 /* Flag for nonzero P_1 or L  */
#define GRENZE 1.0e-3 /* minimium amount to print out */

int MAIN(){
  int i,iist,ist,k,npi/*,hlist[HINDEX]*/,vv;
  size_t j;
  sectors x={0};
  /*           1 */
  clock_t t1;
  static struct {int n; partype *s; element *p;} shapes[NPMAX];
  partype *temp,shape[NPMAX],p;
  element beta,structure[KKMAX],*pvec,rescale2,rest;
#define NOUT 100 /* Output file name length */
  char out[NOUT]="structure.out",*pout;
  FILE *fp;
#if HINDEX==1
#if MOM
  int   ht[HINDEX]={0},np=6,kt=6,cyclic=1,o=1,multi=1,multi2=2,nfig=8;
  element kmax=3.0,ll=3.0;
#else /*  1             0     24    8       1      8   */
  int   o=1,ht[HINDEX]={0},kt=6,np=4,multi=0,nfig=12,cyclic=1;
#endif
#if 1
  element couplings[NPARAMS]={0.1,2.0,-0.2,-0.23,20.0,0.0};
#else
  element couplings[NPARAMS]={0.0324919696232906307,1.0,0.0,0.0,0.0,-2.0};
#endif

#elif HINDEX==2

#if MOM
  /*                              8 26        */
  int     ht[HINDEX]={0,0},np=6,kt=26,cyclic=1,o=1,multi=8,multi2=9,nfig=8;
  element kmax=4.0,ll=4.0;
#else  /*   1                    24    8       1      8   */
  int     ht[HINDEX]={0,0},np=6,kt=36,cyclic=1,o=1,multi=3,nfig=9;
#endif
#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 0 /* couplings used in first 3+1 paper */
  element rescale=8.099820052, couplings[NPARAMS]={0.0517766953,1,0.7584261311,
	    -0.07623059914,-0.1458991891,8.889796764,0.1608104465,
            1.657995607,-3.268403913};
#elif 1  /* couplings from best fit in Spring 1999  */
  element rescale=11.4617271, couplings[NPARAMS]={0.001967542671,1,0.3909839522,
	  -0.01250282378,-0.06421231123,353.2492529,0.07049127522,160.6273835,
          -0.7123520719,-0.6519355869};
#endif
#endif
  integer nval=5;
  doublereal *val,*vec=NULL;
#if MOM
  struct sectors y={0};
  int ipt[HINDEX];
  element transmom[HINDEX];
#endif

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

#if MOM
  for(i=0;i<HINDEX;i++)ipt[i]=0;
#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;
  }

  /* Main loop of program */
  for(;;){

    printf("Input number of eigenvalues (0=exit):  ");
    scanf("%li",&nval);
    if(nval<1)break;
    
    
    /******  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 MOM
    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("%f",&multi2);
    }
#endif
#endif
    
    /*  construct basis, any previously constructed basis in x is
	deallocated  */
    printf("Constructing basis");
    t1=clock();
    initbasis(&x,ht,np,kt,cyclic,o,multi);
    makebasis(&x);
#if MOM
    if(o==0){
       fprintf(stderr,__FILE__":  Need nonzero o, exititing\n");
      exit(10);
    }
    if(kt<0&&!cyclic)
      initbasis(&y,ht,np,kt,cyclic,-o,multi);
    else
      initbasis(&y,ht,np,kt,cyclic,-o,multi2);
    makebasis(&y);
#endif
    printf(" in %f CPU seconds\n",
	   (float)(clock()-t1)/(float) CLOCKS_PER_SEC);
    
#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
    printf("Setting up tables");
    t1=clock();
    beta=mapcouplings(kt,couplings);
    
#if MOM
    if(x.kt<0 && !x.cyclic){
      maketenslist(ll,x.kt,kmax,beta); 
    } else {    
      /*  check for nonzero winding number */
      if(!iswinding(x.ht)){
	printf("Input transverse momenta (%i): ",HINDEX);
	for(i=0; i<HINDEX; i++)scanf("%lg",transmom+i);
	/*printf("P_1 a = %e\n",transmom);*/
      } else {
	if(multi!=0) fprintf(stderr,__FILE__":  Need multi=0 for nonzero winding\n");
	printf("Input integer transverse momenta; to (");
	for(i=0; i<HINDEX; i++)printf("%i,",ht[i]-1);
	printf("):");
	for(i=0; i<HINDEX; i++)scanf("%i",ipt+i);
	for(i=0; i<HINDEX; i++)
	  transmom[i]=2.0*(atan(1.) * 4.0)*(element) ipt[i]/(element) ht[i];
      }
      rebasis(&x,ipt,transmom);
      rebasis(&y,ipt,transmom);
      makephaselist(transmom,kt);
    }
#endif
    printf(" in %f CPU seconds\n",
	   (float)(clock()-t1)/(float) CLOCKS_PER_SEC);
    
    
    printf("Starting calculation of spectrum.\n");
    t1=clock();
#if MOM
    spectrum(&x,&y,&val,&vec,nval,couplings,nfig);
#else
    spectrum(&x,&x,&val,&vec,nval,couplings,nfig);
#endif
    printf("%f CPU seconds to calculate Hamiltonian and run Lanczos\n\n",
	   (float)(clock()-t1)/(float) CLOCKS_PER_SEC);
    
    /* use right resale for nonzero L */
    
    rescale2=cyclic?rescale:sqrt(rescale);
    
    for(vv=0; vv<nval; vv++)
      printf("Eigenvalue %21.14e\n",rescale2*val[vv]);
    
    /*  Print out parameters for basis and couplings */
    
#if HINDEX==1
    fprintf(fp,"***** Eigenfunctions 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,"****** Eigenfunctions 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 && kt<0)
      fprintf(fp,"\nHeavy sources, K_max = %f\n",kmax);
#else
    fprintf(fp," with multiplet=%i and orientation symmetry %i.\n",multi,o);
#endif
    fprintf(fp,"Transverse shapes are distinguished to within symmetries\n"
	    " of the lattice subgroup chosen.\n"
	    "Configurations with probability greater than %f are shown.\n\n",
	    GRENZE);
    
#if 1 && !MOM /* Calculate structure functions, etc. */
    for(vv=0; vv<nval; vv++){
      fprintf(fp,"Eigenvalue %21.14e\n",rescale2*val[vv]);
      pvec=vec+x.length*vv;
      for(i=0; i<KKMAX; i++)structure[i]=0.0;
      beta=0.0;
      for(i=0, iist=0; i<x.nsectors; i++){
	temp=x.basis[i]->locate;
	shapes[i].n=0;
	npi=x.np[i];
	for(ist=0; ist<x.basis[i]->length; ist++, temp+=npi)
	  if(x.inner[i][ist]!=0){
	    /* Find contribution of current state to structure function */
	    for(j=0; j<npi; j++)
	      structure[temp[j]&KMASK]+=pow(pvec[iist],2);
	    /* Find transverse structure of current state */ 
	    for(j=0; j<npi;j++)
	      shape[j]=temp[j]&(~KMASK);
	    wellorder(shape,npi,x.cyclic,x.o,subgroup(x.multi));
	    for(j=0; j<shapes[i].n; j++){
	      for(k=0; k<npi; k++)
		if(shape[k]!=shapes[i].s[npi*j+k])break;
	      if(k==npi){
		shapes[i].p[j]+=pow(pvec[iist],2);
		break;
	      }
	    }
	    if(j==shapes[i].n){
	      shapes[i].n++;
	      shapes[i].s=realloc(shapes[i].s,
				  npi*shapes[i].n*sizeof(partype));
	      shapes[i].p=realloc(shapes[i].p,shapes[i].n*sizeof(element));
	      for(k=0; k<npi; k++)
		shapes[i].s[npi*j+k]=shape[k];
	      shapes[i].p[j]=pow(pvec[iist],2);
	    }
	    /* increment current state (beta is a debugging check) */
	    beta+=pow(pvec[iist],2);
	    iist++;
	  }
      }
      
      /* Print out results for wavefunction */
      
#if 1  /* debugging print to check probability adds to 1 */
      fprintf(fp,"Probability = %f\n",beta);
#endif
      fprintf(fp,"structure function:\n");
      for(i=0; 2*i+1<=x.kt; i++)fprintf(fp,"  %i/2 %f\n",2*i+1,structure[i]);
      fprintf(fp,"Transverse shape distribution:\n");
      rest=0.0;
      for(i=0; i<x.nsectors; i++){
	npi=x.np[i];
	for(j=0; j<shapes[i].n; j++)
	  if(shapes[i].p[j]>GRENZE){
	    fprintf(fp,"  %f ",shapes[i].p[j]);
	    for(k=0; k<npi; k++){
	      p=shapes[i].s[j*npi+k]>>HSHIFT;
	      /*          This expression should agree with the output of gluedecode() */
	      fprintf(fp," %i",(p&1)?+((p>>1)+1):-((p>>1)+1));
	    }
	    fprintf(fp,"\n");
	  } else
	    rest += shapes[i].p[j];
      }
      fprintf(fp,"  %f  other shapes\n\n",rest);
    }
#if 0 /* print out helicity map */
    fprintf(fp,"Helicity map:\n");
    for(i=0; i < 2*HINDEX; i++){
      p=i<<HSHIFT;
      gluedecode(p,hlist);
      fprintf(fp,"  %i: ",i);
      for(j=0; j<HINDEX; j++) fprintf(fp,"%i ",hlist[j]);
      fprintf(fp,"\n");
    }
#endif
    fprintf(fp,"\n");
#endif /* end Calculate structure function */

    freebasis(&x);
#if MOM
    freebasis(&y);
#endif

  }   /* End main loop */
    
cleanup:
  fclose(fp);
  return 0;
}