#include #include #include #include #include #include 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_levels.tmp"); 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 << ie << " " << pair_rho[ie] << " " << endl; } levels_file.close(); /* Part four: Write to a series of files the energy and the number of ways to distribute the energy. */ for(a=2;a<=250;a++){ levels_file.open(form("a%d.tmp",a)); levels_file << form("!a=%d\n",a); dele=PI*PI*LEVELPAR/double(3*a); if(dele>dele_min){ nlevels=1+int((10.0+8.0*double(a))/dele); if(nlevels>NEMAX) nlevels=NEMAX; levels_file << form("0 %d\n",nlevels); for(ie=0;ieNEMAX) nlevels=NEMAX; levels_file << form("0 %d\n",nlevels); for(ie=0;ie