// This function reads in the energies, degeneracies and binding energy
// for the nuclide n,z. It stores the energies and degeneracies in arrays
// called energy_tmp and degen_tmp. The arrays are sent to omega finder,
// where they are used for calculating the single-fragment partition functions
// In omegafinder, they are also cropped to the right size and stored 
// as energy and degen respectively for later use.
void nucl_info_class::levelinit(){
  double ek0,e,weight0,levelprob,dtemp,dtemp1,dtemp2;
  double levelsep,dele;
  double energy_tmp[MAXLEVELS],degen_tmp[MAXLEVELS];
  static double *generic_degen;
  static int nfakelevels;
  int ie1,ie2,a,ifake,ilevel;
  char filename[30];
  char dummy[100];
  ifstream real_levels_file;
  ifstream fake_levels_file;
  nlevels=0;
  a=n+z;
  if(a<A_MAKEBELIEVE){
    real_levels_file.open(form("real_levels/a%dz%d.tmp",a,z));
      real_levels_file.getline(dummy,80);
      if(real_levels_file){
	real_levels_file >> ek0 >> nlevels;
	if(energy || degen){
	  cout << form("z=%d a=%d, energy or degen already exists??\n",z,a);
	  //exit(1);
	}
	for(ilevel=0;ilevel<nlevels;ilevel++){
	  real_levels_file >> e >> dtemp;
	  if(ilevel==0) be0=e;
	  energy_tmp[ilevel]=e;
	  degen_tmp[ilevel]=dtemp;
	}
      }
      real_levels_file.close();
  }
  else{
    if(!generic_degen){
      fake_levels_file.open("fake_levels/generic.tmp");
      fake_levels_file >> nfakelevels;
      generic_degen=new double[nfakelevels];
      for(ifake=0;ifake<nfakelevels;ifake++)
	fake_levels_file >> generic_degen[ifake];
      fake_levels_file.close();

    }

    frldm(); //finds binding energy be0
    levelsep=PI*PI*LEVELPAR/double(3*a);

    if(levelsep>ERESOLUTION){
      dele=levelsep;
      nlevels=MAXLEVELS;
      if(nfakelevels<MAXLEVELS) nlevels=nfakelevels;
      for(ilevel=0;ilevel<nlevels;ilevel++){
	degen_tmp[ilevel]=generic_degen[ilevel];
	energy_tmp[ilevel]=be0+ilevel*dele;
      }
    }
    else{
      dele=ERESOLUTION;
      for(ifake=0;ifake<MAXLEVELS;ifake++){
	degen_tmp[ifake]=energy_tmp[ifake]=0.0;
      }
      degen_tmp[0]=generic_degen[0];
      ifake=1;
      do{
	dtemp=generic_degen[ifake];
	ie1=int(double(ifake)*levelsep/dele);
	ie2=ie1+1;
	dtemp1=dtemp*(double(ie2)*dele-ifake*levelsep)/dele;
	dtemp2=dtemp-dtemp1;
	if(ie1<MAXLEVELS){
	  degen_tmp[ie1] += dtemp1;
	}
	if(ie2<MAXLEVELS){
	  degen_tmp[ie2] += dtemp2;
	}
	ifake+=1;
      } while(ie1<=MAXLEVELS && ifake<nfakelevels);
      nlevels=ie1;
      if(nlevels>MAXLEVELS) nlevels=MAXLEVELS;
      for(ilevel=0;ilevel<nlevels;ilevel++){
	energy_tmp[ilevel]=be0+dele*ilevel;
      }
    }
  }

  omegafinder(degen_tmp,energy_tmp);
  if(nlevels<0) exit(1);

}

#define MUA -6 /*Divides out of all observables, but helps reduce
		 problems with numerical overflow 
		 Do not put it below -9, or you might be
		 throwing away good parts of the ensemble. */
