// This function runs the entire decay process from start to finish. It
// begins by filling the energy levels of the smallest nuclei, those that
// will break off from the nuclei during their decays. Then, it starts from
// the largest nuclei, running them through a decay process. Once it has 
// finished with a nucleus, it might print the populations. Next, it 
// deletes the memory hungry arrays, leaving the population of the nucleus, 
// and then moves on to the next nucleus. Once all nuclei have been decayed, 
// it calls print_summary(), which writes all the results to disk.  
void decay_main(int NSYSTEM, int ZSYSTEM, int cutoff){
  int nmother,zmother,amother,ibaby,zmin,zmax;
  double testfinalsum=0.0,testinitialsum=0.0;
  nucl_info_class *mother,*baby;  

  for(ibaby=0;ibaby<NBABIES;ibaby++){
    baby=nucl_info[NBABY[ibaby]][ZBABY[ibaby]];
    baby->levelinit();
    baby->levelfiller();
  }
  for(amother=cutoff;amother>0;amother--){
    zmax=ZSYSTEM;
    zmin=0;
    if(zmax>amother) zmax=amother;
    if(amother>NSYSTEM) zmin=amother-NSYSTEM;
    for(zmother=zmax;zmother>=zmin;zmother--){
      nmother=amother-zmother;
      mother=nucl_info[nmother][zmother];
      decay(mother);
      if(mother->nlevels>0) {
	testinitialsum+=amother*mother->initialsum;
	testfinalsum+=amother*mother->finalsum;
#ifdef WRITE_POPS
	mother->printinitialpops();
	mother->printfinalpops();
#endif
      }
      mother->collapse();
    }
    cout << form("finished decaying a=%d nuclides\n",amother);
  }
  cout << form("testfinalsum=%g testinitialsum=%g\n",
	       testfinalsum,testinitialsum);
}

// This function takes a value for n and z and decays the nucleus with those 
// values. First, it checks to see if the energy levels of the nucleus are
// filled (nucl_info[n][z]->filled==1). If not, it fills those levels. Next,
// it checks to make sure that the daughter nucleus has filled levels, and 
// creates the levels if not. Next, it calls calc_decay_weights to determine
// the additional population of each of the daughters. Finally, it calls 
// add_decay_weights to add the new population to the current population.
void decay(nucl_info_class *mother){
  int ndaughter,zdaughter,ibaby,motherlevel,nmother,zmother,amother;
  double weight_norm;
  static double decayweight[NBABIES][MAXLEVELS];
  double ndecaylevels[NBABIES];
  nucl_info_class *daughter,*baby;  

  nmother=mother->n;
  zmother=mother->z;
  amother=nmother+zmother;

  if(mother->filled==0){
    mother->levelinit();
    mother->levelfiller();
  }

  for(ibaby=0;ibaby<NBABIES;ibaby++){
    ndaughter=nmother-NBABY[ibaby];
    zdaughter=zmother-ZBABY[ibaby];
    if(ndaughter >= 0 && zdaughter >= 0 && ndaughter+zdaughter>0){
      daughter=nucl_info[ndaughter][zdaughter];
      if(daughter->filled==0 && daughter->nlevels>0) {
	daughter->levelinit();
	daughter->levelfiller();
      }
    }
  }
  if (mother->nlevels>0){
    mother->finalsum=0.0;
    mother->nstablelevels=0;
    for(motherlevel=0;motherlevel<mother->nlevels;motherlevel++){
      weight_norm=0.0;
      for(ibaby=0;ibaby<NBABIES;ibaby++){
	ndaughter=nmother-NBABY[ibaby];
	zdaughter=zmother-ZBABY[ibaby];
	if(ndaughter >= 0 && zdaughter >= 0 && ndaughter+zdaughter>0){
	  daughter=nucl_info[ndaughter][zdaughter];
	  if(daughter->nlevels>0){
	    calc_decay_weights(mother,motherlevel,ibaby,
			       decayweight,ndecaylevels, &weight_norm);
	  }
	}
      }
      if(weight_norm>1.0E-50){
	for(ibaby=0;ibaby<NBABIES;ibaby++){
	  ndaughter=nmother-NBABY[ibaby];
	  zdaughter=zmother-ZBABY[ibaby];
	  if(ndaughter>=0 && zdaughter>=0 && ndaughter+zdaughter>0){
	    daughter=nucl_info[ndaughter][zdaughter];

	    if(daughter->nlevels>0){
	      add_decay_weights(mother,motherlevel,ibaby,decayweight,
				ndecaylevels,weight_norm);
	    }
	  }
	}
	mother->finalprob[motherlevel]=0.0;
      }
      else{
	mother->nstablelevels=motherlevel+1;
      }
      mother->finalsum+=mother->finalprob[motherlevel];
    }
  }
}

