```c
else
    s_star = s_hi;
end
[Xi_d1, Xi_d2] = scaled(Xi_d1, Xi_d2, s_star);
[~, Efiss_MeV] = sum_and_Efiss(Xi_d1, Xi_d2);
ok = isfinite(Efiss_MeV) && Efiss_MeV > 0;
fprintf(' [FISSION] Final Efiss = %.3f MeV (target %.1f) | (unbracketed
fallback)\n', ...;
    Efiss_MeV, Efiss_target_MeV);
return;
end

% Bisection to tolerance
tol = 0.05; % MeV
for it = 1:50
    sm = sqrt(s_lo * s_hi); % geometric mid to keep positivity well-
    behaved
        [X1m, X2m] = scaled(Xi_d1, Xi_d2, sm);
        [~, Ef_m] = sum_and_Efiss(X1m, X2m);
    if abs(Ef_m - Efiss_target_MeV) <= tol
        Xi_d1 = X1m; Xi_d2 = X2m; Efiss_MeV = Ef_m;
        ok = true;
        fprintf(' [FISSION] Final Efiss = %.3f MeV (target %.1f) | EO=%.3f, (E1+
E2)=%.3f\n', ...;
            Efiss_MeV, Efiss_target_MeV, EO_MeV, EO_MeV - Ef_m);
        return;
    end
    if Ef_m > Efiss_target_MeV
        % need more daughter energy -> increase s: move lo up
        s_lo = sm; Ef_lo = Ef_m;
    else
        % too much daughter energy -> decrease s: move hi down
        s_hi = sm; Ef_hi = Ef_m;
    end
end

% If here, didn't meet tol; return best midpoint
sm = sqrt(s_lo * s_hi);
[Xi_d1, Xi_d2] = scaled(Xi_d1, Xi_d2, sm);
[~, Efiss_MeV] = sum_and_Efiss(Xi_d1, Xi_d2);
ok = isfinite(Efiss_MeV) && Efiss_MeV > 0;
fprintf(' [FISSION] Final Efiss = %.3f MeV (target %.1f) | (max iters)
\n',
    Efiss_MeV, Efiss_target_MeV);
end

function [Xi_d1, Xi_d2] = try_build_scission(bname, G, C, Xiparent, EO_J, tf)
    try, [Xi_d1, Xi_d2] = feval(bname, G, C, Xiparent, EO_J, tf); return; catch, end
    try, [Xi_d1, Xi_d2] = feval(bname, G, C, Xiparent, EO_J);     return; catch, end
    try, [Xi_d1, Xi_d2] = feval(bname, G, C, EO_J, tf);           return; catch, end
    try, [Xi_d1, Xi_d2] = feval(bname, G, C, Xiparent);          return; catch, end
    [Xi_d1, Xi_d2] = feval(bname, G, C); % last attempt
end

function w = simpson_weights(N)
    w = ones(1,N);
```
```mat
if N>2, w(2:2:end-1)=4; w(3:2:end-2)=2; end
w([1 end]) = 1; w = (1/3)*w;
end

function assert_file(fname)
    if exist(fname, 'file') ~= 2
        error('FISSION:missingFile', 'Required file missing: %s', fname);
    end
end

function s = ternaryLTR技術a, b)
    if cond, s=a; else, s=b; end
end
```

<caption>Listing 28: Fission summary driver (U-235; calibrated scission and energy scale)</caption>

```
% FILE: TSRT_Calib_Fission_ScissWth_U235.m
% PURPOSE: Sweep/tune the smooth scission width so the Efiss range is well-behaved;
%          persist C.fission.scission_w_fm in tsrtconstants Calibrated.mat.
% OUTPUTS: tsrtconstants Calibrated.mat with C.fission.scission_w_fm
function TSRT_Calib_Fission_ScissWth_U235()
% TSRT_Calib_Fission_ScissWth_U235
% Robust sweep-based calibration of C.fission.scission_w_fm for U-235(n,f).
% Compatible with scission builders that require any of:
% (G,C, Xiparent, E0_J, target_fraction)
% (G,C, Xiparent, E0_J)
% (G,C, E0_J, target_fraction)
% (G,C, E0_J)
% (G,C, Xiparent)
% (G,C)
%
% Writes: tsrtconstants Calibrated.mat with C.fission.scission_w_fm

fprintf(' [CAL:SCISION] U-235 scission-width calibration (sweep based) 
');

% ---- constants & merge calibrated MAT ----
C = tsrtconstants();
if exist('tsrtconstants Calibrated.mat', 'file') == 2
    S = load('tsrtconstants Calibrated.mat');
    if isstruct(S)
        f = fieldnames(S);
        for i=1:numel(f), C.(f{i}) = S.(f{i}); end
        disp(' [CAL:SCISION] Merged tsrtconstants Calibrated.mat into C.');
    end
end

% Ensure global normalization exists
if ~isfield(C,'norm_energy') || ~isscalar(C.norm_energy) || ~isfinite(C.
    norm_energy) || C.norm_energy<=0
    error(' [CAL:SCISION] C.norm_energy missing/invalid. Run
        TSRT_Calib_NormEnDeuteron first.');
end

% Provide reasonable defaults for fission sub-struct
if ~isfield(C,'fission') || ~isstruct(C.fission), C.fission = struct(); end
if ~isfield(C.fission,'target_fraction') || ~isfinite(C.fission.target_fraction)
    C.fission.target_fraction = 0.98; % used by tuned builders
end
```
end

% ---- Local Simpson spherical grid (with sin theta) folded into theta weights)

---

if ~isfield(C,'N_r') || isempty(C.N_r), C.N_r = 96; end
if ~isfield(C,'N_theta') || isempty(C.N_theta), C.N_theta = 64; end
if ~isfield(C,'N phi') || isempty(C.N phi), C.N phi = 64; end
if ~isfield(C,'r_max') || isempty(C.r_max), C.r_max = 12e-15; end

Nr = max(3, C.N_r);
Nt = max(3, C.N_theta);
Np = max(3, C.N phi);

r = linspace(0, C.r_max, Nr);
theta = linspace(0, pi, Nt);
phi = linspace(0, 2.0*pi, Np);

wr = simpson weights (Nr) * (C.r_max/(Nr-1));
wt = simpson weights (Nt) * (pi/(Nt-1));
wp = simpson weights (Np) * ((2*pi)/(Np-1));
wt = wt .* sin theta;

G = struct();

G.r = r; G theta = theta; G = phi;
G.w_r = wr; G.w_theta = wt; G.w phi = wp;

% legacy aliases:

G.wr = wr; G.wrad = wr;
G.w_th = wt; G.wtheta = wt; G.wth = wt;
G.w_ph = wp; G.wphi = wp; G.wph = wp;

disp(' [CAL:SCISION] Using local Simpson-weighted grid. ');

% ---- Parent once and E0 ----

assert_file('tsrt_buildVO235.m');
XiVO parent = tsrt_buildVO235(G, C);
KO = tsrtcurvaturetensor(G, XiVO parent, C);
EO_J = tsrtenergyfromcurvature(G, KO, C);
EO_MeV = EO_J / 1.602176634e-13;

% ---- Choose scission builder ----

use_tuned = (exist('tsrt_build_daughters_scission_tuned.m', 'file') == 2);
if use_tuned
    bname = 'tsrt_build_daughters_scission_tuned';
    fprintf(' [CAL:SCISION] Using tuned scission builder: %s
', bname);
else
    assert_file('tsrt_build_daughters_scission.m');
    bname = 'tsrt_build_daughters_scission';
    fprintf(' [CAL:SCISION] Using scission builder: %s
', bname);
end

% ---- Sweep widths (fm) ----

w FM grid = [0.3 0.5 0.7 1.0 1.3 1.6 2.0 2.5 3.0 3.5 4.0];
Ef_mev = nan(size(w FM grid));

for i = 1: numel(w FM grid)

w FM = w FM grid (i);
C.fission.scission w FM = w FM; % let builders read from C
% Build daughters with robust signature attempts

[Xi_d1, Xi_d2] = try_build_scission(bname, G, C, Xi_parent, EO_J);

% Energies and Efiss

K1 = tsrtcurvaturetensor(G, Xi_d1, C);
E1 = tsrtenergyfromcurvature(G, K1, C);
K2 = tsrtcurvaturetensor(G, Xi_d2, C);
E2 = tsrtenergyfromcurvature(G, K2, C);

Ef_J = EO_J - (E1 + E2);
Ef_mev(i) = Ef_J / 1.602176634e-13;

end

% ---- Print sweep table ----

fprintf('\n[CAL:SCISION] width (fm) Efiss (MeV) EO (MeV)\n');
for i = 1:numel(w_fm_grid)
    fprintf(' %8.3f %10.4f %10.4f\n', w_fm_grid(i), Ef_mev(i),
        EO_MeV);
end
fprintf('\n');

% ---- Choose width ----

target = 170;          % MeV
band   = [120, 240];    % MeV
cand = find(isfinite(Ef_mev) & Ef_mev >= band(1) & Ef_mev <= band(2));
if ~isempty(cand)
    [~,j] = min(abs(Ef_mev(cand) - target));
    w_best = w_fm_grid(cand(j));
    why = sprintf('closest to %g MeV in [%g,%g]', target, band(1), band(2));
else
    cand2 = find(isfinite(Ef_mev) & Ef_mev > 0);
    if ~isempty(cand2)
        [~,j] = max(Ef_mev(cand2));
        w_best = w_fm_grid(cand2(j));
        why = 'maximum positive Efiss';
    else
        [~,j] = max(Ef_mev);
        w_best = w_fm_grid(j);
        why = 'largest Efiss (no positive values found)';
    end
end

Ef_best = Ef_mev(w_fm_grid==w_best);
fprintf(' [CAL:SCISION] Chosen width = %.3f fm (%s) ; Efiss = %.4f MeV\n' , w_best,
why, Ef_best);

% ---- Persist calibrated width ----

C_out = struct();
if exist('tsrtconstantsCalibrated.mat','file') == 2
S2 = load('tsrtconstantsCalibrated.mat');
if isstruct(S2)
    f2 = fieldnames(S2);
    for i=1:numel(f2), C_out.(f2{i}) = S2.(f2{i}); end
end
end
C_out.fission.scission_w_fm = w_best;
save('tsrtconstantsCalibrated.mat','-struct','C_out');
fprintf(' [CAL:SCISION] Saved C.fission.scission_w_fm = %.3f into
```matlab
tsrtconstantsCalibrated.mat\n', w_best);

fprintf('\Paste this into tsrtconstants.m for it to be hard-coded:\n');
fprintf(' C.fission.scission_w_fm = %.3f; %% calibrated scission smoothness (fm)\n', w_best);
fprintf('\nDone.\n');
end

% ---- helpers ----
function [Xi_d1, Xi_d2] = try_build_scission(bname, G, C, Xiparent, EO_J)
tf = C.fission.target; % Most informative -> least:
try
    [Xi_d1, Xi_d2] = feval(bname, G, C, Xiparent, EO_J, tf); % (G, C, Xiparent, EO_J, target_fra
    return;
catch, end
try
    [Xi_d1, Xi_d2] = feval(bname, G, C, Xiparent, EO_J); % (G, C, Xiparent, EO_J)
    return;
catch, end
try
    [Xi_d1, Xi_d2] = feval(bname, G, C, EO_J, tf); % (G, C, EO_J, target_fra
    return;
catch, end
try
    [Xi_d1, Xi_d2] = feval(bname, G, C, EO_J); % (G, C, EO_J)
    return;
catch, end
try
    [Xi_d1, Xi_d2] = feval(bname, G, C, Xiparent); % (G, C, Xiparent)
    return;
catch, end
% last attempt
[Xi_d1, Xi_d2] = feval(bname, G, C); % (G, C)
end

function assert_file(fname)
    if exist(fname, 'file') ~= 2
        error('[CAL:SCISS] Required file missing: %s', fname);
    end
end

function w = simpson_weights(N)
w = ones(1,N);
if N > 2
    w(2:2:end-1) = 4;
    w(3:2:end-2) = 2;
end
w([1 end]) = 1;
w = (1/3) * w;
end
```

Listing 29: Calibrate scission smoothness for $^{235}$U (sets C.fission.scission_w_fm)

tsrt_build_parent_U235.m
```matlab
% FILE: tsrt_buildVO235.m
% PURPOSE: Construct the compact U-235 parent curvature/field proxy on grid G.
function XiVO235 = tsrt_buildVO235(G, C)
%TSRT_BUILD_PARENT_U235 Composite trembling field for U-235 (Z=92, A=235).
% Strategy:
% - Generate Z proton and N neutron centers inside a nuclear radius
% R = r0 * A^(1/3) (fm) with a soft-surface (Woods-Saxon-like) radial bias.
% - Assemble via tsrtnucleonfieldshifted using modest per-nucleon amplitudes.
% % Dependencies: tsrtnucleonfieldshifted, tsrtnucleonfield, tsrtcompositefield helpers.
Z = 92; A = 235; N = A - Z;

% Geometry scale (meters)
r0_fm = getfield_def(C, 'r0_fm', 1.20);
R_fm = r0_fm * A^(1/3);
R_m = R_fm * 1e-15;

% Nucleon radii / amplitudes (tweakable but fixed across code)
r0_p = 0.85e-15; % m
r0_n = 0.90e-15; % m
Aamp = 1.0;      % unit amplitude per nucleon
anis = 0.0;      % isotropic for baseline

% Draw positions with soft edge (Woods-Saxon radial PDF)
Rp = sample positions woods saxon(Z, R_m);
Rn = sample positions woods saxon(N, R_m);

% Pack into the 'nuc' struct expected by tsrtcompositefield loop
types = [repmat {'proton'}, 1, Z), repmat {'neutron'}, 1, N) ];
R = [Rp, Rn]; % 3xA

% Optional per-nucleon parameters (scalar apply to all)
nuc = struct();
nuc.types = types;
nuc.R = R;
nuc.r0 = [r0_p * ones(1,Z), r0_n * ones(1,N)];
nuc.Aamp = Aamp;
nuc.anis = anis;
nuc.phi = 0.0;

% Assemble composite field (time argument unused; pass 0)
XiVO235 = tsrtcompositefield(G, nuc, 0);

% ---- helpers ----
function val = getfield_def(S, k, d)
    if isfield(S,k) && ~isempty(S.(k)), val = S.(k); else, val = d; end
end
end

function R = sample positions woods saxon(M, R_m)
% Draw M positions with an approximate Woods-Saxon radial bias: p(r) ~ r^2 / (1+exp((r-R)/a))
% using simple rejection; place isotropically on the sphere.

a = 0.6e-15;      % surface thickness (m)
Rmax = 1.2 * R_m;  % loose outer bound
R = zeros(3, M);
```
```c
i = 1;
while i <= M
    r = Rmax * rand;      % trial radius
    x = 2*rand - 1;      % for acceptance
    pdf = (r^2) / (1 + exp((r - R_m)/a));
    pdf_max = (R_m^2);   % safe crude bound near interior
    if x^2 <= pdf / pdf_max % accept (cheap acceptance trick)
        % isotropic direction
        u = 2*rand-1;      phi = 2*pi*rand;
        s = sqrt(1-u^2);
        R(:,i) = r * [s*cos phi; s*sin phi; u];
        i = i + 1;
    end
end
end
```

<caption>Listing 30: U-235 parent field builder</caption>

```c
% FILE: tsrt_build_daughters_scission_tuned.m
% PURPOSE: Build daughters from the parent using a tuned scission logic; optionally
%          accepts (G, C, Xi_parent, E0_J, target fraction) signatures.
function [Xi_d1, Xi_d2] = tsrt_build_daughters_scission_tuned(G, C, Xiparent, E0_J,
    target_fraction)
%TSRT_BUILD_DAUGHTERS_SCISSION_TUNED Asymmetric scission daughters with energy cap.
% [Xi_d1, Xi_d2] = tsrt_build_daughters_scission_tuned(G, C, Xiparent, E0_J,
    target_fraction)
%
% Steps:
% 1) Build two daughter composites at ±d/2 on z, with (Z1,A1)~(37,95), (Z2,A2)~
    (55,140).
% 2) Compute E1+E2; rescale BOTH daughters by s = sqrt( (tf*E0)/(E1+E2) ).
% This enforces E1+E2 = tf * E0 (<= E0), guaranteeing Efiss >= 0.

if nargin < 5 || isempty(target_fraction), target_fraction = getf(C,'fission','
    target_fraction',0.90); end
if target_fraction > 1, target_fraction = 1; end
if target_fraction <= 0, target_fraction = 0.85; end

% Daughter splits (typical U-235(n,f))
Z1=37; A1=95; N1 = A1 - Z1;
Z2=55; A2=140; N2 = A2 - Z2;

% Scission separation (meters)
R1_fm = getf(C,'r0_fm',1.20) * A1^(1/3);
R2_fm = getf(C,'r0_fm',1.20) * A2^(1/3);
gap_fm = getf(C,'fission','gap_fm',1.0); % small gap at scission
    d FM = R1 FM + R2 FM + gap FM;    % center-to-center
    d M = d FM * 1e-15;

% Nucleon radii / amplitude
r0_p = 0.85e-15; r0_n = 0.90e-15; Aamp = 1.0; anis = 0.0;

% Daughter 1 positions around z = -d/2
R1 = sample positions woods saxon (A1, (getf(C,'r0_fm',1.20)*A1^(1/3))*1e-15);
R1(3,:) = R1(3atz) - d_m/2;
types1 = [ repmat {'proton'}, 1, Z1), repmat {'neutron'}, 1, N1) ];
nuc1 = struct('types', {types1}, 'R', R1, ...)
```
'r0', [r0_p*ones(1, Z1), r0_n*ones(1, N1)], ...
'Aamp', Aamp, 'anis', anis, 'phi', 0.0);
Xi_d1 = tsrtcompositefield(G, nuc1, 0);

% Daughter 2 positions around z = +d/2
R2 = sample positions woods saxon (A2, (getf(C,'r0_fm',1.20)*A2^(1/3))*1e-15);
R2(3, :) = R2(3, :) + d_m/2;
types2 = [ repmat {'proton'}, 1, Z2), repmat {'neutron'}, 1, N2) ];
nuc2 = struct('types', {types2}. 'R', R2, ...
'r0', [r0_p*ones(1, Z2), r0_n*ones(1, N2)], ...
'Aamp', Aamp, 'anis', anis, 'phi', 0.0);
Xi_d2 = tsrtcompositefield(G, nuc2, 0);

% Energies
K1 = tsrtcurvaturetensor(G, Xi_d1, C); E1 = tsrtenergyfromcurvature(G, K1, C);
K2 = tsrtcurvaturetensor(G, Xi_d2, C); E2 = tsrtenergyfromcurvature(G, K2, C);
Esum = E1 + E2;

% Rescale daughters to enforce E1+E2 = target_fraction * E0
if Esum > 0 && isfinite(E0_J)
    tf = min(max(target_fraction, 0), 1);
    s = sqrt(tf * E0_J) / Esum;
    Xi_d1.profile = s * Xi_d1.profile;
    Xi_d2.profile = s * Xi_d2.profile;
end
end

% --- small helpers ---
function val = getf(C, k1, k2, d)
    if nargin==3, d=k2; k2=''; end
    if isempty(k2)
        if isfield(C,k1)&&~isempty(C.(k1)), val=C.(k1); else, val=d; end
    else
        if isfield(C,k1)&&isstruct(C.(k1))&&isfield(C.(k1),k2)&&~isempty(C.(k1).(k2))
            val=C.(k1).(k2);
        else
            val=d;
        end
    end
end

function R = sample positions woods saxon (M, R_m)
    a = 0.6e-15; Rmax = 1.2 * R_m;
    R = zeros(3, M);
    i = 1;
    while i <= M
        r = Rmax * rand; u = 2*rand-1; phi = 2*pi*rand;
        pdf = (r^2) / (1 + exp((r - R_m)/a));
        pdf_max = (R_m^2);
        if (rand * pdf_max) <= pdf
            s = sqrt(1-u^2);
            R(:,i) = r * [s*cos phi; s*sin phi; u];
            i = i + 1;
        end
    end
end
end

<caption>Listing 31: Scission daughters (tuned builder)</caption>
tsrt_build_daughters_partition.m

% FILE: tsrt_build_daughters_partition.m
% PURPOSE: Build daughters by a smooth partition-of-unity split along z with width w_m

function [Xi_d1, Xi_d2, tag] = tsrt_build_daughters_partition(G, Xiparent, z0, w)
%TSRT_BUILD_DAUGHTERS_PARTITION
% Smoothly split the parent field into two daughters along z with a
% partition-of-unity w1(z) + w2(z) = 1 and small overlap. This removes
% cross terms in a quadratic energy, giving E1+E2 <= E0 (=> Efiss >= 0).

% Inputs:
% G : grid struct with fields r,theta,phi (and weights)
% Xiparent : struct with .profile (Nr x Nt x Np)
% z0 : split center (m), typically 0
% w : transition width (m), e.g. 1 fm = 1e-15

% Outputs:
% Xi_d1, Xi_d2 : daughter field structs with .profile
% tag : short descriptor for logging

if nargin < 3 || isempty(z0), z0 = 0.0; end
if nargin < 4 || isempty(w), w = 1.0e-15; end

% Coordinates
[RR,TT,PP] = ndgrid(G.r, G theta, G.phi);
[x,y,z] = sph2cart(PP, pi/2 - TT, RR);

% Smooth logistic windows along z
% w1 ~ 1 for z > z0, ~0 for z < z0; w2 = 1 - w1
s = (z - z0) ./ max(w, eps);
w1 = 1 ./ (1 + exp(-6*s)); % sharper transition with factor 6
w2 = 1 - w1;

% Optional: reduce overlap further by softening the mid-plane

mid = abs(z - z0) <= 0.5*w;
if any(mid(:))
    % taper both slightly at the midplane to suppress double counting
    w1(mid) = 0.5 * w1(mid);
    w2(mid) = 0.5 * w2(mid);
    % Renormalize to keep w1+w2 ~ 1 outside midplane, <=1 on the plane
    % (This only helps ensure E1+E2 strictly <= E0 numerically.)
end

% Apply windows (pointwise product)
P = Xiparent.profile;
Xi_d1 = struct('profile', w1 .* P);
Xi_d2 = struct('profile', w2 .* P);

tag = sprintf('z0=%2g fm, w=%2g fm', z0/1e-15, w/1e-15);
end

<caption>Listing 32: Scission daughters via smooth partition (fallback)</caption>

tsrtfissionenergy.m

% FILE: tsrtfissionenergy.m
% PURPOSE: Compute fission-barrier heights and semiclassical action integrals
% for TSRT curvature-derived fields.
```mat
% INPUTS: tsrtconstants.m, tsrtcurvaturetensor.m, tsrtenergyfromcurvature.m
% OUTPUTS: Efiss (J)
function Efiss = tsrtfissionenergy(G, Xi_init, Xi_d1, Xi_d2, C)
K0 = tsrtcurvaturetensor(G, Xi_init, C);
E0 = tsrtenergyfromcurvature(G, K0, C);
K1 = tsrtcurvaturetensor(G, Xi_d1, C);
E1 = tsrtenergyfromcurvature(G, K1, C);
K2 = tsrtcurvaturetensor(G, Xi_d2, C);
E2 = tsrtenergyfromcurvature(G, K2, C);
Efiss = (E0 - (E1 + E2));
end
```

<caption>Listing 33: TSRT fission energetics core routine</caption>

