function results = CAS_MASTER_RUN()
% CAS_MASTER_RUN
% Master integration harness: runs CAS_F0000_VERIFY ... CAS_F0031_VERIFY sequentially
% in the BASE workspace (so each block's clear; clc; is self-contained),
% then performs cross-block constant consistency checks.
%
% F0000 = geometric scaling lemmas (axiomatic root)
% F0001–F0031 = physics corollaries
%
% Usage:
%   results = CAS_MASTER_RUN();
%
% Expected location:
%   Place this file in the same folder as CAS_F0000_VERIFY.m ... CAS_F0031_VERIFY.m

cd(fileparts(mfilename('fullpath')));
clc;
fprintf('=== VMS CAS MASTER RUN ===\n');
tStart = tic;

blockFiles = {
  'CAS_F0000_VERIFY.m'
  'CAS_F0001_VERIFY.m'
  'CAS_F0002_VERIFY.m'
  'CAS_F0003_VERIFY.m'
  'CAS_F0004_VERIFY.m'
  'CAS_F0005_VERIFY.m'
  'CAS_F0006_VERIFY.m'
  'CAS_F0007_VERIFY.m'
  'CAS_F0008_VERIFY.m'
  'CAS_F0009_VERIFY.m'
  'CAS_F0010_VERIFY.m'
  'CAS_F0011_VERIFY.m'
  'CAS_F0012_VERIFY.m'
  'CAS_F0013_VERIFY.m'
  'CAS_F0014_VERIFY.m'
  'CAS_F0015_VERIFY.m'
  'CAS_F0016_VERIFY.m'
  'CAS_F0017_VERIFY.m'
  'CAS_F0018_VERIFY.m'
  'CAS_F0019_VERIFY.m'
  'CAS_F0020_VERIFY.m'
  'CAS_F0021_VERIFY.m'
  'CAS_F0022_VERIFY.m'
  'CAS_F0023_VERIFY.m'
  'CAS_F0024_VERIFY.m'
  'CAS_F0025_VERIFY.m'
  'CAS_F0026_VERIFY.m'
  'CAS_F0027_VERIFY.m'
  'CAS_F0028_VERIFY.m'
  'CAS_F0029_VERIFY.m'
  'CAS_F0030_VERIFY.m'
  'CAS_F0031_VERIFY.m'
};

nBlocks = numel(blockFiles);

% Run results
results = struct();
results.blocks = repmat(struct( ...
  'name','', ...
  'status','', ...
  'err','', ...
  'constants',struct() ...
), nBlocks, 1);

% Constant collection across blocks
allConst = struct(); % canonicalName -> struct('values',[],'sources',{{}})
canonList = {'c','eps0','mu0','h','hbar','kB','NA','sigma_sb','G','e','me'};

for i = 1:nBlocks
  f = blockFiles{i};
  results.blocks(i).name = f;

  fprintf('[%02d/%02d] %s ... ', i, nBlocks, f);
  try
    % Run in BASE so the block's clear; clc; doesn't nuke this function workspace.
    evalin('base', sprintf("run('%s');", f));
    results.blocks(i).status = 'PASS';
    fprintf('PASS\n');
  catch ME
    results.blocks(i).status = 'FAIL';
    results.blocks(i).err = getReport(ME, 'basic', 'hyperlinks', 'off');
    fprintf('FAIL\n');
  end

  % Collect constants from BASE after the block ran (even if it failed mid-way,
  % some constants may still exist, but we only use them for diagnostics).
  consts = collect_constants_from_base();
  results.blocks(i).constants = consts;

  % Aggregate constants
  for k = 1:numel(canonList)
    cn = canonList{k};
    if isfield(consts, cn) && ~isempty(consts.(cn)) && isnumeric(consts.(cn))
      if ~isfield(allConst, cn)
        allConst.(cn) = struct('values', [], 'sources', {{}} );
      end
      allConst.(cn).values(end+1,1) = double(consts.(cn)); %#ok<AGROW>
      allConst.(cn).sources{end+1,1} = f; %#ok<AGROW>
    end
  end
end

% -------------------------
% Global consistency checks
% -------------------------
fprintf('\n=== GLOBAL CONSISTENCY CHECKS ===\n');

globalChecks = {};
globalPass = true;

% Helper: compare a constant across blocks
tolRelDefault = 1e-6;
tolRelTight   = 1e-9;

% Per-constant tolerances (keep conservative; many blocks use rounded CODATA)
tolMap = containers.Map();
tolMap('c')        = tolRelDefault;
tolMap('eps0')     = tolRelDefault;
tolMap('mu0')      = tolRelDefault;
tolMap('h')        = tolRelDefault;
tolMap('hbar')     = tolRelDefault;
tolMap('kB')       = tolRelDefault;
tolMap('NA')       = tolRelDefault;
tolMap('sigma_sb') = tolRelDefault;
tolMap('G')        = 1e-4;  % often rounded
tolMap('e')        = tolRelDefault;
tolMap('me')       = 1e-6;

for k = 1:numel(canonList)
  cn = canonList{k};
  if isfield(allConst, cn)
    vals = allConst.(cn).values;
    srcs = allConst.(cn).sources;
    if numel(vals) >= 2
      ref = median(vals);
      relErr = max(abs(vals - ref) ./ max(abs(ref), eps));
      tol = tolRelDefault;
      if tolMap.isKey(cn), tol = tolMap(cn); end
      pass = relErr <= tol;
      globalPass = globalPass && pass;

      fprintf('%-10s : %s  (n=%d, max rel dev=%.3g, tol=%.3g)\n', cn, tern(pass,'PASS','FAIL'), numel(vals), relErr, tol);
      if ~pass
        fprintf('  Values:\n');
        for j = 1:numel(vals)
          fprintf('    %-22s  %.15g\n', srcs{j}, vals(j));
        end
      end
    else
      fprintf('%-10s : SKIP (found in %d block)\n', cn, numel(vals));
    end
  else
    fprintf('%-10s : SKIP (not found)\n', cn);
  end
end

% Identity checks per-block where possible
% 1) c^2 * eps0 * mu0 == 1
id1Pass = true;
for i = 1:nBlocks
  cst = results.blocks(i).constants;
  if isfield(cst,'c') && isfield(cst,'eps0') && isfield(cst,'mu0') && ...
     ~isempty(cst.c) && ~isempty(cst.eps0) && ~isempty(cst.mu0)
    lhs = double(cst.c)^2 * double(cst.eps0) * double(cst.mu0);
    rel = abs(lhs - 1) / max(1, abs(1));
    pass = rel <= tolRelDefault;
    id1Pass = id1Pass && pass;
  end
