void recursive(){
  double **logpartfunc;
  double partfunc,prob,ncheck,zcheck;
  int n,z,a,ntot,ztot,atot;
  int **existence;

  logpartfunc=new double *[NSYSTEM+1];
  existence=new int *[NSYSTEM+1];
  for(n=0;n<=NSYSTEM;n++){
    logpartfunc[n]=new double[ZSYSTEM+1];
    existence[n]=new int[ZSYSTEM+1];
  }
  for(ntot=0;ntot<=NSYSTEM;ntot++){
    for(ztot=0;ztot<=ZSYSTEM;ztot++) existence[ntot][ztot]=0;
  }

  /* First calculate partition functions */

  logpartfunc[0][0]=0.0;
  existence[0][0]=1;
  nucl_info[0][0]->nlevels=1;
  for(ntot=0;ntot<=NSYSTEM;ntot++){
    for(ztot=0;ztot<=ZSYSTEM;ztot++){
      atot=ntot+ztot;
      if(atot>0){
	partfunc=0.0;
	for(n=0;n<=ntot;n++){
	  for(z=0;z<=ztot;z++){
	    a=n+z;
	    if(a>0 && nucl_info[n][z]->nlevels>0
	       && existence[ntot-n][ztot-z]==1){
	      partfunc=partfunc
		+double(a)*exp(nucl_info[n][z]->logomega
			       +logpartfunc[ntot-n][ztot-z]);
	      existence[ntot][ztot]=1;
	    }
	  }
	}
	if(existence[ntot][ztot]==1){
	  logpartfunc[ntot][ztot]=log(partfunc/double(atot));
	}
      }
    }
  }
  
  /* Now calculate probability of specific nuclei */
  cout << form("Beginning calc. of mass dist.s \n");
  ncheck=0.0;
  zcheck=0.0;
  for(n=0;n<=NSYSTEM;n++){
    for(z=0;z<=ZSYSTEM;z++){
      a=n+z;
      if(a>0 && nucl_info[n][z]->nlevels>0){
	prob=0.0;
	if(existence[n][z]==1 && existence[NSYSTEM-n][ZSYSTEM-z]==1)
	  prob=exp(nucl_info[n][z]->logomega+logpartfunc[NSYSTEM-n][ZSYSTEM-z]
		   -logpartfunc[NSYSTEM][ZSYSTEM]);
	ncheck=ncheck+n*prob;
	zcheck=zcheck+z*prob;
	nucl_info[n][z]->initialsum=prob;
      }
    }
  }

  // The net numbers of neutrons and protons should equal the number initially 
  // entered. This is a simple error check.
  cout << form("CHECK: net neutrons=%g  net protons=%g\n",ncheck,zcheck); 
  if(!(fabs(ncheck-double(NSYSTEM))<0.001 
       && fabs(zcheck-double(ZSYSTEM))<0.001)){
    cout << form("Does not Add up!\n");
    exit(1);
  }
}  

// This function writes the population of each nuclide before/after decays.
void print_summary(int cutoff){
  int n,z,a;
  ofstream summary_file;
  summary_file.open(form("results/summary_%s.dat",RESULTSFILEPREFIX));
  summary_file << "! z   a   finalsum   initialsum\n";
  for(z=0;z<=ZSYSTEM;z++){
    for(n=0;n<=NSYSTEM;n++){
      if(n+z>0 && nucl_info[n][z]){
	a = n+z;
	if(a<=cutoff){
	  // Only print out summary of nuclides with pops > 1.0E-50.
	  if(nucl_info[n][z]->finalsum > LOWERLIMIT
	     || nucl_info[n][z]->initialsum > LOWERLIMIT)
	    summary_file << setw(4) << z << setw(4) << a
			 << setw(15) << nucl_info[n][z]->finalsum
			 << setw(15) << nucl_info[n][z]->initialsum << endl;
	}
      }
    }
  }
  summary_file.close();
}