```mat
% FILE: fissiontablemake.m
% PURPOSE: Format fission_results.mat into a LaTeX table (tab_fission_summary.tex)
%          for direct inclusion in the main article (Table\ref{tab:fission-summary-235U}).
% INPUTS: fission_results.mat (struct array F)
% OUTPUTS: paper-assets-tex/tab_fission_summary.tex (or ./tab_fission_summary.tex)
function fissiontablemake(out_txt)
% FISSIONTABLEMAKE Write a minimal fission summary table from fission_results.mat.
% Usage:
% fissiontablemake;                 % writes ./tab_fission_summary.tex
% fissiontablemake('paper-assets-tex/tab_fission_summary.tex');

if length < 1 || isempty(out_txt)
    out_txt = 'tab_fission_summary.tex';
end

% --- Load results (struct array F expected) ---
if exist('fission_results.mat', 'file') =~ 2
    error('fissiontablemake:missingFile', 'fission_results.mat not found.');
end
S = load('fission_results.mat');
if ~isfield(S, 'F') || ~isstruct(S.F)
    error('fissiontablemake:badShape', 'Expected struct array F in fission_results.mat');
end
F = S.F;

need = {'nucleus', 'B_TSRT_MeV', 'S_action', 'Efiss_MeV'};
for k = 1:numel(need)
    if ~isfield(F, need{k})
        error('fissiontablemake:missingField', 'Missing field "%s" in F.', need{k});
    end
end

% --- Ensure directory OR fallback to current folder if not writable ---
target = out_txt;
[outDir, base, ext] = fileparts(out_txt);
if isempty(ext), ext = '.tex'; end
if ~isempty(outDir)
    if exist(outDir, 'dir') =~ 7
        ok = mkdir(outDir);
        if ~ok
```
```
warning('fissiontablemake:mkdir', 'Cannot create directory: %s. Falling
back to current folder.', outDir);
        target = [base ext];
    end
end
end

% --- Open file (with fallback) ---
fid = fopen(target, 'w');
if fid < 0
    warning('fissiontablemake:fopen', 'Cannot open for write: %s. Falling back to
./tab_fission_summary.tex', target);
    target = 'tab_fission_summary.tex';
    fid = fopen(target, 'w');
    if fid < 0
        error('fissiontablemake:fopenFinal', 'Cannot open for write: %s', target);
    end
end

% --- Emit LaTeX tabular (no caption/label; paper wraps it) ---
fprintf fid, '%%% Auto-generated by fissiontablemake.m\n');
fprintf fid, '\\begin{tabular}{lccc}\n');
fprintf fid, '\\hline\n');
fprintf fid, 'System & $B_{\\rm TSRT}$ (MeV) & $\mathcal{S}$ (arb.) & $E_{\\rm rm}$
fiss)$ (MeV)\\\\\n');
fprintf fid, '\\hline\n');

for i = 1:numel(F)
    nuc = aschar(F(i).nucleus);

    % allow char or numeric-like strings for these:
    BMeV = aschar(field_or_dash(F(i),'B_TSRT_MeV'));
    Sact = aschar(field_or_dash(F(i),'S_action'));
    EfMeV= aschar(field_or_dash(F(i),'Efiss_MeV'));

    fprintf fid, '%s & %s & %s & %s\\\\\n', nuc, BMeV, Sact, EfMeV);
end

 fprintf fid, '\\hline\n');
 fprintf fid, '\\end{tabular}\n');
 fclose fid;

 fprintf('[FISSION] Table written: %s\n', target);
 fprintf('[FISSION] In LaTeX, use: \\input{%s}\n', target);
end

function s = aschar(x)
    if ischar(x) || isstring(x), s = char(x); else, s = '<unnamed>'; end
end

function s = asnum(x)
    if isempty(x) || ~isfinite(x), s = '---'; else, s = sprintf('%-.6g', x); end
end

function s = field_or_dash(S, key)
    if ~isfield(S, key) || isempty(S.(key))
        s = '---'; return;
    end
```
```c
v = S.(key);
if ischar(v) || isstring(v)
    s = char(v);
else if isnumeric(v) && isscalar(v) && isfinite(v)
    s = sprintf('%.*g', v);
else
    s = '---';
end
end
```

<caption>Listing 34: Fission summary-table generator (feeds Table 42 (p. 257))</caption>

*Usage.* This script generates `tab_fission_summary.tex`, which is included verbatim in the main article as Table 42 (p. 257). It requires no manual editing and ensures numerical identity between the fission energetics discussed in Section 6.2 and those tabulated in Appendix M.
## Y.16 Fission Path Correlations (optional)

Dependencies: Listings 33 (p. 270), 34 (p. 271).

Purpose: Quantify and visualize correlations between computed fission-barrier indicators and equilibrium binding energies.

This module is diagnostic and optional: it does not feed into the main pipeline but confirms geometric consistency across the fission-energy landscape.

Used in the paper: Figures in the supplementary diagnostics (Appendix Y.10) and cross-validation plots comparing curvature-derived and empirical barrier trends.

Inputs: `fission_results.mat` from Listing 33 (p. 270), and `fission_tables.xlsx` from Listing 34 (p. 271).

Outputs: Correlation plots `fission_path_corr.pdf` and numerical correlation coefficients saved in `fission_corr.mat`.

Dependencies: Requires prior execution of the fission energetics (Listing 33 (p. 270)) and fission table assembler (Listing 34 (p. 271)). No additional toolboxes or external data are needed.

`fission_path.m`

```matlab
% FILE: fission_path.m
% PURPOSE: Plot correlation between computed fission-barrier indicators and
% equilibrium binding energies or other observables.
% INPUTS: fission_results.mat, fission_tables.xlsx
% OUTPUTS: fission_path_corr.pdf, fission_corr.mat
% NOTES: Optional diagnostic. Run after fissiontablemake.m to validate trends.
function fission_path()
% FISSION_PATH Correlate TSRT fission outputs and make a small diagnostic plot.
% Inputs (optional/if exist):
% - fission_results.mat (struct array F with fields: nucleus, B_TSRT_MeV, S_action, Efiss_MeV, note)
% - fission_tables.xlsx (optional tabular, any columns; used if present)
% - exp_bindings.csv, tsrt_bindingscan_cln.csv, tsrt_bindingscan_q.csv (optional, if present)
% Outputs:
% - fission_corr.mat (struct with .data and .corr Pearson correlations)
% - fission_path_corr.pdf (scatter or histogram diagnostic)
% Notes:
% - Safe to run anytime after make_fission_summary_now and/or fissiontablemake.
% - All optional files are consumed only if present; otherwise ignored.
% - No toolboxes required.

clc;

%% ---------- Load core results ----------
if exist('fission_results.mat', 'file') ~= 2
    error('fission_path:missingCore', 'fission_results.mat not found. Run
           make_fission_summary_now first.');
end
S = load('fission_results.mat');
if ~isfield(S, 'F') || ~isstruct(S.F)
    error('fission_path:badShape', 'Expected struct array F in fission_results.mat');
end
F = S.F();

% Coerce to numeric where possible
Efiss = to_num_field(F, 'Efiss_MeV');
Btsrt = to_num_field(F, 'B_TSRT_MeV');
Saction = to_num_field(F, 'S_action');
```
% Build core table
core = table(...
    string(field_or_default(F, 'nucleus', '<unnamed>')),
    ...,
    Efiss, Btsrt, Saction, ...
    'VariableNames', {'System', 'Efiss_MeV', 'B_TSRT_MeV', 'S_action'});
%% ---------- Optionally load fission_tables.xlsx ----------
Txls = table();
if exist('fission_tables.xlsx', 'file') == 2
    try
        Txls = readtable('fission_tables.xlsx', 'FileType', 'spreadsheet');
    catch ME
        warning('fission_path:xlsx', 'Could not read fission_tables.xlsx (%s).
Continuing without it.', ME.message);
        Txls = table();
    end
end

%% ---------- Optionally load additional CSV diagnostics ----------
Texp = try_read_csv('exp bindings.csv');
Tscan = try_read_csv('tsrt bindingscan cln.csv');
if isempty(Tscan); Tscan = try_read_csv('tsrt bindingscan_q.csv'); end

%% ---------- Merge what we can (lightweight, robust) ----------
% We don't assume a common key beyond 'System' as a string.
% If xlsx/CSV use a different key (e.g., Z,A), we keep them separate.

data = struct();
data.core = core;
if ~isempty(Txls), data.xlsx = Txls; end
if ~isempty(Texp), data.exp = Texp; end
if ~isempty(Tscan), data.scan = Tscan; end

%% ---------- Compute simple correlations (available columns only) ----------
corrS = struct();
corrS.available = {};

% Helper to compute Pearson if both vectors are finite and length>2
function c = pcorr(a,b)
    m = isfinite(a) & isfinite(b);
    if nnz(m) > 2
        c = corr(a(m),b(m),'type','Pearson');
    else
        c = NaN;
    end
end

% Core correlations (when present)
corrS.Efiss_vs_B = pcorr(core.Efiss_MeV, core.B_TSRT_MeV);
corrS.Efiss_vs_S = pcorr(core.Efiss_MeV, core.S_action);
corrS.B_vs_S = pcorr(core.B_TSRT_MeV, core.S_action);

% List what we computed
corrS available = {'Efiss vs B', 'Efiss vs S', 'B vs S'};

%% ---------- Plot: If B_TSRT_MeV exists, scatter Efiss vs B; else histogram Efiss ----------
pdf_out = 'fission_path_corr.pdf';
try
    if any(isfinite(core.B_TSRT_MeV))
        figure('Color', 'w');
        scatter(core.B_TSRT_MeV, core.Efiss_MeV, 36, 'filled'); grid on;
       xlabel('B_{TSRT} (MeV)');
       ylabel('E_{fiss} (MeV)');
        title sprintf('Efiss vs Barrier (Pearson=%3f)', corrs.Efiss_vs_B));
        set(gcf, 'PaperPositionMode', 'auto');
        print '-dpdf', '-r300', pdf_out);
    else %
        histogram(core.Efiss_MeV, max(10, min(50, min(50, numel(core.Efiss_MeV)))); grid on;
       xlabel('E_{fiss} (MeV)'); ylabel('Count');
        title('Distribution of E_{fiss}');
        set(gcf, 'PaperPositionMode', 'auto');
        print '-dpdf', '-r300', pdf_out);
    end
    close(gcf);
catch ME
    warning 'fission_path:plot', 'Could not write %s (%s). Skipping plot.', pdf_out, ME
    .message);
    pdf_out = '';
end

%% Save correlation bundle
save('fission_corr.mat', 'data', 'corrs', 'pdf_out');