end
fprintf('c^2*eps0*mu0=1 : %s\n', tern(id1Pass,'PASS','FAIL'));
globalPass = globalPass && id1Pass;

% 2) hbar == h/(2*pi)
id2Pass = true;
for i = 1:nBlocks
  cst = results.blocks(i).constants;
  if isfield(cst,'h') && isfield(cst,'hbar') && ~isempty(cst.h) && ~isempty(cst.hbar)
    lhs = double(cst.hbar);
    rhs = double(cst.h)/(2*pi);
    rel = abs(lhs - rhs) / max(abs(rhs), eps);
    pass = rel <= tolRelDefault;
    id2Pass = id2Pass && pass;
  end
end
fprintf('hbar=h/(2pi)   : %s\n', tern(id2Pass,'PASS','FAIL'));
globalPass = globalPass && id2Pass;

% 3) sigma_sb formula (optional): sigma = (pi^2 * kB^4) / (60 * hbar^3 * c^2)
id3Pass = true;
id3Any  = false;
for i = 1:nBlocks
  cst = results.blocks(i).constants;
  need = {'sigma_sb','kB','hbar','c'};
  ok = true;
  for j = 1:numel(need)
    if ~isfield(cst, need{j}) || isempty(cst.(need{j}))
      ok = false; break;
    end
  end
  if ok
    id3Any = true;
    sig = double(cst.sigma_sb);
    kB  = double(cst.kB);
    hb  = double(cst.hbar);
    c0  = double(cst.c);
    sig_formula = (pi^2 * kB^4) / (60 * hb^3 * c0^2);
    rel = abs(sig - sig_formula) / max(abs(sig_formula), eps);
    pass = rel <= 5e-4; % allow looser: many scripts use rounded constants
    id3Pass = id3Pass && pass;
  end
