#include <stream.h>
#include <iostream.h>
#include <iomanip.h>
#include <fstream.h>
#include <stdlib.h>
#include <math.h>

const double PI=3.14159265358979323844;

long IDUM=-12345; double ran2(void);
const int NMAX=11;
const int ZMAX=11;
const int AMAX=11;
// LEVELPAR has default value 10.0
const double LEVELPAR=10.0;
const double WELLDEPTH=8.0;
const int NEMAX=10000;
#define EXCLUDE_CONTINUUM

void main(){
  int n,a,ie,pe,ne,iiee,sign,nlevels;
  static double rho[AMAX+1][NEMAX+1]={0};
  static double pair_rho[NEMAX+1]={0};
  static double netrho[NEMAX+1]={0};
  static double netrho_prime[NEMAX+1]={0};
  double dele,guess,dele_min=1.0;
  ofstream levels_file;

  /* Part the firste: find the number of ways to distribute a integers to 
     sum to a total equal to ie */
  rho[0][0]=1;
  for(a=1;a<=AMAX;a++){
    for(ie=0;ie<=100;ie++){
      sign=-1;
      for(n=1;n<=a;n++){
	sign=-sign;
	iiee=0;
	while(ie-n*iiee>=0){
	  rho[a][ie]=rho[a][ie]+sign*rho[a-n][ie-n*iiee];
	  iiee=iiee+1;
	}
      }
      rho[a][ie]=rho[a][ie]/a;
    }
  }

  /* Part two: find the number of ways to distibute an infinite number of
     integers with "activation energy" ie. */
  levels_file.open("generic.tmp");
  levels_file << NEMAX << endl;
  for(ie=0;ie<=NEMAX;ie++){
    if(ie<=100){
      for(a=0;a<=AMAX;a++){
	for(iiee=a;iiee<=ie;iiee++){
	  netrho[ie]=netrho[ie]+rho[a][iiee-a]*rho[a][ie-iiee];
	}
      }
    }
    /* Part three: approximate the number of ways as "guess", write to 
       a file the activation energy, number of ways, and ratio between 
       guess and the number of ways. For larger activaton energies, just
       write the approximation. */
    guess=exp(2.0*sqrt(PI*PI*double(ie)/6.0));
    guess=guess/(sqrt(48.0*double(ie*ie)));
    guess=guess/sqrt(1.0+1.0/sqrt(double(ie)));
    guess=guess/0.99704;
    if(ie>100) netrho[ie]=guess;
    for(pe=0; pe<=ie; pe++){
      ne = ie - pe;
      pair_rho[ie]+=netrho[pe]*netrho[ne];
    }
    levels_file << pair_rho[ie] << endl;
  }
  levels_file.close();
  
}