fprintf(' [PATH] Wrote fission_corr.mat
');
if ~isempty(pdf_out)
    fprintf(' [PATH] Wrote %s
', pdf_out);
end

end % fission_path

%% helpers
function v = to_num_field(F, key)
    v = nan(numel(F),1);
    for i = 1:numel(F)
        if isfield(F(i), key)
            v(i) = asnum(F(i).(key));
        end
    end
end

function T = try_read_csv(fname)
    T = [];
    if exist(fname, 'file') == 2
        try
            T = readtable(fname);
        catch ME
            warning('fission_path:csv', 'Could not read %s (%s). Ignoring.', fname, ME.
                message);
            T = [];
        end
    end
end

function s = field_or_default(S, key, def)
```matlab
s = repmat(string(def), numel(S), 1);
for i = 1:numel(S)
    if isfield(S(i), key)
        val = S(i).(key);
        if isstring(val) || ischar(val), s(i) = string(val); end
    end
end

function x = asnum(val)
    if isempty(val)
        x = NaN; return;
    end
    if ischar(val) || isstring(val)
        str = strtrim(char(val));
        ifstrcmp(str,'---')
            x = NaN;
        else
            y = str2double(str);
            if isfinite(y), x = y; else, x = NaN; end
        end
    elseif isnumeric(val) && isscalar(val)
        if isfinite(val), x = double(val); else, x = NaN; end
    else
        x = NaN;
    end
end
```

<caption>Listing 35: Fission-path correlation plot generator</caption>

*Usage.* Run after Listing 28 (p. 257). This produces `fission_corr.mat` and (if possible) `fission_path_corr.pdf`. Optional CSV/XLSX inputs are consumed only if present.

*Figures produced correspond to a supporting diagnostic (not shown in the main text).*
## Y.17 Parameter Scan Utilities

Dependencies: Listings 2 (p. 213).

**Purpose:** Provide a general orchestration framework for automated sensitivity studies, parameter sweeps, and large-scale diagnostic runs. Used for cross-checking robustness of fits and exploring parameter influence on curvature-derived observables.

**Used in the paper:** All parameter-sensitivity tests discussed in Appendix Y.11 and in the convergence/robustness scans (Appendix Y.3).

**Inputs:** `tsrtconstants.mat` (global constants) and the relevant TSRT module scripts (binding, fission, lifetime, etc.).

**Outputs:** Parameter-response maps and tabulated scan results stored under `results/scan/`.

**Dependencies:** All core routines (Listing 33 (p. 270), 16 (p. 238), 8 (p. 221), etc.). No special MATLAB toolboxes required.

### `scan_parameters.m`

```matlab
% FILE: scan_parameters.m
% PURPOSE: Execute batch parameter scans for sensitivity and robustness tests.
% INPUTS: tsrtconstants.mat (optional), TSRT module builders & energy routines
% OUTPUTS: results/scan/*.mat, results/scan/*.csv, optional diagnostic plot
%
% Usage examples:
% % Scan scission width (fm) for Efiss(U-235)
% scan_parameters('param', 'fission.scission_w_fm', ...
% 'values', linspace(0.2,3.0,12), ...
% 'task', 'Efiss_U235');
%
% % Same, but also compute a single-shot rescaled Efiss to a target (diagnostic)
% scan_parameters('param', 'fission.scission_w_fm', ...
% 'values', [0.3 0.5 1.0 2.0], ...
% 'task', 'Efiss_U235', ...
% 'force_target', true, 'Efiss_target_MeV', 170.0);
```

```matlab
function scan_parameters(varargin)
% ---------- Parse inputs ----------
p = inputParser;
addParameter(p, 'param', 'fission.scission_w_fm',@(s)ischar(s)||isstring(s));
addParameter(p, 'values', linspace(0.3,2.0,10),@isnumeric);
addParameter(p, 'task', 'Efiss_U235',@(s)ischar(s)||isstring(s));
addParameter(p, 'force_target',false,@islogical);
addParameter(p, 'Efiss_target_MeV',170.0,@(x)isnumeric(x)&@isscalar(x));
addParameter(p, 'note', '',@(s)ischar(s)||isstring(s));
parse(p, varargin{:});
P = p.Result;
param_path = char(P.param);
values = P.values(:).'; % row vector
```

```matlab
% ---------- Ensure path & output folder ----------
if isempty(strfind(path, pwd)), addpath(genpathascwd); end
outdir = fullfile('results', 'scan');
if exist(outdir, 'dir')~=7, mkdir(outdir); end
```

```matlab
% ---------- Load constants & merge calibrated ----------
C = tsrtconstants();
if exist('tsrtconstants Calibrated.mat', 'file')==2
    S = load('tsrtconstants Calibrated.mat');
    C = deepmerge(C,S);
end
```
```c
must_have(C, 'norm_energy', 'C.norm_energy missing/invalid; run
TSRT_Calib_NormEn_Deuteron.');
if ~isfield(C, 'fission') || ~isstruct(C.fission), C.fission = struct(); end
if ~isfield(C.fission, 'scission_w_fm') || ~isfinite(C.fission.scission_w_fm)
    C.fission.scission_w_fm = 1.0; % fallback
end
if ~isfield(C.fission, 'target fraction') || ~isfinite(C.fission.target fraction)
    C.fission.target fraction = 0.98;
end
```

```c
% ---------- Build compact spherical grid (no toolboxes) ----------
if ~isfield(C, 'N_r') || isempty(C.N_r), C.N_r=96; end
if ~isfield(C, 'N_theta') || isempty(C.N_theta), C.N_theta=64; end
if ~isfield(C, 'N') || isempty(C.N Phi), C.N Phi=64; end
if ~isfield(C, 'r_max') || isempty(C.r_max), C.r_max=12e-15; end
[G, grid note] = local_simpson_grid(C);
fprintf(' [SCAN] Using %s\n', grid note);
```

```c
% ---------- Prepare task ----------
switch upper(char(P task))
    case 'EFISS_U235'
        taskfun =@(C) task_Efiss_U235(G,C,P.force_target,P.Efiss_target_MeV);
        y_name = 'Efiss_MeV';
    otherwise
        error('Unknown task: %s', P task);
end
```

```c
% ---------- Scan loop ----------
N = numel(values);
Y_raw = NaN(1,N);
Y_forced= NaN(1,N); % filled if force_target=true
ok_vec = false(1,N);
```

```c
CO = C;
for i = 1:N
    Ci = CO;
    Ci = setfield_deep(Ci, param_path, values(i));
    [y_raw, y_forced, ok, meta] = taskfun(Ci);
    Y_raw(i) = y_raw;
    Y_forced(i) = y_forced;
    ok_vec(i) = ok;
    fprintf(' [SCAN] %s = %.6g -> %s = %.6g\s\n', ...
        param_path, values(i), y_name, y_raw, tern(P.force_target, ...
        sprintf(' (forced: %.6g)', y_forced), ''));
end
```

```c
% ---------- Save MAT & CSV ----------
ts = datestr(now, 'yyyymmdd_HHMMSS');
base = sprintf('scan_%s_%s_%s', y_name, param_path_to_tag(param_path), ts);
matfile = fullfile(outdir, [base '.mat)];
csvfile = fullfile(outdir, [base '.csv]);
S = struct();
S.param = param_path;
S.values = values;
S.y_name = y_name;
S.Y_raw = Y_raw;
S.Y_forced = Y_forced;
S.ok = ok_vec;
```
```R
S task = P task;
S force_target = P force_target;
S Efiss_target_MeV = P Efiss_target_MeV;
S note = char P note;
save (matfile, '-struct', 'S');

fid = fopen (csvfile, 'w');
fprintf (fid, 'param, value, %s, forced_%s, ok
', y_name, y_name);
for i=1:N
    fprintf (fid, '%s, %12g, %12g, %12g, %d
', param_path, values (i), Y_raw (i),
           Y_forced (i), ok_vec (i));
end
fclose (fid);

% ---------- Minimal plot (optional; no toolboxes) ----------
try
    fig = figure('visible', 'off');
    plot (values, Y_raw, '-o'); xlabel (strrep (param_path, '.', '\_'));
    ylabel (y_name)
;
    title (sprintf ('Scan: %s vs %s', y_name, param_path), 'interpreter', 'none');
    pngfile = fullfile (outdir, [base '.png']);
    print (fig, pngfile, '-dpng', '-r150'); close (fig);
catch
    % plotting is optional
end

fprintf ('[SCAN] Wrote: 
 %s 
 %s 
', matfile, csvfile);
end

% ======== Task: Efiss for thermal U-235 (raw + optional single-shot target-rescale) ========
function [Efiss_raw_MeV, Efiss_forced_MeV, ok, meta] = task_Efiss_U235(G,C,
force_target, Efiss_target_MeV)
meta = struct();
% Parent
assert_file('tsrt_buildVOparent_U235.m');
XiVOparent = tsrt_buildVOparent_U235(G, C);
KO = tsrtcurvaturetensor(G, XiVOparent, C);
E0_J = tsrtenergyfromcurvature(G, KO, C);
E0_MeV= E0_J / 1.602176634e-13;

% Daughters (prefer tuned builder)
bname = '';
if exist('tsrt_build_daughters_scission_tuned.m', 'file')==2
    bname = 'tsrt_build_daughters_scission_tuned';
else if exist('tsrt_build_daughters_scission.m', 'file')==2
    bname = 'tsrt_build_daughters_scission';
end

% Try tuned builder (mild sweep on targetrovraction), else partition fallback
Xd1 = []; Xd2 = [];
tf_list = [C.fission.targetrovaction,
0.96,0.94,0.92,0.90,0.88,0.86,0.84,0.82,0.80];
if ~isempty (bname)
    for tf = unique (tf_list)
        [Xd1, Xd2] = try_build_scission (bname, G, C, XiVOparent, E0_J, tf);
        if ~isempty (Xd1), break; end
    end
```
end

if isempty(Xd1)
    z0 = 0.0; w_m = C.fission.scission_w_fm * 1e-15;
    [Xd1, Xd2] = tsrt_build_daughters_partition(G, Xiparent, z0, w_m);
end

% Raw Efiss
[Efiss_raw_MeV, ok_raw] = eval_Efiss(G, C, EO_J, EO_MeV, Xd1, Xd2, false,
Efiss_target_MeV);

% Optional forced single-shot to target
if force_target
    [Efiss_forced_MeV, ok_forced, Xd1, Xd2] = eval_Efiss(G, C, EO_J, EO_MeV, Xd1,
Xd2, true, Efiss_target_MeV);
else
    Efiss_forced_MeV = NaN; ok_forced = true;
end

ok = ok_raw && ok_forced;
meta.EO_MeV = EO_MeV;
end

% ====== Energy eval with optional one-shot target forcing ======
function [Efiss_MeV, ok, Xi_d1, Xi_d2] = eval_Efiss(G, C, EO_J, EO_MeV, Xi_d1, Xi_d2,
do_force, Efiss_target_MeV)
    K1 = tsrtcurvaturetensor(G, Xi_d1, C);
    E1_J = tsrtenergyfromcurvature(G, K1, C);
    K2 = tsrtcurvaturetensor(G, Xi_d2, C);
    E2_J = tsrtenergyfromcurvature(G, K2, C);

    Efiss_J = EO_J - (E1_J + E2_J);
    Efiss_MeV = Efiss_J / 1.602176634e-13;

    if do_force
        E_target_sum_MeV = max(EO_MeV - Efiss_target_MeV, 0.0);
        E_sum_MeV = (E1_J + E2_J) / 1.602176634e-13;
        if isfinite(E_sum_MeV) && E_sum_MeV > 0
            s = sqrt(E_target_sum_MeV / E_sum_MeV);
            if isfield(Xi_d1,'profile') && isfield(Xi_d2,'profile')
                Xi_d1.profile = s * Xi_d1.profile;
                Xi_d2.profile = s * Xi_d2.profile;
            % recompute
            K1 = tsrtcurvaturetensor(G, Xi_d1, C); E1_J =
            tsrtenergyfromcurvature(G, K1, C);
            K2 = tsrtcurvaturetensor(G, Xi_d2, C); E2_J =
            tsrtenergyfromcurvature(G, K2, C);
            Efiss_J = EO_J - (E1_J + E2_J);
            Efiss_MeV = Efiss_J / 1.602176634e-13;
        end
    end
end

ok = isfinite(Efiss_MeV);
end

% ====== Helpers (no toolboxes) ======
function [Xi_d1, Xi_d2] = try_build_scission(bname, G, C, Xiparent, EO_J, tf)
    try, [Xi_d1, Xi_d2] = feval(bname, G, C, Xiparent, EO_J, tf); return; catch, end
    try, [Xi_d1, Xi_d2] = feval(bname, G, C, Xiparent, EO_J); return; catch, end
```matlab
try, [Xi_d1, Xi_d2] = feval(bname, G, C, Xiparent);
return; catch, end

try, [Xi_d1, Xi_d2] = feval(bname, G, C);
return; catch, end

Xi_d1 = []; Xi_d2 = [];
end

function [G,note] = local_simpson_grid(C)
    Nr = max(3, C.N_r); Nt = max(3, C.N_theta); Np = max(3, C.N);
    r = linspace(0, C.r_max, Nr);
    th = linspace(0, pi, Nt);
    ph = linspace(0, 2*pi, Np);
    wr = simpson weights (Nr) * (C.r_max / (Nr - 1));
    wt = simpson weights (Nt) * (pi / (Nt - 1));
    wp = simpson weights (Np) * ((2*pi) / (Np - 1));
    wt = wt .* sin(th);
    G = struct('r', r, 'theta', th, 'phi', ph, ...
        'w_r', wr, 'w_theta', wt, 'w', wp, ...
        'wr', wr, 'wrad', wr, 'w_th', wt, 'wtheta', wt, 'wth', wt, 'w ph', wp, 'wphi', wp, '
        wph', wp);
    note = sprintf('Simpson spherical grid: Nr=%d, Nt=%d, Np=%d, r_max=%d.3g m', Nr, Nt,
        Np, C.r_max);
end

function w = simpson weights (N)
    w = ones(1, N);
    if N > 2, w(2:2:end-1) = 4; w(3:2:end-2) = 2; end
    w([1 end]) = 1; w = (1/3) * w;
end

function must_have(C, field, msg)
    ok = isfield(C, 'norm_energy') && isscalar(C.norm_energy) && isfinite(C.norm_energy)
        && C.norm_energy > 0;
    if ~ok, error('SCAN:missing', msg); end
end

function tag = param_path_to_tag(s)
    tag = regexp(s, '[^A-Za-z0-9]+', '_');
end

function C = deepmerge(C, S)
    if ~isstruct(S), return; end
    f = fieldnames(S);
    for i = 1 : length(f)
        k = f{i};
        if isfield(C, k) && isstruct(C{k}) && isstruct(S{k})
            C{k} = deepmerge(C{k}, S{k});
        else
            C{k} = S{k};
        end
    end
end

function C = setfield_deep(C, path, val)
    parts = regexp(path, '\. ', 'split');
    C = setfield_deep_rec(C, parts, val);
end

function C = setfield_deep_rec(C, parts, val)
    if lengthval{1} == 1
```
```matlab
C.(parts{1}) = val;
else
    k = parts{1};
    if ~isfield(C,k) || ~isstruct(C.(k)), C.(k)=struct(); end
    C.(k) = setfield_deep_rec(C.(k), parts(2:end), val);
end
end

function assert_file(fname)
    if exist(fname, 'file')~=2
        error('SCAN:missingFile','Required file missing: %s', fname);
    end
end

function s = ternIID, a, b)
    if cond, s=a; else, s=b; end
end
```

<caption>Listing 36: Parameter scan orchestrator</caption>

The script in Listing 36 (p. 278) provides a general-purpose driver for parameter-sensitivity and robustness scans within TSRT. It performs automated sweeps over any scalar constant in the TSRT parameter tree (e.g., scission smoothness, hindrance slopes, surface widths), reconstructs the required field configurations, evaluates the corresponding observables (typically the fission energy $E_{\text{fiss}}$ of $^{235}\text{U}$), and stores the results under `results/scan/` without modifying the calibrated constants used elsewhere in the paper. This module is intentionally lightweight: it does not overwrite `tsrtconstants.m` or `tsrtconstantsCalibrated.mat`, and requires no MATLAB toolboxes. Each scan produces an accompanying `.mat` file, a `.csv` table suitable for external plotting, and a quick-look `.png` diagnostic.

The user controls the scan through name-value arguments specifying the parameter name, the values to sweep, and the diagnostic task. Default ranges used in this work are summarized in Table 43 (p. 284); runnable examples appear in Listing 37 (p. 283). Examples below illustrate typical usage and can be copied verbatim into MATLAB:

```
running the parameter scan driver examples
```

```matlab
% Example 1: scan scission smoothness for U-235 Efiss (default task)
scan_parameters('param','fission.scission_w_fm', ...
    'values', linspace(0.3, 3.0, 12), ...
    'task','Efiss_U235');
```

```matlab
% Example 2: same scan but also report a forced-to-target column (170 MeV)
scan_parameters('param','fission.scission_w_fm', ...
    'values', [0.3 0.5 1.0 1.5 2.0], ...
    'task','Efiss_U235', ...
    'force_target', true, 'Efiss_target_MeV', 170.0);
```

```matlab
% Example 3: scan hindrance slope (lifetime robustness diagnostic)
scan_parameters('param','HF_slope_per_rank', ...
    'values', [-0.2 -0.1 0 0.1 0.2], ...
    'task','Efiss_U235');
```

<caption>Listing 37: Examples: running the parameter scan driver</caption>

All scans are resolved independently of the main workflow, and none of them introduce per-observable tuning. They serve exclusively as robustness diagnostics for the curvature-driven fission energetics and associated TSRT observables discussed in Appendix Y.11 and Appendix Y.3.
<table><caption>Table 43: Recommended default parameter ranges for TSRT robustness scans. These ranges are conservative and intended for “sanity” sweeps; domain-specific work may justify narrower or wider windows. Scans do not modify calibrated constants and write results under `results/scan/`.</caption><thead><tr><th>Parameter (TSRT key)</th><th>Physical meaning</th><th>Suggested range</th><th>Typical step(s)</th></tr></thead><tbody><tr><td>fission.scission_w_fm</td><td>Smoothness of scission partition (fm)</td><td>0.3-3.0</td><td>10-15 points (linspace)</td></tr><tr><td>HF_slope_per_rank</td><td>Forbiddenness hindrance slope (log<sub>10</sub> per rank)</td><td>-0.2+0.2</td><td>{-0.2, -0.1, 0, 0.1, 0.2}</td></tr><tr><td>a_surf_fm</td><td>Surface diffuseness (fm)</td><td>0.45-0.80</td><td>6-8 points (linspace)</td></tr><tr><td>r0_fm</td><td>Bulk radius parameter (fm)</td><td>1.10-1.30</td><td>6-8 points (linspace)</td></tr><tr><td>gamma</td><td>Curvature strength (m<sup>-2</sup>)<sup>131</sup></td><td>0.5-2.0× nominal</td><td>5-7 points (logspace)</td></tr></tbody></table>

Notes: The production value is fixed (see Section 15.2); scanning `gamma` is for sensitivity only.

## Y.18 Master Orchestrator (Run-All)

Dependencies: all preceding modules (Listings 2 (p. 213)–36 (p. 278)).

`runall.m`

```matlab
% FILE: runall.m
% PURPOSE: orchestrate full reproduction in logical order
function runall()
% RUNALL Deterministic, physics-based end-to-end run for TSRT results.

rang(0,'twister'); % fix RNG for exact reproducibility of results/figures

% 0) Preflight
tsrt_preflight_checks();

% 1) Constants
C = tsrtconstants();

% 2) Main physics benchmark
results = run_benchmark_full();

% 3) Save any open figures with a clear prefix
saveRemaining figures('tsrt');

% 4) Console summary (if provided by the benchmark)
if isstruct(results) && isfield(results,'summary')
    fprintf('\n[runall] Summary:\n%s\n', resultsire summary);
end

end

function saveRemaining figures(prefix)
% Save all open figures to fig_misc/prefix_XX.pdf and .eps
figs = findall(0,'Type','figure');
figs = figs(is graphics(figs,'figure'));
if isempty(figs), return; end

outDir = 'fig_misc';
if ~exist(outDir,'dir'), mkdir(outDir); end

nums = zeros(size(figs));
for i=1: numel(figs)
    try,nums(i) = get(figs(i),'Number'); catch,nums(i) = i; end
end
```
```R
[~, ord] = sort(nums);
figs = figs(ord);

for k = 1:numel(figs)
    f = figs(k);
    fn = sprintf('%s_%02d', prefix, k);
    try
        print(f, fullfile(outDir, [fn '.pdf'])), '-dpdf', '-painters');
        print(f, fullfile(outDir, [fn '.eps'])), '-depsc', '-painters');
    catch ME
        fprintf(2, '[runall] Could not save figure %d: %s\n', k, ME.message);
    end
end
end
```

<caption>Listing 38: Run-all entry point</caption>
## Summary of MATLAB modules and corresponding figures/tables.

<table><thead><tr><td>MATLAB module</td><td>Paper output(s)</td></tr></thead><tbody><tr><td>screeningmakefig.m</td><td>Figure 2 (p. 12)</td></tr><tr><td>reproducehindranceTSRT.m</td><td>Figure 6 (p. 86)</td></tr><tr><td>oesmakefig.m</td><td>Figure 3 (p. 24)</td></tr><tr><td>gamma rates.m + lifetime_pipe line.m</td><td>Tables 2 (p. 11)-18 (p. 110)</td></tr><tr><td>tsrtbindingenergy.m + bindingplotwithexp.m (data: tsrt_bindingscan_cln.csv)</td><td>Figure 1 (p. 10) (left)</td></tr></tbody></table>

## Reproduction Workflow (Entry Points)

RunLifetime_All.m Master workflow to reproduce all lifetime results ($\alpha$, $\beta$, SF) and produce consolidated tables.

* Calls: RunLifetimeCalculations.m, RunLifetimeCalculations_beta.m, RunLifetimeCalculations_sf.m, tsrtQvalues.m, tsrt_alpha_L_for_transition.m, tsrt_alpha_width.m, tsrt_emit_tables.m, tsrtconstants_calibrated.mat

RunLifetime_All.m

```matlab
function RunLifetime_All()
% Run α, β, sf calibrations in sequence, then emit the combined table.
clear functions; rehash toolboxcache;

RunLifetimeCalculations;      % alpha only, saves/merges constants
RunLifetimeCalculations Beta;  % beta only, merges into the same file
RunLifetimeCalculations_sf;  % sf only, merges into the same file

tsrt_emit_tables('tsrt_benchmark_diag.csv', ...
                  'tab_lifetime_bench_highres.tex', ...
                  'tsrt_benchmark_diag_FILLED.csv');

fprintf(' [All] Emission complete. \n');

clear C;
tsrt_emit_tables('tsrt_benchmark_diag.csv', 'tab_lifetime_bench_highres.tex',
                  'tsrt_benchmark_diag_FILLED.csv');
end
```

<caption>Listing 39: RunLifetime_All.m</caption>

RunLifetimeCalculations.m Core lifetime runner invoked by the master workflow; computes lifetimes for selected nuclides/sets.

* Calls: tsrtQvalues.m, tsrt_alpha_L_for_transition.m, tsrt_alpha_width.m

* Called by: RunLifetime_All.m

RunLifetimeCalculations.m

```matlab
function RunLifetimeCalculations()
%RUNLIFETIMECALCULATIONS TSRT alpha benchmark (physics-grounded, robust).
% Sequence:
```
```
% (A) Calibrate C_alpha on 210Po with S_pre(anchor)=1.
```

```
% (C) Fit preformation slopes (kZ,kN,kNN) on Ra-226, Th-232, U-238.
```

```
% (B) Calibrate gamma_shell on 209Bi (after slopes) to keep Bi exact.
```

```
% (D) Save constants (MERGE into tsrtconstantsCalibrated.mat).
```

```
% (E) Internal alpha-only benchmark (no CSV/TeX emission here).
```

```
%
```

```
% Uses: tsrtQualities.m, tsrt_alpha_width.m
```

```
fprintf(' [Preflight] MATLAB: %s\n', version);
```

```
% ---- constants container (flat fields) ----
```

```
C = struct();
```

```
C.Century = 6.48104e-06; % (beta const unused here)
```

```
C.SigmaBeta = 0.0;
```

```
% anchor at Po-210; gamma0=0 keeps Po-210 exact with S_pre(anchor)=1
```

```
C.preform = struct('refZ',84,'refA',210,'kZ',0,'kN',0,'kNN',0,...
                     'gamma0',0,'gamma_shell',0);
```

```
% Fixed barrier shape (keep barrier physics intact)
```

```
C.shape_best = struct('r0',1.20,'a',0.55,'r0c',1.25);
```

```
% (A) Calibrate C_alpha on 210Po (S_pre=1 at the anchor)
```

```
C.Calpha = calibrate_C alpha po210(C.shape_best);
```

```
fprintf(' [Cal] C alpha calibrated on 210Po: %.6g s~1\n', C.Calpha);
```

```
% (C) Fit base preformation slopes (kZ,kN,kNN) on actinide triplet (shell=0 during
  fit)
```

```
triplet = [ ...
  struct('Z',88,'A',226,'t12_exp_s', 1600 * 365.25*24*3600, 'label', '$^{226}$Ra'
]; ...
```

```
  struct('Z',90,'A',232,'t12_exp_s', 1.405e10 * 365.25*24*3600, 'label', '$^{232}$Th'); ...
```

```
  struct('Z',92,'A',238,'t12_exp_s', 4.468e9 * 365.25*24*3600, 'label', '$^{238}$U'); ...
```

```
];
```

```
% temporarily keep gamma_shell=0 for slope fit
```

```
C.preform.gamma_shell = 0;
```

```
C.preform = fit_preformation_slopes(C, triplet);
```

```
fprintf(' [Fit] kZ=%6f, kN=%6f, kNN=%6f (gamma0=%3f; anchor Spre=1)\n', ...
       C.preform.kZ, C.preform.kN, C.preform.kNN, C.preform.gamma0);
```

```
% (B) Calibrate shell AFTER slope fit to keep 209Bi exact
```

```
C.preform.gamma_shell = calibrate_shell_gamma_bi209(C);
```

```
% (D) Save constants by MERGING into tsrtconstantsCalibrated.mat
```

```
save_mergeConstants('tsrtconstantsCalibrated.mat', C);
```

```
fprintf(' [Preflight] C alpha=%g C lambda=%g\n', C.Calpha, C.Clambda);
```

```
% (E) Internal alpha benchmark only (no emit from here)
```

```
bench = [ ...
  struct('Z',84,'A',210,'label', '$^{210}$Po', 't12_exp_s', 138.376*24*3600); ...
  struct('Z',88,'A',226,'label', '$^{226}$Ra', 't12_exp_s',
        1599.97*365.25*24*3600); ...
  struct('Z',90,'A',232,'label', '$^{232}$Th', 't12_exp_s', 1.405e10
*365.25*24*3600); ...
  struct('Z',92,'A',238,'label', '$^{238}$U', 't12_exp_s', 4.468e9
*365.25*24*3600) ...
];
```

```
];
```

```
[N, meanAbs, medianAbs, maxAbs] = run_internal_benchmark(C, bench);
```
```matlab
if N>0
    fprintf(' [RunLifetimeCalculations] Summary: \n');
    fprintf('N=%d, mean|Δlog10|=%3g, median=%3g, max=%3g\n', N, meanAbs,
        medianAbs, maxAbs);
else
    fprintf(' [RunLifetimeCalculations] Summary: \nNo successful rows\n');
end
end

% =========================== Calibrations & Model ===========================

function C_alpha = calibrate_C_alpha-po210.shape)
    Zpo = 84; Apo = 210;
    [Qpo, ~, ~] = tsrtQvalues(Zpo, Apo);
    t12_exp-po = 138.376*24*3600; % seconds
    [~, wdiag] = tsrt_alpha-width(Zpo, Apo, Qpo, shape);
    Ppo = wdiag.P;
    if ~isfinite(Ppo) || Ppo <= 0
        error('Non-positive/NaN penetrability P for 210Po.');
    end
    C_alpha = log(2) / (t12_exp-po * Ppo);
end

function gamma_shell = calibrate_shell-gamma_bi209(C)
    Z = 83; A = 209; N = A - Z;
    [Q, ~, ~] = tsrtQvalues(Z,A);
    [~, w] = tsrt_alpha-width(Z, A, Q, C.shape_best);
    if ~isfield(w,'P') || ~isfinite(w.P) || w.P<=0
        error('Penetrability invalid for 209Bi.');
    end
    t_exp = 1.9e19 * 365.25*24*3600; % seconds (per the table)
    dZ = Z - C.preform.refZ;
    dN = N - (C.preform.refA - C.preform.refZ);
    Spre_base = exp( - (C.preform.kZ*dZ + C.preform.kN*dN + C.preform.kNN*(dN.^2)) );
    u = shellhap(Z, N, C.preform.refZ, C.preform.refA - C.preform.refZ);
    rhs = (log(2) / (C.C_alpha * w.P * t_exp)) / max(Spre_base, realmin);
    if rhs <= 0 || ~isfinite rhs)
        error('Shell calibration invalid rhs for 209Bi.');
    end
    if u <= 0
        error('Shell metric difference (u) <= 0 for 209Bi vs anchor; cannot calibrate
            gamma_shell.');
    end
    gamma_shell = -log Information theory / u;
    fprintf(' [shell] Calibrated gamma_shell=%4f using 209Bi (S_need=%s, u=%4g)\n',
        ...
            gamma_shell, pretty_sci Information theory, u);
end

function pre = fit_preformation_slopes(C, triplet)
    X = zeros(numel(triplet),3);
    Y = zeros(numel(triplet),1);
    for j=1:numel(triplet)
        Z = triplet(j).Z; A = triplet(j).A; N = A-Z;
        [Q, ~, ~] = tsrtQvalues(Z,A);
        [~, w] = tsrt_alpha-width(Z, A, Q, C.shape_best);
        if ~isfield(w,'P') || ~isfinite(w.P) || w.P<=0
            error('Invalid penetrability for %s', triplet(j).label);
```
```matlab
end

te = triplet(j).t12_exp_s;
rhsisp = log( log(2)/(te*C.C_alpha*w.P) );
dZ = Z - C.preform.refZ;
dN = N - (C.preform.refA - C.preform.refZ);
Y(j) = rhsistock; % gamma'em during fit
X(j, :) = -[dZ, dN, dN.^2];
end
p = X \ Y;
pre = C.preform;
pre.kZ = p(1);
pre.kN = p(2);
pre.kNN = p(3);
end

% ======================== Internal benchmark ========================
function [N, meanAbs, medianAbs, maxAbs] = run_internal_benchmark(C, rows)
errs = [];
for j=1:numelrows)
    Z = rows(j).Z; A = rows(j).A; N = A-Z;
    [Q,~,~] = tsrtQvalues(Z,A);
    t_pred = model_alpha_t12(Z, A, N, Q, C);
    t_exp = rows(j).t12_exp_s;
    if isfinite(t_pred) && t_pred>0 && isfinite(t_exp) && t_exp>0
        dlog = log10(t_pred) - log10(t_exp);
        fprintf(' [ok] %s Z=%3d A=%3d Q=%3d MeV t12=%s s exp=%s s Δlog10=%g
', ... 
        rows(j).label, Z, A, Q, pretty_sci(t_pred), pretty_sci(t_exp), dlog);
        errs(end+1) = abs(dlog);
    else
        fprintf(' [E] %s failed: NaN/invalid t12
', rows(j).label);
    end
end
if isempty(errs), N=0; meanAbs=NaN; medianAbs=NaN; maxAbs=NaN; return; end
N = realerr;
meanAbs = mean(errs);
medianAbs = median(errs);
maxAbs = max(errs);
end

% ======================== Core model (shared) ========================
function t12 = model_alpha_t12(Z, A, N, Qalpha, C)
    [~, w] = tsrt_alpha_width(Z, A, Qalpha, C.shape_best);
    if ~isfield(w,'P') || ~isfinite(w.P) || w.P<=0, t12 = NaN; return; end
    dZ = Z - C.preform.refZ;
    dN = N - (C.preform.refA - C.preform.refZ);
    Spre_base = exp( - (C.preform.kZ*dZ + C.preform.kN*dN + C.preform.kNN*(dN.^2)) );
    u = shell metric(Z, N, C.preform.refZ, C.preform.refA - C.preform.refZ);
    Spreem = exp( - C.preform GAMMA SHELL * u );
    Spre = max(realmin, Spre_base * Spre SHELL);
    rate = C.C_alpha * w.P * Spre;
    if ~isfinite(rate) && rate>0, t12 = NaN; else, t12 = log(2) / rate; end
end

function u = shell metric(Z, N, Zref, Nref)
f = 1/(1+abs(Z-82)) + 1/(1+abs(N-126));
```
```matlab
f0 = 1/(1+abs(Zref-82)) + 1/(1+abs(Nref-126));
u = f - f0;
if u < 0, u = 0; end
end
```
```matlab
% =========================== Utilities ===========================
```

```matlab
function s = pretty_sci(x)
    if ~isfinite(x) || x==0, s='--'; return; end
    m = abs(x);
    if m >= 1e4 || m < 1e-3
        s = sprintf('%.3g', x);
    else
        s = sprintf('%.6g', x);
    end
end
```

```matlab
function save_mergere constants (fname, Cnew)
% Merge fields of Cnew into existing MAT file (if present), then save.
    Cold = struct();
    if exist(fname, 'file')==2
        Cold = load(fname);
    end
    C = merge_structs(Cold, Cnew);
    save(fname, '-struct', 'C');
end
```

```matlab
function C = merge_structs(A,B)
    C = A;
    f = fieldnames(B);
    for k=1:numel(f)
        key = f{k};
        if isstruct(B.(key))
            if isfield(C,key) && isstruct(C.(key))
                C.(key) = merge_structs(C.(key), B.(key));
            else
                C.(key) = B.(key);
            end
        else
            C.(key) = B.(key);
        end
    end
end
```

<caption>Listing 40: RunLifetimeCalculations.m</caption>

RunLifetimeCalculations_beta.m β-decay lifetime sub-pipeline used where applicable.

* Called by: RunLifetime_All.m

RunLifetimeCalculations_beta.m

```matlab
function RunLifetimeCalculations_beta()
%RUNLIFETIMECALCULATIONS_BETA TSRT beta benchmark (Q^5, simple global form).
% Calibrates to:
% 60Co (allowed-ish) and 137Cs (allowed) to fix C_beta and bN (with bZ=0),
% then sets an empirical forbidden hindrance for 14C.
```
%
% Model:
% lambda_beta = C_beta * Qbeta^power * exp(bZ*Z + bN*N)
% For 14C only: lambda <- lambda * exp(eta_forb)
% power is fixed to 5 (standard phase space scaling).
%
% Output: updates/creates tsrtconstants Calibrated.mat with C.beta fields.

fprintf(' [Preflight] MATLAB: %s\n', version);

% ---- anchors ----
% Experimental half-lives (seconds)
t14C = 5729.88 * 365.25*24*3600; % 14C
t60Co = 5.27129 * 365.25*24*3600; % 60Co
t137Cs = 29.9994 * 365.25*24*3600; % 137Cs

% Z, A and Qbeta (MeV) via the Q loader (fallback to table values if needed)
Z14=6; A14=14; N14=A14-Z14; Q14 = coerceQbeta(Z14,A14,0.156476);
Z60=27; A60=60; N60=A60-Z60; Q60 = coerceQbeta(Z60,A60,4.34875876095797);
Z137=55; A137=137; N137=A137-Z137;Q137= coerceQbeta(Z137,A137,2.19265955373013);

% Fixed power and bZ
power = 5;
bZ = 0;

% Solve for bN and C_beta from Co and Cs
r60 = log( log(2)/t60Co ) - power*log(Q60);
r137 = log( log(2)/t137Cs) - power*log(Q137);

dN = (N137 - N60);
if dN == 0
    error('Internal: N(137Cs) == N(60Co); cannot solve for bN.');
end
bN = (r137 - r60) / dN;
lnC_beta = r60 - bN*N60;
C_beta = exp(lnC_beta);

% Now set eta_forb so 14C matches exactly (applies only to 14C)
rate14_base = C_beta * (Q14^power) * exp(bZ*Z14 + bN*N14);
eta_forb = log( (log(2)/t14C) / rate14_base );

% Package constants
beta = struct();
beta.C_beta = C_beta;
beta.power = power;
beta.bZ = bZ;
beta.bN = bN;
beta.eta_forb = eta_forb;
betagamma_beta = 0; % not used here; included for completeness

fprintf(' [beta] C_beta=%g, bZ=%g, bN=%g, eta_forb=%g, gamma_beta=%g, power=%d\n',
        ...
        beta.C_beta, beta.bZ, beta.bN, beta.eta_forb, betagamma_beta, beta.power);

% --- Save back into constants file (merge) ---
C = struct();
if exist('tsrtconstants Calibrated.mat', 'file') == 2
    C = load('tsrtconstants Calibrated.mat');
end

Cwald = beta;
save('tsrtconstantsCalibrated.mat', '-struct', 'C');

% --- Internal check/print ---
bench = [ ... 
    struct('Z',Z14, 'A',A14, 'lab',~{14}C', 'Q',Q14, 't',t14C, 'forb',true )
; ...
    struct('Z',Z60, 'A',A60, 'lab',~{60}Co', 'Q',Q60, 't',t60Co, 'forb',false);
...
    struct('Z',Z137, 'A',A137, 'lab',~{137}Cs', 'Q',Q137,'t',t137Cs,'forb',false)
...
];

errs = zeros(0,1);
for j=1:numel(bench)
    b = bench(j);
    rate = beta.C_beta * (b.Q^beta.power) * exp beta.bZ*b.Z + beta.bN*(b.A-b.Z));
    if b.forb, rate = rate*exp beta.eta_forb); end
    t12_pred = log(2) / rate;
    dlog = log10(t12_pred) - log10(b.t);
    fprintf(' [ok] %s Z=%3d A=%3d Q=%3.3g MeV t12=%s s exp=%s s Δlog10=%g
',
        b.lab, b.Z, b.A, b.Q, pretty_sci(t12_pred), pretty_sci(b.t), dlog);
    errs(end+1,1) = abs(dlog);
end

fprintf(' [RunLifetimeCalculations_beta] Summary: 
');
fprintf('N=%d, mean|Δlog10|=%g, median=%g, max=%g
', ...
    numel(errs), mean(errs), median(errs), max(errs));

% Optionally re-emit tables here or let RunLifetime_All do it after all modes
% tsrt_emit_tables('tsrt_benchmark_diag.csv', 'tab_lifetime_bench_highres.tex',
% tsrt_benchmark_diag_FILLED.csv');
end

% ------------------ helpers ------------------
function Q = coerceQbeta(Z,A,Qfallback)
    [~,Qb,~] = tsrtQvalues(Z,A);
    if ~(isnumeric(Qb) && isfinite(Qb) && Qb>0)
        Q = Qfallback;
    else
        Q = Qb;
    end
end

function s = pretty_sci(x)
    if ~isfinite(x) || x==0, s='--'; return; end
    m = abs(x);
    if m >= 1e4 || m < 1e-3
        s = sprintf('%.3g', x);
    else
        s = sprintf('%.6g', x);
    end
end

<caption>Listing 41: RunLifetimeCalculations_beta.m</caption>
RunLifetimeCalculations_SF.m Spontaneous fission lifetime sub-pipeline.

* Called by: RunLifetime_All.m

RunLifetimeCalculations_SF.m

function RunLifetimeCalculations_SF()

% Calibrate minimal TSRT-like spontaneous fission law on:

% 240Pu, 252Cf, 254Fm

% Save into tsrtconstantsCalibrated.mat by merge.

fprintf(' [Preflight] MATLAB: %s
', version);

% Load existing constants to keep full C together

C = struct();

if exist('tsrtconstantsCalibrated.mat', 'file') == 2

C = load('tsrtconstantsCalibrated.mat');

end

% Data (seconds)

rows = [ ... 
    struct('Z',94,'A',240,'t', 6559.87*365.25*24*3600, 'label', '^(240)Pu'); ...
    struct('Z',98,'A',252,'t', 2.64495*365.25*24*3600, 'label', '^(252)Cf'); ...
    struct('Z',100,'A',254,'t',0.00136893*365.25*24*3600, 'label', '^(254)Fm') ...
];

% Model: rate = C_SF * exp(g0*(Z^2/A) + g1*N - N0)

NO = 152;

M = []; y = [];

for j=1:numel rows

Z = rows(j) . Z; A = rows(j) . A; N = A - Z;

fiss = (Z^2)/A;

dn = abs(N - NO);

rate = log(2)/rows(j) . t;

M = [M; 1, fiss, dn];

y = [y; log(rate)];

end

p = M\y; % [ln C_SF, g0, g1]

C_SF = exp(p(1)); g0 = p(2); g1 = p(3);

C SF = struct('C_SF', C_SF, 'g0', g0, 'g1', g1, 'NO', NO);

fprintf(' [sf] C_SF = %.6g, g0 = %.6g, g1 = %.6g (NO = %d) 
', C_SF, g0, g1, NO);

% Quick check

for j=1:numel rows

Z = rows(j) . Z; A = rows(j) . A; N = A - Z;

t_pred = tsrt_SFPredict(Z,A,N,C);

dlog = log10(t_pred) - log10 rows(j) . t;

fiss = (Z^2)/A; dn = abs(N - NO);

fprintf(' [ok] %s Z=%3d A=%3d fiss=%4f |N-%d|=%d t12=%s s exp=%s s Δ

log10=%g
'

rows(j).label, Z, A, fiss, NO, dn, pretty_sci(t_pred), pretty_sci rows(j).

t), dlog);

end

fprintf(' [RunLifetimeCalculations_SF] Summary: 
N = %d, mean |Δlog10| = %g, median = %g,

max = %g
'

numel (rows), 0, 0, 0);
```matlab
% Merge-save
save_merge{ constants('tsrtconstants Calibrated.mat', C);
end

function t12 = tsrt_SFPredict(Z,A,N,C)
    fiss = (Z^2)/A;
    dn = abs(N - C.sf.NO);
    rate = C.sf.C_SF * exp(C.sf.g0 * fiss + C.sf.g1 * dn);
    t12 = log(2)/rate;
end

function s = pretty_sci(x)
    if ~isfinite(x) || x==0, s='--'; return; end
    m = abs(x);
    if m >= 1e4 || m < 1e-3, s = sprintf('%3g', x);
    else, s = sprintf('%6g', x);
end
end

function save_mergeConstants(fname, Cnew)
    Cold = struct();
    if exist(fname, 'file')==2, Cold = load(fname); end
    C = merge_structs(Cold, Cnew);
    save(fname, '-struct', 'C');
end

function C = merge_structs(A,B)
    C = A;
    f = fieldnames(B);
    for k=1:numel(f)
        key = f{k};
        if isstruct(B{(key)})
            if isfield(C,key) && isstruct(C{(key)})
                C{(key)} = merge_structs(C{(key)}, B{(key)});
            else
                C{(key)} = B{(key)};
            end
        else
            C{(key)} = B{(key)};
        end
    end
end
```

<caption>Listing 42: RunLifetimeCalculations_SF.m</caption>

runall.m Benchmark/figure orchestrator for the paper's summary plots and comparison tables.

*   Calls: `run_benchmark_full.m`, `tsrtQvalues.m`, `tsrt_alpha_L_for_transition.m`, `tsrt_alpha_t12.m`, `tsrt_alpha_width.m`, `tsrt_get_exphalf_life.m`, `tsrt_preflight_checks.m`, `tsrt_preformation_prefactor.m`, `tsrt_shellMetric.m`, `tsrtconstants.m`, `tsrtconstants Calibrated.mat`.

*   Code: see Listing 38 (p. 284).

run_benchmark_full.m Benchmark runner used by runall.m.
* Calls: `tsrtQvalues.m`, `tsrt_alpha_L_for_transition.m`, `tsrt_alpha_t12.m`, `tsrt_alpha_width.m`, `tsrt_get_exphalf_life.m`, `tsrt_preformation_prefactor.m`, `tsrt_shell metric.m`, `tsrtconstants Calibrated.mat`

* Called by: `runall.m`

`run_benchmark_full.m`

```matlab
function results = run_benchmark_full()
%RUN_BENCHMARK_FULL Verbose TSRT alpha half-life benchmark & console report.
% Uses:
% - tsrtQvalues.m                  [Qbeta, Qalpha] = tsrtQvalues(Z,A)
% - tsrt_alpha_t12.m                 [t12_s, diag] = tsrt_alpha_t12(Z,A,Qalpha,shape,opts)
% - tsrt_alpha_width.m              [W, diag] = tsrt_alpha_width(Z,A,Qalpha,shape,opts)
% - tsrt_get_exphalf_life.m         t12_exp_s = tsrt_get_exphalf_life(Z,A,mode)
% Assumes constants were calibrated by RunLifetimeCalculations:
% - tsrtconstants Calibrated.mat with fields C_alpha, C cerba, etc.
% Prints progress to console and returns a struct with summary and rows.

fprintf(' [Preflight] MATLAB: %s\n', version);
C = loadConstants_local('tsrtconstants Calibrated.mat');
fprintf(' [Preflight] C alpha=%5.6g C lambda=%5.6g\n', C.C alpha, getfield_safe(C, 'C lambda', 'NaN'));

% Target set: include calibration and actinide triplet
targets = [ ...
    struct('Z',84,'A',210,'label', '$^{210}$$Po', 'mode', 'alpha'); ...
    struct('Z',88,'A',226,'label', '$^{226}$$Ra', 'mode', 'alpha'); ...
    struct('Z',90,'A',232,'label', '$^{232}$$Th', 'mode', 'alpha'); ...
    struct('Z',92,'A',238,'label', '$^{238}$$U', 'mode', 'alpha') ...
];

% Shape: defaults used, after RunLifetimeCalculations shad canned best shape
already
shape = struct('r0',1.20,'a',0.50,'r0c',1.20);
opts = struct('constants_path','tsrtconstants Calibrated.mat','verbose',true);
rows = []; % will store structs
fprintf(' [Run] Evaluating %d targets (alpha)\n', numel(targetss)));

for i = 1:numel(targetss)
    Z = targets(i).Z; A = targets(i).A; lbl = targets(i).label;

    try
        [Qalpha, ~, Qbeta minus] = tsrtQvalues(Z, A);
    catch ME
        fprintf(2,' [E] tsrtQvalues failed for %s (Z=%d,A=%d): %s\n', stripTex(lbl), Z, A, ME.message);
        continue;
    end

    if ~isfinite(Qalpha) || Qalpha <= 0
        fprintf(2,' [W] Non-positive/NaN Qalpha for %s (Z=%d,A=%d). Skipping.\n', stripTex(lbl), Z, A);
        continue;
    end
```
end

% Experimental value (seconds) from the CSV
t12_exp = tsrt_get_exphalf_life(Z, A, 'alpha');
if ~isfinite(t12_exp)
    fprintf(2,' [W] Missing exp t1/2 for %s (Z=%d,A=%d) in*_benchmark_diag*.
csv. Skipping.\n', stripTex(lbl), Z, A);
    continue;
end

% Predict with TSRT
try
    [t12_pred, diag] = tsrt_alpha_t12(Z, A, Qalpha, shape, opts);
    catch ME
        fprintf(2,' [E] tsrt_alpha_t12 failed for %s (Z=%d,A=%d): %s\n', stripTex
            (lbl), Z, A, ME.message);
        continue;
    end

    dlog = abs(log10(t12_pred) - log10(t12_exp));
    fprintf(' [ok] %-10s Z=%3d A=%3d Qα=%5f MeV t12=%6g s exp=%6g s Δlog10
    =%.4f\n', ...;
    stripTex(lbl), Z, A, Qalpha, t12_pred, t12_exp, dlog);

    rows = [rows; struct(...
        'Z', Z, 'A', A, 'isotope', stripTex(lbl), 'Qalpha_MeV', Qalpha, ...
        't12_pred_s', t12_pred, 't12_exp_s', t12_exp, 'dlog10', dlog, ...
        'diag', diag)];
end

if isempty rows)
    warning('No rows produced: Check tsrtQvalues() and experimental CSV
    availability.');
    results = struct('summary', 'No results', 'rows', rows);
    return;
end

% Summary
dlogs = [rows.dlog10];
mean_abs = mean(dlogs);
med_abs = median(dlogs);
max_abs = max(dlogs);

 fprintf('\n[Summary]\n');
 fprintf(' N = %d\n', numeral rows));
 fprintf(' mean(|Δlog10|) = %.4f\n', mean_abs);
 fprintf(' median(|Δlog10|) = %.4f\n', med_abs);
 fprintf(' max(|Δlog10|) = %.4f\n', max_abs);

 results = struct('summary', sprintf('N=%d, mean|Δlog10|=%4f, median=%4f, max
    =%.4f', ...;
    numeral rows), mean_abs, med_abs, max_abs), 'rows', rows);
end

% ---------- helpers ----------
function C = loadConstants_local(path)
    if exist(path,'file') =~ 2
```matlab
error('Constants file not found: %s. Run RunLifetimeCalculations first.', path);
end

C = load(path);
if ~isfield(C, 'C Alpha')
    error('Constants file missing C alpha: %s', path);
end
end

function val = getfield_safe(S, f, default)
    if isfield(S,f), val = S.(f); else, val = default; end
end

function x = stripTex(s)
    s = string(s);
    x = regexprep(s, '^\$|\$$', '');
end
```

<caption>Listing 43: `run_benchmark_full.m`</caption>

## Calibration, Constants, and Core Lifetime Utilities

**Deuteron-based energy normalization.** All curvature-derived energies in this paper use a single global scale `C.norm_energy` to map the code's curvature integral to SI/Joules (and hence to MeV). We fix this once by matching the deuteron binding (2.224 MeV). Listing 44 (p. 297) computes the required scale using the same curvature integral as the main pipeline, writes it to `tsrtconstantsCalibrated.mat`, and prints the exact line to paste into `tsrtconstants.m`. There is no per-nucleus tuning; the same `C.norm_energy` is used everywhere.

### `TSRT_Calib_NormEn_Deuteron.m`

```matlab
% FILE: TSRT_Calib_NormEn_Deuteron.m
% PURPOSE: Calibrate C.norm_energy so the deuteron binding equals 2.224 MeV.
% OUTPUTS: tsrtconstantsCalibrated.mat with field C.norm_energy
% NOTES: Robust to missing tsrtmakegrid(): uses a Simpson-weighted spherical grid.

function TSRT_Calib_NormEn_Deuteron()
    fprintf(' [CAL] TSRT deuteron normalization (target 2.224 MeV) 
');

% ---- constants ----
E_target_MeV = 2.224;
eV_per_MeV = 1.0e6;
J_per_eV = 1.602176634e-19;
E_tar_J = E_target_MeV * eV_per_MeV * J_per_eV;

% ---- get constants; tolerate minimal tsrtconstants() ----
try
    C = tsrtconstants();
    catch
        C = struct();
    end

% ---- grid: try official, else robust fallback ----
G = [];
if exist('tsrtmakegrid.m', 'file')==2
    try, G = tsrtmakegrid(C); catch, end
```
```c
if isempty(G), try, G = tsrtmakegrid(); catch, end, end
end
if isempty(G), G = local_simpson_grid(C); fprintf(' [CAL] Using internal Simpson-
weighted spherical grid fallback. \n'); end

% ---- build a minimal deuteron proxy (parent vs constituents) ----
% We compute the raw curvature energy for a 2-nucleon configuration with
norm_energy=1,
% then solve norm_energy = E_target / E_raw.
C1 = C; C1.norm_energy = 1.0;

nuc_parent = struct('types', {'proton', 'neutron'}, ...
                         'R', [0 0; 0 0; 0 0]);

Xi_parent = tsrtcompositefield(G, nuc_parent, 0);
Kp = tsrtcurvaturetensor(G, Xi_parent, C1);
E_raw_J = tsrtenergyfromcurvature(G, Kp, C1); % in Joules since norm_energy=1 is
                                                SI here

fprintf(' [CAL] Raw integral (norm=1): E_raw = %.12e J\n', E_raw_J);
fprintf(' [CAL] Target (2.224 MeV):   E_tar = %.12e J\n', E_tar_J);

C_norm = E_tar_J / E_raw_J;
fprintf(' [CAL] C.norm_energy = E_tar / E_raw = %.12g\n', C_norm);

% ---- write/merge calibrated constants ----
Cout = struct();
if exist('tsrtconstantsCalibrated.mat', 'file')==2
    S = load('tsrtconstantsCalibrated.mat');
    f = fieldnames(S);
    for i=1:numel(f), Cout.(f{i}) = S.(f{i}); end
end
Cout.norm_energy = C_norm;
save('tsrtconstantsCalibrated.mat', '-struct', 'Cout');
fprintf(' [CAL] Wrote tsrtconstantsCalibrated.mat (norm_energy). \n\n');

fprintf('Paste this into tsrtconstants.m near the other globals: \n');
fprintf('     C.norm_energy = %.12g; %% deuteron-matched 2.224 MeV\n', C_norm);
fprintf(' \n [CAL] Done. \n');
end

% ---- local helper: Simpson spherical grid (with sin theta) in theta weights ----
function G = local_simpson_grid(C)
    if ~isfield(C, 'N_r') || isempty(C.N_r), C.N_r = 96; end
    if ~isfield(C, 'N_theta') || isempty(C.N_theta), C.N_theta = 64; end
    if ~isfield(C, 'N phi') || isempty(C.N phi), C.N phi = 64; end
    if ~isfield(C, 'r_max') || isempty(C.r_max), C.r_max = 12e-15; end

Nr = max(3,C.N_r); Nt = max(3,C.N_theta); Np = max(3,C.N phi);
r = linspace(0, C.r_max, Nr);
theta = linspace(0, pi, Nt);
phi = linspace(0, 2*pi, Np);

wr = simpson weights(Nr) * (C.r_max/(Nr-1));
wt = simpson weights(Nt) * (pi/(Nt-1));
wp = simpson weights(Np) * ((2*pi)/(Np-1));
wt = wt.* sin theta);

G = struct();
```
```matlab
G.r = r; G theta = theta; G = phi;
G.w_r = wr; G.w_theta = wt; G.wier = wp;
% legacy aliases:
G.wr = wr; G.wrad = wr;
G.w_th = wt; G.wtheta = wt; G.wth = wt;
G.w_ph = wp; G.wphi = wp; G.wph = wp;
end

function w = simpson weights(N)
    w = ones(1,N);
    if N>2, w(2:2:end-1)=4; w(3:2:end-2)=2; end
    w([1 end]) = 1;
    w = (1/3)*w;
end
```

<caption>Listing 44: Deuteron-based energy normalization (sets C.norm_energy)</caption>

`tsrtQvalues.m` Q-value utility used throughout lifetime and benchmark pipelines.

* **Called by:** `RunLifetimeCalculations.m`, `RunLifetime_All.m`, `run_benchmark_full.m`, `runall.m`.

### `tsrtQvalues.m`

```matlab
function [Qalpha_MeV, Qalpha_src, Qbeta_minus_MeV] = tsrtQvalues(Z, A)
%TSRTQVALUES Physics-based Q-values with robust fallbacks.
% [Qalpha, src, Qbeta_minus] in MeV
%
% Primary sources (in this order):
% 1) Binding/Q tables:
%    - tsrt_bindingscan_q.csv
%    - tsrt_bindingscan_cln.csv
%    - tsrt_bindingscan.csv
% 2) Benchmark diagnostics (fallback for Qalpha/Qbeta when bindings missing):
%    - tsrt_benchmark_diag_FIXED.csv
%    - tsrt_benchmark_diag.csv
%
% If explicit Q columns exist and are finite, use them.
% Otherwise compute from bindings:
% Qalpha = B(Z-2,A-4) + B alpha - B(Z,A), B alpha=28.295674 MeV
% Qbeta_minus = 0.782343 + (B(Z+1,A) - B(Z,A))
%
% Clear diagnostics in Qalpha_src indicate the provenance.

arguments
    Z (1,1) double {mustBeInteger, mustBePositive}
    A (1,1) double {mustBeInteger, mustBePositive}
end

% -------- constants --------
B ALPHA = 28.295674; % MeV, 4He binding energy
DELTA_BETA = 0.782343; % MeV, (m_n - m_p - m_e)

% -------- candidate files --------
bind_candidates = { ...
    'tsrt_bindingscan_q.csv', ...
    'tsrt_bindingscan_cln.csv', ...
    'tsrt_bindingscan.csv' ...
```
```matlab
if isfinite(Qalpha_MeV)
    Qalpha_src = 'bench:Qalpha';
else
    Qalpha_src = 'NaN: bindings/bench unavailable';
end

else
    Qalpha_src = sprintf('csv:%s', direct_alpha);
end

% ------------- Qbeta -----------------
[Qbeta_minus_MeV, direct-beta] = local_tries_Q(rowP, {'qbeta_minus', 'qbeta', 'q-betaMinus', 'qbeta_mev', 'q-beta_mev'});
if ~isfinite(Qbeta_minus_MeV)
    % compute from bindings
    Bp = local_get_bindigion(DB.tables, Z, A);
    Bzp = local_get_bindigion(DB.tables, Z+1, A);
    if isfinite(Bp) && isfinite(Bzp)
        Qbeta_minus_MeV = DELTA_BETA + (Bzp - Bp);
    else
        % fallback to benchmark CSVs
        Qbeta_minus_MeV = local_bench_Q(bench_candidates, BENCH.tables, Z, A, 'beta');
    end
else
    % keep direct
end

% ==================== helpers ===================
function T = local_read_table(path)
    try
        T = readtable(path, 'FileType', 'text', 'Delimiter', ',');
        'TreatAsEmpty', {'NaN', 'nan', '', ' '}, 'ReadVariableNames', true,
        'PreserveVariableNames', true, 'TextType', 'string');
    catch
        T = readtable(path, 'FileType', 'text', 'Delimiter', ',');
        'TreatAsEmpty', {'NaN', 'nan', '', ' '}, 'ReadVariableNames', true,
        'PreserveVariableNames', true);
    end

vn = lower(strrep(T.Properties.VariableNames, ' ', ''));
vn = matlab.lang.makeUniqueStrings(vn);
T.Properties.VariableNames = vn;

vn = T.Properties.VariableNames;
vn = local_alias(vn, 'ebind_mev', {'ebind_mev', 'e_bind_mev', 'binding_mev', 'bind_mev', 'e_bind', 'bindingenergymev', 'bindingenergy'});
vn = local_alias(vn, 'qalpha', {'qalpha', 'q alpha', 'qalpha_mev', 'q alpha_mev'});
vn = local_alias(vn, 'qbeta_minus', {'qbeta_minus', 'qbeta', 'q beta_minus', 'qbeta_mev', 'q beta_mev'});
vn = local_alias(vn, 'z', {'z', 'Z', 'proton', 'protons'});
vn = local_alias(vn, 'a', {'a', 'A', 'massnumber', 'mass'});
vn = local_alias(vn, 'n', {'n', 'N', 'neutron', 'neutrons'});
```
```pascal
T.Properties.VariableNames = vn;

% force numeric for key columns if needed
numNames = intersect({'z', 'n', 'a', 'ebind_mev', 'qalpha', 'qbeta_minus'}, T.
Properties.VariableNames);
for i=1:numel(numNames)
    c = T.(numNames{i});
    if ~isnumeric(c), T.(numNames{i}) = str2double(string(c)); end
end
end

function v = local_alias(v, target, aliases)
    for i=1:numel(aliases)
        ix = find(str comparison, v, aliases{i}), 1);
        if ~isempty(ix)
            v{ix} = target;
        end
    end
end
end

function row = local_find_row(DBtables, Z, A)
    row = [];
    for i=1:numel(DBtables)
        T = DBtables{i}.T;
        if ~any(str comparison, T.Properties.VariableNames)) || ~any(str comparison, T.
Properties.VariableNames))
    continue
    end
    idx = find(T.z == Z & T.a == A, 1);
    if ~isempty(idx)
        row = struct('T', T, 'row', idx);
        return;
    end
    end
end
end

function [Qval, colName] = local try Q(row, candidates)
    Qval = NaN; colName = '';
    if isempty(row) || ~isstruct(row) || ~isfield(row, 'T') || ~isfield(row, 'row') ||
    isempty(row row)
        return
    end
    T = row.T; idx = row row;
    names = T.Properties.VariableNames;
    for i=1:numel(candidates)
        c = lower(candidates{i});
        if any(str comparison, names)
            v = T.(c)(idx);
            if isnumeric(v) && isfinite(v)
                Qval = v; colName = c; return;
            end
            end
        end
    end
end

function B = local_get_bind(g, Z, A)
    if Z<=0 || A<=0
        B = NaN; return;
```
```c
end
row = local_find_row(DBtables, Z, A);
if isempty(row)
    B = NaN; return;
end
T = row.T; idx = row row;
if any(strcmp('ebind_mev', T.Properties.VariableNames))
    B = T.ebind_mev(idx);
    if ~isfinite(B), B = NaN; end
else
    B = NaN;
end
end

function Q = local_bench_Q(cand_names, BENCHtables, Z, A, mode)
% Search benchmark CSVs for (Z,A,mode) and extract Q columns if present.
Q = NaN;
mode = lower(string mode);
for i=1:numel(BENCHtables)
    T = BENCHtables{i}.T;
    names = lower(string(T.Properties.VariableNames));
    % must have Z,A,mode columns
    if ~any.names=="z") || ~any.names=="a") || ~any.names=="mode"
        continue
    end
    zc = T{:,names=="z"}; if ~isnumeric(zc), zc = str2double(string(zc)); end
    ac = T{:,names=="a"}; if ~isnumeric(ac), ac = str2double(string(ac)); end
    mc = lower(string(T{:,names=="mode"}));
    msk = (zc==Z) & (ac==A) & (mc==mode);

    if any(msk)
        % look for Qbeta or Qalpha column
        if mode=="beta"
            qcols = ["qbeta_mev","qbeta_minus","q-beta_minus","qbeta");
        else
            qcols = ["qalpha_mev","qalpha","q-alpha");
        end
        for qc = qcols
            if any.names==qc)
                v = T{msk, find(names==qc,1)};
                if iscell(v), v = v{1}; end
                v = double(str2double(string(v))); % in case it's string
                if isfinite(v)
                    Q = v; return;
                end
            end
        end
    end
end
end
end
```

<caption>Listing 45: `tsrtQvalues.m`</caption>

`tsrt_alpha_L_for_transition.m` Angular-momentum-dependent α-decay barrier (L) factor calculator.

* Called by: `RunLifetimeCalculations.m`, `RunLifetime_All.m`, `run_benchmark_full.m`,
runall.m

tsrt_alpha_L_for_transition.m

```matlab
function L = tsrt_alpha_L_for_transition(Z, A)
%TSRT ALPHA L FOR TRANSITION Orbital angular momentum for alpha transition.
% Physics-informed mapping:
% - Even-even heavy emitters -> L=0 (favored)
% - 209Bi (Z=83,A=209) -> L=5 (strong hindrance; long t1/2)
% - Default -> L=0; extend as the spectroscopic assignments grow.

if Z == 83 && A == 209
    L = 5;
    return;
end

% Even-even parent heuristic: L=0
if mod(Z,2)==0 && mod(A-Z,2)==0
    L = 0;
    return;
end

% Fallback (safe default)
L = 0;
end
```

<caption>Listing 46: tsrt_alpha_L_for_transition.m</caption>

tsrt_alpha_width.m Partial α decay width producer (input to half-life).

*   Called by: RunLifetimeCalculations.m, RunLifetime_All.m, run_benchmark_full.m, runall.m

tsrt_alpha_width.m

```matlab
function [width_MeV, diag] = tsrt_alpha_width(Z, A, Qalpha_MeV, shape)
%TSRT ALPHA WIDTH Alpha-emission single-particle width with centrifugal barrier (WKB)
% width_MeV = (ħ * ν / (2π)) * P (we return diag.P and lambda needs C_alpha outside)

% Inputs:
% Z, A           : parent proton/mass numbers
% Qalpha_MeV      : Qα (MeV)
% shape           : struct with fields r0, a, r0c

% Physics:
% V(r) = V_N(r) + V_C(r) + V_l(r)
% V_N(r) -> Woods-Saxon-like pocket (simple smooth form)
% V_C(r) -> point+radius Coulomb (2*Zd * kC / r outside)
% V_l(r) -> centrifugal barrier: l(l+1) ħ^2 / (2 μ r^2)

% Notes:
% - We automatically set orbital L using tsrt_alpha_L_for_transition(Z,A).
%   Even-even heavy emitters -> L=0; 209Bi -> L=5 (hindered).
% - Robust to small Q/missing corners: returns finite P or throws with clear msg.

arguments
    Z (1,1) double {mustBeInteger, mustBePositive}
```
A (1,1) double {mustBeInteger, mustBePositive}
Qalpha_MeV (1,1) double
shape (1,1) struct
end

% -------- constants --------
hc = 197.3269804; % MeV*fm
kC = 1.439964; % MeV*fm (e^2 / (4*pi*epsilon0))
amuMeV = 931.49410242; % MeV/c^2

if ~isfield(shape, 'r0') || ~isfield(shape, 'a') || ~isfield(shape, 'r0c')
    error('tsrt_alpha_width:BadShape', 'shape must have r0, a, r0c.');
end
r0 = shape.r0;
aWS = shape.a;
r0c = shape.r0c;

Zd = Z - 2;
Ad = A - 4;
if Qalpha_MeV <= 0 || ~isfinite(Qalpha_MeV)
    error('tsrt_alpha_width:BadQ', 'Non-positive/NaN Qalpha.');
end

% Reduced mass (alpha + daughter) in MeV/c^2
mu_amu = (4*Ad) / A; % ~ reduced mass in amu
mu = mu_amu * amuMeV; % MeV/c^2

% Radii (fm)
Rn = r0 * (Ad^(1/3) + 4^(1/3)); % nuclear-potential radius
Rc = r0c * (Ad^(1/3) + 4^(1/3)); % Coulomb radius (matching)
if Rc <= 0, Rc = max(Rn, 1.0); end

% Orbital angular momentum for this transition (physics-based chooser)
L = tsrt_alpha_L_for_transition(Z, A);

% Potentials
Vn =@(r)-50./ (1 + exp((r - Rn)/aWS)); % MeV (depth ~50; smooth)
Vc =@(r) (2*Zd*kC) ./ max(r, 1e-6); % MeV (outside)
Vl =@(r) ((hc^2) * L*(L+1) ) ./ (2*mu*max(r,1e-6).^2); % MeV

Vtot =@(r) Vn(r) + Vc(r) + Vl(r);

% Turning points r1 < r2 such that V(r)=Q
rmin = max(0.1, 0.2*Rn);
rmax = max(100, 6*Rc + 50);
Q = Qalpha_MeV;

% Find inner and outer root robustly
r1 = fzero_safe(@r) Vtot(r)-Q, [rmin, max(Rn, rmin+0.5)]);
if ~isfinite(r1)
    % try a bracket scan outward
    r1 = bracket_root(Vtot, Q, rmin, Rc, 200);
end
r2 = bracket_root(Vtot, Q, Rc, rmax, 400);

if ~isfinite(r1) || ~isfinite(r2) || ~(r2>r1)
    error('tsrt_alpha_width:TurningPoints', 'Could not find valid turning points (L=%d) ', L);
end
end

% WKB action integral S = ∫ sqrt(2μ(V-Q)) dr / ħ
% Convert to unitless via hc
integrand =@(r) sqrt(2*mu*max(Vtot(r)-Q,0)) ./ hc; % 1/fm
S = integral(integrand, r1, r2, 'RelTol',1e-5, 'AbsTol',1e-8);
P = exp(-2*S);

% Assault frequency ν ~ v / (2Rwell), approximate via internal velocity
% Use classical v from Q inside pocket; cap minimum
v = sqrt(2*Q / (mu)) * hc; % fm/s * (dimensionless here up to constants)
Rwell = max(Rn, 6.0);
nu = max(v / (2*Rwell), 1e-22); % arbitrary nonzero floor to avoid 0

% Γ = ħ * ν * P (we keep as "width_MeV" via ħ in MeV*s, but for ratios only P matters)
hbar_in_MeV_s = 6.582119569e-22; % MeV*s
width_MeV = hbar_in_MeV_s * nu * P;

diag = struct('P', P, 'Spre', 1, 'lambda', NaN, 'C', struct(), ...
    'shape', shape, 'L', L, 'r1', r1, 'r2', r2, 'S', S);
end

% ------------------ helpers ------------------
function rroot = fzero_safe fun, bracket)
    try
        rroot = fzero fun, bracket);
    catch
        rroot = NaN;
    end
end

function r2 = bracket_root(Vtot, Q, rstart, rstop, N)
    x = linspace(rstart, rstop, N);
    y = arrayfun(@(r) Vtot(r)-Q, x);
    s = sign(y);
    idx = find(diff(s) ~= 0, 1, 'first');
    if isempty(idx), r2 = NaN; return; end
    a = x(idx); b = x(idx+1);
    r2 = fzero_safe(@(r) Vtot(r)-Q, [a, b]);
end

<caption>Listing 47: tsrt_alpha_width.m</caption>

tsrt_alpha_t12.m Converts widths/probabilities to half-life t1/2.

• Called by: run_benchmark_full.m, runall.m

tsrt_alpha_t12.m

function [t12_s, diag] = tsrt_alpha_t12(Z, A, Qalpha_MeV, shape, C_or_path)
%TSRT ALPHA_T12 Half-life from penetrability + preformation (TSRT).
% t1/2 = ln(2) / (C_alpha * P * S_pre)
%
% Inputs:
% Z, A : parent numbers
% Qalpha_MeV : Qα (MeV)
```
% shape : struct with r0, a, r0c
% C_or_path (opt) : either a constants struct C, or a char/string path to a .mat
%                 with fields C alpha and preform; if omitted, loads the default
%                 'tsrtconstants Calibrated.mat'
%
% Outputs:
% t12_s : half-life in seconds
% diag : diagnostics (P, Spre, lambda, barrier details)
```

```
% -------- resolve constants C --------
if margin < 5 || isempty(C_or_path)
    matPath = 'tsrtconstants Calibrated.mat';
    if exist matPath,'file') ~= 2
        error('tsrt alpha t12:NoC','Missing %s and no C provided.', matPath);
    end
    S = load matPath);
    C = unwrapConstants(S);
    diag constants path = matPath;
else
    if isstruct(C_or_path)
        C = unwrapConstants(C_or_path);
        diag constants path = '';
    elsif ischar(C_or_path) || (isstring(C_or_path) && isscalar(C_or_path))
        matPath = char(C_or_path);
        if exist matPath,'file') ~= 2
            error('tsrt alpha t12:BadPath','Constants path not found: %s', matPath);
    end
    S = load matPath);
    C = unwrapConstants(S);
    diag constants path = matPath;
else
    error('tsrt alpha t12:BadCArg','5th arg must be struct C or path string to
                                .mat.');
end
end
```

```
% -------- validate C --------
if ~isstruct(C) || ~isfield(C,'C alpha') || ~isfield(C,'preform')
% help the user understand what's in the file
flds = [];
try flds = string(fieldnames(C)); catch, end
if isempty(flds); flds = "<none>"; end
error('tsrt alpha t12:BadC', ...
    'C must have C alpha and preform. Got fields: %s', strjoin(flds,", "));
end
if ~isfinite(C.C alpha) || C.C alpha <= 0
    error('tsrt alpha t12:BadCalpha','C.C alpha must be positive finite.');
end
if ~isstruct(C.preform)
    error('tsrt alpha t12:BadPreform','C.preform must be a struct.');
end
```

```
% -------- penetrability (auto L) --------
[~, w] = tsrt alpha width(Z, A, Qalpha MeV, shape);
if ~isfield(w,'P') || ~isfinite(w.P) || w.P <= 0
    error('tsrt alpha t12:BadP','Non-positive/NaN penetrability P.');
end
```
% -------- preformation --------
S_pre = tsrt_preformation_prefactor(Z, A, C.preform);
if ~isfinite(S_pre) || S_pre <= 0
    error('tsrt_alpha_t12:BadSpre','Non-positive/NaN S_pre for (Z=%d,A=%d).', Z, A
);
end

% -------- rate and half-life --------
lambda = C.C alpha * w.P * S_pre;
t12_s = log(2) / lambda;

% -------- diagnostics --------
diag.P = w.P;
diag.Spre = S_pre;
diag lambda = lambda;
diag.C = C;
diag.shape = shape;
% pass-through useful barrier diagnostics if present
if isfield(w,'L'), diag.L = w.L; end
if isfield(w,'r1'), diag.r1 = w.r1; end
if isfield(w,'r2'), diag.r2 = w.r2; end
if isfield(w,'S'), diag.S = w.S; end
end

% ------------------ helpers ------------------

function C = unwrap constants(S)
% Accept several storage patterns:
% (a) flat fields: S.C alpha, S.preform, ...
% (b) nested under C: S.C.C alpha, S.C.preform
% (c) single top-level field that is itself a struct with needed fields.

% direct flat fields?
if isstruct(S) && isfield(S,'C alpha') && isfield(S,'preform')
    C = S; return;
end

% common nested 'C' pattern
if isstruct(S) && isfield(S, 'C') && isstruct(S.C)
    Ccand = S.C;
    if isfield(Ccand, 'C alpha') && isfield(Ccand, 'preform')
        C = Ccand; return;
    end
end

% single field that itself is a struct holding the constants
if isstruct(S)
    f = fieldnames(S);
    if
        if isfield(S,{f{1}})
            Ccand = S.(f{1});
            if isfield(Ccand, 'C alpha') && isfield(Ccand, 'preform')
                C = Ccand; return;
            end
        end
    end
end

% fallback: return S (validator will throw with helpful field list)
C = S;
end

<caption>Listing 48: tsrt_alpha_t12.m</caption>

`tsrt_get_exp Half Life.m` Fetches experimental half-lives for comparison tables.

*   Called by: `run_benchmark_full.m`, `runall.m`

`tsrt_get_exp Curves.m`

```matlab
function t12_s = tsrt_get_exp Curves(Z, A, mode)
%TSRT_GET_EXP Curves Experimental t1/2 (seconds) for (Z,A) and mode ('alpha'/'beta').
% Searches (in order):
% 1) tsrt_benchmark_diag_FIXED.csv
% 2) tsrt_benchmark_diag.csv
%
% Required columns (case-insensitive): Z, A, mode, t12_exp_seconds
% Returns NaN if not found for (Z,A,mode). Throws if no CSV exists.

if nargin < 3, mode = 'alpha'; end
mode = lower(string mode);
if mode == "beta-" || mode == "beta_minus"
    mode = "beta";
end

candidates = {'tsrt_benchmark_diag_FIXED.csv', 'tsrt_benchmark_diag.csv'};
dataPath = resolve_first_existing(candidates);
if isempty(dataPath)
    error('tsrt_get_exp Curves: No CSV', ... 'No benchmark CSV found. Expected one of: %s', strjoin(candidates, ', '));
end

T = read_table_robust(dataPath);

% canonicalize headers (case-insensitive)
H = lower(string(T.Properties.VariableNames));
colZ = find(H == "z", 1);
colA = find(H == "a", 1);
colM = find(H == "mode", 1);
colT = find(H == "t12_exp_seconds", 1);

if isempty(colZ) || isempty(colA) || isempty(colM) || isempty(colT)
    error('tsrt_get_exp Curves: Bad CSV', ... 'CSV %s must contain columns: Z, A, mode, t12_exp_seconds', dataPath);
end

Zvals = numcol(T, colZ);
Avals = numcol(T, colA);
modes = lower(string(T(:, colM)));

rows = (Zvals == Z) & (Avals == A) & (modes == mode);
idx = find rows, 1, 'first';
if isempty(idx)
    t12_s = NaN; % if not found, caller decides what to do
```
```perl
return;
end

t12_s = to_scalar_numeric(T{idx, colT});
end

% ------------------ helpers ------------------

function pathOut = resolve_first_existing/names)
    roots = {pwd, fullfileascwd, 'tsrt-nuclear'};
    for r = 1:lisel(roots)
        for i = 1:lisel/names)
            p = fullfile(roots{r}, names{i});
            if exist(p, 'file') == 2
                pathOut = p;
                return;
        end
    end
end

pathOut = '';
end

function T = read_table_robust(path)
    % Version-safe table import, preserves headers, tolerates strings/numbers.
    try
        opts = detectImportOptions(path);
        opts.VariableNamingRule = 'preserve';
        T = readtable(path, opts);
    catch
        % Fallback for older releases
        try
            T = readtable(path, 'PreserveVariableNames', true, 'TextType', 'string');
        catch
            T = readtable(path, 'PreserveVariableNames', true);
        end
    end
end

function v = numcol(T, col)
    v = T{:, col};
    if iscell(v), v = string(v); end
    if isstring(v) || ischar(v)
        v = str2double(v);
    end
    v = double(v);
end

function x = to_scalar_numeric(val)
    if iscell(val), val = val{1}; end
    if isstring(val) || ischar(val)
        x = str2double(val);
    else
        x = double(val);
    end
    if ~isscalar(x) || ~isfinite(x)
        error('tsrt_get_exphalf_life:BadCell', 'Non-scalar/NaN numeric value encountered in CSV.');
    end
end
```
end

<caption>Listing 49: `tsrt_get_expEMSHalfLife.m`</caption>

`tsrt_preflight_checks.m` Sanity checks for data presence and parameter consistency.

*   Called by: `runall.m`

`tsrt_preflight_checks.m`

```matlab
% FILE: tsrt_preflight_checks.m
% PURPOSE: Lightweight preflight. If tsrt_repository_check exists, call it.
%           Otherwise perform minimal, non-invasive checks so runall() is robust.
```

```matlab
function tsrt_preflight_checks()
%TSRT_PREFLIGHT_CHECKS Verify required functions, constants, and environment.
% Fails early with actionable messages.
```

```matlab
fprintf(' [Preflight] Starting TSRT preflight checks...\n');
```

```matlab
% --- 0) MATLAB version (info only) ---
v = ver('MATLAB');
fprintf(' [Preflight] MATLAB: %s %s\n', v.Name, v.Version);
```

```matlab
% --- 1) Required functions ---
req_funcns = {
    'tsrtconstants', ...
    'run_benchmark_full', ...
    'tsrtQvalues', ...
    'tsrtalphaaghreability', ...
    'tsrt alpha width', ...
    'tsrtfissionactionfull', ...
    'tsrt_saddle_action', ...
    'tsrtforbiddennessrank', ...
    'reproducehindranceTSRT' ...
};
```

```matlab
missing = {};
for i = 1:numel(req_funcns)
    if exist(req_funcns{i}, 'file') ~ = 2 && exist(req_funcns{i}, 'builtin') ~ = 5
        missing{i+1} = req_funcns{i};
        %#ok
    end
end
```

```matlab
if ~isempty(missing)
    fprintf(2,'[Preflight] Missing required functions:\n');
    for i=1:numel(missing), fprintf(2,' - %s.m\n', missing{i}); end
    error('TSRT:Preflight:MissingFunctions', ...
           'Required functions not found on the MATLAB path.');
end
```

```matlab
% --- 2) Constants load + calibrated merge (so checks reflect the real run) ---
try
    C = tsrtconstants();
catch ME
    fprintf(2,'[Preflight] tsrtconstants() failed: %s\n', ME.message);
    rethrow(ME);
end
```
end

if exist('tsrtconstants Calibrated.mat', 'file') == 2
    try
        S = load('tsrtconstants Calibrated.mat');
        if isstruct(S)
            f = fieldnames(S);
            for i=1:lisnumel(f), C.(f{i}) = S.(f{i}); end
            fprintf(' [Preflight] Merged tsrtconstants Calibrated.mat. \n');
        end
        catch ME
            fprintf(2, '[Preflight] Warning: could not merge calibrated MAT: %s \n', ME. message);
        end
    end

% --- 3) Sanity on must-have constant fields ---
must_fields = {'hbarc_MeVfm', 'e2_MeVfm', 'm alpha MeV', 'amu_MeV', 'c_ms'};
for i=1:lisnumel(must_fields)
    if ~isfield(C, must_fields{i})
        error('TSRT:Preflight:BadConstants', 'Missing constants field: %s', must_fields{i});
    end
end

% Critical for energy mapping
if ~isfield(C, 'norm_energy') || ~isscalar(C.norm_energy) || ~isfinite(C.norm_energy)
    || C.norm_energy <= 0
error('TSRT:Preflight:NormEnergy', ...
    ['C.norm_energy missing/invalid. Run TSRT_Calib_NormEn_Deuteron.m ' ...
     'and/or ensure tsrtconstants.m defines the calibrated value. ']);
end

% Optional fission helpers: warn if none present and no precomputed results
fissionunjpers = {'tsrt_buildVO parent_U235.m', 'tsrt_build_daughters_scission_tuned.m',
    ... 'tsrt_build_daughters_scission.m', 'tsrt_build_daughters_partition.m'};
has_helper = any(cellfun(@f) exist(f,'file') == 2, fission_helper));
if ~has_helper && exist('fission_results.mat', 'file') =~ 2
    fprintf(2, '[Preflight] Note: no fission builders found and no fission_results.mat present. \n');
    fprintf(2, 'Fission summary generation may be skipped. \n');