end
if id3Any
  fprintf('sigma formula  : %s\n', tern(id3Pass,'PASS','FAIL'));
  globalPass = globalPass && id3Pass;
else
  fprintf('sigma formula  : SKIP (sigma_sb/kB/hbar/c not co-present)\n');
end

% -------------------------
% Final summary
% -------------------------
nPass = sum(strcmp({results.blocks.status}, 'PASS'));
nFail = nBlocks - nPass;

fprintf('\n=== SUMMARY ===\n');
fprintf('Blocks: %d PASS, %d FAIL\n', nPass, nFail);
fprintf('Global: %s\n', tern(globalPass,'PASS','FAIL'));
fprintf('Elapsed: %.2fs\n', toc(tStart));

if nFail ~= 0 || ~globalPass
  fprintf('\nFailures detected. See results.blocks(i).err for details.\n');
else
  fprintf('\nSTACK CLOSED: 32 green lights (F0000–F0031) + global closure.\n');
end

% ── Write JSON audit report ──
report = struct();
report.audit = 'VMS CAS verification';
report.engine = 'MATLAB Symbolic Math Toolbox';
report.timestamp = char(datetime('now','TimeZone','UTC','Format','yyyy-MM-dd''T''HH:mm:ss''Z'''));
report.matlab_version = version;
report.blocks_run = nBlocks;
report.blocks_passed = nPass;
report.blocks_failed = nFail;
report.global_consistency = tern(globalPass, 'PASS', 'FAIL');
report.elapsed_seconds = round(toc(tStart), 2);

% Build results map
for i = 1:nBlocks
  fname = strrep(results.blocks(i).name, '.m', '');
  fname = strrep(fname, 'CAS_', '');
  fname = strrep(fname, '_VERIFY', '');
  report.results.(fname) = results.blocks(i).status;
end

% Constants registered
cnames = fieldnames(allConst);
report.constants_registered = cnames;

reportJson = jsonencode(report, 'PrettyPrint', true);
fid = fopen('audit_report.json', 'w');
fprintf(fid, '%s', reportJson);
fclose(fid);
fprintf('\nReport written to: audit_report.json\n');

end % function CAS_MASTER_RUN

% -------------------------
% Local helpers
% -------------------------
function consts = collect_constants_from_base()
% Pull known constants from BASE workspace using multiple candidate names.

consts = struct();

% Query names in base
W = evalin('base','whos');
names = {W.name};

% Canonical -> candidates (ordered)
candidates = struct();
candidates.c        = {'c_val','c0','c','c_light','c_computed','speed_of_light'};
candidates.eps0     = {'eps0_val','eps0','epsilon0','e0'};
candidates.mu0      = {'mu0_val','mu0','mu_0','u0'};
candidates.h        = {'h_val','h','planck_h'};
candidates.hbar     = {'hbar_val','hbar','h_bar','h_bar_val','hbar0'};
candidates.kB       = {'kB_val','kB','kb','k_B'};
candidates.NA       = {'NA_val','NA','N_A','avogadro','NAv'};
candidates.sigma_sb = {'sigma_sb_val','sigma_sb','sigma','sigmaSB','sigma_num'};
candidates.G        = {'G_val','G','gravG'};
candidates.e        = {'e_val','e','q_e','qe'};
candidates.me       = {'m_e_val','me_val','m_e','me','electron_mass'};

canon = fieldnames(candidates);
for i = 1:numel(canon)
  cn = canon{i};
  cand = candidates.(cn);
  val = [];
  for j = 1:numel(cand)
    nm = cand{j};
    if any(strcmp(names, nm))
      try
        v = evalin('base', nm);
        if isnumeric(v) && isscalar(v)
          val = double(v);
          break;
        end
      catch
        % ignore
      end
    end
  end
  if ~isempty(val)
    consts.(cn) = val;
  end
end
end

function out = tern(cond, a, b)
if cond, out = a; else, out = b; end
end