Bethe-Weizsäcker - Radioactive Decay

bw-radio-decay.cpp

// Bethe-Weizsäcker Radioactive Decay

#include <iostream>
#include <cmath>

using namespace std;

const long double M_PROTON = 1.67262192369e-27; // kg
const long double M_NEUTRON = 1.67492749804e-27; // kg
const long double M_ELECTRON = 9.1093837015e-31; // kg
const long double CONST_U = 1.66053906660e-27; // kg
const long double CONST_E = 1.602176634e-19; // J
const long double CONST_C = 299792458; // m/s

// Spektrum.DE Values + University
/*
const double av = 15.85; // MeV
const double as = 18.34; // MeV
const double ac = 0.71; // MeV
const double aa = 92.86/4.; // MeV
const double ap = 11.46; // MeV
*/

// Springer Values, A>50, https://link.springer.com/article/10.1007/s41365-019-0718-8
const double av = 14.6433;
const double as = 14.0788;
const double ac = 0.6442;
const double aa = 21.0680;
const double ap = 11.5398;

// Returns Nucleus Mass In kg
long double get_mass(unsigned int A, unsigned int Z){
	if(A == 4 && Z == 2){
		return 6.644657230e-27; // Alpha Particle Mass, Computation Inaccurate
	}

	long double energy = av * A - as * pow(A, 2./3) - ac * Z*Z / pow(A, 1./3) - aa * pow(A-2*Z, 2) / A;

        if(Z%2 == 0 && A%2 == 0){ // gg-Kern
        	energy += ap / sqrt(A);
      	}
        if(Z%2 != 0 && A%2 == 0){ // uu-Kern
                energy -= ap / sqrt(A);
        }
        energy *= 1e6 * CONST_E; // Joule

	return (A-Z)*M_NEUTRON + Z*M_PROTON - energy/(CONST_C * CONST_C);
}

// True if Beta+ Decay Possible
bool beta_plus(unsigned int A, unsigned int Z){
	return get_mass(A,Z) > get_mass(A,Z-1) + 2 * M_ELECTRON;
}

// True if Beta- Decay Possible
bool beta_minus(unsigned int A, unsigned int Z){
	return get_mass(A,Z) > get_mass(A,Z+1);
}

// True if Alpha Decay Possible
bool alpha(unsigned int A, unsigned int Z){
	return get_mass(A,Z) > get_mass(4,2) + get_mass(A-4,Z-2);
}

int main(){
	unsigned int A,Z;
	cout << "A: ";
	cin >> A;
	cout << "Z: ";
	cin >> Z;

	cout << "Mass (Nucleus): " << get_mass(A,Z)/CONST_U << " u" << endl;
	cout << "Alpha Decay: " << (get_mass(A,Z) -get_mass(A-4,Z-2) -get_mass(4,2))*pow(CONST_C,2)/CONST_E/1e6 << " MeV" << endl;
	cout << "EC: " << (get_mass(A,Z) -get_mass(A,Z-1) +M_ELECTRON)*pow(CONST_C,2)/CONST_E/1e6 << " MeV" << endl;
	cout << "Beta+ Decay: " << (get_mass(A,Z) -get_mass(A,Z-1) -M_ELECTRON)*pow(CONST_C,2)/CONST_E/1e6 << " MeV" << endl;
	cout << "Beta- Decay: " << (get_mass(A,Z) -get_mass(A,Z+1) -M_ELECTRON)*pow(CONST_C,2)/CONST_E/1e6 << " MeV" << endl;
	//cout << "Alpha Decay Possible: " << alpha(A,Z) << endl;
	//cout << "Beta+ Decay Possible: " << beta_plus(A,Z) << endl;
	//cout << "Beta- Decay Possible: " << beta_minus(A,Z) << endl;
	return 0;
}