// This function finds the single-fragment partition function, crops
// the energy_tmp and degen_tmp arrays and stores them as energy and degen.
void nucl_info_class::omegafinder(double *degen_tmp,double *energy_tmp){
  double ek0,e,weight0,relweight0,relweight;
  double volume,internalpart,exparg;
  int a,ilevel,nlevels0,relcheck;
  nlevels0=nlevels;
  if(nlevels>0){
    a=n+z;
    internalpart=0.0;
    ilevel=0;
    relweight0=0.0;
    do{
      relcheck=1;
      e=energy_tmp[ilevel]+coulombcorr();
      exparg=-(e-MUA*a)/TEMPERATURE;
      if(exparg>600){
	cout << form("n=%d z=%d ilevel=%d exparg=%g\n",n,z,ilevel,exparg);
	cout << form("e/a=%g\n",e/a);
	cout << form("Possible Overflow might result from large arguments\n");
	cout << form("in exponentials, try reducing MUA.\n");
	cout << form("This problem can appear when running at low T\n");
	cout << form("It is not reccommended to run EASY for T<1.5 MeV\n");
	exit(1);
      }
      if(exparg>-600){
	relweight=degen_tmp[ilevel]*exp(exparg);
	internalpart+=relweight;
	if(relweight>relweight0) relweight0=relweight;
	ilevel+=1;
	if(relweight/relweight0<LOWERLIMIT && degen_tmp[ilevel]>1.0E-5)
	  relcheck=0;
      }
      else{
	// This attempts to go around underflow problems. But,
	// be careful. If MUA is set below -9 MeV, you might be throwing
	// away good parts of the ensemble!
	relcheck=0;
      }
    } while(ilevel<nlevels0 && relcheck==1);
    nlevels=ilevel;

    energy=new double[nlevels];
    degen=new double[nlevels];
    for(ilevel=0;ilevel<nlevels;ilevel++){
      energy[ilevel]=energy_tmp[ilevel];
      degen[ilevel]=degen_tmp[ilevel];
    }

    volume=(NSYSTEM+ZSYSTEM)*AVAILABLE_VOL_FRACTION*RHO_RATIO/RHO0;
    weight0=volume*pow(931.5*a*TEMPERATURE
		       /(2.0*PI*HBARC*HBARC),1.5);
    omega=internalpart*weight0;
    logomega=log(omega);
  }
}


// This function finds the initial weight of each level of each level.
// It also sets the flag nucl_info[n][z]->filled to 1.
void nucl_info_class::levelfiller(){
  double ek0,e,weight0,levelprob,levelcheck,volume;
  int ie,a;

  if (nlevels>0){
    a=n+z;
    volume=(NSYSTEM+ZSYSTEM)*AVAILABLE_VOL_FRACTION*RHO_RATIO/RHO0;
    weight0=volume*pow(931.5*a*TEMPERATURE
		     /(2.0*PI*HBARC*HBARC),1.5);
    finalprob= new double[nlevels];
    initialprob= new double[nlevels];
    // When created, the energy, weight, etc. pointers are all null. This 
    // creates something for them to point at.
    levelcheck=0.0;
    for(ie=0;ie<nlevels;ie++){
      e=energy[ie]+coulombcorr();
      levelprob=degen[ie]*weight0
	*exp((a*MUA-e)/TEMPERATURE-logomega);
      levelcheck=levelcheck+levelprob;
      levelprob=levelprob*initialsum;
      initialprob[ie] = levelprob;
      finalprob[ie] = levelprob;
    }
    if(fabs(levelcheck-1.0)>1.0E-5 && nlevels>0){
      cout << form("does not add up in levelfiller, n=%d, z=%d",n,z);
      cout << form("initialsum=%g,levelcheck=%g,omega=%g\n",
		   initialsum,levelcheck,omega);
      exit(1);
    }
    filled=1;
  }
}
#undef MUA

/* This function finds the binding energy of a nucleus with N neutrons 
   and Z protons. It is only used with "cheap_nucl_levels" */