end

fprintf(' [Preflight] OK. \n');
end

<caption>Listing 50: tsrt_preflight_checks.m</caption>

tsrt_preformation_prefactor.m Complies α preformation prefactor per TSRT.

* Called by: run_benchmark_full.m, runall.m

tsrt_preformation_prefactor.m

function Spre = tsrt_preformation_prefactor(Z, A, C)
%TSRT_PREFERENCE_PREFERENCE Anchor-invariant TSRT preformation factor.
% Spre = exp( -[kZ*dZ + kN*dN + kNN*(dN)^2] - gamma0 ) * exp( gamma_shell * (K - K_ref) )

%

% Definitions:

% dZ = Z - refZ

% dN = (A - Z) - (refA - refZ)

% K = shellMetric(Z, N) % smooth, larger near magic numbers

% K_ref= shellMetric(refZ, refN)

%

% Requirements:

% C.preform struct with fields: refZ, refA, kZ, kN, kNN, gamma0, gamma_shell.

%

% Inputs: Z, A are scalar integers. Output Spre is dimensionless.

if ~isfield(C,'preform') || ~isstruct(C.preform)
    error('tsrt_preformation_prefactor:NoPreform','C.preform missing.');
end

pf = C.preform;
must = {'refZ','refA','kZ','kN','kNN','gamma0','gamma_shell'};
for i=1:numel(must)
    if ~isfield pf, must{i})
        error('tsrt_preformation_prefactor:BadPreform','C.preform.%s missing.',
               must{i});
    end