/* This function calculates decay weights into every possible level
   of a daughter nucleus determined by the properties of the mother
   and the baby. Weights are not added at the time as they cannot be
   normalized until they have been calculated for all possible babies */
void calc_decay_weights(nucl_info_class *mother,int motherlevel,int ibaby,
			double decayweight[NBABIES][MAXLEVELS],
			double *ndecaylevels,double *weight_norm){
  nucl_info_class *daughter,*baby;  
    // Mother decays into daughter and baby.
  double alpha,radius,V,rmass,Ek,sum,dw,dwprime;
  int j;
  int nmother,zmother,ndaughter,zdaughter,adaughter,ababy,amother;
  nmother=mother->n;
  zmother=mother->z;
  ndaughter=nmother-NBABY[ibaby];
  zdaughter=zmother-ZBABY[ibaby];
  daughter=nucl_info[ndaughter][zdaughter];
  baby=nucl_info[NBABY[ibaby]][ZBABY[ibaby]];
  ababy=NBABY[ibaby]+ZBABY[ibaby];
  adaughter=ndaughter+zdaughter;
  amother=nmother+zmother;
  alpha=1.44*double(ZBABY[ibaby])*double(zdaughter);
  radius=1.2*(pow(double(ababy),1.0/3.0)+pow(double(adaughter),1.0/3.0));
  V=alpha/radius;
  rmass=931.5*adaughter*ababy/double(ababy+adaughter);

  sum=0.0;
  ndecaylevels[ibaby]=0;
  j=0;
  if(daughter->filled==0){
    daughter->levelfiller();
  }
  Ek=mother->energy[motherlevel] - daughter->energy[j] - baby->energy[0];
  while(Ek > 1.0E-20 && j < daughter->nlevels){ 
    // While mother still has energy for every energy level of the daughter
    
    if(Ek<V){
            dw=daughter->degen[j]*PI*HBARC*HBARC/
      2*exp((2*sqrt(2*rmass)/HBARC)*(radius*sqrt(V-Ek)
      -(alpha/sqrt(Ek))*atan(sqrt((V-Ek)/Ek))));	      
    }
    else{
      dwprime=daughter->degen[j]*PI*HBARC*HBARC/2;
      dw=daughter->degen[j]*PI*radius*radius*(Ek-V)
	*rmass;
      if(dwprime>dw) {
	dw=dwprime;
      }
    }
    decayweight[ibaby][j]=dw;
    sum=sum+dw;
    j=j+1;
    ndecaylevels[ibaby]=j;
    if(j<daughter->nlevels)
      Ek=mother->energy[motherlevel]-daughter->energy[j]-baby->energy[0];
  }
  *weight_norm+=sum;
}

// This function takes an n and z value, an energy level, an integer 
// corresponding to the the product, and the weight normalization parameter.
// The decay_weights array is filled in calc_decay_weights.
// It modifies the weight of the mother, daughter, and baby.
void add_decay_weights(nucl_info_class *mother,int motherlevel,int ibaby,
		       double decayweight[NBABIES][MAXLEVELS],
		       double *ndecaylevels,double weight_norm){
  nucl_info_class *daughter,*baby;
  int j,nmother,zmother;
  double dw;

  nmother=mother->n;
  zmother=mother->z;
  daughter=nucl_info[nmother-NBABY[ibaby]][zmother-ZBABY[ibaby]];
  baby=nucl_info[NBABY[ibaby]][ZBABY[ibaby]];

  for(j=0;j<ndecaylevels[ibaby];j++){
    dw=decayweight[ibaby][j]*mother->finalprob[motherlevel]/weight_norm;
    daughter->finalprob[j]+=dw;
    baby->finalprob[0]+=dw;
  }
}