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

/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
/* PARAMETERS YOU MIGHT BE INTERESTED IN ADJUSTING */
const int NBABIES=8;
const int NBABY[NBABIES]={1,0,1,2,0,2,1,2};
const int ZBABY[NBABIES]={0,1,1,0,2,1,2,2};
const double ERESOLUTION=1.0;
/* For calculating nuclear decays*/
const int MAXLEVELS=1000; //Maximum number of levels for an individual nucleus
/* Sets limit on neglecting low populations */
const double LOWERLIMIT=1.0E-20;
/* These are used to calculate partition functions of individual nuclei */
const double RHO_RATIO=6.0; // gives breakup density (rho0/rho)
const double AVAILABLE_VOL_FRACTION=1.0-2.0/RHO_RATIO;
const int A_MAKEBELIEVE=6;
/* If a>=A_MAKEBELIEVE levels are generated from levels in FAKEPATH dir.
   If a<A_MAKEBELIEVE levels are read from data in REALPATH dir. */
const double LEVELPAR=10.0; /* Level density parameter in MeV
//#define WRITE_POPS
/* Define this if you wish populations of individual levels to be written */

/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
/* THINGS YOU WOULD NEVER CHANGE */
const double PI=3.14159265358979323844;
const double HBARC=197.326;
const double ONETHIRD=1.0/3.0;
const double TWOTHIRDS=2.0/3.0;
const double RHO0=0.15; // nuclear saturation density in fm^-3

/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
/* Global variables */
int ZSYSTEM,NSYSTEM,ASYSTEM; /* These and the temperature are user entered and 
				refered to throughout the program. */
double TEMPERATURE;
/* This is only a string used by several subroutines to label output files */
char RESULTSFILEPREFIX[100];

/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
/* Class Definitions and prototypes */
class nucl_info_class{
public:
  nucl_info_class(int aenter,int zenter);
  ~nucl_info_class();
  int n,z;
  int nlevels;
  int nstablelevels;
  int filled;
  double finalsum,initialsum; // weights summed over all levels
  double be0;// binding energy without Coulomb screening.
  double omega; // Found by omegafinder
  double logomega; // Filled by omegafinder
  double *energy;  // Filled by levelinit
  double *degen;   // Also filled by levelinit
  double *finalprob,*initialprob; // levelfiller
  void printlevels();
  void levelinit();
  void omegafinder(double *degen_tmp,double *energy_tmp);
  void levelfiller();
  void frldm();
  void printinitialpops();
  void printfinalpops();
  void collapse();
  double coulombcorr();
};

nucl_info_class::nucl_info_class(int n_enter,int z_enter) {
  n=n_enter;
  z=z_enter;
  nlevels=0;
  nstablelevels=0;
  filled = 0;
  initialsum=0;
  energy=0;
  degen=0;
  finalprob=0;
  initialprob=0;
  omega=logomega=-1.0;
}

nucl_info_class::~nucl_info_class() {
  if(energy) delete [] energy;
  if(degen) delete [] degen;
  if(finalprob) delete [] finalprob;
  if(initialprob) delete [] initialprob;
  energy = finalprob = initialprob = degen = 0; 
  // This sets the pointers to NULL
}

nucl_info_class ***nucl_info;

void calc_decay_weights(nucl_info_class *mother,int momlevel,
			int ibaby,double decayweight[NBABIES][MAXLEVELS],
			double *ndecaylevels,double *weight_norm);
void add_decay_weights(nucl_info_class *mother,int momlevel,
		       int ibaby,double decayweight[NBABIES][MAXLEVELS],
		       double *ndecaylevels,double weight_norm);
void recursive();
void decay_main(int n, int z, int cutoff);
void decay(nucl_info_class *mother);
void print_summary(int cutoff);