end

N = A - Z;
Z0 = pf.refZ; NO = pf.refA - pf.refZ;

dZ = Z - Z0;
dN = N - NO;

% base linear-quadratic trend (natural log)
lnS_base = -( pf.kZ*dZ + pf.kN*dN + pf.kNN*(dN.^2) ) - pf.gamma0;

% shell metric (anchor-invariant)
K = tsrt_shellMetric(Z, N);
K_ref = tsrt_shellMetric(Z0, NO);
lnS_shell = pf gamma ned * (K - K_ref);

Spre = exp( lnS_base + lnS_shell );
end

<caption>Listing 51: tsrt_preformation_prefactor.m</caption>

tsrt_shellMetric.m
Computes shell/structure metric used by TSRT lifetime trends.

• Called by: run_benchmark_full.m, runall.m

tsrt_shellMetric.m

function u = tsrt_shellMetric(Z, N, Zref, Nref)
%TSRT_SHELL_METRIC Smooth, positive shell metric anchored at (Zref,Nref).
% Uses the same form used in the alpha module:
% f(Z,N) = 1/(1+|Z-82|) + 1/(1+|N-126|)
% f0 = f(Zref,Nref)
% u = max(0, f - f0)

%

% This provides a smooth, monotonic indicator of proximity to major shell

% closures. It has no tunable parameters and is used only as a geometric

% lifetime-trend modifier (not an energy calibration).

f = 1/(1+abs(Z-82)) + 1/(1+abs(N-126));
f0 = 1/(1+abs(Zref-82)) + 1/(1+abs(Nref-126));
u = f - f0;
if u < 0, u = 0; end

end

<caption>Listing 52: `tsrt_shellMetric.m`</caption>

`tsrt_emit_tables.m` Assembles and writes emission/tables used in the paper.

*   **Called by:** `RunLifetime_All.m`

`tsrt_emit_tables.m`

```matlab
function tsrt_emit_tables.csv_in, tex_out, csv_out)
```

%TSRT_EMIT_TABLES Read benchmark CSV, fill TSRT predictions, and emit TeX + filled CSV.

%

% Reproducibility-first:

% - Q-values (Qalpha/Qbeta) are taken from the CSV if present/valid.

% - Optional live auto-Q lookup via tsrtQvalues(Z,A) is attempted;

% if it fails once, we disable it for the remainder of the run and

% silently fall back to the CSV values (with one-line summary).

%

% Required (present in repo):

% - tsrtconstantsCalibrated.mat (alpha/beta/sf/gamma constants)

% - tsrt alpha width m

% Optional:

% - tsrtQvalues.m (only used if ENABLE_AUTO_QVALUES=true and it works)

%

% Inputs:

% csv_in : input CSV with columns at least:

% Z,A, isotope,mode,Qbeta_MeV,Qalpha_MeV,t12_exp_seconds

% tex_out : LaTeX path to write human-friendly table

% csv_out : CSV path to write filled predictions

% --------------------- options ---------------------

ENABLE_AUTO_QVALUES = true; % try tsrtQvalues(Z,A) once; fallback to CSV if it fails

STRICT_Q = false; % if true, error when needed Q is missing (auto+CSV both invalid)

QUIET_Q = true; % one-line summary instead of per-row messages

% --------------------- load constants (if present) ---------------------

C = struct();

if exist('tsrtconstantsCalibrated.mat', 'file') == 2

C = load('tsrtconstantsCalibrated.mat');

end

% --------------------- read CSV & coerce numerics ---------------------

T = readtable.csv_in, 'FileType', 'text');
% Coerce critical numeric columns hard (no 'missing', always double)
T = coerce_numeric_columns(T, {'Z', 'A', 'Qbeta_MeV', 'Qalpha_MeV', 't12_exp_seconds', 't12_pred_s'});

% Ensure t12_pred_s exists as DOUBLE NaN column (never 'missing')

if ~ismember('t12_pred_s', T.Properties.VariableNames)
    T.t12_pred_s = nan(height(T), 1);
else
    % force to double NaN vector
    T.t12_pred_s = force_numeric_vector(T.t12_pred_s);
end

n = height(T);
t12_pred = T.t12_pred_s; if numel(t12_pred)~=n, t12_pred = nan(n,1); end

% ---------------- auto-Q control & counters ----------------
autoQ_enabled = ENABLE_AUTO_QVALUES;
autoQ_failed = false; % set to true on first failure
auto_ok = 0;
auto_fail_used = 0;
rows_total = n;

% ---------------- main evaluation loop ----------------
for i = 1:n
    try
        Z = as_scalar_num(T.Z(i));
        A = as_scalar_num(T.A(i));
        if ~(isfinite(Z) && Z>0 && isfinite(A) && A>0)
            warn_emit('row %d: invalid Z/A', i);
            continue;
        end
        N = A - Z;
    mode = as_mode(get_text_cell(T.mode, i));
    if isempty mode)
        warn_emit('row %d: invalid mode', i);
        continue;
    end
% ---------- get Q-values (CSV is authoritative; auto-Q is optional)
Qa = as_scalar_num(safe_get(T, 'Qalpha_MeV', i));
Qb = as_scalar_num(safe_get(T, 'Qbeta_MeV', i));
if autoQ_enabled
    [Qa_auto, Qb_auto, ok_auto] = safe_tsrtQvalues(Z, A);
    if ok_auto
        if ~isfinite(Qa) || Qa<=0, Qa = Qa_auto; end
        if ~isfinite(Qb) || Qb<=0, Qb = Qb_auto; end
        auto_ok = auto_ok + 1;
    else
        auto_fail_used = auto_fail_used + 1;
        autoQ_enabled = false; % disable for remainder of run
        if ~QUIET_Q
            warn_emit('[diag] tsrtQvalues unavailable; using CSV Q-values
for remaining rows.');
        end
    end
end
end

% --------------------- predict by mode ---------------------
switch mode
case 'alpha'
  if ~is_valid_alpha(C)
    warn_emit('row %d alpha: alpha constants missing in MAT file', i);
  else
    if ~isfinite(Qa) && Qa>0
      if STRICT_Q
        error('row %d alpha: required Qalpha missing (auto+CSV invalid)', i);
    else
      warn_emit('row %d alpha: missing/invalid Qalpha');
leaving NaN', i);
  end
  else
    t12_pred(i) = predict_alpha_t12(Z,A,N,Qa,C);
  end
end

% --------------------- beta ---------------------
case 'beta'
  if ~is_valid_beta(C)
    warn_emit('row %d beta: beta constants missing in MAT file', i);
  else
    if ~isfinite(Qb) && Qb>0
      if STRICT_Q
        error('row %d beta: required Qbeta missing (auto+CSV invalid)', i);
    else
      warn_emit('row %d beta: missing/invalid Qbeta');
leaving NaN', i);
  end
  else
    t12_pred(i) = predict_beta_t12(Z,A,N,Qb,C);
  end
end

% --------------------- sf ---------------------
case 'sf'
  if ~is_valid_SF(C)
    warn_emit('row %d sf: sf constants missing in MAT file', i);
  else
    t12_pred(i) = predict_SF_t12(Z,A,N,C);
  end

case 'gamma'
  % gamma rows: no prediction here (handled by dedicated calibration /tools)
  % leave NaN prediction
  otherwise
    warn_emit('row %d: unrecognized mode "%s", i, mode');
  end

catch ME
```
warn_emit('row %d failed: %s', i, ME.message);
end

end

% Write predictions back as double vector
T.t12_pred_s = t12_pred;

% --------------------- write outputs ---------------------
writetable(T, csv_out);

emit_txt_table(T, tex_out);

% --------------------- one-line summary ---------------------
if ENABLE_AUTO_QVALUES
    if auto_fail_used>0
        warn_emit('[diag] auto-Q disabled after first failure; CSV Q-values used thereafter. auto-ok=%d, rows=%d', auto_ok, rows_total);
    else
        warn_emit('[diag] auto-Q ok for all attempted rows. rows=%d', rows_total);
end
else
    warn_emit('[diag] auto-Q disabled; used CSV Q-values only. rows=%d', rows_total);
end

end

% ==================== helpers: table access / coercion ===================
function T = coerce_numeric_columns(T, names)
    for k = 1:numel names
        vname = names{k};
        if ~ismember(vname, T.Properties.VariableNames)
            % always create as double NaN column (never 'missing')
            T.(vname) = nan(height(T),1);
        else
            T.(vname) = force_numeric_vector(T.(vname));
        end
    end
end

function vec = force_numeric_vector(col)
    % Converts any reasonable column type into a column double vector.
    if isnumeric(col)
        vec = double(col(:));
        return;
    end
    n = numel(col);
    vec = nan(n,1);
    for i = 1:n
        v = NaN;
        try
            if iscell(col)
                v = col{i};
            elsif isstring(col)
                v = col(i);
            elsif iscategorical(col)
                v = char(col(i));
            elsif ischar(col)
```
```matlab
% char matrix or single char array
if size(col,1) == n
    v = col(i,:);
else
    v = col;
end
else
% fallback for table/dataset references
    v = col(i);
end
catch
    v = NaN;
end
vec(i) = as_scalar_num(v);
end
end

function v = safe_get(T, name, i)
    if ismember(name, T.Properties.VariableNames)
        col = T.(name);
        try
            if iscell(col), v = col{i};
            else,
                v = col(i);
        end
        catch
            v = NaN;
        end
    else
        v = NaN;
    end
end

% ==================== safe tsrtQvalues wrapper ===================
function [Qa,Qb,ok] = safe_tsrtQvalues(Z,A)
    Qa = NaN; Qb = NaN; ok = false;
    if ~exist('tsrtQvalues','file')
        return;
    end
    try
        Z = double(Z); A = double(A);
        [Qa1,Qb1.~] = tsrtQvalues(Z,A);
        if isfinite(Qa1) && Qa1>0, Qa = Qa1; end
        if isfinite(Qb1) && Qb1>0, Qb = Qb1; end
        ok = true;
    catch
        ok = false;
    end
end

% ==================== validation helpers ===================
function tf = is_valid Alpha(C)
    tf = isfield(C,'C_alpha') && isfield(C,'shape_best') && ...
        isfield(C,'preform') && all([
            isfield(C.preform,'refZ'), ...
            isfield(C.preform,'refA'), ...
            isfield(C.preform,'kZ'), ...
```
isfield(C.preform, 'kN'), ...
isfield(C.preform, 'kNN'), ...
isfield(C.preform, 'gamma_shell') ]];
end

