#!/usr/bin/env python3
"""
2-Photon Momentum & Geometric Plot Analysis
Produces two figures:
- Figure 1: master_transition_diagnostics.png (3 panels: A, B, D)
- Figure 2: normalized_consistency_check.png (2 panels: barycenter vs n, residuals)

Copyright Malcolm J. Macleod (malcolm@theprogrammergod.com)

4. Geometrical origins of quantization in H atom electron transitions
https://www.doi.org/10.2139/ssrn.3703266

 * This software is free to use, modify, and distribute under creative commons
 * https://creativecommons.org/licenses/by-nc-sa/4.0/ (CC CC BY-NC-SA 4.0).
 * Any derivative works or redistributions must include appropriate credit to
 * Malcolm Macleod and this permission notice shall
 * be included in all copies or substantial portions of the Software.

"""

import numpy as np
import sys
from scipy.interpolate import interp1d, UnivariateSpline
import matplotlib.pyplot as plt


# -------------------------
# Settings
# -------------------------
data_file = "simulation-data-m66-n7.txt"
max_nshell = 7
total_mass = 66.0

print("\n" + "=" * 80)
print("GEOMETRIC FRAMEWORK ANALYSIS")
print("=" * 80 + "\n")

# -------------------------
# Physical Constants
# -------------------------
c = 299792458  # m/s
lambda_e = 2.42631023538e-12
lambda_p = 1.32140985360e-15
H = 4.0 * np.pi * c / (lambda_e + lambda_p)

# -------------------------
# Load Simulation Data
# -------------------------
try:
    data = np.loadtxt(data_file)
    if data.ndim == 1:
        data = data.reshape(1, -1)
    if data.shape[1] < 9:
        raise ValueError(f"Expected at least 9 columns, found {data.shape[1]}")
    n_values = data[:, 0].astype(float)
    r_radius = data[:, 1].astype(float)
    xe = data[:, 2].astype(float)
    ye = data[:, 3].astype(float)
    xp = data[:, 4].astype(float)
    yp = data[:, 5].astype(float)
    angle_deg = data[:, 6].astype(float)
    path_length = data[:, 7].astype(float)
    step_counter = data[:, 8].astype(float)
    if len(n_values) < 3:
        raise ValueError("Need at least 3 saved samples")
    if np.any(step_counter <= 0):
        raise ValueError("step_counter must be strictly positive")
    print(f"Loaded {len(n_values)} data points")
    print(f"n range: {n_values[0]:.3f} -> {n_values[-1]:.3f}")
except Exception as e:
    print(f"ERROR: Cannot load {data_file}. Reason: {e}")
    sys.exit(1)

# -------------------------
# Experimental Frequency Data
# -------------------------
fh_exp = {  
    2: 2466061413187035.0,      # Parthey et al., 2011
    3: 2922743278665790.0,      # Grinin et al., 2020
    4: 3082581563822629.0,      # NIST (estimated)
    5: 3156563118433117.0,      # NIST (estimated)
    6: 3196751090001330.0,      # NIST (estimated)
    7: 3220983339500000.0,      # NIST (estimated)
}
fh_ionised = 3288086857127600.0  # Hz

# -------------------------
# Exact transition anchors from raw data (CODATA-consistent reduction)
# -------------------------
raw_crossings = {
    1: {"n_sq": 1.0,  "step": 973196},
    2: {"n_sq": 4.0000002380,  "step": 3892791},
    3: {"n_sq": 8.9999952028,  "step": 8758769},
    4: {"n_sq": 15.9999867061, "step": 15571134},
    5: {"n_sq": 24.9999745134, "step": 24329887},
    6: {"n_sq": 35.9999552306, "step": 35035025},
    7: {"n_sq": 48.9999287512, "step": 47686548},
}

missing_levels = [k for k in range(2, max_nshell + 1) if k not in raw_crossings or k not in fh_exp]
if missing_levels:
    raise ValueError(f"Missing anchor or experimental data for levels: {missing_levels}")

# -------------------------
# Frequency calibration
# -------------------------
lam = 2.0 * total_mass**2 / (total_mass - 1.0)**2
freq = (n_values**2 - 1.0) * H * lam / step_counter
fh_anchor = {k: (v["n_sq"] - 1.0) * H * lam / v["step"] for k, v in raw_crossings.items()}

