/*
  This program calculates the amount of various nuclides produced in a 
  canonical ensemble of nucleons. It first populates the various states that 
  would exist immediately after a collision, then follows them through a 
  decay. 
*/

#include "easy.h"

void main(){
  int a,z,n;
  int cutoff;

  cout << "How many protons(Z) and nucleons(A) in the system?\n";
  cin >> ZSYSTEM >> ASYSTEM;
  NSYSTEM=ASYSTEM-ZSYSTEM;

  nucl_info=new nucl_info_class**[NSYSTEM+1];
  for(n=0;n<=NSYSTEM;n++){
    nucl_info[n]=new nucl_info_class*[ZSYSTEM+1];
    for(z=0;z<=ZSYSTEM;z++) nucl_info[n][z]=new nucl_info_class(n,z);
  }

  cout << "What is the TEMPERATURE? (in MeV)\n";
  cin >> TEMPERATURE;

  strcpy(RESULTSFILEPREFIX,
	 form("Z%d_A%d_T%g",ZSYSTEM,ZSYSTEM+NSYSTEM,TEMPERATURE));


  cout << "What is the decay cutoff?" << endl;
  cout << "(Will only decay masses <= cutoff to save CPU)" << endl;
  cin >> cutoff;
  if(cutoff>NSYSTEM+ZSYSTEM) cutoff=NSYSTEM+ZSYSTEM;

  /* Now, initialize levels which means finding the energy levels
     and obtaining partition functions of individual nuclei.
     The collapse function erases the energy and degen arrays after they
     have been used to obtain the partition functions */
  for(n=0;n<=NSYSTEM;n++){
    for(z=0;z<=ZSYSTEM;z++){
      if(n+z>0){
	nucl_info[n][z]->levelinit();
	nucl_info[n][z]->collapse();
      }
    }
  }

  /* Using the partition functions calculated above, recursive methods
     are now employed to generate partition functions for the entire 
     system and obtain the probabilities of specific nuclei */
  cout << form("Beginning recursive calc.s of all part. functions \n");
  recursive();

  /* Nuclei with masses <= cutoff are decayed */
  cout << "Beginning Decay Process\n";
  decay_main(NSYSTEM,ZSYSTEM,cutoff); 

  /* Prints the Summary */
  print_summary(cutoff);

}

#include "easysubs.cc"
#include "decaysubs.cc"
#include "nucl_info.cc"