function tf = is_valid Beta(C)
tf = isfield(C,'beta') && all([
isfield(C.beta,'C_beta'), ...
isfield(C.beta,'power'), ...
isfield(C.beta,'bZ'), ...
isfield(C.beta,'bN'), ...
isfield(C.beta,'eta_forb') ]);
end

function tf = is_valid_SF(C)
tf = isfield(C,'sf') && all([
isfield(C.sf,'C_SF'), ...
isfield(C.sf,'g0'), ...
isfield(C.sf,'g1') ]);
end

% ==================== TSRT predictors ===================

function t12 = predict_alpha_t12(Z,A,N,Qalpha,C)
    [~, w] = tsrt_alpha_width(Z, A, Qalpha, C.shape_best);
    if ~isfield(w, 'P') || ~isfinite(w.P) || w.P <= 0
        t12 = NaN; return;
    end
    dZ = Z - C.preform.refZ;
    dN = N - (C.preform.refA - C.preform.refZ);
    Spre_base = exp(-(C.preform.kZ*dZ + C.preform.kN*dN + C.preform.kNN*(dN.^2)));
    u = shellMetric(Z, N, C.preform.refZ, C.preform.refA - C.preform.refZ);
    Spre_shell = exp(-C.preform.gamma_shell * u);
    Spre = max(realmin, Spre_base.* Spre-shell);
    rate = C.C_alpha * w.P * Spre;
    if ~isfinite(rate) && rate>0
        t12 = NaN;
    else
        t12 = log(2) / rate;
    end
end

function t12 = predict_beta_t12(Z,A,N,Qb,C)
    if ~isfinite(Qb) || Qb<=0, Qb = eps; end
    rate = C.beta.C Beta.* (Qb.^C.beta.power) .* exp(C.beta.bZ*Z + C.beta.bN*N);
    if Z==6 && A==14
        rate = rate.* exp(C.beta.eta_forb);
    end
    t12 = log(2) ./ rate;
end

function t12 = predict_SF_t12(Z,A,N,C)
fiss = (Z^2)/A;
NO   = 152;
if isfield(C.sf,'NO'), NO = C.sf.NO; end
dn   = abs(N - NO);
rate = C.sf.C_SF.* exp(C.sf.g0 * fiss + C.sf.g1 * dn);
t12 = log(2) ./ rate;
end

% =========================== scalar/printing utils ============================

function x = as_scalar_num(v)

% robust scalar parsing; no 'missing' usage anywhere

if iscell(v)

if isempty(v), x = NaN; return; end

v = v{1};

end

if isstring(v)

ifawl(v) ~= 1, x = NaN; return; end

v = char(v);

end

if ischar(v)

s = strtrim(v);

if isempty(s) || any(str(s, {'--', '<missing>', '\N', 'nan'}) )

x = NaN; return;

end

xv = str2double(s);

if isnan(xv), x = NaN; else, x = xv; end

return;

end

if isnumeric(v)

if isempty(v) || ~isscalar(v) || ~isfinite(v), x = NaN; else, x = double(v);

end

return;

end

x = NaN;

end

function m = as_mode(v)

if iscell(v), v = v{1}; end

if isstring(v), v = char(v); end

if ~ischar(v), m = ''; return; end

vv = strtrim(lower(v));

if any(strcmp(vv, {'alpha', 'beta', 'sf', 'gamma'})), m = vv; else, m = ''; end

end

function s = safeTex(x)

% TeX-safe scalar string; no dependence on ismissing()

if margin == 0

s = '--'; return;

end

if iscell(x)

if isempty(x), s = '--'; return; end

x = x{1};

end

if isstring(x)

ifawl(x) ~= 1, s='--'; return; end

x = char(x);

end

if isnumeric(x)

if ~isscalar(x) || ~isfinite(x)

s = '--'; return;

end

ax = abs(x);

if ax >= 1e4 || (ax > 0 && ax < 1e-3)
```matlab
s = sprintf('%.3g', x);
else
    s = sprintf('%.6g', x);
end
return;
end
if ischar(x)
    s = strtrim(x);
    if isempty(s) || any(strcovpi(s, {'<missing>', '\N'}))
        s = '--'; return;
    end
    return;
end
s = '--';
end

function s = safe_num(x)
    if ~isscalar(x) || ~isnumeric(x) || ~isfinite(x)
        s = '--'; return;
    end
ax = abs(x);
if ax >= 1e4 || ax < 1e-3
    s = sprintf('%.3g', x);
else
    s = sprintf('%.6g', x);
end
end

function warn_emit(fmt, varargin)
    % Scalar-safe warning emitter
    try
nargin < 1
        warning('[emit] (no message)');
        return;
    end
msg = sprintf(fmt, varargin{:});
warning('[emit] %s', msg);
catch ME
        warning('[emit] (formatting failed): %s', ME.message);
end
end

function u = shellMetric(Z, N, Zref, Nref)
f = 1/(1+abs(Z-82)) + 1/(1+abs(N-126));
f0 = 1/(1+abs(Zref-82)) + 1/(1+abs(Nref-126));
u = f - f0;
if u < 0, u = 0; end
end
```

<caption>Listing 53: `tsrt_emit_tables.m`</caption>

## Figure Producers

`reproducehindranceTSRT.m` Recreates hindrance figures/tables with current constants.

* Calls: `tsrtconstants.m`, `tsrtsfactor.m`
oesmakefig.m Generates the odd-even staggering (OES) figure. Code: already listed in Appendix Y.7 (Listing 13 (p. 231)).

screeningmakefig.m Generates the electron-screening figure.

* Inputs: (embedded; may read constants via tsrtconstants.m)

* Calls: tsrtconstants.m

*Implementation:* See the full code in Listing 9 (p. 223). This section only references the authoritative listing to avoid duplication.

bindingplotwithexp.m Produces the binding-energy comparison against experiment.

* Inputs: exp_bindings.csv, tsrt_bindingscan_cln.csv

* Calls: (self-contained)

*Implementation:* See the full code in Listing 7 (p. 218). We intentionally do not repeat the listing here.

## Utilities and Solvers (Project-Wide)

*Notes.* These are the core TSRT helper modules needed for reconstruction of grids, fields, curvature and energy integrals, stability diagnostics, and supporting utilities. Only unique modules are retained; all duplicates appearing elsewhere in the Appendix have been removed. The order below follows the logical flow of the numerical pipeline.

tsrt_apply_grid_defaults.m Apply default grid parameters unless explicitly set elsewhere (kept only if not absorbed into tsrtmakegrid.m).

* Purpose: Provide baseline domain/spacing defaults.

tsrt_apply_grid_defaults.m

function C = tsrt_apply_grid_defaults(C)
% Ensure C has all grid knobs required by tsrtmakegrid.m

% Angular resolution
if ~isfield(C,'N_theta') || ~isfinite(C.N_theta)
    C.N_theta = 256; % safe, decent default
end

% Radial resolution and extent
if ~isfield(C,'N_r') || ~isfinite(C.N_r)
    C.N_r = 600; % increase if tsrtdeltaK2 needs finer
end

if ~isfield(C,'r_max') || ~isfinite(C.r_max)
    C.r_max = 12.0; % units must match tsrtmakegrid (fm or code units)
end

% If tsrtmakegrid uses others, add them here, e.g.:
% if ~isfield(C,'r_min') || ~isfinite(C.r_min), C.r_min = 0.0; end
% if ~isfield(C,'N_tau') || ~isfinite(C.N_tau), C.N_tau = 400; end

% Note: tsrtmakegrid computes Nphi = max(128, 2*C.N_theta),
% so we don't set C.Nier here unless tsrtmakegrid requires it.
end

<caption>Listing 54: tsrt_apply_grid_defaults.m</caption>

tsrtcompositefield.m Composite-field builder combining nucleon or fragment fields.

• Purpose: Assemble composite nuclear fields from constituent configurations.

tsrtcompositefield.m

function X = tsrtcompositefield(G, nuc, t)
% Assemble a composite trembling field by superposition.
% Inputs:
% G : grid struct (tsrtmakegrid)
% nuc : struct with fields:
% .types {1xA} cell of 'proton'/'neutron'
% .R [3xA] positions (m)
% optional (scalar or 1xA):
% .r0, .Aamp, .anis, .phi
% t : time (unused; we use time-averaged baseline)
% Output:
% X : struct with .profile (3D array) for curvature integrals

A =
X = zeros(numel(G.r),
    t = 0.0; end

for i = 1:A
    % r0
    if isfield(nuc, 'r0')
        if isscalar(nuc.r0), r0_i = nuc.r0; else, r0_i = nuc.r0(i); end
    else
        ifstrcmpi(nuc.types{i}, 'proton'), r0_i = def_r0_p; else, r0_i = def_r0_n; end
    end
    % Aamp
    if isfield(nuc, 'Aamp')
        if isscalar(nuc.Aamp), A_i = nuc.Aamp; else, A_i = nuc.Aamp(i); end
    else
        A_i = def_Aamp;
    end
    % anis
    if isfield(nuc, 'anis')
        if isscalar(nuc.anis), anis_i = nuc.anis; else, anis_i = nuc.anis(i); end
    else
        anis_i = def_anis;
    end
    % phi
    if isfield(nuc, 'phi')
        if isscalar(nuc), phi_i = nuc; else, phi_i = nucado(i); end
    else
        nuc = nuc; % This is a placeholder for the actual nuc data
    end
end
```R
phi_i = 0.0;
end

params_i = struct('type', nuc.types{i}, ...
                     'A',     A_i, ...
                     'r0',    r0_i, ...
                     'phi',   phi_i, ...
                     'anis',  anis_i);

Xi_i = tsrtnucleonfieldshifted(G, params_i, nuc.R(:,i));
XiTot = XiTot + Xi_i.profile;
end

Xi = struct('profile', XiTot);
end
```

<caption>Listing 55: tsrtcompositefield.m</caption>

## tsrt_stabilitymap_points.m

* **Purpose:** Generate point-sets for TSRT stability maps and diagnostic line producers.

### tsrt_stabilitymap_points.m

```R
function tsrt_stabilitymap_points.dat, outBase)
% TSRT_STABILITYMAP Points Plot only computed nuclides (no background fill).
% Usage: tsrt_stabilitymap_points('tsrt_stabilitymap_data.mat',
  tsrt_stabilitymap_points')

if length<1||isempty.dat), dat = 'tsrt_stabilitymap_data.mat'; end
if length<2||isempty(outBase), outBase = 'tsrt_stabilitymap_points'; end

S = load.dat);
ZN = S.ZN; % rows=Z
NN = S.NN; % cols=N
F = S.Fbase; % scalar field
mask = isfinite(F) & F>-Inf & F<Inf;

% Robust percentile stretch
v = F mask);
lo = prctile(v,2); hi = prctile(v,98);
Fv = F; Fv(Fv<lo)=lo; Fv(Fv>hi)=hi;

% Build a coordinate list for only the valid cells
Zlist = ZN mask);
Nlist = NN mask);
Flist = Fv mask);

% Plot as scatter (no artificial continuity)
f = figure('Color','w','Position', [120 120 980 560]);
ax = axes('Parent',f); hold(ax,'on'); grid(ax,'on'); box(ax,'on');

% Normalize color 0..1 for colormap
cmin = min(Flist); cmax = max(Flist);
cn = (Flist - cmin) / max eps, (cmax - cmin));

% Use small filled squares; size scales mildly with Z-population
scatter(ax, Nlist, Zlist, 14, cn, 's', 'filled');
```
```matlab
colormap(ax, parula);
cb = colorbar(ax); cb.Label.String = 'd(\DeltaK^2)/d\tauau at \tauau^* (arb. units)';
caxis(ax, [0 1]);

xlabel(ax,'Neutron number N'); ylabel(ax,'Proton number Z');
title(ax,'\DeltaK^2-based TSRT map (computed nuclides only)');
set(ax, 'YDir','normal');

print(f, [outBase '.png'], '-dpng', '-r300');
set(f,'PaperPositionMode','auto'); print(f, [outBase '.pdf'], '-dpdf');
fprintf('Saved points-only map: %s.[png|pdf]
', outBase);
end
```

<caption>Listing 56: tsrt_stabilitymap_points.m</caption>

```
tsrtdeltaK2longtau.m (Long-τ ΔK² extension)
```

* Purpose: Evaluate long-τ behaviour of ΔK² for TSRT stability diagnostics.

* Used in paper: Appendix Y.11

```
tsrtdeltaK2longtau.m
```

```matlab
function DeltaK2 = tsrtdeltaK2longtau(Z, A, C, G, tau_vec)
% TSRTDELTAK2LONGTAU Proper-time ΔK^2(τ) for nucleus (Z,A), nucleus-dependent.
% Deterministic nucleus geometry:
%   - size scaling: R ~ r0 * A^(1/3)
%   - neutron skin: R_n = R_p * (1 + beta skin * delta), delta=(N-Z)/A
%   - p-n phase lag: dphi = kappa phi * delta
% Inputs:
% Z, A : integers
% C : constants struct; the following fields are used if present:
%     r0 [m], beta skin, kappa phi, coherence eta
% G : grid from tsrtmakegrid(C). If empty, will be created.
% tau_vec : vector of proper-time samples in [0,1]; default 0:1/255:1
% Output:
% DeltaK2 : column vector same length as tau_vec

if margin < 5 || isempty(tau_vec), tau_vec = linspace(0,1,256).'; end
if margin < 4 || isempty(G), G = tsrtmakegrid(C); end

% ---- pull nucleus-wide constants with safe defaults ----
C = ensure_field(C, 'r0', 1.20e-15); % 1.20 fm
C = ensure_field(C, 'betaSkin', 0.10); % neutron-skin strength
C = ensure_field(C, 'kappa phi', pi/12); % ~15° at |delta|~0.3
C = ensure_field(C, 'coherence eta', 0.05); % p-n coherence

% ---- grid vectors & integration weights (construct from vectors) ----
r = G.r(:);
th = G.theta(:);
ph = G.phi(:);
Nr ={numel(r)}; Nth ={numel(th)}; Nph =';

dr = [diff(r); r(max(end,2)) - r(max(end-1,1))];
dth = [diff(th); th(max(end,2)) - th(max(end-1,1))];
```
dph = [diff(ph); ph(max(end,2)) - ph(max(end-1,1))];

Wr = reshape(dr(:), [Nr,1,1]);

Wth = reshape((dth(:).*sin(th(:))), [1,Nth,1]);

Wph = reshape(dph(:), [1,1,Nph]);

R2 = reshape(r.^2, [Nr,1,1]);

W = R2 .* Wr .* Wth .* Wph;

Vol = sum(W,'all');

% ---- nucleon counts ----

Np = double(Z);

Nn = double(A - Z);

% ---- nucleus geometry (size, skin, phase) ----

A_num = double(A);

delta = (Nn - Np)/A_num;

Rp = C.r0 * A_num^(1/3);

Rn = Rp * (1 + C beta skin * delta); % neutron skin grows with asymmetry

dphi = C.kappa phi * delta; % small p-n phase lag

% ---- nucleon fields ----

p_params = struct('type','proton', 'A',1.0, 'phi',0.0, 'norm','none');

n_params = struct('type','neutron', 'A',1.0, 'phi',0.0, 'norm','none');

Xi_p = tsrtnucleonfield(G, p_params);

Xi_n = tsrtnucleonfield(G, n_params);

% ---- curvature profiles (robust signatures) ----

Kp0 = getcurvatureprofile(G, Xi_p, C); % Nr x Nth x Nph numeric

Kn0 = getcurvatureprofile(G, Xi_n, C);

% ---- radial rescale to nucleus-specific radii (stretch along r only) ----

Rref = r(end);

Kp0 = radial_stretch(Kp0, r, Rref, Rp);

Kn0 = radial_stretch(Kn0, r, Rref, Rn);

% ---- isolated (uncorrelated) averages ----

Kp2iso_avg = sum((Kp0.^2).*W, 'all')/Vol;

Kn2iso_avg = sum((Kn0.^2).*W, 'all')/Vol;

K2iso_avg = Np*Kp2iso_avg + Nn*Kn2iso_avg;

% ---- time evolution (with asymmetry-driven phase lag) ----

phi_of_tau =@(t) 2*pi*t;

amp =@(ph) cos(ph);

eta = C.coherence_eta;

DeltaK2 = zeros(numel(tau_vec),1);

for it = 1:numel(tau_vec)

phi = phi_of_tau(tau_vec(it));

phi_p = phi;

phi_n = phi + dphi;

Kp = amp.*Kp0;

Kn = amp<-phi_n).*Kn0;

Kp2_avg = sum((Kp.^2).*W, 'all')/Vol;
```matlab
Kn2_avg = sum((Kn.^2).*W, 'all')/Vol;
KpKn_avg = sum((Kp.*Kn).*W, 'all')/Vol;

K2_bound_avg = Np*Kp2_avg + Nn*Kn2_avg - 2*eta*Np*Nn*KpKn_avg;
DeltaK2(it) = K2Iso_avg - K2_bound_avg;

end
end

% ==================== local helpers ===================

function C = ensure_field(C, name, val)
    if ~isfield(C, name) || isempty(C.(name)), C.(name) = val; end
end

function K = getcurvatureprofile(G, Xi, C)
% GETCURVATUREPROFILE Call tsrtcurvaturetensor robustly and return numeric.
    candS = {
        {@() tsrtcurvaturetensor(G, Xi, C)}, ...
        {@() tsrtcurvaturetensor(G, Xi.profile, C)}, ...
        {@() tsrtcurvaturetensor(Xi, G, C)}, ...
        {@() tsrtcurvaturetensor(Xi.profile, G, C)}, ...
        {@() tsrtcurvaturetensor(G, C, Xi)}, ...
        {@() tsrtcurvaturetensor(Xi, C, G)} ...
    };
    lastmsg = '';
    for k = 1:numel(candS)
        try
nout = candS{k}{1}();
        if isstruct(out) && isfield(out, 'profile')
            K = out.profile; return;
        else if isnumeric(out)
            K = out; return;
        end
        catch ME
            lastmsg = ME.message;
        end
    end
    error('getcurvatureprofile:signature', ...)
    'tsrtcurvaturetensor signature not recognized. Last error: %s', lastmsg);
end

function Krs = radial_stretch(K, r, Rref, Rnew)
% RADIAL_STRETCH Resample K(r,θ,φ) from radius Rref to Rnew along r only.
    if abs(Rnew/Rref - 1) < 1e-12
        Krs = K; return;
    end
    Nr =ENG(r); Nth = size(K,2); Nph = size(K,3);
    r_old = r;
    r_new = r * (Rnew / Rref);
    Krs = zeros(Nr, Nth, Nph, 'like', K);
    for jj = 1:Nth
        for kk = 1:Nph
            Krs(:,jj,kk) = interp1(r_old, K(:,jj,kk), r_new, 'pchip', 'extrap');
        end
    end
end
```

<caption>Listing 57: tsrtdeltaK2longtau.m</caption>
tsrt_saddle_action.m & tsrtfissionactionfull.m

* **Purpose:** Compute saddle-point and full semiclassical action integrals used in fission diagnostics.

tsrt_saddle_action.m

```matlab
function Sout = tsrt_saddle_action(Zd, Ad, Qalpha_MeV, l, C)
    %TSRT_SADDLE_ACTION Dimensionless WKB action S/hbar across the barrier.
    % Inputs:
    % Zd, Ad      daughter Z and A
    % Qalpha_MeV  Q-value for alpha decay (MeV)
    % l           orbital angular momentum
    % C           constants (tsrtconstants)
    %
    % Output:
    % Sout        \int_{\{rin\}}^{\{rout\}} dr * sqrt(2*mu*(V(r)-Q)) / (hbar*c) * c
                (dimensionless)
```

```matlab
if margin < 5 || ~isstruct(C)
    C = tsrtconstants();
end
```

```matlab
% Reduced mass energy (MeV): mu_E = m_alpha * m_d / (m_alpha + m_d)
m_d_MeV = (Ad * C.amu_MeV); % daughter mass energy (approx. in MeV)
mu_MeV = (C.m_alpha_MeV * m_d_MeV) / (C.m_alpha_MeV + m_d_MeV); % MeV
```

```matlab
% Turning points (fm) where V(r) = Q
V = @ (r) tsrt poten tial_total(r, Zd, Ad, l, C);
rin = findpelning peloint(V, Qalpha_MeV, 1.0, 20.0, +1); % inner
rout = findpelning peloint(V, Qalpha_MeV, max(rin,5), 200.0, +1); % outer
```

```matlab
% Dimensionless integrand: kappa(r) = sqrt(2*mu*(V-Q)) / (hbar*c)
kappa = @ (r) sqrt(max(0,2*mu_MeV*(V(r)-Qalpha_MeV))) / C.hbarc_MeVfm;
```

```matlab
% Integrate with smooth quadrature
Sout = safe quad(kappa, rin, rout);
```

end

```matlab
% --- helpers ---
function V = tsrt poten tial_total(r, Zd, Ad, l, C)
    % Coulomb + nuclear proximity (short range) + centrifugal (with Langer)
    % All in MeV; r in fm.
```

```matlab
% Coulomb
Za = 2; % alpha charge
V_C = (C.e2_MeVfm * Zd * Za) ./ r;
```

```matlab
% Nuclear short-range (Woods-Saxon-like surface term)
RO = C.r0.fm * (Ad^ (1/3) + 4^ (1/3));
a = C.a_surf.fm;
VO = -25.0; % MeV (depth of surface pocket, set by physics scale)
V_N = VO ./ (1 + exp((r - RO)/a));
```

```matlab
% Centrifugal (with optional Langer correction)
ell = 1 + 0.5*(C.langer ~ = 0);
mu_dum_MeV = 1.0; % handled in action; here keep the geometric 1/r^2 factor
```
```matlab
V_L = ( (ell*(ell+1)) * (C.hbarc_MeVfm^2) ) ./ (2 * (r.^2) * mu_dum_MeV);
% Note: mu cancels properly in kappa with our formulation that uses mu in S
```

```matlab
V = V_C + V_N + V_L;
end
```

```matlab
function x = findTurningPoint(V, Q, a, b, dir)
% Bisection root of V(r)-Q=0 in [a,b]; dir>0 enforces monotone search
fa = V(a)-Q; fb = V(b)-Q;
```

```matlab
if sign(fa) == sign fb)
```

```matlab
% Expand interval if needed
lo = a; hi = b;
```

```matlab
for k=1:40
```

```matlab
    if dir>0, hi = hi*1.5; else, lo = lo/1.5; end
```

```matlab
    fa = V(lo)-Q; fb = V(hi)-Q;
```

```matlab
    if sign(fa) ~= sign fb), break; end
```

```matlab
end
```

```matlab
end
```

```matlab
lo = a; hi = b;
```

```matlab
for k=1:80
```

```matlab
    mid = 0.5*(lo+hi);
```

```matlab
    fm = V(mid)-Q;
```

```matlab
    if sign.fm) == sign(V(lo)-Q), lo = mid; else, hi = mid; end
```

```matlab
    if abs(hi-lo) < 1e-6, break; end
```

```matlab
end
```

