/************************************************************************/
/*                                                                      */
/*  Program    :     t  H o o f t . c p p                               */
/*                                                                      */
/*  Creation   :     9 VIII 2001                                        */
/*                                                                      */
/*  Class      :     Main program                                       */
/*                                                                      */
/*  Programmer :     Uwe Trittmann                                      */
/*                                                                      */
/*  Purpose    :     PGR for calculating the 't Hooft model spectrum in */
/*                   DLCQ                                               */
/*                                                                      */
/************************************************************************/

#include <malloc.h>
#include <stdio.h>
#include <math.h>
#include <unistd.h>
#include <signal.h>
#include <iostream>
#include "NAGdiagII.h"         /* diagonalization routines */

#define PI 3.141592653589793

double *Hamiltonian_Matrix;    /* Hamiltonian matrix */
Complex *r,*v;	               /* evs, efs for NAGLIB diagonalization */
int nos;                       /* number of basis states */

//-----------------------------------------------------------------------

void generate_Hamiltonian(int K, double mu, int test)
// generates Hamiltonian according to old "Two_particles.c" code 
// at harmonic resolution <K>, <test>=1 displays Hamiltonian
{
  int i,j,k,m;
  double sum;
  
  if (K%2==0) nos=K/2+1;
  else nos=(K-1)/2+1;
  printf("\nNumber of basis states: %3d\n",nos);

  Hamiltonian_Matrix=dvector(0,(nos-1)*(nos-1)-1);

  m=0;
  for (i=0; i<nos-1; ++i)
    for (j=0; j<nos-1; ++j)
      {
	Hamiltonian_Matrix[m]=0.0;
	++m;
      }
    
  for (i=1; i<nos; ++i)
    {
      sum=0.0;
      for (m=1; m<i; ++m)   sum+=1.0/m/m;
      for (m=1; m<K-i; ++m) sum+=1.0/m/m;

      sum+=mu*(1.0/i+1.0/(K-i));    // kinetic term, 't Hooft mass <mu>
      
      Hamiltonian_Matrix[(i-1)*(nos-1)+i-1]+=sum;

      for (m=1; m<i; ++m) 
	{
	  sum=1.0/m/m;
	  k=K-i+m;
	  if (k>=nos) k=i-m;
	  Hamiltonian_Matrix[(i-1)*(nos-1)+k-1]-=sum;      
	}
      for (m=1; m<K-i; ++m) 
	{
	  sum=1.0/m/m;
	  k=i+m;
	  if (k>=nos) k=K-i-m;
	  Hamiltonian_Matrix[(i-1)*(nos-1)+k-1]-=sum;      
	}
    }
  if (test==1)
    {
      m=0;
      for (i=0; i<nos-1; ++i)
	{
	  for (j=0; j<nos-1; ++j)
	    {
	      printf("%5.3f ",Hamiltonian_Matrix[m]);
	      ++m;
	    }
	  printf("\n");
	}
    }
}

void write_results(int K, double mu)
// writes the efs and evs to file <Spectrum_tHooft_K<K>.dat>
{
  int i,j;
  char name[40];
  FILE *fp;

  sprintf(name,"Spectrum_tHooft_K%1dm%3.2f.dat",K,mu);
  fp = fopen(name,"w");

  fprintf(fp,"%3d\n",nos);
  for (i=0; i<nos-1; ++i) fprintf(fp,"%12.8f %12.8f\n",K*r[i].re,K*r[i].im);
  for (i=0; i<nos-1; ++i) 
    {
      fprintf(fp,"%2d  ",i);  
      for (j=0; j<nos-1; ++j) fprintf(fp,"%8.5f ",v[i*(nos-1)+j]);
      fprintf(fp,"\n",i);  
    }
  fclose(fp);

  cout<<endl<<"Results in "<<name<<"."<<endl<<endl;
}


void main()
{
  int i;
  int K=1000;      // harmonic resolution 
  double mu;       // fermion mass 

  cout<<endl;
  cout<<"tHooft.cpp"<<endl;
  cout<<"---------------------------------------------------------------"<<endl;
  cout<<"Calculates the spectrum, i.e. eigenvalues and -functions "<<endl;
  cout<<"of QCD_2 at large N_c."<<endl<<endl;

  cout<<"Harmonic resolution  K="; cin>>K;
  cout<<"Fermion mass        mu="; cin>>mu;  

  generate_Hamiltonian(K,mu,0);
  diagonalize_asymmetric(&Hamiltonian_Matrix,nos-1,&r,&v);

  cout<<endl<<"Eigenvalues:"<<endl;
  for (i=0; i<nos-1; ++i) 
    {
      printf("%2d EV = %12.8f+ i %12.8f\n",i,K*r[i].re,K*r[i].im);
    }
  write_results(K,mu);
}