# -------------------------
# Core observables
# -------------------------
m_electron = 1.0
m_nucleus = total_mass - 1.0
xb = (m_electron * xe + m_nucleus * xp) / total_mass
yb = (m_electron * ye + m_nucleus * yp) / total_mass
barycenter_distance = np.sqrt(xb**2 + yb**2)

xe_rel, ye_rel = xe - xb, ye - yb
xp_rel, yp_rel = xp - xb, yp - yb
vx_e = np.gradient(xe_rel)
vy_e = np.gradient(ye_rel)
vx_p = np.gradient(xp_rel)
vy_p = np.gradient(yp_rel)
Lz_electron = xe_rel * vy_e - ye_rel * vx_e
Lz_nucleus = xp_rel * vy_p - yp_rel * vx_p
Lz_total = Lz_electron + Lz_nucleus

bc_interp = interp1d(n_values, barycenter_distance, kind='cubic', bounds_error=False, fill_value='extrapolate')

# -------------------------
# Compute collapse curves (for Figure 1B)
# -------------------------
saved_crossing_info = {}
n_sq_series = n_values**2
for k, anchor in raw_crossings.items():
    target_n_sq = anchor["n_sq"]
    cross_idx = np.where((n_sq_series[:-1] - target_n_sq) * (n_sq_series[1:] - target_n_sq) <= 0)[0]
    if len(cross_idx) == 0:
        saved_crossing_info[k] = None
        continue
    i = cross_idx[0]
    n0, n1 = n_sq_series[i], n_sq_series[i+1]
    s0, s1 = step_counter[i], step_counter[i+1]
    if n1 == n0:
        saved_crossing_info[k] = None
        continue
    frac = (target_n_sq - n0) / (n1 - n0)
    step_interp = s0 + frac * (s1 - s0)
    f_anchor = fh_anchor[k]
    f_saved_interp = (target_n_sq - 1.0) * H * lam / step_interp
    f0 = (target_n_sq - 1.0) * H * lam / s0
    f1 = (target_n_sq - 1.0) * H * lam / s1
    f_low, f_high = min(f0, f1), max(f0, f1)
    saved_crossing_info[k] = {
        "f_exact": f_anchor,
        "f_saved_interp": f_saved_interp,
        "f_low": f_low,
        "f_high": f_high,
        "step_interp": step_interp,
    }