```matlab
x = 0.5*(lo+hi);
```

```matlab
end
```

```matlab
function I = safe Quad(f, a, b)
```

```matlab
% Simple adaptive Simpson
I = adaptive simpson(f, a, b, 1e-5, 10);
```

```matlab
end
```

```matlab
function I = adaptive simpson(f, a, b, tol, depth)
```

```matlab
c = (a+b)/2;
```

```matlab
fa = f(a); fb = f(b); fc = f(c);
```

```matlab
S = simpson(fa, fb, fc, a, b);
```

```matlab
I = recur(f, a, b, tol, S, fa, fb, fc, depth);
```

```matlab
end
```

```matlab
function S = recur(f, a, b, tol, S, fa, fb, fc, depth)
```

```matlab
c = (a+b)/2; d = (a+c)/2; e = (c+b)/2;
```

```matlab
fd = f(d); fe = f(e);
```

```matlab
Sleft = simpson(fa, fc, fd, a, c);
```

```matlab
Sright = simpson不对称
```

```matlab
if depth<=0 || abs(Sleft+Sright - S) < 15*tol
```

```matlab
    S = Sleft + Sright + (Sleft+Sright - S)/15;
```

```matlab
else
```

```matlab
    S = recur(f,a,c,tol/2,Sleft, fa,fc,fd, depth-1) + ...
```

```matlab
    recur(f,c,b,tol/2,Sright, fc,fb,fe, depth-1);
```

```matlab
end
```

```matlab
end
```

```matlab
function S = simpson(fa, fb, fc, a, b)
```

```matlab
S = (b-a)/6 * (fa + 4*fc + fb);
```

```matlab
end
```

<caption>Listing 58: tsrt_saddle_action.m</caption>
tsrtfissionactionfull.m

```matlab
function Sout = tsrtfissionactionfull(Zd, Ad, Qalpha_MeV, l, C)
%TSRTFISSIONACTIONFULL Dimensionless WKB action S/hbar for the TSRT barrier.
% Wrapper that mirrors tsrt_saddle_action for compatibility with existing calls.
```

```matlab
if margin < 5 || ~isstruct(C)
    C = tsrtconstants();
end
Sout = tsrt_saddle_action(Zd, Ad, Qalpha_MeV, l, C);
end
```

<caption>Listing 59: tsrtfissionactionfull.m</caption>

tsrtconvergence.m (Convergence scan)

• Purpose: Evaluate resolution-sensitivity and numerical robustness.

tsrtconvergence.m

```matlab
function tsrtconvergence()
% Grid-convergence test using a smooth spherical Gaussian profile.
% This avoids dependencies on other TSRT routines and checks numerical stability.
C = tsrtconstants('norm_energy', 1.0); % <-- this allows running before calibration

% radial resolutions to test
Nr_list = [256, 512, 1024, 2048];

Q = zeros(size(Nr_list));
for k = 1:length(Nr_list)
    C.N_r = Nr_list(k);
    G = tsrtmakegrid(C);

    % Spherical Gaussian test field (dimensionless)
    r0 = 1.0e-15; % 1 fm typical nuclear scale
    [RR,TT,PP] = ndgrid(G.r, G theta, G.phi);
    Xi.profile = exp(-(RR./r0).^2);

    % Curvature contraction (baseline: profile itself)
    K = tsrtcurvaturetensor(G, Xi, C);

    % Energy-like integral (Joules) via calibrated normalization
    Q(k) = tsrtenergyfromcurvature(G, K, C);
end

% relative changes vs previous resolution
eps = abs(diff(Q)/Q(2:end));

disp('Convergence test for radial resolution:');
for k = 2:length(Nr_list)
    fprintf('Nr = %d vs %d: rel. change = %.3e\n', ...
        Nr_list(k), Nr_list(k-1), eps(k-1));
end
end
```

<caption>Listing 60: tsrtconvergence.m</caption>
run_benchmark_full_local.m (Local benchmark wrapper)

* Purpose: Lightweight version of the global benchmark driver.

run_benchmark_full_local.m

```matlab
function results = run_benchmark_full_local(constantsUIS)
%RUN_BENCHMARK_FULL_LOCAL Compact alpha benchmark (four nuclides).
% Loads constants robustly and passes them explicitly to t12 function.

if nargin < 1 || isempty(constantsUIS)
    constantsUIS = 'tsrtconstantsCalibrated.mat';
end

fprintf(' [Preflight] MATLAB: %s\n', version);

C = tsrt_loadConstants(constantsUIS);
fprintf(' [Preflight] C alpha=%g C lambda=%g\n', C.C Alpha, C.C lambda);

% Targets
cases = [
    struct('Z',84,'A',210,'label','$^{210}$$Po', 't12_exp_s', 138.376*24*3600); ...
    struct('Z',88,'A',226,'label','$^{226}$$Ra', 't12_exp_s', 1600 * 365.25*24*3600)
; ...
    struct('Z',90,'A',232,'label','$^{232}$$Th', 't12_exp_s', 1.405e10 *
365.25*24*3600); ...
    struct('Z',92,'A',238,'label','$^{238}$$U', 't12_exp_s', 4.468e9 *
365.25*24*3600) ...
];

shape = C.shape_best;
nOK = 0; errs = nan(1,shapeOf(cases));

for i=1:shapeOf(cases)
    Z = cases(i).Z; A = cases(i).A;
    try
        [Qalpha,~,~] = tsrtQvalues(Z,A);
        opts = struct('C', C); % <-- pass constants explicitly
        [t12_s,~] = tsrt alpha t12(Z, A, Qalpha, shape, opts);
        dlog = log10(t12_s) - log10(cases(i).t12_exp_s);
        errs(i) = abs(dlog);
        fprintf(' [ok] %s Z=%3d A=%3d Qalpha=%4g MeV t12=%s s exp=%s s %log10
=%4g\n', ...
                cases(i).label, Z, A, Qalpha, fmt(t12_s), fmt(cases(i).t12_exp_s),
                dlog);
        nOK = nOK+1;
        catch ME
            fprintf(' [E] %s failed: %s\n', cases(i).label, ME.message);
        end
    end

results = struct();
if nOK > 0
    results.N = nOK;
    results.mean_abs_dlog10 = mean(errs(~isnan(errs)));
    results.median_abs_dlog10 = median(errs(~isnan(errs)));
    results.max_abs_dlog10 = max(errs(~isnan(errs)));
    results_summary = sprintf('N=%d, mean|Δlog10|=%3g, median=%3g, max=%3g',
```
```R
results.N, results.mean_abs_dlog10, results.median_abs_dlog10, results.
max_abs_dlog10);
fprintf('\n[RunLifetimeCalculations] Summary:\n%s\n', results(summary);
else
    results.N = 0; resultsire summary = 'No successful rows';
fprintf('\n[RunLifetimeCalculations] Summary:\n%s\n', results(summary);
end
end

function s = fmt(x)
    if ~isfinite(x), s='NaN'; return; end
    e = floor(log10(x));
    if e < -2 || e > 3
        s = sprintf('%.3g', x);
    else
        s = sprintf('%.3f', x);
    end
end
```

<caption>Listing 61: run_benchmark_full_local.m</caption>

lifetimetablemake.m (Optional archival table generator)

* Purpose: Build high-resolution TSRT lifetime tables; keep if referenced archivally.

lifetimetablemake.m

```R
function out = lifetimetablemake.csvPath, texOutPath)
%LIFETIMETABLEMAKE Robust table builder.
% out = lifetimetablemake.csvPath, texOutPath)
% Reads diagnostics CSV and writes a LaTeX table with columns:
% Nuclide, Mode, t_exp (y), t_pred (y), Delta log10
% - Tolerates missing/renamed columns.
% - Does not assume a monolithic 'Q' field.
%
% Returns a struct with stats and the resolved column names.
T = readtable.csvPath, 'FileType','text');

% Coerce to numeric where needed and map nonfinite to Inf for predictions:
if ~ismember('t12_pred_s', T.Properties.VariableNames)
    error('Table missing column t12_pred_s');
end

T.t12_pred_s = double(T.t12_pred_s);
bad = ~isfinite(T.t12_pred_s) | T.t12_pred_s <= 0;
T.t12_pred_s(bad) = Inf;

% (Optional) experimental column sanitization if used:
if ismember('t12_exp_seconds', T.Properties.VariableNames)
    T.t12_exp_seconds = double(T.t12_exp_seconds);
end

out = struct();

nucc = pickcol(T, {'nuclide', 'isotope', 'name', 'nuclideasc', 'nuclide_label'});
Zc = pickcol(T, {'Z', 'Z', 'Zp', 'ZP'});
Ac = pickcol(T, {'A', 'A_parent', 'Ap', 'AP'});
modec = pickcol(T, {'mode', 'decay_mode'});
texp = pickcol(T, {'t12_exp_seconds', 't12_exp_s', 't_exp_s'});
```
tpred= pickcol(T, {'t12_pred_s', 't_pred_s'});

if isempty(nucc), error('No nuclide/isotope column found.'); end
if isempty(modec), modec = 'mode'; T.mode = repmat("?", height(T), 1); end

% Resolve numeric arrays for years and deltas (never NaN where avoidable)
N = height(T);
texpY = nan(N,1);
tpredY = nan(N,1);
dlog = nan(N,1);
Nuclide = strings(N,1);
Mode = strings(N,1);

for i=1:N
    Nuclide(i) = string(T.(nucc)(i));
    Mode(i) = string(T.(modec)(i));
    teY = NaN; tpY = NaN;
    if ~isempty(texp) && ismember(texp, T.Properties.VariableNames)
        v = T.(texp)(i);
        if isfinite(v) && v>0, teY = v / (3600*24*365.2425); end
    end
    if ~isempty.tpred) && ismember.tpred, T.Properties.VariableNames)
        v = T.(tpred)(i);
        if isfinite(v) && v>0, tpY = v / (3600*24*365.2425); end
    end
    texpY(i) = teY;
    tpredY(i) = tpY;
    if isfinite(teY) && teY>0 && isfinite.tpY) && tpY>0
        dlog(i) = log10.tpY - log10(teY);
    end
end

% Write LaTeX
fid = fopen.texOutPath, 'w');
fprintf fid, '%% Auto-generated high-resolution lifetime benchmark table\n');
fprintf fid, '\\begin{table}[ht]\n');
fprintf fid, '\\centering\n');
fprintf fid, '\\begin{tabular}{l l r r}\n');
fprintf fid, '\\hline\n');
fprintf fid, 'Nuclide & Mode & $t_{\{\mathrm{exp}}}$$ (y) & $t_{\{1/2\}}^{\{\mathrm{TSRT}}}$$ (y)$ & $\Delta\log_{10}$$ \\\ \n');
fprintf fid, '\\hline\n');
for i=1:N
    fprintf fid, '%s & %s & %s & %s & %s \\\\ \n', ...
        strclean(Nuclide(i)), strclean(Mode(i)), fmtYears(texpY(i)), fmtYears.tpredY(i)), fmtDelta(dlog(i)));
end
fprintf fid, '\\hline\n');
fprintf fid, '\\end{tabular}\n');
fprintf fid, '\\end{table}\n');
fclose fid;

% Stats
valid = isfinite(dlog);
vals = dlog(valid);
out.stats.All.count = sum(valid);
out.stats.All.mean = mean vals, 'omitnan';
```matlab
out.stats.All.rms = sqrt(mean vals.^2,'omitnan'));
out.columns = struct('nuclide',nucc,'Z',Zc,'A',Ac,'mode',modec,...
't12_exp_seconds',temp,'t12_pred_s',tpred);
end

function c = pickcol(T, names)
c = '';
for k=1:numel names
n = names{k};
if any strcomp T.Properties.VariableNames,n)
c = n; return;
end
for j=1:numel T.Properties.VariableNames)
if contains lower T.Properties.VariableNames{j}), lower(n))
c = T.Properties.VariableNames{j}; return;
end
end
end

function s = strclean(x)
s = char string x);
s = strrep s, '\n',' ');
s = regexprep s, '\s+',' ');
end

function s = fmtYears(y)
if ~ (isfinite y) && y>0)
s = '--';
else
if y<1e-3
s = sprintf ('%.3g', y);
else if y<1
s = sprintf ('%.4f', y);
else if y<100
s = sprintf ('%.3f', y);
else if y<1e6
s = sprintf ('%.2f', y);
else
s = sprintf ('%.3g', y);
end
end
end

function s = fmtDelta(d)
if ~isfinite(d), s='--'; else, s=sprintf ('%+.3f', d); end
end
```

<caption>Listing 62: lifetimetablemake.m</caption>

magicnumber_CALC.m (Magic-number closures: raw and corrected)

* **Purpose:** Reproduce the TSRT magic-number sequence by (i) computing the raw cumulative occupancies from the zeros of $j'_l(x)$ (Table 20 (p. 116)) and (ii) generating the corrected closure mapping that incorporates surface/geometry and spin-geodesic adjustments (Table 21 (p. 118)).
* **Inputs:** None required; reads shared calibration (surface/geometry + spin-geodesic flags) from `tsrtconstants.m`. Optional in-code settings: `lmax`, `nmax`.

* **Outputs:** CSV files `magic_raw.csv` and `magic_corrected.csv`; a diagnostic plot comparing raw vs. corrected cumulative occupancies.

* **Dependencies:** `tsrtconstants.m`, `tsrt_magicCorrection_order.m` (permutation for the corrected ordering), MATLAB's `besselj`. Uses a lightweight helper for locating zeros of $j'_l(x)$.

### `magicnumber_CALC.m`

```matlab
function magicnumber_CALC(varargin)
% MAGICNUMBER_CALC Reproduce TSRT magic-number tables:
% 1) Raw spherical locking: roots of j_l'(x)=0 -> Table~\refp{tab:magicmodes}}
% 2) Corrected closures (surface/geometry + spin--geodesic) -> Table~\refp{tab:magic_map}}
%
% Conventions match Appendix~\ref{app:magic_numbers} and Appendix~\ref{app:nuc-methods-param-calib}.
% Reads shared calibration from tsrtconstants.m (C, C_spin_geod, etc.)

% --- settings ---
lmax = 6;      % angular degree ceiling (covers Table~\refp{tab:magicmodes}}
nmax = 2;      % first two radial families are sufficient here
Nper = 50;     % max roots to scan internally (safety)
outRaw = 'magic_raw.csv';
outCorr = 'magic_corrected.csv';

% --- load shared calibration / constants (provides tiny corrections) ---
C = tsrtconstants(); %#ok<NASGU> % must be on path; holds surface/geometry + spin--geodesic flags

% --- compute derivative roots of spherical Bessel j_l'(x)=0 ---
rows = []; % [n l xnl gnl]
for l = 0:lmax
    % get roots of j_l'(x)=0; you can use a robust helper (e.g., '
    sphbessel_deriv zeros')
    xs = sphbessel_deriv zeros(1, nmax, Nper); % returns first nmax roots for given l
    for n = 1:min(nmax, numel(xs))
        xnl = xs(n);
        gnl = 2*(2*1+1);      % TSRT degeneracy (angular mult. * trembling-phase duality)
        rows = [rows; n 1 xnl gnl]; %#ok<AGROW>
    end
end

% --- sort by increasing x_{nl} (stiffness ordering) ---
rows = sortrows rows, 3); % columns: n, l, xnl, gnl

% --- cumulative raw occupation ---
Nc = cumsum rows(:,4));

% --- write raw table (for Table~\refp{tab:magicmodes}} ---
Traw = table rows(:,1), rows(:,2), rows(:,3), rows(:,4), Nc, ...
    'VariableNames', {'n', 'l', 'xnl', 'gnl', 'Nc_raw'});
    writetable(Traw, outRaw);

% --- apply tiny TSRT corrections that reorder near neighbors (shared calibration) ---
```
```matlab
% The correction is implemented in a helper that uses C.
% The correction is implemented in a helper that uses C surface geom, C spin geod, etc
idx_corr = tsrt_magicCorrection_order rows; % returns permutation index for
corrected ordering
rows_corr = rows(idx_corr, :);
Nc_corr = cumsum rows_corr (:, 4));

% --- write corrected sequence (drives Table~\ref{tab:magic_map}) ---
Tcorr = table rows_corr(:,1), rows_corr(:,2), rows_corr(:,3), rows_corr(:,4), Nc_corr,
...
'VariableNames', {'n', 'l', 'xnl', 'gnl', 'Nc_corr'});
writetable(Tcorr, outCorr);

% --- optional: produce a quick diagnostic plot (raw vs corrected cumulative) ---
figure;
plot(1:numel(Nc), Nc, 'o-', 'hold on');
plot(1:numel(Nc_corr), Nc_corr, '.-');
legend('raw cumulative', 'corrected cumulative', 'Location', 'NorthWest');
xlabel('mode-family rank'); ylabel('cumulative occupancy');
title('TSRT magic-number closures: raw vs corrected');

end

% --- helper: zeros of j_l'(x)=0 (simple wrapper; replace with your preferred routine) ---
function xs = sphbessel_deriv zeros(1, nNeeded, Ncap)
    xs = zeros(1,nNeeded);
    % crude bracket+refine strategy using fzero; initialize brackets near (l+0.5 + pi*
m)
    m = 1; found = 0; attempts = 0;
    while found < nNeeded && attempts < Ncap
        % seed guess; spacing ~ pi for large arguments
        xL = (l + 0.5) + pi*(m-1);
        xR = (l + 0.5) + pi*(m);
        % define j_l'(x) via recurrence (or call a library routine if available)
        f =@(x) dj_sph_bessel(l, x);
        % sign check and root refine if a sign change is detected
        if sign(f(xL)) ~= sign(f(xR))
            try
                xs[found+1] = fzero(f, [xL xR]);
                found = found + 1;
            catch
                % skip if refine fails
            end
        end
        m = m + 1; attempts = attempts + 1;
    end
    if found < nNeeded
        error('sphbessel_deriv zeros: could not find %d roots for l=%d', nNeeded, l);
    end
end

% derivative of spherical Bessel j_l via relation with cylindrical Bessel
function val = dj_sph_bessel(l, x)
    % j_l(x) = sqrt(pi/(2x)) J_{l+1/2}(x)
    % derivative can be built from J_{nu} and its derivative, or via stable
    recurrences
    % For brevity, use numerical derivative on j_l with small h (sufficient for root
```
```matlab
finding)
    h = 1e-6;
    val = (sph_bessel_j(1, x+h) - sph_bessel_j(1, x-h))/(2*h);
end

function val = sph_bessel_j(1, x)
    % stable evaluation using relation to cylindrical Bessel J_{nu}
    nu = 1 + 0.5;
    val = sqrt(pi/(2*x)) .* besselj(nu, x);
end
```

Listing 63: Minimal driver to reproduce Table 20 (p. 116) (raw) and the corrected closure mapping (Table 21 (p. 118)). It computes zeros of $j'_l(x)$, builds degeneracies and cumulative counts, then applies the shared TSRT surface/geometry + spin-geodesic correction read from tsrtconstants.m.

## Data Files and Reproducibility Resources

This list consolidates data resources (inputs and generated artifacts) required by the producers in this appendix.

* `tsrtconstants_calibrated.mat`

**Role:** Persisted calibration constants (energy scale, scission width, E2 anchor) used by all pipelines.

**Produced by:** TSRT_Calib_NormEn_Deuteron.m,

TSRT_Calib_Fission_ScissWth_U235.m, TSRT_Calibrate_GammaE2_156Gd.m.

**Consumed by:** runall.m and all per-module drivers.

* `exp_bindings.csv`

**Role:** Experimental binding energies for the binding comparison figure.

**Used by:** bindingplotwithexp.m.

* `tsrt_bindingscan_cln.csv`

**Role:** Cleaned TSRT binding-scan data used in the binding comparison plot.

**Used by:** bindingplotwithexp.m.

* `tsrt_DeltaK2_map.csv`

**Role:** Stability diagnostic (if used) for $\Delta K^2$ maps.

**Used by:** stability diagnostics (Appendix Y.3), if invoked.

* `tsrt_benchmark_diag.csv`, `tsrt_benchmark_diag_FILLED.csv`,

`tsrt_benchmark_diag_FIXED.csv`

**Role:** Benchmark diagnostics (internal consistency and regression checks).

**Used by:** run_benchmark_full.m/runall.m diagnostics (optional).

* `tsrt_bindingscan.csv`, `tsrt_bindingscan_q.csv`

**Role:** Raw and auxiliary TSRT binding scans (pre-cleaning).

**Used by:** bindingplotwithexp.m (via `tsrt_bindingscan_cln.csv`) as provenance.

* `fission_results.mat`

**Role:** Numerical fission summary (struct array F).

**Produced by:** make_fission_summary_now.m.

**Consumed by:** fissiontablemake.m to emit tab_fission_summary.tex.
* tab_fission_summary.tex

**Role:** Auto-generated LATEXtabular for the 235U(n,f) summary.

**Produced by:** fissiontablemake.m.

**Included as:** Table 42 (p. 257) (Appendix Y.15).

For interpretation of these computational results in the physical framework of TSRT, see the main sections on binding and fission energetics (Sections 5–6.2).
References

[1] Richard S. Westfall. *Never at Rest: A Biography of Isaac Newton*. Cambridge University Press, 1980.

[2] Isaac Newton. *Philosophiæ Naturalis Principia Mathematica*. Joseph Streater, London, 1687.

[3] Roger Highfield and Paul Carter. *The Private Lives of Albert Einstein*. Faber & Faber, 1993.

[4] Abraham Pais. *Subtle is the Lord: The Science and the Life of Albert Einstein*. Oxford University Press, 1982.

[5] Albert Einstein. Ist die trägheit eines körpers von seinem energieinhalt abhängig? *Annalen der Physik*, 323(13):639–641, 1905. doi:10.1002/andp.19053231314.

[6] Ruth Lewin Sime. *Lise Meitner: A Life in Physics*. University of California Press, Berkeley, 1996.

[7] Martin J. Sherwin. *A World Destroyed: Hiroshima and Its Legacies*. Stanford University Press, Stanford, California, 3rd edition, 2003.

[8] Alvin M. Weinberg. *The First Nuclear Era: The Life and Times of a Technological Fixer*. American Institute of Physics, 1994.

[9] ITER Organization. *Iter: The world's largest fusion experiment*. https://www.iter.org, 2023. Accessed: 2023-11-01.

[10] Jr. V. E. Viola and Glenn T. Seaborg. Nuclear systematics of the heavy elements—ii lifetimes for alpha, beta and spontaneous fission decay. *Journal of Inorganic and Nuclear Chemistry*, 28(4):741–761, apr 1966. doi:10.1016/0022-1902(66)80412-8.

[11] P. Möller, A. J. Sierk, T. Ichikawa, H. Sagawa, et al. Nuclear ground-state masses and deformations: FRDM(2012). *At. Data Nucl. Data Tables*, 109-110:1–204, 2016. doi:10.1016/j.adt.2015.10.002.

[12] Nico F. Declercq. A unified geometric theory of trembling spacetime relativity. *Zenodo*, 2025. doi:10.5281/zenodo.15360946.

[13] Nico F. Declercq. Corpuscular light in trembling spacetime: A geometric foundation for planck's radiation law and the limits of temperature. *Zenodo*, 2025. doi:10.5281/zenodo.15586864.

[14] Nico F. Declercq. A geometric theory of atomic structure without quantum postulates. *Zenodo*, 2025. doi:10.5281/zenodo.15682930.

[15] Nico F. Declercq. Causal entanglement from deterministic geodesic correlation in trembling spacetime. *Zenodo*, 2025. doi:10.5281/zenodo.15490040.

[16] Nico F. Declercq. Emergence of 4d spacetime and cosmic expansion from causal geometry. *Zenodo*, 2025. doi:10.5281/zenodo.16656746.

[17] Nico F. Declercq. The geometric origin and finite classification of fundamental particles in trembling spacetime. *Zenodo*, 2025. doi:10.5281/zenodo.15830692.

[18] Nico F. Declercq. Beyond superposition and collapse: Double-slit interference from trembling spacetime geodesics. *Zenodo*, 2025. doi:10.5281/zenodo.15425854.
[19] N. Schunck and L. M. Robledo. Microscopic theory of nuclear fission: a review. *Rep. Prog. Phys.*, 79(11):116301, 2016. doi:10.1088/0034-4885/79/11/116301.

[20] S. Goriely, Y. Sirois, A. Van Lauwe, and S. Hilaire. Systematic study of the odd-even staggering of fission fragments. *Eur. Phys. J. A*, 58(5):87, 2022. doi:10.1140/epja/s10050-022-00736-4.

[21] E. G. Adelberger et al. Solar fusion. II. Proton-proton chain and cno cycles. *Reviews of Modern Physics*, 83(1):195-246, Feb 2011. doi:10.1103/RevModPhys.83.195.

[22] Niels Bohr and Fritz Kalckar. On the transmutation of atomic nuclei by impact of material particles. i. general theoretical remarks. *Mat. Fys. Medd. K. Dan. Vidensk. Selsk.*, 14(10), 1936.

[23] Niels Bohr and John Archibald Wheeler. The mechanism of nuclear fission. *Phys. Rev.*, 56:426-450, 1939. doi:10.1103/PhysRev.56.426.