void nucl_info_class::frldm(){
  const double ael=1.433E-5,rp=0.80,r0=1.16,a=0.68,aden=0.70;
  const double rmac=4.80,h=6.6,W=30.0,av=16.00126,kv=1.92240,as=21.18466;
  const double ks=2.345,a0=2.615,ca=0.10289;
  const double e2=1.4399764,Mn=8.071431,MH=7.289034,Mel=0.511;

  int A;
  double x0,y0,B1,B3,kf,I,f,Deltabar_n,Deltabar_p,delta_np,Bs,Bc,c1,c4;
  double coulomb;

  double Nonethird,Zonethird,Aonethird,deltasymm;
  int noddeven,zoddeven;

  double E;
  A=z+n;
  Nonethird=pow(double(n),1.0/3.0)+1.0E-6;
  Zonethird=pow(double(z),1.0/3.0)+1.0E-6;
  Aonethird=pow(double(A),1.0/3.0);
    
  I=double(n-z)/double(A);
  x0=r0*Aonethird/a;
  y0=r0*Aonethird/aden;
  B1=1.0-3.0/(x0*x0)+(1.0+x0)*(2.0+(3.0/x0)+3.0/(x0*x0))*exp(-2.0*x0);
  B3=1.0-(5.0/(y0*y0))*(1.0-15.0/(8.0*y0)+21.0/(8.0*y0*y0*y0)
			-0.75*(1.0+9.0/(2.0*y0)+7.0/(y0*y0)+7.0/(2.0*y0*y0))
			*exp(-2.0*y0));
  kf=pow(9.0*PI*double(z)/(4.0*double(A)),1.0/3.0)/r0;
  Bs=B1;Bc=B3;
  c1=(3.0/5.0)*e2/r0;
  c4=(5.0/4.0)*c1*pow(3.0/(2.0*PI),2.0/3.0);
  f=-(1.0/8.0)*rp*rp*e2/(r0*r0*r0);
  f=f*(145.0/48.0-(327.0/2880.0)*kf*rp*kf*rp
       +(1527.0/1209600.0)*kf*rp*kf*rp*kf*rp*kf*rp);
  Deltabar_n=rmac*Bs/Nonethird;
  Deltabar_p=rmac*Bs/Zonethird;
  delta_np=h/(Bs*Aonethird*Aonethird);
  E=0.0;
  //E=MH*double(z)+Mn*double(n);
  E=E-av*(1.0-kv*I*I)*double(A)
    +as*(1.0-ks*I*I)*B1*(Aonethird*Aonethird);
  E=E+a0;
  coulomb=c1*double(z*z)*B3/Aonethird-c4*double(z)*Zonethird/Aonethird
    +f*double(z*z)/double(A);
  E=E+coulomb;
  E=E-ca*double(n-z);
  noddeven=zoddeven=-1;
  if(2*(n/2) == n) noddeven=1;
  if(2*(z/2) == z) zoddeven=1;
  if(noddeven==-1 && zoddeven==-1) deltasymm=Deltabar_p+Deltabar_n-delta_np;
  if(noddeven==-1 && zoddeven==1) deltasymm=Deltabar_n;
  if(noddeven==1 && zoddeven==-1) deltasymm=Deltabar_p;
  if(noddeven==1 && zoddeven==1) deltasymm=0.0;
  E=E+deltasymm;
  E=E+W*fabs(I);
  if(noddeven=-1 && zoddeven==-1 && n==z) E=E+W/double(A);
  //E=E+ael*pow(double(z),2.39);
  //cout << form("E=%g, E/A=%g\n",E,E/double(A));
  be0=E;
}

// This function writes energies, degeneracies, and population of nuclei 
// to a file called a#z#.tmp, where the two #'s are n+z and z respectively.
void nucl_info_class::printfinalpops(){
  int a,ilevel;
  ofstream finalpops_file;
  a=n+z;
  finalpops_file.open(form("results/finalpops/%s_z%da%d.tmp",
			   RESULTSFILEPREFIX,z,a));
  finalpops_file << form("!nlevels=%d\n",nstablelevels);
  finalpops_file << nstablelevels << endl;
  for(ilevel=0;ilevel<nstablelevels;ilevel++){
    finalpops_file << energy[ilevel]
		 << " " << degen[ilevel]
		 << " " << finalprob[ilevel] << endl;
  }
  finalpops_file.close();
}
void nucl_info_class::printinitialpops(){
  int a,ilevel;
  ofstream initialpops_file;
  a=n+z;
  initialpops_file.open(form("results/initialpops/%s_z%da%d.tmp",
			     RESULTSFILEPREFIX,z,a));
  initialpops_file << form("!nlevels=%d\n",nstablelevels);
  initialpops_file << nlevels << endl;
  for(ilevel=0;ilevel<nlevels;ilevel++){
    initialpops_file << energy[ilevel]
		 << " " << degen[ilevel]
		 << " " << initialprob[ilevel] << endl;
  }
  initialpops_file.close();
}

void nucl_info_class::printlevels(){
  int ilevel,a;
  a=z+n;
  cout << "----------PRINTING LEVELS-----------------\n";
  cout << form("n=%d  z=%d\n",n,z);
  for(ilevel=0;ilevel<nlevels;ilevel++){
    cout << form("%d  %g   %g\n",ilevel,energy[ilevel],degen[ilevel]);
  }
}

void nucl_info_class::collapse() {
  if(energy) delete [] energy;
  if(finalprob) delete [] finalprob;
  if(initialprob) delete [] initialprob;
  if(degen) delete [] degen;
  energy = finalprob = initialprob = degen = 0;
  filled=0;
}

double nucl_info_class::coulombcorr(){
  int a;
  double coulomb,answer;
  a=n+z;
  coulomb=0.7*double(z*z)/pow(a,ONETHIRD);
  answer=-coulomb*pow(RHO_RATIO,-ONETHIRD);
  return answer;
}