collapse_curves = {}
for k in sorted(raw_crossings):
    if k == 1:
        continue
    info = saved_crossing_info.get(k)
    if info is None:
        continue
    f0 = info["f_exact"]
    if k == 2:
        left_scale = fh_anchor[3] - fh_anchor[2]
    else:
        left_scale = fh_anchor[k] - fh_anchor[k-1]
    if k == max_nshell:
        right_scale = fh_anchor[k] - fh_anchor[k-1]
    else:
        right_scale = fh_anchor[k+1] - fh_anchor[k]
    half_width = 0.20 * min(left_scale, right_scale)
    mask = (freq >= f0 - half_width) & (freq <= f0 + half_width)
    if np.count_nonzero(mask) < 20:
        continue
    f_local = freq[mask]
    L_local = Lz_total[mask]
    order = np.argsort(f_local)
    f_local = f_local[order]
    L_local = L_local[order]
    x_norm = (f_local - f0) / half_width
    edge_n = max(5, len(L_local)//10)
    y_left = np.median(L_local[:edge_n])
    y_right = np.median(L_local[-edge_n:])
    amp = y_right - y_left
    if abs(amp) < 1e-12:
        continue
    y_norm = (L_local - y_left) / amp
    collapse_curves[k] = (x_norm, y_norm)

# -------------------------
# Residuals for models (Figure 2C)
# -------------------------
def ppm_residual_dict(model_dict, exp_dict):
    return {n: 1e6 * (model_dict[n] - exp_dict[n]) / exp_dict[n] for n in sorted(exp_dict) if n in model_dict}

levels_cmp = np.array(sorted(fh_exp))
fh_bohr = {n: fh_ionised * (1.0 - 1.0/n**2) for n in levels_cmp}
# Reduced Dirac (simplified)
reduced_compton_freq = c / (lambda_e + lambda_p)
alpha_sq_eff = 2.0 * fh_ionised / reduced_compton_freq
alpha_eff = np.sqrt(alpha_sq_eff) if alpha_sq_eff > 0 else 0.007297
def dirac_binding_freq(n, j, alpha, nu_c):
    gamma = np.sqrt((j+0.5)**2 - alpha**2)
    n_eff = n - (j+0.5) + gamma
    return nu_c * (1.0 - 1.0/np.sqrt(1.0 + (alpha/n_eff)**2))
f_bind_1s_dirac = dirac_binding_freq(1, 0.5, alpha_eff, reduced_compton_freq)
fh_dirac = {n: f_bind_1s_dirac - dirac_binding_freq(n, 0.5, alpha_eff, reduced_compton_freq) for n in levels_cmp}
ppm_anchor = ppm_residual_dict(fh_anchor, fh_exp)
ppm_bohr = ppm_residual_dict(fh_bohr, fh_exp)
ppm_dirac = ppm_residual_dict(fh_dirac, fh_exp)

# Compute statistics for the printout
def residual_stats(ppm_dict):
    vals = np.array(list(ppm_dict.values()))
    mae = np.mean(np.abs(vals))
    rms = np.sqrt(np.mean(vals**2))
    maxabs = np.max(np.abs(vals))
    return mae, rms, maxabs

anchor_mae, anchor_rms, anchor_max = residual_stats(ppm_anchor)
bohr_mae, bohr_rms, bohr_max = residual_stats(ppm_bohr)
dirac_mae, dirac_rms, dirac_max = residual_stats(ppm_dirac)

ppm_anchor_mean = np.mean(list(ppm_anchor.values()))
ppm_bohr_mean = np.mean(list(ppm_bohr.values()))
ppm_dirac_mean = np.mean(list(ppm_dirac.values()))

print("\n" + "=" * 80)
print("MODEL COMPARISON VS EXPERIMENT")
print("=" * 80)
print("Note: gravity simulation is dimensionless internally; constants are used only for conversion to frequency.")
print("Dirac comparison uses a reduced hydrogenic reference and excludes Lamb/QED corrections.")
print()
print(f"{'n':>3} {'exp (Hz)':>18} {'gravity ppm':>14} {'Bohr ppm':>12} {'Dirac ppm':>12}")
for n in levels_cmp:
    print(f"{n:>3d} {fh_exp[n]:18.3f} {ppm_anchor[n]:14.3f} {ppm_bohr[n]:12.3f} {ppm_dirac[n]:12.3f}")
print()
print(f"{'Model':<28} {'MAE ppm':>12} {'RMS ppm':>12} {'Max |ppm|':>12}")
print(f"{'Gravity simulation anchor':<28} {anchor_mae:12.3f} {anchor_rms:12.3f} {anchor_max:12.3f}")
print(f"{'Bohr (ionization matched)':<28} {bohr_mae:12.3f} {bohr_rms:12.3f} {bohr_max:12.3f}")
print(f"{'Dirac (reduced baseline)':<28} {dirac_mae:12.3f} {dirac_rms:12.3f} {dirac_max:12.3f}")
print()
print(f"Gravity mean residual: {ppm_anchor_mean:.3f} ppm")
print(f"Bohr mean residual:    {ppm_bohr_mean:.3f} ppm")
print(f"Dirac mean residual:   {ppm_dirac_mean:.3f} ppm")



# -------------------------
# Diagnostic: local max-curvature finding in 1.2 < n < 2.1
# -------------------------
def local_spline_feature(x, y, window=(1.2, 2.1), smooth_frac=1e-5, grid_size=1200):
    """
    Fit a local smoothing spline inside the chosen window and extract:
    - minimum location
    - steepest descent location
    - maximum curvature location

    smooth_frac controls smoothing:
      0      -> interpolating spline
      1e-6   -> very light smoothing
      1e-5   -> mild smoothing
      1e-4   -> stronger smoothing
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)

    mask = (x >= window[0]) & (x <= window[1])
    xw = x[mask]
    yw = y[mask]

    if len(xw) < 8:
        raise ValueError("Not enough points in the local spline window")

    order = np.argsort(xw)
    xw = xw[order]
    yw = yw[order]

    # scale x for numerical stability
    x_mean = np.mean(xw)
    x_std = np.std(xw)
    xs = (xw - x_mean) / x_std

    s_val = smooth_frac * len(yw) * np.var(yw)
    spl = UnivariateSpline(xs, yw, s=s_val, k=4)

    x_fit = np.linspace(xw[0], xw[-1], grid_size)
    x_fit_s = (x_fit - x_mean) / x_std

    y_fit = spl(x_fit_s)
    d1 = spl.derivative(1)(x_fit_s) / x_std
    d2 = spl.derivative(2)(x_fit_s) / (x_std**2)

    # Scale-invariant (visual) curvature
    # Normalize d1 and d2 by data extents so that curvature shape is independent of physical units
    x_range = np.max(xw) - np.min(xw)
    y_range = np.max(yw) - np.min(yw)
    if x_range == 0: x_range = 1.0
    if y_range == 0: y_range = 1.0
    
    d1_norm = d1 * (x_range / y_range)
    d2_norm = d2 * ((x_range**2) / y_range)

    # geometric curvature (normalized)
    curvature = np.abs(d2_norm) / (1.0 + d1_norm**2)**1.5

    idx_min = np.argmin(y_fit)
    idx_steep = np.argmin(d1)
    idx_curv = np.argmax(curvature)

    return {
        "x_data": xw,
        "y_data": yw,
        "x_fit": x_fit,
        "y_fit": y_fit,
        "d1_fit": d1,
        "d2_fit": d2,
        "curvature_fit": curvature,
        "x_min": float(x_fit[idx_min]),
        "y_min": float(y_fit[idx_min]),
        "x_steepest": float(x_fit[idx_steep]),
        "slope_steepest": float(d1[idx_steep]),
        "x_curvature": float(x_fit[idx_curv]),
        "curvature_max": float(curvature[idx_curv]),
    }


# Restrict to 1 < n < 2 first
mask_n12 = (n_values > 1.2) & (n_values < 2.1)
n_12 = n_values[mask_n12]
bc_12 = barycenter_distance[mask_n12]

baseline_sim = 471964.337
norm_sim_12 = step_counter[mask_n12] / (n_12**2 * lam)
diff_sim_12 = baseline_sim - norm_sim_12

# Analyze only the visually relevant region
feature_window = (1.2, 2.1)

bc_feat = local_spline_feature(
    n_12, bc_12,
    window=feature_window,
    smooth_frac=1e-5
)

diff_feat = local_spline_feature(
    n_12, diff_sim_12,
    window=feature_window,
    smooth_frac=1e-5
)

print("\n" + "=" * 80)
print("LOCAL FEATURE DIAGNOSTICS (1.4 < n < 1.9)")
print("=" * 80)
print(f"Barycenter minimum:          n = {bc_feat['x_min']:.5f}, value = {bc_feat['y_min']:.5f}")
print(f"Barycenter max curvature:    n = {bc_feat['x_curvature']:.5f}, curvature = {bc_feat['curvature_max']:.5f}")
print(f"Baseline minimum:            n = {diff_feat['x_min']:.5f}, value = {diff_feat['y_min']:.5f}")
print(f"Baseline steepest descent:   n = {diff_feat['x_steepest']:.5f}, slope = {diff_feat['slope_steepest']:.5f}")
print(f"Baseline max curvature:      n = {diff_feat['x_curvature']:.5f}, curvature = {diff_feat['curvature_max']:.5f}")
print(f"Offset: barycenter min vs baseline max curvature = {bc_feat['x_min'] - diff_feat['x_curvature']:+.5f}")
print(f"Offset: barycenter min vs baseline steepest descent = {bc_feat['x_min'] - diff_feat['x_steepest']:+.5f}")











# ============================================================================
# Figure 1: Three panels in one row (A, B, D)
# ============================================================================
plt.style.use('default')
fig1, axs1 = plt.subplots(1, 3, figsize=(18, 5))
fig1.suptitle("Broad Diagnostics: Continuous Hydrogen-Transition Structure (n = 2 ... )", fontsize=16, fontweight='bold')

# Panel A: calibrated frequency vs continuous n (swapped axes)
ax = axs1[0]
stride = max(1, len(n_values)//20000)
sort_idx = np.argsort(freq)
freq_sorted_all = freq[sort_idx]
n_sorted = n_values[sort_idx]
ax.plot(freq_sorted_all[::stride], n_sorted[::stride], linewidth=1.4, alpha=0.85, label='Saved samples')
anchor_levels = np.array(sorted(raw_crossings))
anchor_n_exact = np.array([np.sqrt(raw_crossings[k]["n_sq"]) for k in anchor_levels])
anchor_freqs = np.array([fh_anchor[k] for k in anchor_levels])
exp_levels = np.array(sorted(fh_exp))
exp_freqs = np.array([fh_exp[k] for k in exp_levels])
ax.plot(anchor_freqs, anchor_n_exact, 'o', markersize=7, label='Simulation anchors')
ax.plot(exp_freqs, exp_levels, 'x', markersize=7, mew=1.5, label='Experiment')
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Continuous state coordinate n")
ax.set_title("A. Calibrated Frequency vs Continuous n")
ax.grid(True, alpha=0.3)
ax.legend(fontsize=9)

# Panel B: normalized transition-collapse of Lz (exclude n = max_nshell)
ax = axs1[1]
for k in sorted(collapse_curves):
    if k <= max_nshell - 1:
        x_norm, y_norm = collapse_curves[k]
        ax.plot(x_norm, y_norm, linewidth=1.8, alpha=0.95, label=f'n={k}')
ax.axvline(0.0, linestyle='--', alpha=0.7)
ax.set_xlabel(r"Normalized detuning $(f-f_k)/\Delta f_k$")
ax.set_ylabel(r"Normalized $L_z$")
ax.set_title(r"B. Normalized Transition Collapse of $L_z$")
ax.grid(True, alpha=0.3)
ax.legend(fontsize=9)

# Panel C: save-interval residuals in ppm
ax = axs1[2]
n_list = []
ppm_residual = []
ppm_low = []
ppm_high = []
for k in sorted(saved_crossing_info):
    info = saved_crossing_info[k]
    if info is None:
        continue
    f_exact = info["f_exact"]
    f_est = info["f_saved_interp"]
    f_low = info["f_low"]
    f_high = info["f_high"]
    ppm = 1e6 * (f_est - f_exact) / f_exact
    ppm_lo = 1e6 * (f_est - f_low) / f_exact
    ppm_hi = 1e6 * (f_high - f_est) / f_exact
    n_list.append(k)
    ppm_residual.append(ppm)
    ppm_low.append(ppm_lo)
    ppm_high.append(ppm_hi)
n_arr = np.array(n_list)
ppm_arr = np.array(ppm_residual)
ppm_err = np.vstack([ppm_low, ppm_high])
ax.errorbar(n_arr, ppm_arr, yerr=ppm_err, fmt='o', capsize=4, linewidth=1.2)
ax.axhline(0.0, linestyle='--', alpha=0.7)
ax.set_xlabel("Transition n")
ax.set_ylabel("Residual (ppm)")
ax.set_title("C. Save-Interval Residuals in ppm")
ax.grid(True, alpha=0.3)
for x, y in zip(n_arr, ppm_arr):
    ax.annotate(f"{y:.2f}", (x, y), textcoords="offset points", xytext=(5, 5), fontsize=8)

plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.savefig("master_transition_diagnostics.png", dpi=300, bbox_inches='tight')
print("\n[OK] Saved: master_transition_diagnostics.png")




# ============================================================================
# Figure 2: Three panels in one row (barycenter, baseline-subtracted, residuals)
# ============================================================================
fig2, axs2 = plt.subplots(1, 3, figsize=(18, 5.5))
fig2.suptitle("Geometric Origin, Baseline Consistency & Agreement with Experiment", fontsize=15, fontweight='bold')

# ----------------------------------------------------------------------------
# Panel A: Barycenter distance vs continuous n (n from ~1 to max_nshell)
# ----------------------------------------------------------------------------
ax = axs2[0]
mask_n = n_values <= max_nshell
n_bc = n_values[mask_n]
bc_plot = barycenter_distance[mask_n]
ax.plot(n_bc, bc_plot, linewidth=1.8, color='C0', label='Simulation (barycenter distance)')
for n_int in range(2, max_nshell+1):
    ax.axvline(n_int, linestyle='--', alpha=0.5, color='gray', linewidth=0.8)
anchor_n = [np.sqrt(raw_crossings[n]["n_sq"]) for n in range(2, max_nshell+1)]
anchor_bc = bc_interp(anchor_n)
ax.plot(anchor_n, anchor_bc, 's', markersize=6, color='C0', label='Simulation anchors')

# Quadratic fit (from diagnostic)
ax.plot(bc_feat["x_fit"], bc_feat["y_fit"], color='purple', linewidth=1.4, alpha=0.9, label='Local spline fit')
ax.axvline(bc_feat["x_min"], color='purple', linestyle=':', linewidth=1.4, label='Barycenter minimum')
ax.axvline(bc_feat["x_curvature"], color='magenta', linestyle='--', linewidth=1.4, label='Max curvature')

ax.set_xlabel('State coordinate n')
ax.set_ylabel('Barycenter distance (relative)')
ax.set_title('A. Barycenter distance vs n (geometric signature)')
ax.grid(True, alpha=0.3)
ax.legend(fontsize=8, ncol=2)
ax.set_xlim(0.9, max_nshell + 0.1)

# ----------------------------------------------------------------------------
# Panel B: Baseline-subtracted consistency check (simulation vs experiment)
# ----------------------------------------------------------------------------
ax = axs2[1]
# Use already defined baseline_sim and lam (from earlier)
# baseline_sim = 471964.337, lam already computed
norm_sim = step_counter / (n_values**2 * lam)
diff_sim = baseline_sim - norm_sim
sort_idx_n = np.argsort(n_values)
n_sorted = n_values[sort_idx_n]
diff_sim_sorted = diff_sim[sort_idx_n]

ax.plot(n_sorted, diff_sim_sorted, linewidth=1.0, alpha=0.35, color='C0', label='Simulation (raw)')
ax.plot(diff_feat["x_fit"], diff_feat["y_fit"], linewidth=1.8, color='C0', label='Local spline fit')
ax.axvline(diff_feat["x_curvature"], color='magenta', linestyle='--', linewidth=1.4, label='Max curvature')
ax.axvline(diff_feat["x_steepest"], color='purple', linestyle=':', linewidth=1.2, label='Steepest descent')

# Experimental points
exp_n = [1.0, 2.0, 3.0, 4.0]
exp_diff = [0.0, -0.713788744242, -0.095972375167, 0.020709850987]
ax.plot(exp_n, exp_diff, 'ro', markersize=8, label='Experiment (raw)')
ax.plot(exp_n, exp_diff, 'g--', linewidth=1.5, label='Connecting line (experiment)')

exp_diff_shifted = [d - exp_diff[0] for d in exp_diff]
ax.plot(exp_n, exp_diff_shifted, 'bs', markersize=6, label='Experiment (shifted to zero at n=1)')
ax.plot(exp_n, exp_diff_shifted, 'b--', linewidth=1.0, alpha=0.7)

ax.set_xlabel('State coordinate n')
ax.set_ylabel(r'Baseline $-$ normalized quantity')
ax.set_title('B. Baseline-subtracted consistency check')
ax.grid(True, alpha=0.3)
ax.legend(fontsize=8)
ax.set_xlim(0.9, max_nshell + 0.1)

# ----------------------------------------------------------------------------
# Panel C: Residuals vs experiment (ppm)
# ----------------------------------------------------------------------------
ax = axs2[2]
ax.plot(levels_cmp, [ppm_anchor[n] for n in levels_cmp], marker='s', linewidth=1.8, label='Gravity simulation')
ax.plot(levels_cmp, [ppm_bohr[n] for n in levels_cmp], marker='^', linewidth=1.6, label='Bohr')
ax.plot(levels_cmp, [ppm_dirac[n] for n in levels_cmp], marker='d', linewidth=1.6, label='Dirac (reduced)')
ax.axhline(0.0, linestyle='--', alpha=0.7)
ax.set_xlabel('Transition n')
ax.set_ylabel('Residual vs experiment (ppm)')
ax.set_title('C. Residual comparison in ppm')
ax.grid(True, alpha=0.3)
ax.legend(fontsize=9)

plt.tight_layout(rect=[0, 0, 1, 0.94])
plt.savefig("normalized_consistency_check.png", dpi=300, bbox_inches='tight')
print("[OK] Saved: normalized_consistency_check.png")





plt.show()
print("\nAnalysis successful. Two figures produced.")