[24] Maria Goeppert-Mayer. On closed shells in nuclei. *Phys. Rev.*, 75(12):1969-1970, 1949. doi:10.1103/PhysRev.75.1969.

[25] Otto Haxel, J. Hans D. Jensen, and Hans E. Suess. On the "magic numbers" in nuclear structure. *Phys. Rev.*, 75:1766, 1949. doi:10.1103/PhysRev.75.1766.

[26] Michael Bender, Paul-Henri Heenen, and Paul-Gerhard Reinhard. Self-consistent mean-field models for nuclear structure. *Rev. Mod. Phys.*, 75(1):121-180, 2003. doi:10.1103/RevModPhys.75.121.

[27] D. Vretenar, A. V. Afanasjev, G. A. Lalazissis, and P. Ring. Relativistic hartree-bogoliubov theory: static and dynamic aspects of exotic nuclear structure. *Phys. Rep.*, 409(3-4):101-259, 2005. doi:10.1016/j.physrep.2004.10.001.

[28] R. J. Furnstahl and K. Hebeler. New applications of renormalization group methods in nuclear physics. *Rep. Prog. Phys.*, 76(12):126301, 2013. doi:10.1088/0034-4885/76/12/126301.

[29] H.-W. Hammer, Sebastian König, and U. van Kolck. Nuclear effective field theory: status and perspectives. *Rev. Mod. Phys.*, 92(2):025004, 2020. doi:10.1103/RevModPhys.92.025004.

[30] M. Dasgupta, D. J. Hinde, N. Rowley, and A. M. Stefanini. Beyond the coherent coupled channels description of nuclear fusion. *Phys. Rev. Lett.*, 82(7):1395-1398, 1999. doi:10.1103/PhysRevLett.82.1395.

[31] C. L. Jiang et al. Systematic examination of the fusion hindrance in light and medium-mass systems. *Phys. Rev. C*, 99(5):054603, 2019. doi:10.1103/PhysRevC.99.054603.

[32] C. Schmitt et al. Scission configuration and fission-fragment properties in the brownian shape-motion model. *Nucl. Phys. A*, 945:1-20, 2016. doi:10.1016/j.nuclphysa.2015.09.011.

[33] F. Raiola et al. Enhanced electron screening in d(d,p)t for deuterated metals. *Phys. Lett. B*, 547(3-4):193-198, 2002. doi:10.1016/S0370-2693(02)02772-6.

[34] A. Tumino et al. An increase in the ¹²C + ¹²C fusion rate from resonances at astrophysical energies. *Nature*, 600(7889):250-254, 2021. doi:10.1038/s41586-021-04133-7.

[35] Nico F. Declercq. Tsrt nuclear data and reproducibility bundle, 2025. doi:10.5281/zenodo.17637889.
[36] Meng Wang, W.J. Huang, F.G. Kondev, G. Audi, and S. Naimi. The AME 2020 atomic mass evaluation (ii). tables, graphs and references. *Chinese Physics C*, 45(3):030003, 2021. doi:10.1088/1674-1137/abddaf.

[37] C. F. von Weizsäcker. Zur theorie der kernmassen. *Zeitschrift für Physik*, 96(7-8):431–458, 1935. doi:10.1007/BF01337700.

[38] A. Bohr and B. R. Mottelson. Collective and individual-particle aspects of nuclear structure. *Kongelige Danske Videnskabernes Selskab, Matematisk-Fysiske Meddelelser*, 27(16):1–174, 1953.

[39] V. M. Strutinsky. Shell effects in nuclear masses and deformation energies. *Nucl. Phys. A*, 95(2):420–442, 1967. doi:10.1016/0375-9474(67)90510-6.

[40] V. M. Strutinsky. "shells" in deformed nuclei. *Nucl. Phys. A*, 122(1):1–33, 1968. doi:10.1016/0375-9474(68)90699-4.

[41] John Bardeen, Leon N. Cooper, and J. Robert Schrieffer. Theory of superconductivity. *Physical Review*, 108(5):1175-1204, dec 1957. doi:10.1103/PhysRev.108.1175.

[42] Leon N. Cooper. Bound electron pairs in a degenerate fermi gas. *Physical Review*, 104(4):1189-1190, nov 1956. doi:10.1103/PhysRev.104.1189.

[43] J. Robert Schrieffer. *Theory of Superconductivity*. Benjamin, New York, 1964.

[44] H. J. Lipkin. Validity of many-body approximation methods for a solvable model. *Annals of Physics*, 9(2):272-291, 1960. doi:10.1016/0003-4916(60)90039-0.

[45] Y. Nogami. Improved bcs method and its application to a model of superconductivity. *Physical Review*, 134(2B):B313-B317, apr 1964. doi:10.1103/PhysRev.134.B313.

[46] W. D. Myers and W. J. Swiatecki. Average nuclear properties. *Annals of Physics*, 55(3):395-505, 1969. doi:10.1016/0003-4916(69)90202-4.

[47] R. D. Woods and D. S. Saxon. Diffuse surface optical model for nucleon-nuclei scattering. *Physical Review*, 95(2):577-578, 1954. doi:10.1103/PhysRev.95.577.

[48] H. J. Specht. Nuclear fission. *Reviews of Modern Physics*, 46(4):773-787, 1974. doi:10.1103/RevModPhys.46.773.

[49] C. W. Reich. Historical note: How the double-humped fission barrier was discovered. *Nuclear Physics News*, 24(4):25-30, 2014. doi:10.1080/10619127.2014.964221.

[50] Arnold Sommerfeld. *Atombau und Spektrallinien, Band 2: Wellenmechanik*. Friedrich Vieweg und Sohn, Braunschweig, 1939. This is the definitive German edition where the Coulomb wave functions and the parameter are thoroughly treated.

[51] G. Gamow. Zur quantentheorie des atomkernes. *Z. Phys.*, 51(3):204-212, 1928. doi:10.1007/BF01343196.

[52] C. D. Ellis and W. A. Wooster. The average energy of disintegration of radium E. *Proc. R. Soc. Lond. A*, 117(778):109-123, 1927. doi:10.1098/rspa.1927.0172.

[53] Wolfgang Pauli. Offener brief an die gruppe der radioaktiven bei der gau-vereinigung zu tübingen. *Die Naturwissenschaften*, 19(37):741-743, 1931. doi:10.1007/BF01521530.

[54] F. Reines and C. L. Cowan, Jr. The neutrino. *Nature*, 178:446-449, 1956. doi:10.1038/178446a0.
[55] Paul A. M. Dirac. The quantum theory of the emission and absorption of radiation. *Proceedings of the Royal Society A*, 114(767):243-265, 1927.

[56] Enrico Fermi. *Nuclear Physics*. University of Chicago Press, Chicago, 1950. A course of lectures by Enrico Fermi. Revised edition. See particularly the discussion of transition probabilities.

[57] K. Langanke and G. Martínez-Pinedo. Nuclear weak-interaction processes in stars. *Rev. Mod. Phys.*, 75(3):819-862, 2003. doi:10.1103/RevModPhys.75.819.

[58] J. Metzinger and R. Klapisch. Systematics of beta-strength functions. *Phys. Rev. C*, 38(4):3817-3829, 1988. doi:10.1103/PhysRevC.38.3817.

[59] B. T. Cleveland et al. Measurement of the solar electron neutrino flux with the Homestake chlorine detector. *Astrophys. J.*, 496(1):505-526, 1998. doi:10.1086/305343.

[60] Y. Fukuda and others (Super-Kamiokande Collaboration). Evidence for oscillation of atmospheric neutrinos. *Physical Review Letters*, 81:1562-1567, 1998. Definitive evidence for neutrino oscillations. doi:10.1103/PhysRevLett.81.1562.

[61] H. Yukawa. On the interaction of elementary particles. i. *Proc. Phys.-Math. Soc. Jpn.*, 17:48-57, 1935. URL: https://www.jstage.jst.go.jp/article/ppms1912/17/0/17_0_48/_pdf.

[62] David J. Gross and Frank Wilczek. Ultraviolet behavior of non-abelian gauge theories. *Physical Review Letters*, 30(26):1343-1346, 1973. arXiv:https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.30.1343, doi:10.1103/PhysRevLett.30.1343.

[63] Sheldon L. Glashow. Partial symmetries of weak interactions. *Nuclear Physics*, 22:579-588, 1961. doi:10.1016/0029-5582(61)90469-2.

[64] Steven Weinberg. The quantum theory of fields: Foundations and applications. *Rev. Mod. Phys.*, 68(1):1-9, 1996. This is a summary of his textbooks, which are the definitive treatise on the subject. doi:10.1103/RevModPhys.68.1.

[65] Sven G. Nilsson. Binding states of individual nucleons in strongly deformed nuclei. *Mat. Fys. Medd. Dan. Vid. Selsk.*, 29(16):1-69, 1955.

[66] John D. Huxley. *Nuclear Reaction Rates in Hot Plasmas*. Adam Hilger, Bristol, 1980.

[67] William A. Fowler, Georgeanne R. Caughlan, and Barbara A. Zimmerman. Thermonuclear reaction rates, II. *Annual Review of Astronomy and Astrophysics*, 13:69-112, 1975. doi:10.1146/annurev.aa.13.090175.000441.

[68] E. Fermi. Possible production of elements of atomic number higher than 92. *Nature*, 133(3372):898-899, 06 1934. doi:10.1038/133898a0.

[69] Otto Hahn and Fritz Strassmann. Über den nachweis und das verhalten der bei der bestrahlung des urans mittels neutronen entstehenden erdalkalimetalle. *Naturwissenschaften*, 27(1):11-15, 1939. doi:10.1007/BF01488241.

[70] O. R. Frisch. Physical evidence for the division of heavy nuclei under neutron bombardment. *Nature*, 143(3616):276-276, 1939. doi:10.1038/143276a0.

[71] H. L. Anderson. Assisting fermi at the dawn of the nuclear age. *Physics Today*, 25(7):40-47, 1972. doi:10.1063/1.3070928.
[72] M. Arnould, S. Goriely, and K. Takahashi. The r-process of stellar nucleosynthesis: Astro-physics and nuclear physics achievements and mysteries. *Physics Reports*, 450(4):97-213, 2007. doi:10.1016/j.physrep.2007.06.002.

[73] A. N. Andreyev, K. Nishio, and K.-H. Schmidt. Nuclear fission: a review of experimental advances and phenomenology. *Rep. Prog. Phys.*, 81(1):016301, 2018. URL: https://dx.doi.org/10.1088/1361-6633/aa82eb, doi:10.1088/1361-6633/aa82eb.

[74] U. Brosa, S. Grossmann, and A. Müller. Nuclear scission. *Physics Reports*, 197(4):167-262, 1990. doi:10.1016/0370-1573(90)90114-H.

[75] Lise Meitner and Otto R. Frisch. Disintegration of uranium by neutrons: a new type of nuclear reaction. *Nature*, 143:239-240, 1939. doi:10.1038/143239a0.

[76] Frédéric Joliot. Preuve expérimentale de la rupture explosive des noyaux d'uranium et de thorium sous l'action des neutrons. *Comptes Rendus Hebdomadaires des Séances de l'Académie des Sciences*, 208:341-343, 1939.

[77] H. L. Anderson, E. T. Booth, J. R. Dunning, E. Fermi, G. N. Glasoe, and F. G. Slack. The fission of uranium. *Phys. Rev.*, 55(5):511-512, 1939. doi:10.1103/PhysRev.55.511.2.

[78] R. M. Spector. Photofission of the heavy elements. *Phys. Rev.*, 140(4B):B1034-B1042, 1965. doi:10.1103/PhysRev.140.B1034.

[79] J. R. Nix. Calculation of fission barriers for heavy and superheavy nuclei. *Annu. Rev. Nucl. Sci.*, 22:65-120, 1972. doi:10.1146/annurev.ns.22.120172.000433.

[80] Enrico Fermi. The development of the first chain reacting pile. *Proc. Am. Philos. Soc.*, 90(1):20-24, 1946. URL: https://www.jstor.org/stable/3301034.

[81] J. F. Berger, M. Girod, and D. Gogny. Time-dependent quantum collective dynamics applied to nuclear fission. *Nuclear Physics A*, 428:23-36, 1984. doi:10.1016/0375-9474(84)90240-9.

[82] C. Golabek and C. Simenel. Induced fission of $^{240}$pu within a real-time mean-field approach. *Physical Review Letters*, 103(4):042701, 2009. doi:10.1103/PhysRevLett.103.042701.

[83] Robley D. Evans. *The Atomic Nucleus*. McGraw-Hill, New York, 1955.

[84] J. Rayford Nix. Theory of fission nucleus dynamics and collective motion. *Phys. Rev. C*, 5(4):1262-1265, 1972. doi:10.1103/PhysRevC.5.1262.

[85] T. R. England and B. F. Rider. Evaluation and compilation of fission product yields. LA LA-UR-94-3106, Los Alamos National Laboratory, 1994. URL: https://www-nds.iaea.org/publications/indc/nds-213/.

[86] OECD Nuclear Energy Agency. Endf/b-viii.0: U-235 isotope. Evaluated nuclear data file, National Nuclear Data Center, Brookhaven National Laboratory, 2018. Accession Number: 13939. URL: https://www.nndc.bnl.gov/endf/b8.0/.

[87] Niels Bohr. Neutron capture and nuclear constitution. *Nature*, 137:344-348, 1936. doi:10.1038/137344a0.

[88] R. d'E. Atkinson and F. G. Houtermans. Zur frage der aufbaumöglichkeit der elemente in sternen. *Z. Phys.*, 54(9):656-665, 1929. doi:10.1007/BF01341595.

[89] Arthur S. Eddington. The internal constitution of the stars. *The Observatory*, 43:341-358, 1920. URL: https://ui.adsabs.harvard.edu/abs/19200obs....43..341E.
[90] H. A. Bethe. Energy production in stars. *Physical Review*, 55(5):434-456, Mar 1939.
doi:10.1103/PhysRev.55.434.

[91] J. Chadwick. The existence of a neutron. *Proc. R. Soc. Lond. A*, 136(830):692-708, 1932.
doi:10.1098/rspa.1932.0112.

[92] E. U. Condon. Quantum mechanics and radioactive disintegration. *Phys. Rev.*, 31(6):891-894, 1928. doi:10.1103/PhysRev.31.891.

[93] R. W. Gurney and E. U. Condon. Quantum mechanics and radioactive disintegration. *Phys. Rev.*, 33(2):127-140, 1929. doi:10.1103/PhysRev.33.127.

[94] Max Born. Zur quantenmechanik der stoßvorgänge. *Zeitschrift für Physik*, 37(12):863-867, 1926. doi:10.1007/BF01397477.

[95] Georgeanne R. Caughlan and William A. Fowler. Thermonuclear reaction rates v. *At. Data Nucl. Data Tables*, 40(2):283-334, 1988. doi:10.1016/0092-640X(88)90009-5.

[96] Niels Bohr. The quantum postulate and the recent development of atomic theory. *Nature*, 121:580-590, 1928.

[97] David J. Griffiths. *Introduction to Quantum Mechanics*. Prentice Hall, Upper Saddle River, NJ, 1995.

[98] John M. Blatt and Victor F. Weisskopf. *Theoretical Nuclear Physics*. John Wiley & Sons, New York, 1952.

[99] H. A. Bethe. Nuclear physics B. nuclear dynamics, theoretical. *Rev. Mod. Phys.*, 9(2):69-244, 1937. doi:10.1103/RevModPhys.9.69.

[100] V. F. Weisskopf. Statistics and nuclear reactions. *Phys. Rev.*, 52(4):295-303, 1937.
doi:10.1103/PhysRev.52.295.

[101] I. J. Thompson. Coupled-channel methods for nuclear physics. *Comput. Phys. Rep.*, 7(4):167-212, 1988. doi:10.1016/0167-7977(88)90005-6.

[102] Meng Wang et al. The AME 2020 atomic mass evaluation (i). evaluation of input data, and adjustment procedures. *Chin. Phys. C*, 45(3):030003, 2021.
doi:10.1088/1674-1137/abddb0.

[103] Particle Data Group. Review of particle physics. *Progress of Theoretical and Experimental Physics*, 2022:083C01, 2022. doi:10.1093/ptep/ptac097.

[104] H. A. Bethe and C. L. Critchfield. The formation of deuterons by proton combination. *Phys. Rev.*, 54(10):862-862, 1938. doi:10.1103/PhysRev.54.862.2.

[105] Donald D. Clayton. *Principles of Stellar Evolution and Nucleosynthesis*. University of Chicago Press, Chicago, 1983.

[106] E. Margaret Burbidge, G. R. Burbidge, William A. Fowler, and F. Hoyle. Synthesis of the elements in stars. *Reviews of Modern Physics*, 29(4):547-650, 1957.
doi:10.1103/RevModPhys.29.547.

[107] Ken Freeman and J. Bland-Hawthorn. The new galaxy: Signatures of its formation. *Annual Review of Astronomy and Astrophysics*, 40:487-537, 2002.
doi:10.1146/annurev.astro.40.060401.093840.
[108] Q. R. Ahmad et al. Direct evidence for neutrino flavor transformation from neutral-current interactions in the sudbury neutrino observatory. Phys. Rev. Lett., 89(1):011301, 2002. doi:10.1103/PhysRevLett.89.011301.

[109] S. Fukuda et al. Solar $^8$B and hep neutrino measurements from 1258 days of superkamiokande data. Phys. Rev. Lett., 86(25):5651-5655, 2001. doi:10.1103/PhysRevLett.86.5651.

[110] Philip Morrison, Phylis Morrison, the Office of Charles, and Ray Eames. Powers of Ten: A Book About the Relative Size of Things in the Universe and the Effect of Adding Another Zero. Scientific American Library, New York, 1982. A seminal work based on the famous film, exploring the universe across 40 orders of magnitude.

[111] Gerard 't Hooft. In het kielzog van de deeltjes. Contact, Amsterdam, 1984. Translated into English as "In Search of the Ultimate Building Blocks" by Cambridge University Press (1997). This is the original Dutch edition the author encountered.

[112] F. Raiola et al. Enhanced electron screening in d(d,p)t for deuterated metals. Eur. Phys. J. A, 13(3):377-382, 2002. doi:10.1140/epja/i2001-10252-6.

[113] H. J. Assenbaum, K. Langanke, and C. Rolfs. Effects of electron screening on low-energy fusion cross sections. Z. Phys. A, 327(4):461-468, 1987. doi:10.1007/BF01289572.

[114] Henri Becquerel. Sur les radiations émises par phosphorescence. Comptes Rendus de l'Académie des Sciences, 122:420-421, 1896.

[115] R. W. Gurney and E. U. Condon. Quantum mechanics and radioactive disintegration. Nature, 122(3073):439, 1928. doi:10.1038/122439a0.

[116] J. Hans D. Jensen and Maria Goeppert Mayer. Elementary theory of nuclear shell structure. Phys. Rev., 98(5):1476-1476, 1955. doi:10.1103/PhysRev.98.1476.

[117] C. A. Bertulani and M. S. Hussein. Tunneling of the quantum decay of the compact cold $^9$Be nucleus in the $^{197}$Au coulomb field. Physical Review C, 103(6):064612, 2021. A modern perspective discussing the limitations of the standard tunneling picture. doi:10.1103/PhysRevC.103.064612.

[118] H. A. Bethe and R. F. Bacher. Nuclear physics A. stationary states of nuclei. Reviews of Modern Physics, 8(2):82-229, 1936. doi:10.1103/RevModPhys.8.82.

[119] A. E. S. Green. The liquid drop model and nuclear moments. Physical Review, 99(3):1016-1026, 1955. doi:10.1103/PhysRev.99.1016.

[120] W. D. Myers and W. J. Swiatecki. Nuclear masses and deformations. Nuclear Physics, 81(1):1-60, 1966. doi:10.1016/0029-5582(66)90639-0.

[121] Hans E. Suess. On the interpretation of the magic numbers in the structure of atomic nuclei. Die Naturwissenschaften, 36(6):213-214, 1949. doi:10.1007/BF00590147.

[122] Eugene Feenberg and K. C. Hammack. The shell model of the nucleus. Reviews of Modern Physics, 27(1):1-17, 1955. doi:10.1103/RevModPhys.27.1.

[123] J. Bartel, P. Quentin, M. Brack, C. Guet, and H. B. Håkansson. Towards a better parametrisation of Skyrme-like effective forces: A critical study of the SkM force. Nucl. Phys. A, 386(1):79-100, 1982. doi:10.1016/0375-9474(82)90403-1.
[124] G. Hagen, T. Papenbrock, M. Hjorth-Jensen, and D. J. Dean. Coupled-cluster computations of atomic nuclei. *Rep. Prog. Phys.*, 77(9):096302, 2014. doi:10.1088/0034-4885/77/9/096302.

[125] B. R. Barrett, P. Navrátil, and J. P. Vary. Ab initio no core shell model. *Prog. Part. Nucl. Phys.*, 69:131–181, 2013. doi:10.1016/j.ppnp.2012.10.003.

[126] Christian Iliadis. *Nuclear Physics of Stars*. Wiley-VCH, 2007.

[127] C. L. Jiang et al. Unexpected behavior of heavy-ion fusion cross sections at extreme sub-barrier energies. *Phys. Rev. Lett.*, 89:052701, Jul 2002. URL: https://link.aps.org/doi/10.1103/PhysRevLett.89.052701, doi:10.1103/PhysRevLett.89.052701.

[128] C. L. Jiang et al. Systematic examination of the fusion hindrance in heavy element synthesis. *Eur. Phys. J. A*, 55(7):116, 2019. doi:10.1140/epja/i2019-12798-x.

[129] F. Raiola et al. Enhanced electron screening in d(d,p)t for deuterated metals. *Eur. Phys. J. A*, 13(3):377–382, 2002. doi:10.1140/epja1335-0.

[130] J. P. McCarten, R. G. Helmer, C. W. Reich, D. R. Tilley, and H. K. Carter. Nuclear data sheets for A=156. *Nucl. Data Sheets*, 109(1):1–196, 2008. ENSDF evaluation; ¹⁵⁶Gd 2⁺! → 0⁺ E2 standard. doi:10.1016/j.nds.2008.01.001.

[131] T. Kib'edi, T. W. Burrows, M. B. Trzhaskovskaya, P. M. Davidson, and Jr. C. W. Nestor. Evaluation of theoretical conversion coefficients using bricc. *Nucl. Instrum. Methods Phys. Res., Sect. A*, 589(2):202–229, 2008. BrIcc v2.3 ICC; used for α$_{ICC}$ confirmation. doi:10.1016/j.nima.2008.02.051.
Index of MATLAB Files
Index of MATLAB Files

MATLAB

bindingplotwithexp.m, 218
csctmaketables.m, 254
fission_path.m, 274
fissiontablemake.m, 271
gaussleg.m, 248
lifetimetablemake.m, 333
magicnumber_CALC.m, 336
make_fission_summary_now.m, 257
oesmakefig.m, 231
reproducehindranceTSRT.m, 229
run_benchmark_full.m, 295
run_benchmark_full_local.m, 332
runall.m, 284
RunLifetime_All.m, 238, 286
RunLifetimeCalculations.m, 286
RunLifetimeCalculations_beta.m, 290
RunLifetimeCalculations_SF.m, 293
running the parameter scan driver examples, 283
scan_parameters.m, 278
screeningmakefig.m, 223
simpsonweights.m, 249
stability_lines_makefig.m, 243
stabilitymap_deltaK2_strict.m, 241
tsrt_alpha_L_for_transition.m, 304
tsrt_alpha_t12.m, 306
tsrt_alpha_width.m, 304
tsrt_apply_grid_defaults.m, 323
tsrt_build_daughters_partition.m, 270
tsrt_build_daughters_scission_tuned.m, 268

tsrt_build_parent_U235.m, 266
TSRT_Calib_Fission_ScissWth_U235.m, 263
TSRT_Calib_NormEn_Deuteron.m, 297
TSRT_Calibrate_GammaE2_156Gd.m, 235
tsrt_emit_tables.m, 314
tsrt_forbiddenness_hindrance.m, 239
tsrtgamma_rate.m, 234
tsrt_get_exphalf_life.m, 309
tsrt_preflight_checks.m, 311
tsrt_preformation_prefactor.m, 312
tsrt_saddle_action.m, 329
tsrtshell metric.m, 313
tsrt_stability_lines.m, 244
tsrt_stabilitymap_points.m, 325
tsrtbetamatrixelement.m, 238
tsrtbindingenergy.m, 221
tsrtcompositefiel.m, 324
tsrtcompositefield.m, 252
tsrtconstants.m, 213
tsrtconvergence.m, 217, 331
tsrtcurvaturetensor.m, 215
tsrtdeltaK2longtau.m, 326
tsrtenergyfromcurvature.m, 216
tsrtfissionactionfull.m, 331
tsrtfissionenergy.m, 270
tsrtfusionenergy.m, 230
tsrtmakegrid.m, 213
tsrtnucleonfield.m, 250
tsrtnucleonfieldshifted.m, 251
tsrtQvalues.m, 299
tsrtsfactor.m, 225