% ===========================================================================================================
%  On the Observer-Frame Interpretation of Radial Coordinates in Gravitational Potentials
%                                 Brahim Benaissa
% ===========================================================================================================
%  LOAD & PREP   → Read rotation/mass → Enclosed Masse Profiles 
%  MEZZI         → Integrate PG flow → ζ = exp(377*(K-1)) → transform r_obs→r_true, M_obs→M_true
%  VALIDATION    → Rotation curves → compute v_obs vs v_Newton residuals → RMS/R²/χ²  
%  SCALING       → Analyze surface density ratios, Deflection angles, Cross checkings
%  RAR TEST      → Validate g_dyn vs g_bar for C1(Newton), C2(Empirical), C3(Geometric) → log-scatter
%  BTFR FIT      → Observed and True slope fits → compute β-offset vs ζ
% ===========================================================================================================

clear;

% Free Parameter : Vacuum Compliance = 377
Compliance = 377;       


%% ===== SECTION 1: DATA PREPARATION =====
% Read galaxy rotation data (including gas velocity contribution)
[galaxy_names, observed_radius, observed_speed, e_Vobs, Vgas] = read_rotation_data('OBSERVED_Rotation.txt');
% Save the rotation data
save('galaxy_Rotation_data.mat', 'galaxy_names', 'observed_radius', 'observed_speed', 'Vgas');

% Read galaxy mass data  
[gs_galaxy_names, gs_M2L, e_M2L, gs_M2Lb, gs_L36, e_L36, gs_Lbulge, ...
gs_Rdisk, gs_Reff, Distance_to_G, e_D, Mbar, e_Mbar] = read_mass_data('OBSERVED_Mass.txt');

% Gravitational constant in convenient units: kpc * (km/s)^2 / Msun
G_kpc = 4.300917270437e-06;  % G = 4.300917270437 × 10^{-6} kpc (km/s)^2 / Msun

% Initialize storage with error containers
total_baryonic_mass = cell(size(galaxy_names));
total_baryonic_mass_err = cell(size(galaxy_names));
total_baryonic_mass_direct = cell(size(galaxy_names));
total_baryonic_mass_direct_err = cell(size(galaxy_names));
baryonic_mass_ratio = zeros(length(galaxy_names), 1);
total_gas_mass = cell(size(galaxy_names));
total_gas_mass_err = cell(size(galaxy_names));
enclosed_masses = cell(size(galaxy_names));
enclosed_masses_err = cell(size(galaxy_names));
enclosed_gas_masses = cell(size(galaxy_names));
enclosed_gas_masses_err = cell(size(galaxy_names));
radius_uncertainties = cell(size(galaxy_names));

% Calculate the masses for each galaxy in order
for g = 1:length(galaxy_names)
    current_galaxy = galaxy_names{g};
    % Find index of current galaxy in gs_galaxy_names
    idx = find(strcmp(gs_galaxy_names, current_galaxy), 1);
    
    % ===== DIRECT BARYONIC MASS (AUTHORITATIVE) =====
    % Mbar from SPARC ALREADY INCLUDES stars + gas (with helium correction)
    total_baryonic_mass_direct{g} = Mbar(idx) * 1e9;
    total_baryonic_mass_direct_err{g} = e_Mbar(idx) * 1e9;
    
    % Get distance and its uncertainty
    D = Distance_to_G(idx);
    e_D_val = e_D(idx);
    
    % Convert radius uncertainties
    R_kpc = observed_radius{g};
    e_R_kpc = R_kpc * (e_D_val/D);
    radius_uncertainties{g} = e_R_kpc * 3.086e+19;
    
    % Error propagation for disk luminosity
    L_disk = gs_L36(idx) - gs_Lbulge(idx);
    e_L_disk = sqrt(e_L36(idx)^2 + 0^2);
    
    % Disk mass with error propagation
    M_disk = gs_M2L(idx) * L_disk * 1e9;
    e_M_disk = M_disk * (e_M2L(idx)/gs_M2L(idx));
    
    % Bulge mass with error propagation
    M_bulge = gs_M2Lb(idx) * gs_Lbulge(idx) * 1e9;
    e_M_bulge = 0;
    
    % Total stellar mass with error
    total_stellar_mass = M_disk + M_bulge;
    total_stellar_mass_err = sqrt(e_M_disk^2 + e_M_bulge^2);
    
    gs_Rdisk_ordered(g) = gs_Rdisk(idx);
    gs_Reff_ordered(g) = gs_Reff(idx);
    M_disk_ordered(g) = M_disk;
    M_bulge_ordered(g) = M_bulge;
    
    % Process enclosed masses with errors
    radii = observed_radius{g};
    Vgas_vals = Vgas{g};
    enc_mass = zeros(size(radii));
    enc_mass_err = zeros(size(radii));
    enc_gas_mass = zeros(size(radii));
    enc_gas_mass_err = zeros(size(radii));
    fractional_Vgas_error = 0.05;
    
    % Loop over each radius to calculate enclosed mass components
    for r_idx = 1:length(radii)
        r = radii(r_idx);
        Vgas_val = Vgas_vals(r_idx);
        
        % ===== STELLAR COMPONENTS =====
        [enc_disk, e_enc_disk] = compute_enclosed_disk_mass(r, M_disk, e_M_disk, gs_Rdisk(idx));
        [enc_bulge, e_enc_bulge] = compute_enclosed_bulge_mass(r, M_bulge, e_M_bulge, gs_Reff(idx));
        
        % ===== GAS COMPONENT =====
        if Vgas_val > 0
            %  Approximation for disk geometry
            M_enc_gas = (Vgas_val^2 * r) / G_kpc;
            enc_gas_mass(r_idx) = M_enc_gas;
            e_Vgas_val = fractional_Vgas_error * Vgas_val;
            e_M_enc_gas = M_enc_gas * (2 * e_Vgas_val / Vgas_val);
            enc_gas_mass_err(r_idx) = e_M_enc_gas;
        else
            enc_gas_mass(r_idx) = 0;
            enc_gas_mass_err(r_idx) = 0;
        end
        
        % ===== TOTAL ENCLOSED MASS (CALCULATED) =====
        enc_mass(r_idx) = enc_disk + enc_bulge + enc_gas_mass(r_idx);
        enc_mass_err(r_idx) = sqrt(e_enc_disk^2 + e_enc_bulge^2 + enc_gas_mass_err(r_idx)^2);
    end
    
    % Total gas mass = enclosed gas at largest radius
    %   Approximation
    total_gas_mass{g} = enc_gas_mass(end);
    total_gas_mass_err{g} = enc_gas_mass_err(end);
    
    % ===== TOTAL BARYONIC MASS (YOUR CALCULATION - for comparison) =====
    total_baryonic_mass{g} = total_stellar_mass + total_gas_mass{g};
    total_baryonic_mass_err{g} = sqrt(total_stellar_mass_err^2 + total_gas_mass_err{g}^2);
    
    % ===== COMPARISON METRIC =====
    if total_baryonic_mass_direct{g} > 0
        baryonic_mass_ratio(g) = total_baryonic_mass{g} / total_baryonic_mass_direct{g};
    else
        baryonic_mass_ratio(g) = NaN;
    end
    
    % Store the calculated enclosed masses
    enclosed_masses{g} = enc_mass;
    enclosed_masses_err{g} = enc_mass_err;
    enclosed_gas_masses{g} = enc_gas_mass;
    enclosed_gas_masses_err{g} = enc_gas_mass_err;
    
    % Load data with uncertainties
    Robs_kpc = observed_radius{g};
    Vobs_km_s = observed_speed{g};
    e_Vobs_km_s = e_Vobs{g};
    
    % Unit conversions
    Robs_m = Robs_kpc * 3.086e+19;
    Vobs_m_s = Vobs_km_s * 1e3;
    e_Vobs_m_s = e_Vobs_km_s * 1e3;
    
    % Get Mass with uncertainties
    current_enclosed = enclosed_masses{g};
    current_enclosed_err = enclosed_masses_err{g};
    solar_mass_kg = 1.9885e30;
    Enclosed_Mass_kg = current_enclosed * solar_mass_kg;
    e_Enclosed_Mass_kg = current_enclosed_err * solar_mass_kg;
    
    % Get radius uncertainties
    e_Robs_m = radius_uncertainties{g};
    All_Vobs{g} = Vobs_m_s;
    All_e_Vobs{g} = e_Vobs_m_s;
    All_Robs{g} = Robs_m;
    All_e_Robs{g} = e_Robs_m;
    All_Enclosed_Masses{g} = Enclosed_Mass_kg;
    All_e_Enclosed_Masses{g} = e_Enclosed_Mass_kg;
end

% Save results including direct baryonic mass and comparison metrics
save('galaxy_mass_properties.mat', ...
    'galaxy_names', ...
    'total_baryonic_mass', ...
    'total_baryonic_mass_err', ...
    'total_baryonic_mass_direct', ...
    'total_baryonic_mass_direct_err', ...
    'baryonic_mass_ratio', ...
    'total_gas_mass', 'total_gas_mass_err', ...
    'enclosed_masses', 'enclosed_masses_err', ...
    'enclosed_gas_masses', 'enclosed_gas_masses_err', ...
    'Distance_to_G', 'radius_uncertainties', ...
    'gs_Rdisk_ordered', 'gs_Reff_ordered', 'M_disk_ordered', ...
    'M_bulge_ordered', 'All_Vobs', 'All_Robs', ...
    'All_Enclosed_Masses', 'All_e_Vobs', 'All_e_Robs', 'All_e_Enclosed_Masses');



%% ===== SECTION 2: MEZZI ANALYSIS =====
%% --- 2.1 Individual Figures Configuration ---
IndFigures = true;        % true  For Plotting Individual Galaxy Figures

%% --- 2.2 Constants ---
% Fundamental constants
G      = 6.67430e-11;      % Gravitational constant (m^3 kg^-1 s^-2)
c      = 2.99792458e8;     % Speed of light (m s^-1)
arcsec = 1 / 206265;       % 1 arcsecond in radians
Mpc    = 3.085677581e22;   % 1 Megaparsec (m)
pc     = 3.085677581e16;   % 1 parsec (m)
Msun   = 1.9885e30;         % Solar mass (kg)

%% --- 2.3 Load Observational Data ---
fprintf('[LOAD] Loading observational data...\n');
load('galaxy_mass_properties.mat');
numGalaxies = length(All_Robs);
fprintf('[LOAD] Analyzing %d galaxies with Mezzi corrections...\n', numGalaxies);

%% --- 2.4 Main Galaxy Processing ---
Results = struct();
for g = 1:numGalaxies
    fprintf('[Geometric Scaling] Galaxy %d/%d: %s\n', g, numGalaxies, galaxy_names{g});
    Results(g).GalaxyName = galaxy_names{g};
    
    %% --- 2.4.1 Basic Galaxy Parameters ---
    Distance_to_Galaxy = Distance_to_G(g) * Mpc;
    
    % Observed rotation curve data
    r_obs_I = All_Robs{g}(:);
    M_enc_obs = All_Enclosed_Masses{g}(:);
    
    % Quality control
    valid_idx = (r_obs_I > 0) & isfinite(r_obs_I) & isfinite(M_enc_obs);
    r_obs_I = r_obs_I(valid_idx);
    M_enc_obs = M_enc_obs(valid_idx);
    
    % Key scales
    r_max = max(r_obs_I);
    Mtot_obs = M_enc_obs(end);
    
    % Mass profile interpolator
    M_func_obs = @(x) ...
        (x <= r_max) .* interp1(r_obs_I, M_enc_obs, x, 'pchip', 'extrap') + ...
        (x > r_max) .* Mtot_obs;
    
    %% --- 2.4.2 Painlevé-Gullstrand (PG) Flow Accumulation & Mezzi Scale Factor ---
    K_kinematic_compression = ones(size(r_obs_I));
    v_flow = @(x) -sqrt(2 * G * M_func_obs(x) ./ max(x, 1e-10*pc));
    integrand = @(x) v_flow(x) ./ (c .* max(x, 1e-10*pc));
    
    for i = 1:length(r_obs_I)
        r_i = r_obs_I(i);
        if abs(r_i - Distance_to_Galaxy) < 1e-8 * max(r_i, Distance_to_Galaxy)
            A_accumulated_potential = 0;
        else
            r_min_path = min(r_i, Distance_to_Galaxy);
            r_max_path = max(r_i, Distance_to_Galaxy);
            r_path = logspace(log10(max(r_min_path, 1e-3*pc)), ...
                log10(r_max_path), 500);
            A_accumulated_potential = trapz(r_path, integrand(r_path));
        end
        K_kinematic_compression(i) = exp(A_accumulated_potential);
    end
    
    % Mezzi scale factor
    Mezzi_scale_factor_theory = exp(Compliance .* (K_kinematic_compression - 1));
    
    %% --- 2.4.3 Scaling the Mass Profile ---
    [M_true_cum, M_obs_cum, mass_ratio] = MassProfileFromMass(...
        r_obs_I, All_Enclosed_Masses{1, g}, Mezzi_scale_factor_theory);
    M_enc_true = M_true_cum;
    
    % Coordinate transformation
    r_obs = All_Robs{g}(:);
    r_true = r_obs ./ Mezzi_scale_factor_theory;
    
    % True mass profile interpolator
    M_func_true = @(x) ...
        (x <= r_max) .* interp1(r_true, M_enc_true, x, 'pchip', 'extrap') + ...
        (x > r_max) .* M_enc_true(end);
    
    %% --- 2.4.4 Orbital Velocity Recovery Test ---
    % Observed velocities
    v_obs = All_Vobs{g}(:);
    v_obs = v_obs(valid_idx);
    v_true = v_obs;  % Ground truth
    
    % Newtonian prediction in observed frame
    v_Newton_obs = sqrt(G * M_enc_obs ./ r_obs);
    
    % Newtonian prediction in true frame
    v_Newton_true = sqrt(G * M_enc_true ./ r_true);
    v_Newton_true_To_observer_Frame = v_Newton_true;
    
    % Recovery errors
    residual_naive = v_true - v_Newton_obs;
    residual_corrected = v_true - v_Newton_true_To_observer_Frame;
    rms_naive = sqrt(mean(residual_naive.^2));
    rms_corrected = sqrt(mean(residual_corrected.^2));
    
    %% --- 2.4.5 Density Profiles ---
    % Observed frame density
    r_fine_obs = logspace(log10(min(r_obs)), log10(r_max*1.1), 2000);
    M_fine_obs = arrayfun(M_func_obs, r_fine_obs);
    dMdr_obs = gradient(M_fine_obs, r_fine_obs);
    rho_fine_obs = dMdr_obs ./ (4 * pi * r_fine_obs.^2 + 1e-50);
    
    % True frame density
    r_fine_true = logspace(log10(min(r_true)), log10(max(r_true)*1.1), 2000);
    M_fine_true = arrayfun(M_func_true, r_fine_true);
    dMdr_true = gradient(M_fine_true, r_fine_true);
    rho_fine_true = dMdr_true ./ (4 * pi * r_fine_true.^2 + 1e-50);
    
    %% --- 2.4.6 Surface Density ---
    % Observed surface density
    R_grid_obs = r_obs;
    Sigma_obs = zeros(size(R_grid_obs));
    r_upper_obs = max(r_fine_obs);
    
    for j = 1:length(R_grid_obs)
        R = R_grid_obs(j);
        if R <= 0 || R > r_max || isnan(R)
            Sigma_obs(j) = 0;
            continue;
        end
        r_int_obs = linspace(R*1.0001, r_upper_obs, 1000);
        rho_int_obs = interp1(r_fine_obs, rho_fine_obs, r_int_obs, 'pchip', 0);
        integrand_obs = rho_int_obs .* r_int_obs ./ ...
            sqrt(max(r_int_obs.^2 - R^2, 1e-20));
        Sigma_obs(j) = 2 * trapz(r_int_obs, integrand_obs);
    end
    
    % True surface density
    R_grid_true = r_true;
    Sigma_true = zeros(size(R_grid_true));
    r_upper_true = max(r_fine_true);
    
    for j = 1:length(R_grid_true)
        R = R_grid_true(j);
        if R <= 0 || R > max(r_true) || isnan(R)
            Sigma_true(j) = 0;
            continue;
        end
        r_int_true = linspace(R*1.0001, r_upper_true, 1000);
        rho_int_true = interp1(r_fine_true, rho_fine_true, r_int_true, 'pchip', 0);
        integrand_true = rho_int_true .* r_int_true ./ ...
            sqrt(max(r_int_true.^2 - R^2, 1e-20));
        Sigma_true(j) = 2 * trapz(r_int_true, integrand_true);
    end
    
    % Transform true to observed grid
    Sigma_true_on_obs_grid = zeros(size(R_grid_obs));
    for j = 1:length(R_grid_obs)
        R_obs = R_grid_obs(j);
        [~, idx] = min(abs(r_obs - R_obs));
        if idx > 0 && idx <= length(Mezzi_scale_factor_theory)
            Mezzi_at_R = Mezzi_scale_factor_theory(idx);
            R_true = R_obs / Mezzi_at_R;
            Sigma_val = interp1(R_grid_true, Sigma_true, R_true, 'pchip', 0);
            Sigma_true_on_obs_grid(j) = Sigma_val / (Mezzi_at_R^2 + 1e-20);
        end
    end
    
    %% --- 2.4.7 Deflection Angles ---
     % Observed frame
    alpha_rad_obs = zeros(size(R_grid_obs));
    for j = 1:length(R_grid_obs)
        R = R_grid_obs(j);
        if R < 1e-8*pc || isnan(R)
            alpha_rad_obs(j) = 0;
            continue;
        end
        M_within_R_obs = M_func_obs(R);
        alpha_rad_obs(j) = (4 * G / c^2) * (M_within_R_obs / R);
    end
    
    % True frame
    alpha_rad_true = zeros(size(R_grid_true));
    for j = 1:length(R_grid_true)
        R = R_grid_true(j);
        if R < 1e-8*pc || isnan(R) || R > max(r_true)
            alpha_rad_true(j) = 0;
            continue;
        end
        M_within_R_true = M_func_true(R);
        alpha_rad_true(j) = (4 * G / c^2) * (M_within_R_true / R);
    end
    
    %% --- 2.4.8 Store Results ---
    Results(g).galaxy_index = g;
    Results(g).Distance_to_Galaxy = Distance_to_Galaxy;
    Results(g).r_obs = r_obs;
    Results(g).M_enc_obs = M_enc_obs;
    Results(g).v_obs = v_obs;
    Results(g).r_max = r_max;
    Results(g).Mtot_obs = Mtot_obs;
    Results(g).PG = K_kinematic_compression;
    Results(g).Mezzi_th = Mezzi_scale_factor_theory;
    Results(g).C = Compliance;
    Results(g).r_true = r_true;
    Results(g).M_enc_true = M_enc_true;
    Results(g).Mtot_true = M_enc_true(end);
    Results(g).v_true = v_true;
    Results(g).v_Newton_obs = v_Newton_obs;
    Results(g).v_Newton_true = v_Newton_true;
    Results(g).v_Newton_true_To_obs = v_Newton_true_To_observer_Frame;
    Results(g).residual_naive = residual_naive;
    Results(g).residual_corrected = residual_corrected;
    Results(g).rms_naive = rms_naive;
    Results(g).rms_corrected = rms_corrected;
    Results(g).r_fine_obs = r_fine_obs;
    Results(g).M_fine_obs = M_fine_obs;
    Results(g).dMdr_obs = dMdr_obs;
    Results(g).rho_fine_obs = rho_fine_obs;
    Results(g).r_fine_true = r_fine_true;
    Results(g).M_fine_true = M_fine_true;
    Results(g).dMdr_true = dMdr_true;
    Results(g).rho_fine_true = rho_fine_true;
    Results(g).R_grid_obs = R_grid_obs;
    Results(g).Sigma_obs = Sigma_obs;
    Results(g).r_upper_obs = r_upper_obs;
    Results(g).R_grid_true = R_grid_true;
    Results(g).Sigma_true = Sigma_true;
    Results(g).r_upper_true = r_upper_true;
    Results(g).Sigma_true_on_obs_grid = Sigma_true_on_obs_grid;
    Results(g).Sigma_obs_Msun_pc2 = Sigma_obs * pc^2 / Msun;
    Results(g).Sigma_true_Msun_pc2 = Sigma_true * pc^2 / Msun;
    Results(g).alpha_rad_obs = alpha_rad_obs;
    Results(g).alpha_rad_true = alpha_rad_true;
    Results(g).alpha_obs_arcsec = alpha_rad_obs / arcsec;
    Results(g).alpha_true_arcsec = alpha_rad_true / arcsec;
    Results(g).max_alpha_obs_arcsec = max(Results(g).alpha_obs_arcsec);
    Results(g).max_alpha_true_arcsec = max(Results(g).alpha_true_arcsec);
    Results(g).Sigma_ratio = Sigma_true ./ (Sigma_obs + 1e-20);
    Results(g).alpha_ratio_true = alpha_rad_true ./ (alpha_rad_obs + 1e-20);
    Results(g).v_flow_obs = v_flow(r_obs);
    Results(g).integrand_obs = integrand(r_obs);
end

%% --- 2.5 Visualization ---

% Output directory structure
baseDir = 'Results';
folders = {'RotationCurve', 'MassProfile', 'SurfaceDensity', 'Deflection'};
for k = 1:length(folders)
    subDir = fullfile(baseDir, folders{k});
    if ~exist(subDir, 'dir')
        mkdir(subDir);
    end
end

% Process each galaxy
for g = 1:numGalaxies
    %% --- 2.5.1 Rotation Curve Comparison ---
    r_pc_obs = Results(g).r_obs / pc;
    v_obs_kms = Results(g).v_obs / 1e3;
    v_newt_obs_kms = Results(g).v_Newton_obs / 1e3;
    v_true_kms = Results(g).v_Newton_true_To_obs / 1e3;
    
    if IndFigures
        fig = figure('Color', 'white', 'Position', [100, 100, 800, 600]);
        
        % Plot data points
        scatter(r_pc_obs, v_obs_kms, 60, [0.2 0.4 0.8], 'filled', ...
            'DisplayName', 'Observed', 'MarkerFaceAlpha', 0.7);
        hold on;
        
        % Plot model lines
        plot(r_pc_obs, v_newt_obs_kms, 'b--', 'LineWidth', 2, ...
            'DisplayName', 'Newtonian (Obs Frame)');
        plot(r_pc_obs, v_true_kms, 'r-', 'LineWidth', 2.5, ...
            'DisplayName', 'Newton (True Frame)');
        
        % Formatting
        grid on; box on;
        set(gca, 'XScale', 'log', 'FontSize', 12, 'LineWidth', 1.2);
        xlabel('Radius (pc)', 'FontSize', 16, 'FontWeight', 'bold');
        ylabel('Circular Velocity (km/s)', 'FontSize', 16, 'FontWeight', 'bold');
        title(sprintf('Rotation Curve — %s', galaxy_names{g}), ...
            'FontSize', 14, 'FontWeight', 'bold');
        legend('Location', 'best', 'FontSize', 10, 'Box', 'on');
        
        % Performance annotation
        improvement = 100 * (1 - Results(g).rms_corrected/Results(g).rms_naive);
        annotation_text = sprintf(...
            'RMS (naive) = %.1f km/s\nRMS (corr) = %.1f km/s\nImprovement = %.0f%%', ...
            Results(g).rms_naive/1e3, Results(g).rms_corrected/1e3, improvement);
        annotation('textbox', [0.02, 0.70, 0.25, 0.15], ...
            'String', annotation_text, 'Units', 'normalized', ...
            'FontSize', 10, 'BackgroundColor', 'white', ...
            'EdgeColor', [0.6 0.6 0.6], 'LineWidth', 1);
        
        exportgraphics(fig, fullfile(baseDir, 'RotationCurve', ...
        sprintf('Fig_2_5_1_RotationCurve_%d.png', g)), 'Resolution', 300);
        close(fig);
    end
    
    %% --- 2.5.2 Mass Profile Comparison ---
    r_pc_true = Results(g).r_true / pc;
    M_obs_Msun = Results(g).M_enc_obs / Msun;
    M_true_Msun = Results(g).M_enc_true / Msun;
    
    if IndFigures
        fig = figure('Color', 'white', 'Position', [100, 100, 800, 600]);
        
        loglog(r_pc_obs, M_obs_Msun, 'bo-', 'LineWidth', 2, 'MarkerSize', 6, ...
            'DisplayName', 'Observed');
        hold on;
        loglog(r_pc_true, M_true_Msun, 'rs-', 'LineWidth', 2, 'MarkerSize', 6, ...
            'DisplayName', 'Mezzi-Corrected');
        
        grid on; box on;
        set(gca, 'FontSize', 12, 'LineWidth', 1.2);
        xlabel('Radius (pc)', 'FontSize', 16, 'FontWeight', 'bold');
        ylabel('Enclosed Mass (M_\odot)', 'FontSize', 16, 'FontWeight', 'bold');
        title(sprintf('Mass Profile — %s', galaxy_names{g}), ...
            'FontSize', 14, 'FontWeight', 'bold');
        legend('Location', 'best', 'FontSize', 10, 'Box', 'on');
        
        % Mass ratio annotation
        ratio_mass = Results(g).Mtot_true / Results(g).Mtot_obs;
        annotation_text = sprintf(...
            'M_{tot}(obs) = %.2e M_\\odot\nM_{tot}(true) = %.2e M_\\odot\nRatio = %.2f', ...
            Results(g).Mtot_obs/Msun, Results(g).Mtot_true/Msun, ratio_mass);
        annotation('textbox', [0.02, 0.70, 0.25, 0.15], ...
            'String', annotation_text, 'Units', 'normalized', ...
            'FontSize', 10, 'BackgroundColor', 'white', ...
            'EdgeColor', [0.6 0.6 0.6], 'LineWidth', 1);
        
        exportgraphics(fig, fullfile(baseDir, 'MassProfile', ...
            sprintf('Fig_2_5_2_MassProfile_%d.png', g)), 'Resolution', 300);
        close(fig);
    end
    
    %% --- 2.5.3 Surface Density Comparison ---
    R_pc_obs = Results(g).R_grid_obs / pc;
    R_pc_true = Results(g).R_grid_true / pc;
    Sigma_obs_Msun = Results(g).Sigma_obs * pc^2 / Msun;
    Sigma_true_Msun = Results(g).Sigma_true * pc^2 / Msun;
    
    if IndFigures
        fig = figure('Color', 'white', 'Position', [100, 100, 800, 600]);
        
        loglog(R_pc_obs, Sigma_obs_Msun, 'bo-', 'LineWidth', 2, 'MarkerSize', 6, ...
            'DisplayName', 'Observed');
        hold on;
        loglog(R_pc_true, Sigma_true_Msun, 'rs-', 'LineWidth', 2, 'MarkerSize', 6, ...
            'DisplayName', 'Mezzi-Corrected');
        
        grid on; box on;
        set(gca, 'FontSize', 12, 'LineWidth', 1.2);
        xlabel('Projected Radius (pc)', 'FontSize', 16, 'FontWeight', 'bold');
        ylabel('Surface Density (M_\odot/pc^2)', 'FontSize', 16, 'FontWeight', 'bold');
        title(sprintf('Surface Density — %s', galaxy_names{g}), ...
            'FontSize', 14, 'FontWeight', 'bold');
        legend('Location', 'best', 'FontSize', 10, 'Box', 'on');
        
        % Deviation annotation
        core_ratio = Sigma_true_Msun(1)/(Sigma_obs_Msun(1)+1e-20);
        max_dev = max(abs((Sigma_true_Msun - Sigma_obs_Msun)./Sigma_obs_Msun))*100;
        annotation_text = sprintf(...
            '\\Sigma ratio (core) = %.2f\nMax deviation = %.0f%%', ...
            core_ratio, max_dev);
        annotation('textbox', [0.02, 0.70, 0.25, 0.15], ...
            'String', annotation_text, 'Units', 'normalized', ...
            'FontSize', 10, 'BackgroundColor', 'white', ...
            'EdgeColor', [0.6 0.6 0.6], 'LineWidth', 1);
        
        exportgraphics(fig, fullfile(baseDir, 'SurfaceDensity', ...
            sprintf('Fig_2_5_3_SurfaceDensity_%d.png', g)), 'Resolution', 300);
        close(fig);
    end
    
    %% --- 2.5.4 Deflection Angles ---
    alpha_obs_arc = Results(g).alpha_rad_obs / arcsec;
    alpha_true_arc = Results(g).alpha_rad_true / arcsec;
    
    if IndFigures
        fig = figure('Color', 'white', 'Position', [100, 100, 800, 600]);
        
        semilogx(r_pc_obs, alpha_obs_arc, 'bo-', 'LineWidth', 2, 'MarkerSize', 6, ...
            'DisplayName', 'Observed');
        hold on;
        semilogx(r_pc_true, alpha_true_arc, 'rs-', 'LineWidth', 2, 'MarkerSize', 6, ...
            'DisplayName', 'Mezzi-Corrected');
        
        grid on; box on;
        set(gca, 'YScale', 'log', 'FontSize', 12, 'LineWidth', 1.2);
        xlabel('Impact Parameter (pc)', 'FontSize', 16, 'FontWeight', 'bold');
        ylabel('Deflection Angle (arcsec)', 'FontSize', 16, 'FontWeight', 'bold');
        title(sprintf('Gravitational Deflection — %s', galaxy_names{g}), ...
            'FontSize', 14, 'FontWeight', 'bold');
        legend('Location', 'best', 'FontSize', 10, 'Box', 'on');
        
        exportgraphics(fig, fullfile(baseDir, 'Deflection', ...
            sprintf('Fig_2_5_4_Deflection_%d.png', g)), 'Resolution', 300);
        close(fig);
    end
end

%% --- 2.6 Mezzi Scale Factor Analysis ---
% Analysis directory for output files
analysisDir = fullfile(baseDir, 'MezziAnalysis');
if ~exist(analysisDir, 'dir')
    mkdir(analysisDir);
end

% Initialize summary structure to store galaxy-specific metrics
galaxySummary = struct();

% Loop through each galaxy to compute metrics
for g = 1:numGalaxies
    % Extract key quantities from Results structure
    r_obs = Results(g).r_obs;
    r_true = Results(g).r_true;
    M_enc_obs = Results(g).M_enc_obs;
    M_enc_true = Results(g).M_enc_true;
    Mezzi_scale_factor_theory = Results(g).Mezzi_th;
    K_kinematic_compression = Results(g).PG;
    v_obs = Results(g).v_obs;
    v_Newton_obs = Results(g).v_Newton_obs;
    v_Newton_true = Results(g).v_Newton_true_To_obs;
    
    % Calculate scale factor statistics
    PG_mean = mean(K_kinematic_compression);
    PG_range = max(K_kinematic_compression) - min(K_kinematic_compression);
    Mezzi_mean = mean(Mezzi_scale_factor_theory);
    Mezzi_range = max(Mezzi_scale_factor_theory) - min(Mezzi_scale_factor_theory);
    M_total_obs = M_enc_obs(end);
    M_total_true = M_enc_true(end);
    
    % Velocity analysis at characteristic radii (Effective radius approximation)
    r_eff = r_obs(end) / 2;
    [~, idx_eff] = min(abs(r_obs - r_eff));
    v_flat_obs = v_obs(idx_eff);
    v_newt_flat = v_Newton_obs(idx_eff);
    v_true_flat = v_Newton_true(idx_eff);
    
    % Store summary data in structure
    galaxySummary(g).name = galaxy_names{g};
    galaxySummary(g).PG_mean = PG_mean;
    galaxySummary(g).PG_range = PG_range;
    galaxySummary(g).Mezzi_mean = Mezzi_mean;
    galaxySummary(g).Mezzi_range = Mezzi_range;
    galaxySummary(g).M_total_obs_Msun = M_total_obs / Msun;
    galaxySummary(g).M_total_true_Msun = M_total_true / Msun;
    galaxySummary(g).mass_ratio = M_total_true / M_total_obs;
    galaxySummary(g).v_central_obs_kms = v_obs(1) / 1e3;
    galaxySummary(g).v_flat_obs_kms = v_flat_obs / 1e3;
    galaxySummary(g).v_flat_newt_kms = v_newt_flat / 1e3;
    galaxySummary(g).v_flat_true_kms = v_true_flat / 1e3;
    galaxySummary(g).rms_naive_kms = Results(g).rms_naive / 1e3;
    galaxySummary(g).rms_corrected_kms = Results(g).rms_corrected / 1e3;
    galaxySummary(g).improvement_percent = 100 * (1 - Results(g).rms_corrected/Results(g).rms_naive);
end

%% --- 2.6.1 Extract Arrays and Identify Significant Galaxies ---
% Extract arrays for ensemble analysis
mass_ratios = [galaxySummary.mass_ratio];
PG_means = [galaxySummary.PG_mean];
Mezzi_means = [galaxySummary.Mezzi_mean];


%% --- 2.6.2 Figure: Coordinate scaling with Mezzi Scale Factor Inset ---
fig = figure('Color', 'white', 'Position', [100, 100, 800, 600]);
hold on;

% Extract total observed masses for color mapping
galaxy_masses = [galaxySummary.M_total_obs_Msun];
log_masses = log10(galaxy_masses);

% Mass-based colormap
mass_colors = parula(numGalaxies);

% Sort galaxies by mass for color assignment
[mass_sorted, sort_idx] = sort(log_masses, 'ascend');
[~, color_idx] = sort(sort_idx);

% Plot curves colored by mass (main plot)
for g = 1:numGalaxies
    r_obs_pc = Results(g).r_obs / pc;
    Mezzi = Results(g).Mezzi_th;
    this_color = mass_colors(color_idx(g), :);
    semilogx(r_obs_pc, Mezzi, 'LineWidth', 1.5, 'Color', this_color);
end

cbar = colorbar('Location', 'eastoutside');
cbar.Label.String = 'log_{10}(M_{obs}/M_{☉})';
cbar.Label.FontSize = 14;
cbar.Label.FontWeight = 'bold';

% Set colormap and ticks
colormap(parula);
cbar.Ticks = linspace(0, 1, 5);
cbar.TickLabels = arrayfun(@(x) sprintf('%.1f', x), ...
    linspace(min(log_masses), max(log_masses), 5), 'UniformOutput', false);

xlabel('Observed Radius r_{obs} (pc)', 'FontSize', 22, 'FontWeight', 'bold');
ylabel('Mezzi Scale Factor', 'FontSize', 22, 'FontWeight', 'bold');
grid on; box on;
set(gca, 'FontSize', 11, 'LineWidth', 1.2);

insetAxes = axes('Position', [0.45, 0.19, 0.35, 0.30]);
hold(insetAxes, 'on');

% Plot Mezzi scale factor curves in the inset
for g = 1:numGalaxies
    r_norm = Results(g).r_obs / max(Results(g).r_obs);
    Mezzi = Results(g).Mezzi_th;
    this_color = mass_colors(color_idx(g), :);
    plot(insetAxes, r_norm, Mezzi, 'LineWidth', 1.2, 'Color', this_color);
end

% Configure inset axes
set(insetAxes, 'Box', 'on', 'LineWidth', 1.5, 'XGrid', 'on', 'YGrid', 'on', ...
    'XLim', [0, 1], 'FontSize', 8);
xlabel(insetAxes, 'Normalized Radius (r/r_{max})', 'FontSize', 12, 'FontWeight', 'bold');
ylabel(insetAxes, 'Mezzi Scale Factor', 'FontSize', 12, 'FontWeight', 'bold');

% Save Figure
exportgraphics(fig, fullfile(analysisDir, 'Fig_2_6_2_Coordinate_Scaling.png'), 'Resolution', 300);
close(fig);
%% --- 2.6.3 Figure: Total Mass Comparison (Publication Quality) ---
fig = figure('Color', 'white', 'Units', 'inches', 'Position', [1, 1, 7.2, 5.5], ...
    'PaperUnits', 'inches', 'PaperPosition', [0 0 7.2 5.5], ...
    'Renderer', 'painters');

% Main axes - adjusted position to accommodate inset
ax_main = axes('Position', [0.10, 0.12, 0.75, 0.82]);
hold on;

% === DATA PREPARATION ===
M_obs_all = [galaxySummary.M_total_obs_Msun];
M_true_all = [galaxySummary.M_total_true_Msun];
log_M_obs = log10(M_obs_all);
log_M_true = log10(M_true_all);
correction_factor = M_true_all ./ M_obs_all;

% Sort by correction magnitude
[correction_sorted, sort_idx] = sort(correction_factor, 'ascend');
log_M_obs_sorted = log_M_obs(sort_idx);
log_M_true_sorted = log_M_true(sort_idx);
numGalaxies = length(M_obs_all);

% === COLOR MAPPING ===
cmap = turbo(256);
corr_min = min(correction_sorted);
corr_max = max(correction_sorted);
corr_norm = (correction_sorted - corr_min) / (corr_max - corr_min);
color_indices = max(1, min(256, round(corr_norm * 255) + 1));
line_colors = cmap(color_indices, :);

% === PLOT CONNECTING LINES ===
for g = 1:numGalaxies
    if correction_sorted(g) > prctile(correction_sorted, 75)
        plot([g, g], [log_M_obs_sorted(g), log_M_true_sorted(g)], '-', ...
            'Color', [line_colors(g,:) 0.15], 'LineWidth', 6);
    end
end
for g = 1:numGalaxies
    plot([g, g], [log_M_obs_sorted(g), log_M_true_sorted(g)], '-', ...
        'Color', line_colors(g,:), 'LineWidth', 2.5);
end

% === PLOT DATA POINTS ===
marker_size = 70;
edge_color = [0.15 0.15 0.15];
edge_width = 1.2;
h_obs = scatter(1:numGalaxies, log_M_obs_sorted, marker_size, ...
    [0.95 0.95 0.95], 'o', 'filled', ...
    'MarkerEdgeColor', edge_color, 'LineWidth', edge_width, ...
    'DisplayName', 'Observed Mass');
h_true = scatter(1:numGalaxies, log_M_true_sorted, marker_size, ...
    line_colors, '^', 'filled', ...
    'MarkerEdgeColor', edge_color, 'LineWidth', edge_width, ...
    'DisplayName', 'Mezzi-Corrected Mass');

% === COLORBAR ===
colormap(ax_main, cmap);
cbar = colorbar('Location', 'eastoutside', 'Position', [0.88, 0.12, 0.03, 0.82]);
cbar.Label.String = 'M_{true}/M_{obs}';
cbar.Label.FontSize = 20;
cbar.Label.FontWeight = 'bold';
cbar.Label.Interpreter = 'tex';
caxis([corr_min, corr_max]);
if corr_max / corr_min < 10
    cbar.Ticks = linspace(corr_min, corr_max, 5);
else
    cbar.Ticks = [1, 2, 5, 10, 20, 50];
    cbar.Ticks = cbar.Ticks(cbar.Ticks >= corr_min & cbar.Ticks <= corr_max);
end
cbar.TickLabels = arrayfun(@(x) sprintf('%.1f×', x), cbar.Ticks, 'UniformOutput', false);
cbar.FontSize = 10;

% === STATISTICS ===
mean_true = mean(correction_sorted);
median_true = median(correction_sorted);
std_true = std(correction_sorted);
mad_true = median(abs(correction_sorted - median_true));

% === AXIS LABELS ===
xlabel('Galaxy (sorted by M_{true}/M_{obs})', ...
    'FontSize', 22, 'FontWeight', 'bold', 'Interpreter', 'tex');
ylabel('log_{10}(M_{obs}/M_{☉})', ...
    'FontSize', 22, 'FontWeight', 'bold', 'Interpreter', 'tex');

legend([h_obs, h_true], ...
    {'Observed Mass (M_{obs})', 'True (M_{true})'}, ...
    'Location', 'northwest', 'FontSize', 10, 'Box', 'on', ...
    'EdgeColor', [0.3 0.3 0.3], 'LineWidth', 1);

% === AXIS LIMITS ===
all_y_values = [log_M_obs_sorted(:); log_M_true_sorted(:)];
y_min = min(all_y_values);
y_max = max(all_y_values);
y_range = y_max - y_min;
y_pad = y_range * 0.08;
if numGalaxies <= 20
    x_ticks = 1:numGalaxies;
elseif numGalaxies <= 50
    x_ticks = 1:5:numGalaxies;
else
    x_ticks = round(linspace(1, numGalaxies, 10));
end
set(ax_main, 'FontSize', 10, 'LineWidth', 1.2, 'Box', 'on', ...
    'XLim', [-2, numGalaxies+2], ...
    'YLim', [y_min - y_pad, y_max + y_pad], ...
    'XTick', x_ticks, ...
    'XTickLabelRotation', 45, ...
    'YGrid', 'on', 'XGrid', 'off', 'GridAlpha', 0.15, ...
    'Layer', 'top', 'Color', 'white');

% === INSET: DISTRIBUTION ===
ax_inset = axes('Position', [0.52, 0.22, 0.28, 0.20]);
hold on;

% Histogram with more bins for finer granularity
n_bins = 20;
[counts, edges] = histcounts(correction_sorted, n_bins);

% Plot histogram with thinner bars
h_hist = histogram('BinEdges', edges, 'BinCounts', counts, ...
    'FaceColor', [0.5 0.55 0.6], ...
    'EdgeColor', [0.2 0.2 0.2], ...
    'LineWidth', 0.6, ...
    'FaceAlpha', 0.85);

grid on;
set(ax_inset, 'GridLineStyle', '-', 'GridColor', [0.85 0.85 0.85], 'GridAlpha', 0.8);

% === X-AXIS CONFIGURATION ===
xlabel('M_{true}/M_{obs}', 'FontSize', 11, 'Interpreter', 'tex', 'FontWeight', 'bold');

% Determine scale and set MORE x-ticks
if corr_max / corr_min > 5
    % Log scale for large dynamic range
    set(ax_inset, 'XScale', 'log');
    ax_inset.XAxis.Exponent = 0;
    
    % Generate MORE log-spaced ticks
    log_min = floor(log10(corr_min));
    log_max = ceil(log10(corr_max));
    
    tick_values = [];
    for exp = log_min:log_max
        tick_values = [tick_values, [1, 2, 3, 5, 7] * 10^exp];
    end
    
    % Filter to data range
    tick_values = tick_values(tick_values >= corr_min * 0.8 & tick_values <= corr_max * 1.2);
    
    % Limit to ~9 ticks for density but not clutter
    if length(tick_values) > 9
        step = ceil(length(tick_values) / 9);
        tick_values = tick_values(1:step:end);
    end
    ax_inset.XTick = tick_values;
    
    % Format labels cleanly
    ax_inset.XTickLabel = arrayfun(@(x) sprintf('%.0f', x), tick_values, 'UniformOutput', false);
else
    % Linear scale with MORE evenly spaced ticks
    x_range = corr_max - corr_min;
    n_ticks = 8;
    nice_step = x_range / n_ticks;
    
    % Round to nice number
    order = 10^floor(log10(nice_step));
    nice_step = ceil(nice_step / order) * order;
    if nice_step == 0, nice_step = order / 10; end
    x_start = ceil(corr_min / nice_step) * nice_step;
    x_end = floor(corr_max / nice_step) * nice_step;
    ax_inset.XTick = x_start:nice_step:x_end;
    
    if length(ax_inset.XTick) < 6
        fine_step = nice_step / 2;
        ax_inset.XTick = x_start:fine_step:x_end;
    end
end

% === Y-AXIS CONFIGURATION ===
ylabel('Count', 'FontSize', 11, 'FontWeight', 'bold');

% Smart y-ticks based on data
y_max = max(counts);
if y_max <= 5
    y_step = 1;
elseif y_max <= 10
    y_step = 2;
elseif y_max <= 20
    y_step = 5;
else
    y_step = ceil(y_max / 6 / 5) * 5;  % Aim for ~6 ticks
end
ax_inset.YTick = 0:y_step:y_max;
ax_inset.YLim = [0, y_max * 1.15];

% === STYLING ===
set(ax_inset, 'FontSize', 9, ...
    'LineWidth', 0.8, ...
    'Box', 'on', ...
    'Color', 'white', ...
    'XTickLabelRotation', 0, ...
    'TickLength', [0.03 0.025], ...
    'TickDir', 'out');

% === EXPORT ===
exportgraphics(fig, fullfile(analysisDir, 'Fig_2_6_3_Mass_Comparison.png'), 'Resolution', 600);
close(fig);
%% --- 2.6.4 Figure: Mezzi Velocity Residuals ---
% Extract Mezzi-corrected data (using original parameters)
v_flat_obs = [galaxySummary.v_flat_obs_kms];
v_flat_mezzi = [galaxySummary.v_flat_true_kms];

log_masses = log10(galaxy_masses);
numGalaxies = length(galaxy_names);

% Calculate Mezzi residuals
residual_mezzi = v_flat_mezzi - v_flat_obs;
log_v_obs = log10(v_flat_obs);
log_v_mezzi = log10(v_flat_mezzi);
log_residual_mezzi = log_v_mezzi - log_v_obs;

% Calculate dex statistics for ALL DATA
dex_scatter_mezzi = std(log_residual_mezzi);
mean_dex_error_mezzi = mean(abs(log_residual_mezzi));
rms_dex_mezzi = sqrt(mean(log_residual_mezzi.^2));
bias_mezzi = mean(log_residual_mezzi);

%% FIGURE 2.6.4: Mezzi Velocity Residuals (2-Panel, Clean Design)
fig = figure('Color', 'white', 'Position', [100, 100, 1400, 600]);
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact');

% Mass-based coloring
mass_colors = parula(numGalaxies);
[~, sort_idx] = sort(log_masses, 'ascend');
[~, color_idx] = sort(sort_idx);

% Pre-calculate statistics for both panels
% Panel (a) statistics - linear residuals
rms_mezzi = sqrt(mean(residual_mezzi.^2));
std_mezzi = std(residual_mezzi);
mean_bias_mezzi = mean(residual_mezzi);
median_bias_mezzi = median(residual_mezzi);
max_residual = max(abs(residual_mezzi));
min_residual = min(residual_mezzi);
max_residual_idx = find(abs(residual_mezzi) == max_residual, 1);
min_residual_idx = find(residual_mezzi == min_residual, 1);

% Panel (b) statistics - logarithmic residuals
log_v_obs = log10(v_flat_obs);
log_residual_mezzi = log10(v_flat_mezzi ./ v_flat_obs);
rms_dex_mezzi = sqrt(mean(log_residual_mezzi.^2));
dex_scatter_mezzi = std(log_residual_mezzi);
bias_mezzi = mean(log_residual_mezzi);
median_dex_bias = median(log_residual_mezzi);
max_dex_residual = max(abs(log_residual_mezzi));
min_dex_residual = min(log_residual_mezzi);

% Percentage of galaxies within tolerance bands
within_10_percent = sum(abs(log_residual_mezzi) < log10(1.10)) / numGalaxies * 100;
within_20_percent = sum(abs(log_residual_mezzi) < log10(1.20)) / numGalaxies * 100;
within_50_percent = sum(abs(log_residual_mezzi) < log10(1.50)) / numGalaxies * 100;

% Mass-stratified analysis
low_mass_mask = log_masses < 9.0;
mid_mass_mask = (log_masses >= 9.0) & (log_masses < 10.5);
high_mass_mask = log_masses >= 10.5;

% === PANEL (a): Velocity Residuals vs Observed Velocity ===
ax1 = nexttile;
hold on; box on;
for g = 1:numGalaxies
    scatter(v_flat_obs(g), residual_mezzi(g), 80, mass_colors(color_idx(g), :), ...
        'o', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 1);
end
yline(0, 'k--', 'LineWidth', 2, 'Alpha', 0.6);
yline(mean_bias_mezzi, 'r-', 'LineWidth', 2, 'Alpha', 0.9, 'Label', 'Mean');
yline(std_mezzi, 'r:', 'LineWidth', 1.5, 'Alpha', 0.5, 'Label', '+/-1 sigma');
yline(-std_mezzi, 'r:', 'LineWidth', 1.5, 'Alpha', 0.5);
xlabel('v_{obs} (km/s)', 'FontSize', 20, 'FontWeight', 'bold');
ylabel('v_{Newt,True} - v_{obs} (km/s)', 'FontSize', 20, 'FontWeight', 'bold');
% title('(a) Linear Velocity Residuals', 'FontSize', 13, 'FontWeight', 'bold');
set(ax1, 'FontSize', 12, 'LineWidth', 1.2);
grid on;
stats_text_a = {
    'Linear Statistics:'
    sprintf('N = %d galaxies', numGalaxies)
    sprintf('RMS residual = %.1f km/s', rms_mezzi)
    sprintf('Std dev = %.1f km/s', std_mezzi)
    sprintf('Mean bias = %+.1f km/s', mean_bias_mezzi)
    sprintf('Median bias = %+.1f km/s', median_bias_mezzi)
};
annotation('textbox', [0.07, 0.72, 0.22, 0.22], 'String', stats_text_a, ...
    'FontSize', 12, 'BackgroundColor', [1 1 1 0.9], 'EdgeColor', [0.3 0.3 0.3], ...
    'LineWidth', 1, 'FitBoxToText', 'on', ...
    'HorizontalAlignment', 'left', 'VerticalAlignment', 'top');

% === PANEL (b): Log Residuals (Dex) ===
ax2 = nexttile;
hold on; box on;
for g = 1:numGalaxies
    scatter(log_v_obs(g), log_residual_mezzi(g), 80, mass_colors(color_idx(g), :), ...
        'o', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 1);
end
yline(0, 'k--', 'LineWidth', 2, 'Alpha', 0.6);
yline(bias_mezzi, 'b-', 'LineWidth', 2, 'Alpha', 0.9, 'Label', 'Mean');
yline(dex_scatter_mezzi, 'b:', 'LineWidth', 1.5, 'Alpha', 0.5, 'Label', '+/-1 sigma');
yline(-dex_scatter_mezzi, 'b:', 'LineWidth', 1.5, 'Alpha', 0.5);
xlabel('log_{10}(v_{obs})', 'FontSize', 20, 'FontWeight', 'bold');
ylabel('log(v_{Newt,True}/v_{obs}) (dex)', 'FontSize', 20, 'FontWeight', 'bold');
% title('(b) Logarithmic Residuals (Dex)', 'FontSize', 13, 'FontWeight', 'bold');
set(ax2, 'FontSize', 12, 'LineWidth', 1.2);
grid on;
cbar = colorbar('Location', 'eastoutside');
cbar.Label.String = 'log_{10}(M_{obs}/M_{☉})';
cbar.Label.FontSize = 12;
cbar.Label.FontWeight = 'bold';
cbar.Ticks = linspace(0, 1, 5);
cbar.TickLabels = arrayfun(@(x) sprintf('%.1f', x), ...
    linspace(min(log_masses), max(log_masses), 5), 'UniformOutput', false);
stats_text_b = {
    'Logarithmic Statistics:'
    sprintf('N = %d galaxies', numGalaxies)
    sprintf('RMS = %.3f dex (%.1f%%)', rms_dex_mezzi, (10^rms_dex_mezzi - 1)*100)
    sprintf('Scatter = %.3f dex (%.1f%%)', dex_scatter_mezzi, (10^dex_scatter_mezzi - 1)*100)
    sprintf('Mean bias = %+.3f dex (%.1f%%)', bias_mezzi, (10^bias_mezzi - 1)*100)
    sprintf('Median bias = %+.3f dex', median_dex_bias)
    sprintf('Within 10%%: %.1f%% of galaxies', within_10_percent)
    sprintf('Within 20%%: %.1f%% of galaxies', within_20_percent)
};
annotation('textbox', [0.53, 0.18, 0.20, 0.24], 'String', stats_text_b, ...
    'FontSize', 12, 'BackgroundColor', [1 1 1 0.9], 'EdgeColor', [0.3 0.3 0.3], ...
    'LineWidth', 1, 'FitBoxToText', 'on', ...
    'HorizontalAlignment', 'left', 'VerticalAlignment', 'top');
drawnow;
exportgraphics(fig, fullfile(analysisDir, 'Fig_2_6_4_Mezzi_Residuals.png'), 'Resolution', 300);
close(fig);


%% --- 2.6.5 Figure: Mezzi Dex Error Distribution ---
fig = figure('Color', 'white', 'Position', [100, 100, 1600, 700]);
t2 = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact');

% === PANEL (a): Histogram of Log Residuals ===
ax3 = nexttile;
hold on; box on;
bin_edges = linspace(min(log_residual_mezzi) - 0.03, max(log_residual_mezzi) + 0.03, 12);
h = histogram(log_residual_mezzi, bin_edges, 'FaceColor', [0.25 0.55 0.35], ...
    'EdgeColor', 'k', 'LineWidth', 2, 'FaceAlpha', 0.75, 'DisplayName', 'Data');
x_norm = linspace(min(bin_edges), max(bin_edges), 200);
y_norm = normpdf(x_norm, bias_mezzi, dex_scatter_mezzi);
y_norm_scaled = y_norm * sum(h.Values) * h.BinWidth;
plot(x_norm, y_norm_scaled, 'r-', 'LineWidth', 3, 'DisplayName', 'Gaussian fit');
xline(0, 'k--', 'LineWidth', 2.5, 'Alpha', 0.6, 'HandleVisibility', 'off');
xline(bias_mezzi, 'b-', 'LineWidth', 2.5, 'Alpha', 0.9, 'Label', sprintf('Mean: %.3f', bias_mezzi), 'HandleVisibility', 'off');
xline(bias_mezzi + dex_scatter_mezzi, 'b:', 'LineWidth', 2, 'Alpha', 0.5, 'HandleVisibility', 'off');
xline(bias_mezzi - dex_scatter_mezzi, 'b:', 'LineWidth', 2, 'Alpha', 0.5, 'HandleVisibility', 'off');
xlabel('log(v_{Newton,true}​ ) - log(v_{obs})', 'FontSize', 24, 'FontWeight', 'bold');
ylabel('Number of Galaxies', 'FontSize', 24, 'FontWeight', 'bold');
set(ax3, 'FontSize', 14, 'LineWidth', 1.5, 'Layer', 'top');
grid on;
normality_text = {
    sprintf('Distribution:')
    sprintf('N = %d', numGalaxies)
    sprintf('Mean = %.3f dex', bias_mezzi)
    sprintf('Median = %.3f dex', median(log_residual_mezzi))
    sprintf('Skewness = %.2f', skewness(log_residual_mezzi))
};
annotation('textbox', [0.08, 0.77, 0.22, 0.18], 'String', normality_text, ...
    'FontSize', 12, 'BackgroundColor', [1 1 1 0.9], 'EdgeColor', [0.3 0.3 0.3], ...
    'LineWidth', 1.5, 'FitBoxToText', 'on', ...
    'HorizontalAlignment', 'left', 'VerticalAlignment', 'top');

% === PANEL (b): Cumulative Absolute Error ===
ax4 = nexttile;
hold on; box on;
abs_log_residual = abs(log_residual_mezzi);
[sort_abs, ~] = sort(abs_log_residual);
cdf_y = (1:numGalaxies) / numGalaxies;
plot(sort_abs, cdf_y, 'b-', 'LineWidth', 4, 'DisplayName', 'Cumulative');
p50_idx = find(cdf_y >= 0.5, 1, 'first');
p68_idx = find(cdf_y >= 0.68, 1, 'first');
p95_idx = find(cdf_y >= 0.95, 1, 'first');
xline(sort_abs(p50_idx), 'g--', 'LineWidth', 2.5, 'Label', sprintf('50%%: %.3f dex', sort_abs(p50_idx)));
xline(sort_abs(p68_idx), 'r--', 'LineWidth', 2.5, 'Label', sprintf('68%%: %.3f dex', sort_abs(p68_idx)));
xline(sort_abs(p95_idx), 'm--', 'LineWidth', 2.5, 'Label', sprintf('95%%: %.3f dex', sort_abs(p95_idx)));
xlabel('|log(v_{Newton,true}​ ) - log(v_{obs})|', 'FontSize', 24, 'FontWeight', 'bold');
ylabel('Cumulative Fraction', 'FontSize', 24, 'FontWeight', 'bold');
set(ax4, 'FontSize', 14, 'LineWidth', 1.5, 'Layer', 'top');
grid on;
xlim([0 max(sort_abs) * 1.05]);
ylim([0 1.05]);
percentile_text = {
    sprintf('Percentiles:')
    sprintf('N = %d', numGalaxies)
    sprintf('50%% (median) = %.3f dex', sort_abs(p50_idx))
    sprintf('68%% (1 sigma) = %.3f dex', sort_abs(p68_idx))
    sprintf('95%% (2 sigma) = %.3f dex', sort_abs(p95_idx))
    sprintf('Mean |error| = %.3f dex', mean_dex_error_mezzi)
};
annotation('textbox', [0.77, 0.2, 0.22, 0.18], 'String', percentile_text, ...
    'FontSize', 12, 'BackgroundColor', [1 1 1 0.9], 'EdgeColor', [0.3 0.3 0.3], ...
    'LineWidth', 1.5, 'FitBoxToText', 'on', ...
    'HorizontalAlignment', 'left', 'VerticalAlignment', 'top');
drawnow;
exportgraphics(fig, fullfile(analysisDir, 'Fig_2_6_5_Mezzi_DexError.png'), 'Resolution', 300);
close(fig);
%% --- 2.7 Inner Region Analysis (r < 0.3) ---
innerDir = analysisDir;
InnerAnalysis = struct();
core_fraction = 0.3;

for g = 1:numGalaxies
    r_obs = Results(g).r_obs(:);
    v_obs = Results(g).v_obs(:);
    v_newt_obs = Results(g).v_Newton_obs(:);
    v_newt_mezzi = Results(g).v_Newton_true_To_obs(:);
    N = min([length(r_obs), length(v_obs), length(v_newt_obs), length(v_newt_mezzi)]);
    r_obs = r_obs(1:N);
    v_obs = v_obs(1:N);
    v_newt_obs = v_newt_obs(1:N);
    v_newt_mezzi = v_newt_mezzi(1:N);
    valid = (r_obs > 0) & (v_obs > 0) & isfinite(v_obs) & isfinite(v_newt_obs) & isfinite(v_newt_mezzi);
    
    if sum(valid) < 10
        InnerAnalysis(g).valid = false;
        fprintf('[SKIP] Galaxy %d: insufficient valid points\n', g);
        continue;
    end
    
    r = r_obs(valid);
    v = v_obs(valid);
    v_n = v_newt_obs(valid);
    v_m = v_newt_mezzi(valid);
    r_max = max(r);
    r_core = core_fraction * r_max;
    core_mask = r < r_core;
    outer_mask_g = r >= r_core;
    res_naive = v - v_n;
    res_mezzi = v - v_m;
    stats = struct();
    
    if sum(core_mask) > 3
        stats.N_core = sum(core_mask);
        stats.rms_naive_core = sqrt(mean(res_naive(core_mask).^2));
        stats.rms_mezzi_core = sqrt(mean(res_mezzi(core_mask).^2));
        stats.bias_naive_core = mean(res_naive(core_mask));
        stats.bias_mezzi_core = mean(res_mezzi(core_mask));
        stats.std_naive_core = std(res_naive(core_mask));
        stats.std_mezzi_core = std(res_mezzi(core_mask));
    else
        stats.N_core = 0;
    end
    
    if sum(outer_mask_g) > 3
        stats.N_outer = sum(outer_mask_g);
        stats.rms_naive_outer = sqrt(mean(res_naive(outer_mask_g).^2));
        stats.rms_mezzi_outer = sqrt(mean(res_mezzi(outer_mask_g).^2));
        stats.bias_naive_outer = mean(res_naive(outer_mask_g));
        stats.bias_mezzi_outer = mean(res_mezzi(outer_mask_g));
        stats.std_naive_outer = std(res_naive(outer_mask_g));
        stats.std_mezzi_outer = std(res_mezzi(outer_mask_g));
    else
        stats.N_outer = 0;
    end
    
    InnerAnalysis(g).valid = true;
    InnerAnalysis(g).galaxy = galaxy_names{g};
    InnerAnalysis(g).stats = stats;
    InnerAnalysis(g).r = r;
    InnerAnalysis(g).v = v;
    InnerAnalysis(g).v_newt = v_n;
    InnerAnalysis(g).v_mezzi = v_m;
    InnerAnalysis(g).res_naive = res_naive;
    InnerAnalysis(g).res_mezzi = res_mezzi;
    InnerAnalysis(g).core_mask = core_mask;
    InnerAnalysis(g).outer_mask = outer_mask_g;
    InnerAnalysis(g).r_core = r_core;
    
    if stats.N_core > 3
        fprintf('[Rotation] Galaxy %d/%d: %s\n', g, numGalaxies, galaxy_names{g});
%         fprintf('[Results] Galaxy %d RMS Core - Naive: %.2f km/s | Mezzi: %.2f km/s\n', ...
%             g, stats.rms_naive_core/1e3, stats.rms_mezzi_core/1e3);
    end
end

%% Ensemble Inner Region Statistics
valid_galaxies = find(arrayfun(@(x) isfield(x,'valid') && x.valid, InnerAnalysis));
N_valid = length(valid_galaxies);

if N_valid == 0
    fprintf('[ERROR] No valid galaxies for inner region analysis\n');
else
    rms_naive_core = [];
    rms_mezzi_core = [];
    rms_naive_outer = [];
    rms_mezzi_outer = [];
    bias_naive_core = [];
    bias_mezzi_core = [];
    bias_naive_outer = [];
    bias_mezzi_outer = [];
    
    for idx = 1:N_valid
        g = valid_galaxies(idx);
        s = InnerAnalysis(g).stats;
        if s.N_core > 3
            rms_naive_core(end+1) = s.rms_naive_core;
            rms_mezzi_core(end+1) = s.rms_mezzi_core;
            bias_naive_core(end+1) = s.bias_naive_core;
            bias_mezzi_core(end+1) = s.bias_mezzi_core;
        end
        if s.N_outer > 3
            rms_naive_outer(end+1) = s.rms_naive_outer;
            rms_mezzi_outer(end+1) = s.rms_mezzi_outer;
            bias_naive_outer(end+1) = s.bias_naive_outer;
            bias_mezzi_outer(end+1) = s.bias_mezzi_outer;
        end
    end
    
    %% Figure 2.7.1: Bias Distribution (Inner vs Outer, Naive vs Mezzi)
    fig = figure('Color','white','Position',[100 100 700 600], ...
        'Name','Fig_2_7_1_Bias_Comparison');
    ax = axes;
    hold on; box on;
    r_th = core_fraction;
    
    % Convert to km/s and store in persistent variables
    bias_n_core_kms = (bias_naive_core / 1e3).';
    bias_m_core_kms = (bias_mezzi_core / 1e3).';
    bias_n_outer_kms = (bias_naive_outer / 1e3).';
    bias_m_outer_kms = (bias_mezzi_outer / 1e3).';
    
    % Calculate statistics for annotation
    median_n_c = median(bias_n_core_kms);
    median_m_c = median(bias_m_core_kms);
    median_n_o = median(bias_n_outer_kms);
    median_m_o = median(bias_m_outer_kms);
    rms_n_c = mean(rms_naive_core) / 1e3;
    rms_m_c = mean(rms_mezzi_core) / 1e3;
    rms_n_o = mean(rms_naive_outer) / 1e3;
    rms_m_o = mean(rms_mezzi_outer) / 1e3;
    
    % Calculate IQR for display
    iqr_n_c = prctile(bias_n_core_kms,75) - prctile(bias_n_core_kms,25);
    iqr_m_c = prctile(bias_m_core_kms,75) - prctile(bias_m_core_kms,25);
    iqr_n_o = prctile(bias_n_outer_kms,75) - prctile(bias_n_outer_kms,25);
    iqr_m_o = prctile(bias_m_outer_kms,75) - prctile(bias_m_outer_kms,25);
    
    naive_color = [0.85 0.325 0.098];
    mezzi_color = [0.466 0.674 0.188];
    offset = 0.15;
    
    h_n_core = boxplot(bias_n_core_kms,'Positions',1-offset,'Widths',0.22,...
        'Colors','k','Symbol','+');
    h_m_core = boxplot(bias_m_core_kms,'Positions',1+offset,'Widths',0.22,...
        'Colors','k','Symbol','+');
    h_n_outer = boxplot(bias_n_outer_kms,'Positions',2-offset,'Widths',0.22,...
        'Colors','k','Symbol','+');
    h_m_outer = boxplot(bias_m_outer_kms,'Positions',2+offset,'Widths',0.22,...
        'Colors','k','Symbol','+');
    
    positions = [1-offset, 1+offset, 2-offset, 2+offset];
    is_mezzi = [false, true, false, true];
    
    h_boxes = findobj(gca,'Tag','Box');
    h_whiskers = findobj(gca,'Tag','Whisker');
    h_medians = findobj(gca,'Tag','Median');
    [~, sort_idx] = sort(cellfun(@(x) mean(x), get(h_boxes,'XData')));
    h_boxes = h_boxes(sort_idx);
    
    for j = 1:length(h_boxes)
        xmid = mean(get(h_boxes(j),'XData'));
        if is_mezzi(j)
            c = mezzi_color;
        else
            c = naive_color;
        end
        patch(get(h_boxes(j),'XData'),get(h_boxes(j),'YData'), ...
            c,'FaceAlpha',0.35,'EdgeColor',c,'LineWidth',1.5);
        for w = 1:length(h_whiskers)
            wx = mean(get(h_whiskers(w),'XData'));
            if abs(wx - xmid) < 0.1
                set(h_whiskers(w),'Color',c,'LineWidth',1.5);
            end
        end
        for m = 1:length(h_medians)
            mx = mean(get(h_medians(m),'XData'));
            if abs(mx - xmid) < 0.1
                set(h_medians(m),'Color','k','LineWidth',2);
            end
        end
    end
    
    h_outliers = findobj(gca,'Tag','Outliers');
    for k = 1:length(h_outliers)
        xval = mean(get(h_outliers(k),'XData'));
        [~, pos_idx] = min(abs(positions - xval));
        if is_mezzi(pos_idx)
            set(h_outliers(k),'MarkerEdgeColor',mezzi_color, ...
                'MarkerFaceColor',mezzi_color);
        else
            set(h_outliers(k),'MarkerEdgeColor',naive_color, ...
                'MarkerFaceColor',naive_color);
        end
    end
    
    yline(0,'k--','LineWidth',2,'Alpha',0.6);
    set(ax,'XTick',[1 2], ...
        'XTickLabel',{sprintf('Inner (r < %.2f)',r_th), ...
        sprintf('Outer (r > %.2f)',r_th)}, ...
        'FontSize',12,'LineWidth',1.2);
    ylabel('Bias $\langle \mathrm{v}_{obs} - \mathrm{v}_{Newton,True} \rangle$ [km/s]', ...
        'Interpreter','latex','FontSize',14);
    grid on;
    
    h_leg_naive = plot(nan,nan,'s','MarkerFaceColor',naive_color,...
        'MarkerEdgeColor',naive_color,'MarkerSize',10);
    h_leg_mezzi = plot(nan,nan,'s','MarkerFaceColor',mezzi_color,...
        'MarkerEdgeColor',mezzi_color,'MarkerSize',10);
    legend([h_leg_naive, h_leg_mezzi],{'Naive Newtonian','True Frame Newtonian'},...
        'Location','best','FontSize',11);
    
    % Shows MEDIAN (center line) with IQR and RMS
    stats_text = {
        'Median Bias [IQR] (RMS):'
        sprintf('Inner Naive: %+.2f [%.2f] (%.2f) km/s', median_n_c, iqr_n_c, rms_n_c)
        sprintf('Inner Mezzi: %+.2f [%.2f] (%.2f) km/s', median_m_c, iqr_m_c, rms_m_c)
        sprintf('Outer Naive: %+.2f [%.2f] (%.2f) km/s', median_n_o, iqr_n_o, rms_n_o)
        sprintf('Outer Mezzi: %+.2f [%.2f] (%.2f) km/s', median_m_o, iqr_m_o, rms_m_o)
    };
    annotation('textbox',[0.15 0.08 0.38 0.24], ...
        'String',stats_text,'FontSize',10,...
        'BackgroundColor',[1 1 1 0.9],...
        'EdgeColor',[0.3 0.3 0.3],...
        'LineWidth',1,'FitBoxToText','on');
    
    drawnow;
    exportgraphics(fig,fullfile(innerDir,'Fig_2_7_1_Bias_Comparison.png'),'Resolution',300);
    close(fig);
end

%% --- 2.8 Save Summary Data and Final Output ---
save('Galaxy_Summary.mat', 'galaxySummary', 'InnerAnalysis');

%% ===== SECTION 3: COMPREHENSIVE SCALING ANALYSIS =====
ScalingDir = fullfile(baseDir, 'ScalingAnalysis');
if ~exist(ScalingDir, 'dir')
    mkdir(ScalingDir);
end

%% --- 3.1 Rotation Curve Analysis ---
% Initialize arrays to store statistics for all galaxies
RC_stats = struct();
all_rms_naive = zeros(numGalaxies, 1);
all_rms_true = zeros(numGalaxies, 1);
all_R2_naive = zeros(numGalaxies, 1);
all_R2_true = zeros(numGalaxies, 1);

for g = 1:numGalaxies
    % Extract velocity profiles
    v_obs = Results(g).v_obs;
    v_newt = Results(g).v_Newton_obs;
    v_true = Results(g).v_Newton_true_To_obs;
    
    % Calculate R-squared (coefficient of determination)
    SS_tot = sum((v_obs - mean(v_obs)).^2);
    R2_naive = 1 - sum((v_obs - v_newt).^2) / SS_tot;
    R2_true = 1 - sum((v_obs - v_true).^2) / SS_tot;
    
    % Calculate reduced chi-squared (assuming 5% velocity errors)
    sigma_v = 0.05 * v_obs;
    chi2_naive = sum(((v_obs - v_newt) ./ sigma_v).^2) / (length(v_obs) - 1);
    chi2_true = sum(((v_obs - v_true) ./ sigma_v).^2) / (length(v_obs) - 1);
    
    % Store statistics in struct
    RC_stats(g).R2_naive = R2_naive;
    RC_stats(g).R2_true = R2_true;
    RC_stats(g).chi2_naive = chi2_naive;
    RC_stats(g).chi2_true = chi2_true;
    RC_stats(g).v_max_obs = max(v_obs) / 1e3;
    RC_stats(g).v_max_newt = max(v_newt) / 1e3;
    RC_stats(g).v_max_true = max(v_true) / 1e3;
    
    % Store in arrays for ensemble analysis
    all_rms_naive(g) = Results(g).rms_naive / 1e3;
    all_rms_true(g) = Results(g).rms_corrected / 1e3;
    all_R2_naive(g) = R2_naive;
    all_R2_true(g) = R2_true;
end

%% --- 3.2 Enclosed Mass Correction Analysis ---
total_mass_obs = zeros(numGalaxies, 1);
total_mass_true = zeros(numGalaxies, 1);
mass_concentration_obs = zeros(numGalaxies, 1);
mass_concentration_true = zeros(numGalaxies, 1);

for g = 1:numGalaxies
    % Extract radius and mass profiles
    r_obs = Results(g).r_obs;
    r_true = Results(g).r_true;
    M_obs = Results(g).M_enc_obs;
    M_true = Results(g).M_enc_true;
    
    % Calculate mass within half-light radius proxy (r_max/2)
    r_half = r_obs(end) / 2;
    [~, idx_half_obs] = min(abs(r_obs - r_half));
    [~, idx_half_true] = min(abs(r_true - r_half));
    M_half_obs = M_obs(idx_half_obs);
    M_half_true = M_true(idx_half_true);
    
    % Calculate concentration parameter (M_half / M_total)
    C_obs = M_half_obs / M_obs(end);
    C_true = M_half_true / M_true(end);
    
    % Store results
    mass_scaling(g).M_total_obs = M_obs(end) / Msun;
    mass_scaling(g).M_total_true = M_true(end) / Msun;
    mass_scaling(g).M_half_obs = M_half_obs / Msun;
    mass_scaling(g).M_half_true = M_half_true / Msun;
    mass_scaling(g).concentration_obs = C_obs;
    mass_scaling(g).concentration_true = C_true;
    mass_scaling(g).mass_ratio = M_true(end) / M_obs(end);
    total_mass_obs(g) = M_obs(end) / Msun;
    total_mass_true(g) = M_true(end) / Msun;
    mass_concentration_obs(g) = C_obs;
    mass_concentration_true(g) = C_true;
end

%% --- 3.3 Deflection Angle Analysis ---
deflection_stats = struct();
max_alpha_obs = zeros(numGalaxies, 1);
max_alpha_true = zeros(numGalaxies, 1);
 

for g = 1:numGalaxies
    % Extract deflection angles (convert to arcseconds)
    alpha_obs = Results(g).alpha_rad_obs / arcsec;
    alpha_true = Results(g).alpha_rad_true / arcsec;
    
    % Find maximum deflection
    [max_alpha_obs(g), idx_max] = max(alpha_obs);
    max_alpha_true(g) = max(alpha_true);
    
    % Store results
    deflection_stats(g).max_alpha_obs_arcsec = max_alpha_obs(g);
    deflection_stats(g).max_alpha_true_arcsec = max_alpha_true(g);
    deflection_stats(g).alpha_ratio = max_alpha_true(g) / max_alpha_obs(g);
    
%     fprintf('[RESULT] %s: max α_obs=%.3f", max α_true=%.3f"\n', ...
%         galaxy_names{g}, max_alpha_obs(g), max_alpha_true(g));
end


%% --- 3.4 Surface Density Analysis ---
surface_stats = struct();
Sigma_central_obs = zeros(numGalaxies, 1);
Sigma_central_true = zeros(numGalaxies, 1);
Sigma_scale_obs = zeros(numGalaxies, 1);
Sigma_scale_true = zeros(numGalaxies, 1);

for g = 1:numGalaxies
    % Extract surface density profiles
    R_obs = Results(g).R_grid_obs;
    R_true = Results(g).R_grid_true;
    Sigma_obs = Results(g).Sigma_obs * pc^2 / Msun;
    Sigma_true = Results(g).Sigma_true * pc^2 / Msun;
    
    % Central surface density (first valid point)
    Sigma_central_obs(g) = Sigma_obs(1);
    Sigma_central_true(g) = Sigma_true(1);
    
    % Scale radius where Sigma drops to half central
    Sigma_half_obs = Sigma_obs(1) / 2;
    Sigma_half_true = Sigma_true(1) / 2;
    [~, idx_half_obs] = min(abs(Sigma_obs - Sigma_half_obs));
    [~, idx_half_true] = min(abs(Sigma_true - Sigma_half_true));
    R_scale_obs = R_obs(idx_half_obs) / pc;
    R_scale_true = R_true(idx_half_true) / pc;
    Sigma_scale_obs(g) = R_scale_obs;
    Sigma_scale_true(g) = R_scale_true;
    
    % Integrated surface density within effective radius
    R_eff = R_obs(end) / 2;
    [~, idx_eff] = min(abs(R_obs - R_eff));
    Sigma_int_obs = trapz(R_obs(1:idx_eff), Sigma_obs(1:idx_eff)) / 1e6;
    Sigma_int_true = trapz(R_true(1:idx_eff), Sigma_true(1:idx_eff)) / 1e6;
    
    % Store results
    surface_stats(g).Sigma0_obs = Sigma_central_obs(g);
    surface_stats(g).Sigma0_true = Sigma_central_true(g);
    surface_stats(g).R_scale_obs_pc = R_scale_obs;
    surface_stats(g).R_scale_true_pc = R_scale_true;
    surface_stats(g).Sigma_ratio = Sigma_central_true(g) / Sigma_central_obs(g);
    
%     fprintf('[RESULT] %s: Σ0_obs=%.2e, Σ0_true=%.2e Msun/pc², R_scale_obs=%.1f, R_scale_true=%.1f pc\n', ...
%         galaxy_names{g}, Sigma_central_obs(g), Sigma_central_true(g), R_scale_obs, R_scale_true);
end

%% --- 3.5 Individual Publication Figures ---
% Extract improvement percentages from struct array
improvement_percents = arrayfun(@(x) x.improvement_percent, galaxySummary);
[~, sig_idx] = sort(improvement_percents, 'descend');
significant_galaxies = sig_idx(1:min(3, numGalaxies));
sig_names = {galaxy_names{significant_galaxies}};

%% Calculate Outer Region Statistics (r > 0.3 r_max)  
% Initialize outer region arrays
all_rms_true_outer = zeros(numGalaxies, 1);
all_R2_true_outer = zeros(numGalaxies, 1);
chi2_true_outer = zeros(numGalaxies, 1);
outer_fraction = zeros(numGalaxies, 1);
num_outer_points = zeros(numGalaxies, 1);

for g = 1:numGalaxies
    % Get data
    r_norm = Results(g).r_obs / max(Results(g).r_obs);
    v_obs = Results(g).v_obs;
    v_true = Results(g).v_Newton_true_To_obs;
    
    % Find outer region indices (r > 0.3 r_max)
    outer_idx = find(r_norm > 0.3);
    num_outer_points(g) = length(outer_idx);
    outer_fraction(g) = length(outer_idx) / length(r_norm);
    
    if length(outer_idx) >= 3
        % Outer region velocities
        v_obs_outer = v_obs(outer_idx);
        v_true_outer = v_true(outer_idx);
        
        % Calculate RMS for outer region (in km/s)
        residual_true_outer = v_obs_outer - v_true_outer;
        all_rms_true_outer(g) = sqrt(mean(residual_true_outer.^2)) / 1e3;
        
        % R-squared for outer region
        SS_tot_outer = sum((v_obs_outer - mean(v_obs_outer)).^2);
        if SS_tot_outer > 0
            all_R2_true_outer(g) = 1 - sum(residual_true_outer.^2) / SS_tot_outer;
        else
            all_R2_true_outer(g) = NaN;
        end
        
        % Reduced chi-squared
        sigma_v_outer = 0.05 * v_obs_outer;
        chi2_true_outer(g) = sum(((v_obs_outer - v_true_outer) ./ sigma_v_outer).^2) / (length(outer_idx) - 1);
    else
        fprintf('[WARN] %s has only %d points in outer region\n', galaxy_names{g}, length(outer_idx));
        all_rms_true_outer(g) = NaN;
        all_R2_true_outer(g) = NaN;
        chi2_true_outer(g) = NaN;
    end
end

%% Figure 3.5.1: Observed vs Mezzi Deflection (Mass-Colored, Extremes Labeled)
fig = figure('Color','white','Position',[100 100 750 650]);
hold on;

% Mass coloring (same as Figure 2)

log_masses = log10(galaxy_masses);
mass_colors = parula(numGalaxies);
[~, sort_idx] = sort(log_masses, 'ascend');
[~, color_idx] = sort(sort_idx);

% 1:1 Reference line (use min/max of actual data)
min_val = min([min(max_alpha_obs), min(max_alpha_true)]) * 0.95;
max_val = max([max(max_alpha_obs), max(max_alpha_true)]) * 1.05;
plot([min_val max_val], [min_val max_val], 'k--', 'LineWidth', 2, ...
    'DisplayName', '1:1 line');

% Scatter plot
scatter(max_alpha_obs, max_alpha_true, 100, mass_colors(color_idx, :), ...
    'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 0.8, ...
    'DisplayName', 'Galaxies');

% Identify and label extremes
[~, idx_obs_max] = max(max_alpha_obs);
[~, idx_mezzi_max] = max(max_alpha_true);

% Smart label positioning
x_range = range(max_alpha_obs);
y_range = range(max_alpha_true);
offset_x = 0.6 * x_range;
offset_y = 0.03 * y_range;

% Label with background for readability
text(max_alpha_obs(idx_obs_max) + offset_x, max_alpha_true(idx_obs_max), ...
    galaxy_names{idx_obs_max}, 'FontSize', 10, 'FontWeight', 'bold', ...
    'HorizontalAlignment', 'right', 'BackgroundColor', 'none', 'EdgeColor', 'none');
text(max_alpha_obs(idx_mezzi_max) + offset_x, max_alpha_true(idx_mezzi_max) + offset_y, ...
    galaxy_names{idx_mezzi_max}, 'FontSize', 10, 'FontWeight', 'bold', ...
    'HorizontalAlignment', 'right', 'BackgroundColor', 'none', 'EdgeColor', 'none');

% Colorbar
colormap(parula);
cbar = colorbar('Location', 'eastoutside');
cbar.Label.String = 'log_{10}(M_{obs}/M_{☉})';
cbar.Label.FontSize = 12;
cbar.Label.FontWeight = 'bold';
cbar.Ticks = linspace(0, 1, 5);
cbar.TickLabels = arrayfun(@(x) sprintf('%.1f', x), ...
    linspace(min(log_masses), max(log_masses), 5), 'UniformOutput', false);

% Axis formatting
xlabel('Observed $\alpha_{\max}$ (arcsec)','Interpreter','latex','FontSize',20);
ylabel('True $\alpha_{\max}$ (arcsec)','Interpreter','latex','FontSize',20);
xlim([min_val max_val]);
ylim([min_val max_val]);
axis square;
grid on;
box on;
set(gca, 'FontSize', 12, 'LineWidth', 1.3);

% Save Figure
exportgraphics(fig, fullfile(ScalingDir, 'Fig_3_5_1_Deflection_Correlation.png'), 'Resolution', 400);
close(fig);
%% --- 3.6 Scaling Factor Impact on Galaxy Surface Density ---
%% Data Preparation & Derived Quantities
% Extract core data
Sigma0_obs = [surface_stats.Sigma0_obs];
Sigma0_true = [surface_stats.Sigma0_true];
R_s_obs = [surface_stats.R_scale_obs_pc];
R_s_true = [surface_stats.R_scale_true_pc];

log_masses = log10(galaxy_masses);

% Calculate primary ratios
Ratio_Sigma = Sigma0_true ./ Sigma0_obs;
Ratio_Radius = R_s_true ./ R_s_obs;

% Calculate derived physical quantities
Ratio_Sigma_predicted = (R_s_obs ./ R_s_true).^2;
Consistency_residual = Ratio_Sigma ./ Ratio_Sigma_predicted;

% Pre-calculate global statistics
global_stats.mean_Sigma = mean(Ratio_Sigma);
global_stats.std_Sigma = std(Ratio_Sigma);
global_stats.mean_Radius = mean(Ratio_Radius);
global_stats.std_Radius = std(Ratio_Radius);
global_stats.mean_Consistency = mean(Consistency_residual);
global_stats.std_Consistency = std(Consistency_residual);

%% Mass Regime Definition & Stratified Analysis
% Define mass bins based on data distribution
mass_bins.edges = [7.0, 9.0, 10.5, 12.0];
mass_bins.labels = {'Low-Mass', 'Intermediate', 'High-Mass'};
mass_bins.colors = [0.8 0.4 0.2; 0.2 0.6 0.8; 0.4 0.2 0.6];

% Initialize regime statistics
regime_stats = struct();
for r = 1:length(mass_bins.labels)
    % Find galaxies in this mass range
    in_regime = log_masses >= mass_bins.edges(r) & log_masses < mass_bins.edges(r+1);
    idx = find(in_regime);
    if isempty(idx)
        regime_stats(r).count = 0;
        regime_stats(r).mean_Consistency = NaN;
        regime_stats(r).std_Consistency = NaN;
        continue;
    end
    
    % Calculate statistics for this regime
    regime_stats(r).count = length(idx);
    regime_stats(r).mass_range = [min(log_masses(idx)), max(log_masses(idx))];
    
    % Radius statistics
    regime_stats(r).mean_Radius_ratio = mean(Ratio_Radius(idx));
    regime_stats(r).std_Radius_ratio = std(Ratio_Radius(idx));
    regime_stats(r).median_Radius_ratio = median(Ratio_Radius(idx));
    
    % Density statistics
    regime_stats(r).mean_Sigma_ratio = mean(Ratio_Sigma(idx));
    regime_stats(r).std_Sigma_ratio = std(Ratio_Sigma(idx));
    regime_stats(r).median_Sigma_ratio = median(Ratio_Sigma(idx));
    
    % Correlation within regime
    if regime_stats(r).count > 2
        [R_reg, P_reg] = corrcoef(log_masses(idx), Ratio_Radius(idx));
        regime_stats(r).corr_Radius_Mass = R_reg(1,2);
        regime_stats(r).pval_Radius_Mass = P_reg(1,2);
    else
        regime_stats(r).corr_Radius_Mass = NaN;
        regime_stats(r).pval_Radius_Mass = NaN;
    end
    
    % Consistency check within regime
    regime_stats(r).mean_Consistency = mean(Consistency_residual(idx));
    regime_stats(r).std_Consistency = std(Consistency_residual(idx));
    
    % Identify extreme cases in this regime
    regime_stats(r).extreme_radius = idx(Ratio_Radius(idx) > 2 | Ratio_Radius(idx) < 0.5);
    regime_stats(r).extreme_density = idx(Ratio_Sigma(idx) > 2 | Ratio_Sigma(idx) < 0.5);
end

%% --- 3.6.1 Mass-Dependent Radius Expansion ---
%% --- Prepare data ---
x = log_masses;
y = Ratio_Radius;

%% --- Define mass regimes (consistent with Figure 2) ---
mass_edges = [7.0, 9.0, 10.5, 12.0];
mass_labels = {'Low-Mass', 'Intermediate', 'High-Mass'};
regime_colors = [0.8 0.4 0.2;
    0.2 0.6 0.8;
    0.4 0.2 0.6];

%% --- Global regression analysis ---
p_global = polyfit(x, y, 1);
slope_global = p_global(1);
intercept_global = p_global(2);
[R_global, P_global] = corrcoef(x, y);
corr_global = R_global(1,2);
pval_global = P_global(1,2);

% Bootstrap confidence interval for slope
n_boot = 500;
boot_slopes = zeros(n_boot, 1);
for b = 1:n_boot
    idx = datasample(1:numGalaxies, numGalaxies, 'Replace', true);
    pb = polyfit(x(idx), y(idx), 1);
    boot_slopes(b) = pb(1);
end
ci_slope = prctile(boot_slopes, [2.5 97.5]);

%% --- Mass-stratified statistics ---
regime_stats = struct();
for r = 1:3
    idx = x >= mass_edges(r) & x < mass_edges(r+1);
    regime_stats(r).n = sum(idx);
    regime_stats(r).mean_ratio = mean(y(idx));
    regime_stats(r).std_ratio = std(y(idx));
    regime_stats(r).median_ratio = median(y(idx));
    regime_stats(r).min_ratio = min(y(idx));
    regime_stats(r).max_ratio = max(y(idx));
    
    if regime_stats(r).n >= 3
        p_reg = polyfit(x(idx), y(idx), 1);
        regime_stats(r).slope = p_reg(1);
        [R_reg, P_reg] = corrcoef(x(idx), y(idx));
        regime_stats(r).corr = R_reg(1,2);
        regime_stats(r).pval = P_reg(1,2);
    else
        regime_stats(r).slope = NaN;
        regime_stats(r).corr = NaN;
        regime_stats(r).pval = NaN;
    end
end

%% --- 3.6.2 Generate the Figure Mass-regime scatter ---
fig = figure('Color','w','Position',[100 100 850 650],'Name','Fig_3_6_2_MassDependence');
hold on;

h_scatter = gobjects(3,1);
for r = 1:3
    idx = x >= mass_edges(r) & x < mass_edges(r+1);
    if any(idx)
        h_scatter(r) = scatter(x(idx), y(idx), 60, ...
            regime_colors(r,:), ...
            'filled', ...
            'Marker','o', ...
            'MarkerFaceAlpha',0.7, ...
            'MarkerEdgeColor','k', ...
            'MarkerEdgeAlpha',0.4, ...
            'DisplayName', sprintf('%s (N=%d, mean=%.2f)', ...
                mass_labels{r}, regime_stats(r).n, regime_stats(r).mean_ratio));
    end
end

% Reference line: no change (keep this, but exclude from legend)
yline(1, 'k--', 'LineWidth', 1.5, 'HandleVisibility','off');

% regime-specific annotations
for r = 1:3
    idx = x >= mass_edges(r) & x < mass_edges(r+1);
    if any(idx)
        x_anno = mean(x(idx));
        y_anno = mean(y(idx));
        
        x_range = max(x(idx)) - min(x(idx));
        y_range = max(y(idx)) - min(y(idx));
        offset_x = -0.1 * x_range - 0.15;
        offset_y = 0.15 * y_range + 0.3;
        
        text(x_anno + offset_x, y_anno + offset_y, ...
            sprintf('mean = %.2f', regime_stats(r).mean_ratio), ...
            'Color',regime_colors(r,:), 'FontWeight','bold', 'FontSize',9, ...
            'BackgroundColor','w', 'EdgeColor',regime_colors(r,:), 'Margin',3, ...
            'HorizontalAlignment','center', 'VerticalAlignment','bottom');
    end
end

% Legend (only scatter points, no lines) - in box with white background
legend(h_scatter(h_scatter ~= 0), 'Location','northwest', 'Box','on', ...
    'FontSize',9, 'Interpreter','none', 'NumColumns',1, ...
    'Color','w', 'EdgeColor','k');

% Labels and title
xlabel('log_{10}(M_{obs} / M_{☉})', 'FontSize',22, 'FontWeight','bold');
ylabel('R_s^{True} / R_s^{Obs}', 'FontSize',22, 'FontWeight','bold');

% Axis formatting with extended y-axis limits
set(gca, 'FontSize',11, 'LineWidth',1.2, 'GridLineStyle','--', ...
    'XGrid','on', 'YGrid','on', 'YScale','log');
grid on; box on;

ylim([0.6 9.3]);
xlim([7.2 12.1]);

% Comprehensive y-ticks covering full range
yticks([0.6 0.7 0.8 0.9 1 1.2 1.5 2 2.5 3 4 5 6 7 8 9]);

% Save figure
exportgraphics(fig, fullfile(ScalingDir, 'Fig_3_6_2_RadiusRatio_vs_Mass.png'), 'Resolution', 400);
close(fig);


%% --- 3.6.3 Figure 3.6.3 ---
%% --- Prepare log-log data ---
log_R_ratio = log10(Ratio_Radius);
log_Sigma_ratio = log10(Ratio_Sigma);

%% --- Define mass regimes with distinct colors ---
mass_edges = [7.0, 9.0, 10.5, 12.0];
mass_labels = {'Low-Mass', 'Intermediate', 'High-Mass'};
regime_colors = [0.8 0.4 0.2;
    0.2 0.6 0.8;
    0.4 0.2 0.6];
regime_markers = {'o', 's', '^'};

%% --- Perform regression analysis ---
[p_obs, ~, ~, ~, stats] = regress(log_Sigma_ratio', [ones(numGalaxies,1), log_R_ratio']);
observed_slope = p_obs(2);
observed_intercept = p_obs(1);
r_squared = stats(1);
f_stat = stats(2);
p_value = stats(4);
slope_se = stats(3);

% Mass-stratified regressions
regime_slope = zeros(1,3);
regime_n = zeros(1,3);
regime_consistency = zeros(1,3);
for r = 1:3
    idx = log_masses >= mass_edges(r) & log_masses < mass_edges(r+1);
    regime_n(r) = sum(idx);
    if regime_n(r) >= 3
        [p_reg,~,~,~,st_reg] = regress(log_Sigma_ratio(idx)', [ones(sum(idx),1), log_R_ratio(idx)']);
        regime_slope(r) = p_reg(2);
        pred_reg = p_reg(2) * log_R_ratio(idx) + p_reg(1);
        res_reg = log_Sigma_ratio(idx) - pred_reg;
    else
        regime_slope(r) = NaN;
    end
end

%% --- 3.6.3 Generate the Figure with Mass-Regime Colors ---
fig = figure('Color','w','Position',[100 100 850 650],'Name','Fig_3_6_3_MassRegimeConsistency');
hold on;

h_scatter = gobjects(3,1);
for r = 1:3
    idx = log_masses >= mass_edges(r) & log_masses < mass_edges(r+1);
    if any(idx)
        h_scatter(r) = scatter(log_R_ratio(idx), log_Sigma_ratio(idx), 60, ...
            regime_colors(r,:), ...
            'filled', ...
            'Marker','o', ...
            'MarkerFaceAlpha',0.7, ...
            'MarkerEdgeColor','k', ...
            'MarkerEdgeAlpha',0.4, ...
            'DisplayName', sprintf('%s (N=%d, slope=%.2f)', ...
                mass_labels{r}, regime_n(r), regime_slope(r)));
    end
end

% Individual regime slope lines
for r = 1:3
    if regime_n(r) >= 3
        idx = log_masses >= mass_edges(r) & log_masses < mass_edges(r+1);
        x_reg = linspace(min(log_R_ratio(idx)), max(log_R_ratio(idx)), 50);
        intercept_reg = mean(log_Sigma_ratio(idx)) - regime_slope(r) * mean(log_R_ratio(idx));
        y_reg = regime_slope(r) * x_reg + intercept_reg;
        plot(x_reg, y_reg, '--', 'Color', regime_colors(r,:), 'LineWidth', 1.5, 'HandleVisibility','off');
    end
end

% Global observed regression line  
x_line_all = linspace(min(log_R_ratio)-0.3, max(log_R_ratio)+0.3, 100);
y_obs = p_obs(1) + p_obs(2) * x_line_all;
plot(x_line_all, y_obs, 'k-', 'LineWidth', 2, 'HandleVisibility','off');

x_right_line = max(x_line_all);
y_right_line = p_obs(1) + p_obs(2) * x_right_line;
text(x_right_line - 0.4, y_right_line, ...
    sprintf('global slope = %.2f', observed_slope), ...
    'Color','k', 'FontWeight','bold', 'FontSize',9, ...
    'BackgroundColor','w', 'EdgeColor','k', 'Margin',3, ...
    'HorizontalAlignment','left', 'VerticalAlignment','top');

% Legend (only scatter points, no lines) - in box with white background
legend(h_scatter(h_scatter ~= 0), 'Location','northwest', 'Box','on', ...
    'FontSize',12, 'Interpreter','none', 'NumColumns',1, ...
    'Color','w', 'EdgeColor','k');

% Label 
xlabel('log_{10}(R_s^{True} / R_s^{obs})', 'FontSize',22, 'FontWeight','bold');
ylabel('log_{10}(\Sigma_0^{True} / \Sigma_0^{obs})', 'FontSize',22, 'FontWeight','bold');

% Axis formatting with extended limits
set(gca, 'FontSize',11, 'LineWidth',1.2, 'GridLineStyle','--', ...
    'XGrid','on', 'YGrid','on');
grid on; box on;

% Extend axes to ensure all points, lines and tags are visible
xlim([min(log_R_ratio)-0.2 max(log_R_ratio)+0.2]);
ylim([min(log_Sigma_ratio)-0.2 max(log_Sigma_ratio)+0.4]);

% Comprehensive ticks covering full range
xticks(linspace(floor(min(log_R_ratio)), ceil(max(log_R_ratio)), 8));
yticks(linspace(floor(min(log_Sigma_ratio)), ceil(max(log_Sigma_ratio)), 8));

% Save figure
exportgraphics(fig, fullfile(ScalingDir, 'Fig_3_6_3_MassRegime_Consistency.png'), 'Resolution', 400);
close(fig);

%% --- 3.6.4 Figure 3.6.4 ---
%% --- Define mass regimes (consistent with Figures 1 & 2) ---
mass_edges = [7.0, 9.0, 10.5, 12.0];
mass_labels = {'Low-Mass', 'Intermediate', 'High-Mass'};
regime_colors = [0.8 0.4 0.2;
    0.2 0.6 0.8;
    0.4 0.2 0.6];

%% --- Prepare data for boxplots (concatenate with grouping) ---
radius_all = []; sigma_all = [];
group_radius = []; group_sigma = [];
for b = 1:3
    idx = log_masses >= mass_edges(b) & log_masses < mass_edges(b+1);
    if ~isempty(idx)
        radius_all = [radius_all; Ratio_Radius(idx)'];
        sigma_all = [sigma_all; Ratio_Sigma(idx)'];
        group_radius = [group_radius; b * ones(sum(idx), 1)];
        group_sigma = [group_sigma; b * ones(sum(idx), 1)];
    end
end

%% --- Calculate regime statistics ---
bin_stats = struct();
for b = 1:3
    idx = log_masses >= mass_edges(b) & log_masses < mass_edges(b+1);
    bin_stats(b).n = sum(idx);
    if bin_stats(b).n > 0
        % Radius statistics
        bin_stats(b).radius.mean = mean(Ratio_Radius(idx));
        bin_stats(b).radius.std = std(Ratio_Radius(idx));
        bin_stats(b).radius.median = median(Ratio_Radius(idx));
        bin_stats(b).radius.q25 = prctile(Ratio_Radius(idx), 25);
        bin_stats(b).radius.q75 = prctile(Ratio_Radius(idx), 75);
        bin_stats(b).radius.min = min(Ratio_Radius(idx));
        bin_stats(b).radius.max = max(Ratio_Radius(idx));
        
        % Density statistics
        bin_stats(b).sigma.mean = mean(Ratio_Sigma(idx));
        bin_stats(b).sigma.std = std(Ratio_Sigma(idx));
        bin_stats(b).sigma.median = median(Ratio_Sigma(idx));
        bin_stats(b).sigma.q25 = prctile(Ratio_Sigma(idx), 25);
        bin_stats(b).sigma.q75 = prctile(Ratio_Sigma(idx), 75);
        bin_stats(b).sigma.min = min(Ratio_Sigma(idx));
        bin_stats(b).sigma.max = max(Ratio_Sigma(idx));
        
        % Statistical tests: Is mean different from 1.0 (no correction)?
        [h_r, p_r] = ttest(Ratio_Radius(idx), 1.0);
        [h_s, p_s] = ttest(Ratio_Sigma(idx), 1.0);
        bin_stats(b).radius.pval_vs_1 = p_r;
        bin_stats(b).sigma.pval_vs_1 = p_s;
    end
end

%% --- 3.6.4 Generate the Figure ---
fig = figure('Color','white','Position',[100 100 1150 580],'Name','Fig_3_6_4_MassStratified');

% regime colors
regime_colors = [0.8 0.4 0.2;
    0.2 0.6 0.8;
    0.4 0.2 0.6];

% Spaced x-positions for three bins
x_positions = [1 2.5 4];
labels_plot = {'Low-Mass', 'Intermediate', 'High-Mass'};

%% Panel 1: Radius Ratio Boxplots
ax1 = subplot(1,2,1);
hold on; box on;
bp1 = boxplot(radius_all, group_radius, 'Positions', x_positions, ...
    'Colors','k', 'Symbol','+', 'Widths',0.5);
h_boxes1 = findobj(ax1,'Tag','Box');
h_whiskers1 = findobj(ax1,'Tag','Whisker');
h_medians1 = findobj(ax1,'Tag','Median');
h_outliers1 = findobj(ax1,'Tag','Outliers');
[~, sort_idx1] = sort(cellfun(@(x) mean(x), get(h_boxes1,'XData')));
h_boxes1 = h_boxes1(sort_idx1);

for j = 1:length(h_boxes1)
    xmid = mean(get(h_boxes1(j),'XData'));
    c = regime_colors(j,:);
    patch(get(h_boxes1(j),'XData'), get(h_boxes1(j),'YData'), ...
        c, 'FaceAlpha', 0.35, 'EdgeColor', c, 'LineWidth', 1.5);
    delete(h_boxes1(j));
    for w = 1:length(h_whiskers1)
        wx = mean(get(h_whiskers1(w),'XData'));
        if abs(wx - xmid) < 0.1
            set(h_whiskers1(w), 'Color', c, 'LineWidth', 1.5);
        end
    end
    for m = 1:length(h_medians1)
        mx = mean(get(h_medians1(m),'XData'));
        if abs(mx - xmid) < 0.1
            set(h_medians1(m), 'Color', 'k', 'LineWidth', 2);
        end
    end
end

for k = 1:length(h_outliers1)
    xval = mean(get(h_outliers1(k),'XData'));
    [~, bin_idx] = min(abs(x_positions - xval));
    set(h_outliers1(k), 'MarkerEdgeColor', regime_colors(bin_idx,:), ...
        'MarkerFaceColor', regime_colors(bin_idx,:));
end

yline(1, 'k--', 'LineWidth', 2, 'Alpha', 0.6);

% Statistical annotations
for b = 1:3
    if bin_stats(b).n > 0
        y_pos = max(radius_all) * 1;
        text(x_positions(b), y_pos, sprintf('N=%d\n%.2f±%.2f', ...
            bin_stats(b).n, bin_stats(b).radius.mean, bin_stats(b).radius.std), ...
            'HorizontalAlignment','center', 'FontSize',9, ...
            'BackgroundColor','w', 'EdgeColor', regime_colors(b,:), ...
            'LineWidth',1, 'Margin',3);
    end
end

% NORMAL FONT y-axis label
ylabel('R_s^{True} / R_s^{obs}', 'FontSize',14, 'FontWeight','bold');
set(ax1, 'YScale','log', 'FontSize',12, 'LineWidth',1.2, ...
    'XTick',x_positions, 'XTickLabel',labels_plot, ...
    'GridLineStyle','--');
grid on;
ylim(ax1, 'auto');
y_limits1 = ylim(ax1);
log_range1 = log10(y_limits1(2)) - log10(y_limits1(1));
new_min1 = 10^(log10(y_limits1(1)) - 0.15 * log_range1);
new_max1 = 10^(log10(y_limits1(2)) + 0.15 * log_range1);
ylim(ax1, [new_min1, new_max1]);
y_range1 = ylim(ax1);
y_log_min1 = floor(log10(y_range1(1)));
y_log_max1 = ceil(log10(y_range1(2)));
dense_ticks1 = [];
for exp1 = y_log_min1-1:y_log_max1+1
    dense_ticks1 = [dense_ticks1, (1:9) * 10^exp1];
end
yticks1 = dense_ticks1(dense_ticks1 >= y_range1(1)*0.95 & dense_ticks1 <= y_range1(2)*1.05);
set(ax1, 'YTick', yticks1);

%% Panel 2: Density Ratio Boxplots
ax2 = subplot(1,2,2);
hold on; box on;
bp2 = boxplot(sigma_all, group_sigma, 'Positions', x_positions, ...
    'Colors','k', 'Symbol','+', 'Widths',0.5);
h_boxes2 = findobj(ax2,'Tag','Box');
h_whiskers2 = findobj(ax2,'Tag','Whisker');
h_medians2 = findobj(ax2,'Tag','Median');
h_outliers2 = findobj(ax2,'Tag','Outliers');
[~, sort_idx2] = sort(cellfun(@(x) mean(x), get(h_boxes2,'XData')));
h_boxes2 = h_boxes2(sort_idx2);

for j = 1:length(h_boxes2)
    xmid = mean(get(h_boxes2(j),'XData'));
    c = regime_colors(j,:);
    patch(get(h_boxes2(j),'XData'), get(h_boxes2(j),'YData'), ...
        c, 'FaceAlpha', 0.35, 'EdgeColor', c, 'LineWidth', 1.5);
    delete(h_boxes2(j));
    for w = 1:length(h_whiskers2)
        wx = mean(get(h_whiskers2(w),'XData'));
        if abs(wx - xmid) < 0.1
            set(h_whiskers2(w), 'Color', c, 'LineWidth', 1.5);
        end
    end
    for m = 1:length(h_medians2)
        mx = mean(get(h_medians2(m),'XData'));
        if abs(mx - xmid) < 0.1
            set(h_medians2(m), 'Color', 'k', 'LineWidth', 2);
        end
    end
end

for k = 1:length(h_outliers2)
    xval = mean(get(h_outliers2(k),'XData'));
    [~, bin_idx] = min(abs(x_positions - xval));
    set(h_outliers2(k), 'MarkerEdgeColor', regime_colors(bin_idx,:), ...
        'MarkerFaceColor', regime_colors(bin_idx,:));
end

yline(1, 'k--', 'LineWidth', 2, 'Alpha', 0.6);

for b = 1:3
    if bin_stats(b).n > 0
        y_pos = max(sigma_all) * 1.13;
        text(x_positions(b), y_pos, sprintf('N=%d\n%.2f±%.2f', ...
            bin_stats(b).n, bin_stats(b).sigma.mean, bin_stats(b).sigma.std), ...
            'HorizontalAlignment','center', 'FontSize',9, ...
            'BackgroundColor','w', 'EdgeColor', regime_colors(b,:), ...
            'LineWidth',1, 'Margin',3);
    end
end

ylabel('\Sigma_0^{True} / \Sigma_0^{obs}', 'FontSize',14, 'FontWeight','bold');
set(ax2, 'YScale','log', 'FontSize',12, 'LineWidth',1.2, ...
    'XTick',x_positions, 'XTickLabel',labels_plot, ...
    'GridLineStyle','--');
grid on;
ylim(ax2, 'auto');
y_limits2 = ylim(ax2);
log_range2 = log10(y_limits2(2)) - log10(y_limits2(1));
new_min2 = 10^(log10(y_limits2(1)) - 0.15 * log_range2);
new_max2 = 10^(log10(y_limits2(2)) + 0.15 * log_range2);
ylim(ax2, [new_min2, new_max2]);
y_range2 = ylim(ax2);
y_log_min2 = floor(log10(y_range2(1)));
y_log_max2 = ceil(log10(y_range2(2)));
dense_ticks2 = [];
for exp2 = y_log_min2-1:y_log_max2+1
    dense_ticks2 = [dense_ticks2, (1:9) * 10^exp2];
end
yticks2 = dense_ticks2(dense_ticks2 >= y_range2(1)*0.95 & dense_ticks2 <= y_range2(2)*1.05);
set(ax2, 'YTick', yticks2);

%% Final adjustments
set(fig, 'Position', [100 100 1150 580]);

% Save figure
exportgraphics(fig, fullfile(ScalingDir, 'Fig_3_6_4_MassStratified_Summary.png'), 'Resolution', 400);
close(fig);

%% --- 3.7 Lensing-Dynamics Connection Analysis ---
%% 3.7.1 Data Preparation for Lensing-Dynamics Analysis
% Define mass regime edges 
mass_edges = [7.0, 9.0, 10.5, 12.0];
mass_labels = {'Low-Mass', 'Intermediate', 'High-Mass'};
regime_colors = [0.8 0.4 0.2;
    0.2 0.6 0.8;
    0.4 0.2 0.6];

% Extract required quantities from existing data structures
v_flat_obs = [galaxySummary.v_flat_obs_kms];
v_flat_true = [galaxySummary.v_flat_true_kms];
M_baryonic = [galaxySummary.M_total_obs_Msun];
M_true = [galaxySummary.M_total_true_Msun];

% Deflection angles (already computed in Section 3.3)
alpha_max_obs = [deflection_stats.max_alpha_obs_arcsec];
alpha_max_true = [deflection_stats.max_alpha_true_arcsec];

% Calculate baryonic accelerations at characteristic radii
g_bar_obs = zeros(1, numGalaxies);
g_bar_true = zeros(1, numGalaxies);
g_dyn_obs = zeros(1, numGalaxies);
g_dyn_true = zeros(1, numGalaxies);


for g = 1:numGalaxies
    % Use half-light radius or Einstein radius as characteristic scale
    r_char_o = max(Results(g).r_obs) / 2;
    if r_char_o == 0 || isnan(r_char_o)
        r_char_o = Results(g).r_obs(floor(end/2));
    end
    r_char_t = max(Results(g).r_true) / 2;
    if r_char_t == 0 || isnan(r_char_t)
        r_char_t = Results(g).r_obs(floor(end/2));
    end
    
    % Baryonic acceleration: g_bar = G * M_bar / r^2
    g_bar_obs(g) = G * M_baryonic(g) * Msun / (r_char_o)^2;
    g_bar_true(g) = G * M_true(g) * Msun / (r_char_t)^2;
    
    % Dynamical acceleration from rotation curve: g_dyn = v^2 / r
    v_flat_ms = v_flat_obs(g) * 1000;
    g_dyn_obs(g) = v_flat_ms^2 / r_char_o;
    v_true_ms = v_flat_true(g) * 1000;
    g_dyn_true(g) = v_true_ms^2 / r_char_t;
end

% Characteristic acceleration scale from RAR (McGaugh et al. 2016)
g_dagger = 1.20e-10;
g_bar_over_gdagger = g_bar_obs / g_dagger;

%% 3.7.2 Analysis 1: Deflection Angle vs. Rotation Velocity (BTFR Connection)

% Filter valid data for CORRECTED (true) velocities
valid_v_true = isfinite(v_flat_true) & v_flat_true > 0;
valid_alpha_true = isfinite(alpha_max_true) & alpha_max_true > 0;
valid_idx_true = valid_v_true & valid_alpha_true;

% Filter valid data for OBSERVED velocities
valid_v_obs = isfinite(v_flat_obs) & v_flat_obs > 0;
valid_alpha_obs = isfinite(alpha_max_obs) & alpha_max_obs > 0;
valid_idx_obs = valid_v_obs & valid_alpha_obs;

% Analysis for CORRECTED velocities
log_v_true = log10(v_flat_true(valid_idx_true));
log_alpha_true = log10(alpha_max_true(valid_idx_true));
p_alpha_v_true = polyfit(log_v_true, log_alpha_true, 1);
power_law_index_true = p_alpha_v_true(1);
alpha_v_intercept_true = p_alpha_v_true(2);
[R_alpha_v_true, P_alpha_v_true] = corrcoef(v_flat_true(valid_idx_true), alpha_max_true(valid_idx_true));

% Analysis for OBSERVED velocities
log_v_obs = log10(v_flat_obs(valid_idx_obs));
log_alpha_obs = log10(alpha_max_obs(valid_idx_obs));
p_alpha_v_obs = polyfit(log_v_obs, log_alpha_obs, 1);
power_law_index_obs = p_alpha_v_obs(1);
alpha_v_intercept_obs = p_alpha_v_obs(2);
[R_alpha_v_obs, P_alpha_v_obs] = corrcoef(v_flat_obs(valid_idx_obs), alpha_max_obs(valid_idx_obs));

%% Generate Figure 3.7.2: TRUE Deflection vs. OBSERVED Velocity
fig = figure('Color', 'white', 'Position', [100, 100, 800, 650], ...
    'Name', 'Fig_3_7_2_TrueDeflection_vs_ObservedVelocity');
hold on;

% Plot by mass regime - all circles
for r = 1:3
    idx = log_masses >= mass_edges(r) & log_masses < mass_edges(r+1) & ...
        valid_idx_true & valid_idx_obs;
    if any(idx)
        scatter(v_flat_obs(idx), alpha_max_true(idx), 90, regime_colors(r,:), ...
            'filled', 'Marker', 'o', 'MarkerFaceAlpha', 0.75, ...
            'MarkerEdgeColor', 'k', 'LineWidth', 1.2, ...
            'DisplayName', sprintf('%s (N=%d)', mass_labels{r}, sum(idx)));
    end
end

valid_idx_cross = valid_idx_true & valid_idx_obs;

%% Power law fit: alpha_true vs v_obs
% Linear fit in log-log space
log_v_obs = log10(v_flat_obs(valid_idx_cross));
log_alpha_true = log10(alpha_max_true(valid_idx_cross));

% Robust linear regression
p_cross = polyfit(log_v_obs, log_alpha_true, 1);
power_law_index_cross = p_cross(1);
alpha_v_intercept_cross = p_cross(2);

% Calculate fit statistics
[corr_matrix_cross, P_cross] = corrcoef(log_v_obs, log_alpha_true);
R_alpha_v_cross = corr_matrix_cross(1,2);

% Generate fit line
v_theory_cross = linspace(min(v_flat_obs(valid_idx_cross)), ...
    max(v_flat_obs(valid_idx_cross)), 100);
alpha_fit_cross = 10^alpha_v_intercept_cross * v_theory_cross.^power_law_index_cross;

% Black solid line for fit
plot(v_theory_cross, alpha_fit_cross, 'k-', 'LineWidth', 2.5, ...
    'DisplayName', sprintf('Fit: alpha \\propto Obs v^{%.2f}', power_law_index_cross));

%% Axis labels and title
xlabel('$v_{\mathrm{obs}}\;(\mathrm{km\,s^{-1}})$','Interpreter','latex', ...
    'FontSize',30,'FontWeight','bold');
ylabel('$\mathrm{True}\ \alpha_{\max}\;(\mathrm{arcsec})$','Interpreter','latex', ...
    'FontSize',30,'FontWeight','bold');

%% Legend - single column, white box
legend('Location', 'northwest', 'Box', 'on', 'FontSize', 10, 'NumColumns', 1, ...
    'Color', 'white', 'EdgeColor', 'black');
set(gca, 'FontSize', 11, 'LineWidth', 1.2, 'XScale', 'log', 'YScale', 'log', ...
    'GridLineStyle', '--', 'XGrid', 'on', 'YGrid', 'on');
grid on; box on;

ax = gca;
ax.XMinorTick = 'on';
ax.MinorGridLineStyle = ':';
ax.XMinorGrid = 'on';

%% Export figure
exportgraphics(fig, fullfile(ScalingDir, 'Fig_3_7_2_TrueDeflection_vs_ObservedVelocity.png'), ...
    'Resolution', 400);
close(fig);

%% --- 3.8 Rank-Order Preservation Test: ζ Maintains Surface Density Hierarchy ---

%% Data Extraction
% Extract central surface densities from pre-computed structs
Sigma0_obs_Msun_pc2 = [surface_stats.Sigma0_obs];       % Msun/pc^2, observed frame
Sigma0_true_Msun_pc2 = [surface_stats.Sigma0_true];      % Msun/pc^2, true/corrected frame
log_M_obs = log10([mass_scaling.M_total_obs]);           % For mass-regime binning

% Log-transform for linear visualization and correlation analysis
log_Sigma0_obs = log10(Sigma0_obs_Msun_pc2);
log_Sigma0_true = log10(Sigma0_true_Msun_pc2);

%% Mass Regime Definition 
mass_edges = [7.0, 9.0, 10.5, 12.0];
mass_labels = {'Low-Mass', 'Intermediate', 'High-Mass'};
regime_colors = [0.8 0.4 0.2;
                 0.2 0.6 0.8;
                 0.4 0.2 0.6];

%% Per-Regime Statistical Analysis
rank_stats = struct();
for r = 1:3
    idx = log_M_obs >= mass_edges(r) & log_M_obs < mass_edges(r+1);
    if sum(idx) < 3
        rank_stats(r).n = sum(idx);
        rank_stats(r).slope = NaN;
        rank_stats(r).intercept = NaN;
        rank_stats(r).R2 = NaN;
        rank_stats(r).sigma_1to1 = NaN;
        rank_stats(r).pearson_r = NaN;
        continue;
    end
    
    % Extract regime data
    x_reg = log_Sigma0_obs(idx);
    y_reg = log_Sigma0_true(idx);
    
    % Linear regression: log(Σ₀^true) = a + b·log(Σ₀^obs)
    p_reg = polyfit(x_reg, y_reg, 1);
    
    % Goodness-of-fit
    y_fit = polyval(p_reg, x_reg);
    SS_res = sum((y_reg - y_fit).^2);
    SS_tot = sum((y_reg - mean(y_reg)).^2);
    R2_reg = 1 - SS_res / SS_tot;
    
    % Scatter around 1:1 line (dex)
    residuals_1to1 = y_reg - x_reg;
    sigma_1to1_reg = std(residuals_1to1);
    
    % Pearson correlation
    [R_reg, ~] = corrcoef(x_reg, y_reg);
    
    % Store metrics
    rank_stats(r).n = sum(idx);
    rank_stats(r).slope = p_reg(1);
    rank_stats(r).intercept = p_reg(2);
    rank_stats(r).R2 = R2_reg;
    rank_stats(r).sigma_1to1 = sigma_1to1_reg;
    rank_stats(r).pearson_r = R_reg(1,2);
end

%% Global Metrics (for reference)
[p_global, ~, ~, ~, stats_global] = regress(log_Sigma0_true', [ones(numGalaxies,1), log_Sigma0_obs']);
global_slope = p_global(2);
global_R2 = stats_global(1);
global_sigma_1to1 = std(log_Sigma0_true - log_Sigma0_obs);
[global_pearson, ~] = corrcoef(log_Sigma0_obs, log_Sigma0_true);

%% Figure 3.8.1: Observed vs True Central Surface Density  
fig = figure('Color','white','Position',[150 100 900 700],'Name','Fig_3_8_1_RankOrder');
hold on; box on;

% Plot each mass regime with distinct color/marker
h_scatter = gobjects(3,1);
for r = 1:3
    idx = log_M_obs >= mass_edges(r) & log_M_obs < mass_edges(r+1);
    if any(idx)
        if ~isnan(rank_stats(r).slope)
            legendStr = sprintf('%s: slope=%.3f, r=%.3f (N=%d)', mass_labels{r}, rank_stats(r).slope, rank_stats(r).pearson_r, sum(idx));
        else
            legendStr = sprintf('%s (N=%d)', mass_labels{r}, sum(idx));
        end
        h_scatter(r) = scatter(log_Sigma0_obs(idx), log_Sigma0_true(idx), 65, ...
            regime_colors(r,:), ...
            'filled', ...
            'Marker','o', ...
            'MarkerFaceAlpha',0.75, ...
            'MarkerEdgeColor','k', ...
            'MarkerEdgeAlpha',0.5, ...
            'LineWidth',0.8, ...
            'DisplayName', legendStr);
    end
end

% 1:1 reference line (perfect rank preservation)
x_lim = [min(log_Sigma0_obs)-0.15, max(log_Sigma0_obs)+0.15];
plot(x_lim, x_lim, 'k--', 'LineWidth', 2, 'DisplayName', '1:1 line');

% Regime-specific fit lines (dashed, matching regime color)
for r = 1:3
    idx = log_M_obs >= mass_edges(r) & log_M_obs < mass_edges(r+1);
    if sum(idx) >= 3 && ~isnan(rank_stats(r).slope)
        x_fit = linspace(min(log_Sigma0_obs(idx)), max(log_Sigma0_obs(idx)), 50);
        y_fit = rank_stats(r).intercept + rank_stats(r).slope * x_fit;
        plot(x_fit, y_fit, '--', 'Color', regime_colors(r,:), 'LineWidth', 1.5, ...
            'HandleVisibility','off');
    end
end

% Legend: regime labels with slope and correlation, clean box style
legend(h_scatter(h_scatter~=0), 'Location','southeast', 'Box','on', ...
    'FontSize',9, 'Interpreter','none', 'NumColumns',1, ...
    'Color','w', 'EdgeColor','k', 'LineWidth',0.8);

 
% X-axis label (horizontal, normal)
xlabel('$\mathbf{\log_{10}(\Sigma_{0}^{\rm obs} / M_\odot \, {\rm pc}^{-2})}$', ...
       'Interpreter','latex', 'FontSize',30, 'FontWeight','bold');

% Y-axis label (vertical, default)
ylabel('$\mathbf{\log_{10}(\Sigma_{0}^{\rm true} / M_\odot \, {\rm pc}^{-2})}$', ...
       'Interpreter','latex', 'FontSize',30, 'FontWeight','bold');


% Clean axis formatting
set(gca, 'FontSize',11, 'LineWidth',1.2, 'GridLineStyle','--', ...
    'XGrid','on', 'YGrid','on', 'Box','on');
grid on;
axis equal tight;

% Save figure
exportgraphics(fig, fullfile(ScalingDir, 'Fig_3_8_1_RankOrderPreservation.png'), 'Resolution', 400);
close(fig);
 



%% --- 3.8.2.  Mass Scaling Verification  ---

% Extract necessary quantities
R_ratio = [surface_stats.R_scale_true_pc] ./ [surface_stats.R_scale_obs_pc];
M_ratio = [mass_scaling.mass_ratio]; 
log_masses = log10([galaxySummary.M_total_obs_Msun]);

% Calculate deviation from the 1:1 line
log_dev_from_1to1 = log10(M_ratio) - log10(R_ratio);
mean_log_dev = mean(log_dev_from_1to1);
std_log_dev  = std(log_dev_from_1to1);
median_log_dev = median(log_dev_from_1to1);

% Relative deviation in linear space for intuitive percentage reporting
rel_dev = (M_ratio ./ R_ratio) - 1;
mean_rel_dev = mean(rel_dev);

% Power-law fit for displaying fit line
log_R = log10(R_ratio);
log_M = log10(M_ratio);
p_fit = polyfit(log_R, log_M, 1);
measured_exponent = p_fit(1);
fit_intercept = p_fit(2);


fig = figure('Color', 'white', 'Position', [100, 100, 800, 600], 'Name', 'Fig_4_2_ScalingVerification');
ax = axes(fig);
hold(ax, 'on'); 
box(ax, 'on');

% Mass-based coloring  
mass_colors = parula(numGalaxies);
[~, sort_idx] = sort(log_masses, 'ascend');
[~, color_idx] = sort(sort_idx);

% Plot data points  
for g = 1:numGalaxies
    scatter(ax, R_ratio(g), M_ratio(g), 80, mass_colors(color_idx(g), :), ...
        'o', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 1);
end

% Calculate fit line  
R_range = linspace(min(R_ratio)*0.9, max(R_ratio)*1.1, 200);
M_fit = 10.^(p_fit(1) * log10(R_range) + p_fit(2));

% Plot 1:1 reference line  
plot(ax, R_range, R_range, 'k--', 'LineWidth', 2, 'DisplayName', '1:1 Line');

% Plot power-law fit line  
plot(ax, R_range, M_fit, 'r-', 'LineWidth', 2, 'DisplayName', sprintf('Fit (slope=%.2f)', measured_exponent));

 
% Format axes 
set(ax, 'XScale', 'log', 'YScale', 'log', 'FontSize', 12, 'LineWidth', 1.2);
xlabel('R_s^{True} / R_s^{obs}', 'FontSize', 16, 'FontWeight', 'bold', 'Interpreter', 'tex');
ylabel('M_{true} / M_{obs}', 'FontSize', 16, 'FontWeight', 'bold', 'Interpreter', 'tex');
grid(ax, 'on');

% Set axis limits
ax.XLim = [0 10];
ax.YLim = [0 10];

% Set tick values (example: every 1 unit)
ax.XTick = 0:1:12;
ax.YTick = 0:1:12;

% Add colorbar 
cbar = colorbar('Location', 'eastoutside');
cbar.Label.String = 'log_{10}(M_{obs}/M_{☉})';
cbar.Label.FontSize = 12;
cbar.Label.FontWeight = 'bold';
cbar.Ticks = linspace(0, 1, 5);
cbar.TickLabels = arrayfun(@(x) sprintf('%.1f', x), ...
    linspace(min(log_masses), max(log_masses), 5), 'UniformOutput', false);

% Statistics annotation  
stats_text = {
    sprintf('N = %d galaxies', numGalaxies)
    sprintf('Exponent = %.3f', measured_exponent)
    sprintf('Mean log-dev = %.3f dex', mean_log_dev)
    sprintf('Std dev = %.3f dex', std_log_dev)
};
annotation('textbox', [0.13, 0.69, 0.22, 0.22], 'String', stats_text, ...
    'FontSize', 12, 'BackgroundColor', [1 1 1 0.9], 'EdgeColor', [0.3 0.3 0.3], ...
    'LineWidth', 1, 'FitBoxToText', 'on', ...
    'HorizontalAlignment', 'left', 'VerticalAlignment', 'top');


% Export, save, and close
drawnow;
exportgraphics(fig, fullfile(ScalingDir, 'Fig_3_8_2_Scaling_Verification.png'), 'Resolution', 300);
close(fig);

%% --- 3.8.3 Figure: Total Mass Ratio vs Mezzi Scale Factor (ζ) ---
Mean_Mezzi = Mezzi_means(:);
mass_ratio = (M_true_all ./ M_obs_all);
mass_ratio = mass_ratio(:);
lg_all = log10(galaxy_masses(:));

N_data = min([length(Mean_Mezzi), length(mass_ratio), length(lg_all)]);
Mean_Mezzi = Mean_Mezzi(1:N_data);
mass_ratio = mass_ratio(1:N_data);
lg_all = lg_all(1:N_data);

valid_mask = isfinite(Mean_Mezzi) & isfinite(mass_ratio) & (Mean_Mezzi > 0) & (mass_ratio > 0);
z_plot = Mean_Mezzi(valid_mask);
m_plot = mass_ratio(valid_mask);
lg_plot = lg_all(valid_mask);
n_valid = length(z_plot);

[~, sort_idx] = sort(lg_plot, 'ascend');
[~, color_idx] = sort(sort_idx);
valid_colors = parula(n_valid);
point_colors = valid_colors(color_idx, :);

log_z = log10(z_plot);
log_m = log10(m_plot);
p_fit = polyfit(log_z, log_m, 1);
fitted_slope = p_fit(1);
fitted_intercept = p_fit(2);
[R_val, ~] = corrcoef(log_z, log_m);
corr_coef = R_val(1,2);

z_fit = linspace(min(z_plot)*0.9, max(z_plot)*1.1, 200);
m_fit = 10.^(fitted_intercept + fitted_slope * log10(z_fit));

fig = figure('Color', 'white', 'Position', [100, 100, 800, 600], ...
    'Name', 'Fig_3_8_3_Mass_vs_Zeta');
hold on; box on;

h_scatter = scatter(z_plot, m_plot, 80, point_colors, ...
    'o', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 1, ...
    'DisplayName', 'Galaxies');
h_fit = plot(z_fit, m_fit, 'r-', 'LineWidth', 2, ...
    'DisplayName', sprintf('Fit: slope = %.2f (r=%.3f)', fitted_slope, corr_coef));

cbar = colorbar('Location', 'eastoutside');
cbar.Label.String = 'log_{10}(M_{obs}/M_{\odot})';
cbar.Label.FontSize = 14;
cbar.Label.FontWeight = 'bold';
cbar.Ticks = linspace(0, 1, 5);
cbar.TickLabels = arrayfun(@(x) sprintf('%.1f', x), ...
    linspace(min(lg_plot), max(lg_plot), 5), 'UniformOutput', false);

xlabel('Mezzi Scale Factor \zeta', 'FontSize', 22, 'FontWeight', 'bold');
ylabel('M_{true} / M_{obs}', 'FontSize', 22, 'FontWeight', 'bold');
set(gca, 'XScale', 'log', 'YScale', 'log', 'FontSize', 12, 'LineWidth', 1.2, ...
    'GridLineStyle', '--', 'XGrid', 'on', 'YGrid', 'on');
grid on;

legend(h_fit, 'Location', 'southwest', 'Box', 'on', ...
    'FontSize', 10, 'Interpreter', 'none');

stats_text = {
    sprintf('N = %d galaxies', n_valid)
    sprintf('Slope = %.3f (r = %.3f)', fitted_slope, corr_coef)
};
 
drawnow;
exportgraphics(fig, fullfile(ScalingDir, 'Fig_3_8_3_Mass_vs_Zeta.png'), 'Resolution', 400);
close(fig);

%% --- 3.8.4 Figure: Max Deflection Ratio vs Mezzi Scale Factor (ζ) ---
alpha_ratio = (max_alpha_true ./ max_alpha_obs);
alpha_ratio = alpha_ratio(:);

N_data = min([length(Mean_Mezzi), length(alpha_ratio), length(lg_all)]);
Mean_Mezzi = Mean_Mezzi(1:N_data);
alpha_ratio = alpha_ratio(1:N_data);
lg_all = lg_all(1:N_data);

valid_mask = isfinite(Mean_Mezzi) & isfinite(alpha_ratio) & (Mean_Mezzi > 0) & (alpha_ratio > 0);
z_plot = Mean_Mezzi(valid_mask);
a_plot = alpha_ratio(valid_mask);
lg_plot = lg_all(valid_mask);
n_valid = length(z_plot);

[~, sort_idx] = sort(lg_plot, 'ascend');
[~, color_idx] = sort(sort_idx);
valid_colors = parula(n_valid);
point_colors = valid_colors(color_idx, :);

log_z = log10(z_plot);
log_a = log10(a_plot);
p_fit = polyfit(log_z, log_a, 1);
fitted_slope = p_fit(1);
fitted_intercept = p_fit(2);
[R_val, ~] = corrcoef(log_z, log_a);
corr_coef = R_val(1,2);

z_fit = linspace(min(z_plot)*0.9, max(z_plot)*1.1, 200);
a_fit = 10.^(fitted_intercept + fitted_slope * log10(z_fit));

fig = figure('Color', 'white', 'Position', [100, 100, 800, 600], ...
    'Name', 'Fig_3_8_4_Deflection_vs_Zeta');
hold on; box on;

h_scatter = scatter(z_plot, a_plot, 80, point_colors, ...
    'o', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 1, ...
    'DisplayName', 'Galaxies');
h_fit = plot(z_fit, a_fit, 'r-', 'LineWidth', 2, ...
    'DisplayName', sprintf('Fit: slope = %.2f (r=%.3f)', fitted_slope, corr_coef));

cbar = colorbar('Location', 'eastoutside');
cbar.Label.String = 'log_{10}(M_{obs}/M_{\odot})';
cbar.Label.FontSize = 14;
cbar.Label.FontWeight = 'bold';
cbar.Ticks = linspace(0, 1, 5);
cbar.TickLabels = arrayfun(@(x) sprintf('%.1f', x), ...
    linspace(min(lg_plot), max(lg_plot), 5), 'UniformOutput', false);

xlabel('Mezzi Scale Factor \zeta', 'FontSize', 22, 'FontWeight', 'bold');
ylabel('\alpha_{true} / \alpha_{obs}', 'FontSize', 22, 'FontWeight', 'bold');
set(gca, 'XScale', 'log', 'YScale', 'log', 'FontSize', 12, 'LineWidth', 1.2, ...
    'GridLineStyle', '--', 'XGrid', 'on', 'YGrid', 'on');
grid on;

legend( h_fit, 'Location', 'southwest', 'Box', 'on', ...
    'FontSize', 10, 'Interpreter', 'none');

stats_text = {
    sprintf('N = %d galaxies', n_valid)
    sprintf('Slope = %.3f (r = %.3f)', fitted_slope, corr_coef)
};
 
drawnow;
exportgraphics(fig, fullfile(ScalingDir, 'Fig_3_8_4_Deflection_vs_Zeta.png'), 'Resolution', 400);
close(fig);

clear exp
 
%% ===== SECTION 4: RAR VALIDATION ANALYSIS =====
%
% THREE TEST CASES (increasing complexity):
%   C1: Newtonian Baseline    → Validates true-frame self-consistency
%   C2: Empirical Mapping     → Tests if ζ(r) reproduces observed RAR shape
%   C3: Baryonic-Axis Shift   → Tests geometric effect on x-axis alone
%

%% --- 4.1 Compute Accelerations for Three Test Cases ---
% Physical constants
a0_RAR = 1.20e-10;  % Characteristic acceleration scale (m/s²)
g_dagger = a0_RAR;

% Reference curves for RAR analysis
RAR_empirical = @(g) g ./ (1 - exp(-sqrt(g / a0_RAR)));  % McGaugh Empirical fit
RAR_newtonian = @(g_bar) g_bar;                          % True physics: g_dyn = g_bar
RAR_mond = @(g_bar) sqrt(g_bar * g_dagger);              % Deep MOND: g_dyn = √(g_bar·a₀)

% Structure to store results for all galaxies
RAR_Analysis = struct();

for g = 1:numGalaxies
    fprintf('[RAR] Galaxy %d/%d: %s\n', g, numGalaxies, galaxy_names{g});
    
    % Extract data from Results structure
    r_obs_full = Results(g).r_obs(:);
    r_true_full = Results(g).r_true(:);
    M_obs_full = Results(g).M_enc_obs(:);
    M_true_full = Results(g).M_enc_true(:);
    v_obs_full = Results(g).v_obs(:);
    v_Newton_true_full = Results(g).v_Newton_true(:);
    
    N = min([length(r_obs_full), length(r_true_full), length(M_obs_full), ...
        length(M_true_full), length(v_obs_full), length(v_Newton_true_full)]);
    r_obs = r_obs_full(1:N);
    r_true = r_true_full(1:N);
    M_obs = M_obs_full(1:N);
    M_true = M_true_full(1:N);
    v_obs = v_obs_full(1:N);
    v_Newton_true = v_Newton_true_full(1:N);
    
    valid = (r_obs > 1e3) & (M_obs > 0) & (r_true > 1e3) & (M_true > 0) & ...
        (v_obs > 0) & (v_Newton_true > 0) & ...
        isfinite(r_obs) & isfinite(M_obs) & isfinite(r_true) & isfinite(M_true) & ...
        isfinite(v_obs) & isfinite(v_Newton_true);
    
    if sum(valid) < 5
        RAR_Analysis(g).valid = false;
        fprintf('[SKIP] Galaxy %d: insufficient valid points\n', g);
        continue;
    end
    
    % Apply validity mask
    r_o = r_obs(valid);
    r_t = r_true(valid);
    M_o = M_obs(valid);
    M_t = M_true(valid);
    v_o = v_obs(valid);
    v_t_newton = v_Newton_true(valid);
    Mezzi_local = Results(g).Mezzi_th(valid);
    
    r_o = r_o(:); r_t = r_t(:); M_o = M_o(:); M_t = M_t(:);
    v_o = v_o(:); v_t_newton = v_t_newton(:); Mezzi_local = Mezzi_local(:);
    
    %% Compute accelerations for three test cases
    % CASE 1: Newtonian Consistency Check (TRUE-FRAME BASELINE)
    g_bar_C1 = G * M_t ./ r_t.^2;
    g_dyn_C1 = v_t_newton.^2 ./ r_t;
    
    % CASE 2: Empirical Mapping Test (MAIN HYPOTHESIS TEST)
    g_bar_C2 = G * M_t ./ r_t.^2;  % True baryonic field at r_t
    g_dyn_C2 = v_o.^2 ./ r_o;  % Observed dynamics at r_o
    
    % CASE 3: Baryonic-Axis Shift Test (ISOLATED GEOMETRIC EFFECT)
    g_bar_C3 = G * M_t ./ r_o.^2;  % True mass evaluated at compressed radius r_o
    g_dyn_C3 = v_o.^2 ./ r_o;  % Same dynamics as Case 2
    
    % OBSERVED: Standard measurement ( for comparison)
    g_bar_obs = G * M_o ./ r_o.^2;
    g_dyn_obs = v_o.^2 ./ r_o;
    
    g_bar_C1 = g_bar_C1(:); g_dyn_C1 = g_dyn_C1(:);
    g_bar_C2 = g_bar_C2(:); g_dyn_C2 = g_dyn_C2(:);
    g_bar_C3 = g_bar_C3(:); g_dyn_C3 = g_dyn_C3(:);
    g_bar_obs = g_bar_obs(:); g_dyn_obs = g_dyn_obs(:);
    
    %% Compute scatter metrics (log-space, standard for RAR analysis)
    calc_scatter = @(g_dyn, g_pred) std(log10(g_dyn) - log10(g_pred));
    scatter_C1_vs_newt = calc_scatter(g_dyn_C1, g_bar_C1);
    scatter_C2_vs_emp  = calc_scatter(g_dyn_C2, RAR_empirical(g_bar_C2));
    scatter_C3_vs_mond = calc_scatter(g_dyn_C3, RAR_mond(g_bar_C3));
    scatter_obs_vs_emp = calc_scatter(g_dyn_obs, RAR_empirical(g_bar_obs));
    
    %% Store results in structure
    RAR_Analysis(g).valid = true;
    RAR_Analysis(g).galaxy = galaxy_names{g};
    RAR_Analysis(g).N = sum(valid);
    RAR_Analysis(g).g_bar_C1 = g_bar_C1; RAR_Analysis(g).g_dyn_C1 = g_dyn_C1;
    RAR_Analysis(g).g_bar_C2 = g_bar_C2; RAR_Analysis(g).g_dyn_C2 = g_dyn_C2;
    RAR_Analysis(g).g_bar_C3 = g_bar_C3; RAR_Analysis(g).g_dyn_C3 = g_dyn_C3;
    RAR_Analysis(g).g_bar_obs = g_bar_obs; RAR_Analysis(g).g_dyn_obs = g_dyn_obs;
    RAR_Analysis(g).Mezzi = Mezzi_local;
    RAR_Analysis(g).r_norm = r_o / max(r_o);
    
    % Store key scatter metrics only
    RAR_Analysis(g).scatter_C1_newt = scatter_C1_vs_newt;
    RAR_Analysis(g).scatter_C2_emp  = scatter_C2_vs_emp;
    RAR_Analysis(g).scatter_C3_mond = scatter_C3_vs_mond;
    RAR_Analysis(g).scatter_obs_emp = scatter_obs_vs_emp;
end

%% --- 4.2 Global Analysis - Aggregate Results ---
% Collect all valid data from all galaxies
all_gbar_C1 = []; all_gdyn_C1 = [];
all_gbar_C2 = []; all_gdyn_C2 = [];
all_gbar_C3 = []; all_gdyn_C3 = [];
all_gbar_obs = []; all_gdyn_obs = [];

for g = 1:numGalaxies
    if ~isfield(RAR_Analysis(g), 'valid') || ~RAR_Analysis(g).valid
        continue;
    end
    all_gbar_C1 = [all_gbar_C1; RAR_Analysis(g).g_bar_C1];
    all_gdyn_C1 = [all_gdyn_C1; RAR_Analysis(g).g_dyn_C1];
    all_gbar_C2 = [all_gbar_C2; RAR_Analysis(g).g_bar_C2];
    all_gdyn_C2 = [all_gdyn_C2; RAR_Analysis(g).g_dyn_C2];
    all_gbar_C3 = [all_gbar_C3; RAR_Analysis(g).g_bar_C3];
    all_gdyn_C3 = [all_gdyn_C3; RAR_Analysis(g).g_dyn_C3];
    all_gbar_obs = [all_gbar_obs; RAR_Analysis(g).g_bar_obs];
    all_gdyn_obs = [all_gdyn_obs; RAR_Analysis(g).g_dyn_obs];
end

% Clean data (restrict to physically meaningful acceleration range)
valid_all = (all_gbar_C1 > 1e-15) & (all_gbar_C1 < 1e-5) & ...
    (all_gdyn_C1 > 1e-15) & (all_gdyn_C1 < 1e-5);
N_total = sum(valid_all);

%% Global scatter metrics (using relevant reference curve for each case)
calc_scatter_global = @(g_bar, g_dyn, g_pred) std(log10(g_dyn(valid_all)) - log10(g_pred(g_bar(valid_all))));
scatter_global_C1_newt = calc_scatter_global(all_gbar_C1, all_gdyn_C1, RAR_newtonian);
scatter_global_C2_emp  = calc_scatter_global(all_gbar_C2, all_gdyn_C2, RAR_empirical);
scatter_global_C3_mond = calc_scatter_global(all_gbar_C3, all_gdyn_C3, RAR_mond);
scatter_global_obs_emp = calc_scatter_global(all_gbar_obs, all_gdyn_obs, RAR_empirical);

fprintf('[RESULT] Global scatter metrics:\n');
fprintf('  C1 vs Newtonian: %.3f dex\n', scatter_global_C1_newt);
fprintf('  C2 vs Empirical RAR: %.3f dex\n', scatter_global_C2_emp);
fprintf('  C3 vs deep-MOND: %.3f dex\n', scatter_global_C3_mond);
fprintf('  Observed vs Empirical RAR: %.3f dex\n', scatter_global_obs_emp);

%% --- 4.3 Generate Figures (Professional, Informative) ---

rarDir = fullfile(baseDir, 'RAR_Figures');
if ~exist(rarDir, 'dir')
    mkdir(rarDir);
end

% Reference curve for plotting
g_ref = logspace(-2, 2, 100) * 1e-10;

%% Figure 4.3.1: Case 1 - Newtonian Consistency
fig = figure('Color','white','Position',[100 100 1000 600],'Name','Fig_4_3_1_Case1_Newtonian');
subplot('Position',[0.12 0.15 0.82 0.78]); hold on;

% Data
scatter(all_gbar_obs(valid_all)*1e10, all_gdyn_obs(valid_all)*1e10, 15, [0.6 0.6 0.6], ...
    'filled','MarkerFaceAlpha',0.25,'DisplayName','Observed frame');
scatter(all_gbar_C1(valid_all)*1e10, all_gdyn_C1(valid_all)*1e10, 20, [0.2 0.4 0.8], ...
    'filled','MarkerFaceAlpha',0.7,'DisplayName',' True frame');

% Reference curves
plot(g_ref*1e10, g_ref*1e10, 'k--','LineWidth',3,'DisplayName','Newtonian: g_{dyn}=g_{bar}');
plot(g_ref*1e10, RAR_empirical(g_ref)*1e10, 'k-','LineWidth',3,'DisplayName','Empirical RAR');
plot(g_ref*1e10, sqrt(g_ref*1e10)*sqrt(1.2), 'k-.','LineWidth',2.5,'DisplayName','deep-MOND: g_{dyn}=\surd(g_{bar}a_0)');

set(gca,'XScale','log','YScale','log','FontSize',10,'LineWidth',1);
xlabel('g_{bar} [10^{-10} m/s^2]','FontSize',16,'FontWeight','bold');
ylabel('g_{dyn} [10^{-10} m/s^2]','FontSize',16,'FontWeight','bold');
legend('Location','southeast','FontSize',10);
grid on; box on; axis([0.01 100 0.01 100]);

% Info box 
box_pos = [0.15 0.64 0.28 0.28];
box_text = {...
    'True Frame Consistency',...
    '',...
    'DEFINITION:',...
    '  g_{bar} = \color{blue}GM_{True}\color{black} / \color{blue}r_{True}^2\color{black} ',...
    '  g_{dyn} = \color{blue}v_{Newton,True}^2\color{black} / \color{blue}r_{True}\color{black} ',...
    '',...
    sprintf('Scatter vs Newtonian: %.3f dex',scatter_global_C1_newt),...
};
annotation('textbox',box_pos,'String',box_text,...
    'FontSize',11,'FontName','Consolas','BackgroundColor','white',...
    'EdgeColor',[0.2 0.4 0.8],'LineWidth',1.5,'FitBoxToText','on');

exportgraphics(fig, fullfile(rarDir, 'Fig_4_3_1_Case1_Newtonian.png'), 'Resolution', 300);
close(fig);

%% Figure 4.3.2: Case 2 - Empirical Mapping Test

fig = figure('Color','white','Position',[100 100 1000 600],'Name','Fig_4_3_2_Case2_EmpiricalMapping');
subplot('Position',[0.12 0.15 0.82 0.78]); hold on;

scatter(all_gbar_obs(valid_all)*1e10, all_gdyn_obs(valid_all)*1e10, 15, [0.6 0.6 0.6], ...
    'filled','MarkerFaceAlpha',0.25,'DisplayName','Observed frame');
scatter(all_gbar_C2(valid_all)*1e10, all_gdyn_C2(valid_all)*1e10, 20, [0.9 0.5 0.1], ...
    'filled','MarkerFaceAlpha',0.7,'DisplayName','Empirical Mapping');

plot(g_ref*1e10, g_ref*1e10, 'k--','LineWidth',3,'DisplayName','Newtonian: g_{dyn}=g_{bar}');
plot(g_ref*1e10, RAR_empirical(g_ref)*1e10, 'k-','LineWidth',3,'DisplayName','Empirical RAR');
plot(g_ref*1e10, sqrt(g_ref*1e10)*sqrt(1.2), 'k-.','LineWidth',2.5,'DisplayName','deep-MOND: g_{dyn}=\surd(g_{bar}a_0)');

set(gca,'XScale','log','YScale','log','FontSize',10,'LineWidth',1);
xlabel('g_{bar} [10^{-10} m/s^2]','FontSize',16,'FontWeight','bold');
ylabel('g_{dyn} [10^{-10} m/s^2]','FontSize',16,'FontWeight','bold');
legend('Location','southeast','FontSize',10);
grid on; box on; axis([0.01 100 0.01 100]);

box_pos = [0.15 0.64 0.28 0.28];
box_text = {...
    'Empirical Mapping',...
    '',...
    'DEFINITION:',...
    '  g_{bar} = \color{blue}GM_{True}\color{black} / \color{blue}r_{True}^2\color{black} ',...
    '  g_{dyn} = \color{red}v_{obs}^2\color{black} / \color{red}r_{obs}\color{black} ',...
    '',...
    sprintf('Scatter vs RAR: %.3f dex',scatter_global_C2_emp),...
    sprintf('Observed scatter: %.3f dex',scatter_global_obs_emp),...
    sprintf('ratio: %.2f',scatter_global_C2_emp/scatter_global_obs_emp),...
};
annotation('textbox',box_pos,'String',box_text,...
    'FontSize',11,'FontName','Consolas','BackgroundColor','white',...
    'EdgeColor',[0.9 0.5 0.1],'LineWidth',1.5,'FitBoxToText','on');

exportgraphics(fig, fullfile(rarDir, 'Fig_4_3_2_Case2_EmpiricalMapping.png'), 'Resolution', 300);
close(fig);

%% Figure 4.3.3: Case 3 - Baryonic-Scale Shift
fig = figure('Color','white','Position',[100 100 1000 600],'Name','Fig_4_3_3_Case3_BaryonicShift');
subplot('Position',[0.12 0.15 0.82 0.78]); hold on;

scatter(all_gbar_obs(valid_all)*1e10, all_gdyn_obs(valid_all)*1e10, 15, [0.6 0.6 0.6], ...
    'filled','MarkerFaceAlpha',0.25,'DisplayName','Observed frame');
scatter(all_gbar_C3(valid_all)*1e10, all_gdyn_C3(valid_all)*1e10, 20, [0.6 0.2 0.6], ...
    'filled','MarkerFaceAlpha',0.7,'DisplayName',' Baryonic Scale shift ');

plot(g_ref*1e10, g_ref*1e10, 'k--','LineWidth',3,'DisplayName','Newtonian: g_{dyn}=g_{bar}');
plot(g_ref*1e10, RAR_empirical(g_ref)*1e10, 'k-','LineWidth',3,'DisplayName','Empirical RAR');
plot(g_ref*1e10, sqrt(g_ref*1e10)*sqrt(1.2), 'k-.','LineWidth',2.5,'DisplayName','deep-MOND: g_{dyn}=\surd(g_{bar}a_0)');

set(gca,'XScale','log','YScale','log','FontSize',10,'LineWidth',1);
xlabel('g_{bar} [10^{-10} m/s^2]','FontSize',16,'FontWeight','bold');
ylabel('g_{dyn} [10^{-10} m/s^2]','FontSize',16,'FontWeight','bold');
legend('Location','southeast','FontSize',10);
grid on; box on; axis([0.01 100 0.01 100]);

box_pos = [0.15 0.64 0.28 0.28];
box_text = {...
    'Baryonic Scale Shift',...
    '',...
    'DEFINITION:',...
    '  g_{bar} = \color{blue}GM_{True}\color{black} / \color{red}r_{obs}^2\color{black} ',...
    '  g_{dyn} = \color{red}v_{obs}^2\color{black} / \color{red}r_{obs}\color{black} ',...
    '',...
    sprintf('Scatter vs deep-MOND: %.3f dex',scatter_global_C3_mond),...
};
annotation('textbox',box_pos,'String',box_text,...
    'FontSize',11,'FontName','Consolas','BackgroundColor','white',...
    'EdgeColor',[0.6 0.2 0.6],'LineWidth',1.5,'FitBoxToText','on');

exportgraphics(fig, fullfile(rarDir, 'Fig_4_3_3_Case3_BaryonicShift.png'), 'Resolution', 300);
close(fig);

%% --- 4.4 C2_Empirical Residual Analysis ---
% Reference function for C2_Empirical
RAR_empirical = @(g) g ./ (1 - exp(-sqrt(g / a0_RAR)));

% Case configuration
case_name = 'C2_Empirical';
g_bar_field = 'g_bar_C2';
g_dyn_field = 'g_dyn_C2';
ref_func = RAR_empirical;
case_color = [0.9 0.5 0.1];  % Orange

%% Collect all valid data for C2_Empirical (including galaxy mass)
all_g_bar = []; all_g_dyn = []; all_residuals = [];
all_galaxy_mass = []; all_galaxy_idx = [];

for g = 1:numGalaxies
    if ~isfield(RAR_Analysis(g), 'valid') || ~RAR_Analysis(g).valid
        continue;
    end
    if ~isfield(RAR_Analysis(g), g_bar_field) || ~isfield(RAR_Analysis(g), g_dyn_field)
        continue;
    end
    
    % Get accelerations from RAR_Analysis (already filtered)
    g_bar = RAR_Analysis(g).(g_bar_field);
    g_dyn = RAR_Analysis(g).(g_dyn_field);
    
    g_bar = g_bar(:);
    g_dyn = g_dyn(:);
       
    % Get galaxy total mass (use enclosed mass at outermost radius as proxy)
    if isfield(Results(g), 'M_enc_obs') && ~isempty(Results(g).M_enc_obs)
        gal_mass = max(Results(g).M_enc_obs(~isnan(Results(g).M_enc_obs)));
    elseif isfield(Results(g), 'M_enc_true') && ~isempty(Results(g).M_enc_true)
        gal_mass = max(Results(g).M_enc_true(~isnan(Results(g).M_enc_true)));
    else
        gal_mass = NaN;
    end
    
    % Validity mask (only g_bar and g_dyn)
    valid = (g_bar > 1e-15) & (g_bar < 1e-5) & ...
        (g_dyn > 1e-15) & (g_dyn < 1e-5) & ...
        isfinite(g_bar) & isfinite(g_dyn);
    
    if sum(valid) < 5
        continue;
    end
    
    % Compute residuals in log-space vs reference curve
    g_pred = ref_func(g_bar(valid));
    residuals = log10(g_dyn(valid)) - log10(g_pred);
    
    % Store with galaxy metadata
    N_valid = sum(valid);
    all_g_bar = [all_g_bar; g_bar(valid)];
    all_g_dyn = [all_g_dyn; g_dyn(valid)];
    all_residuals = [all_residuals; residuals];
    all_galaxy_mass = [all_galaxy_mass; gal_mass * ones(N_valid, 1)];
    all_galaxy_idx = [all_galaxy_idx; g * ones(N_valid, 1)];
end

N = length(all_residuals);
if N < 50
    error('[ERROR] Insufficient data (%d points) for C2_Empirical analysis\n', N);
end

%% Compute residual statistics
sigma_tot = std(all_residuals);
mean_res = mean(all_residuals);
median_res = median(all_residuals);
mad_res = mad(all_residuals, 1);

%% Collect Radius-Matched Data for C2_Empirical
radius_field_to_use = 'r_obs';
all_g_bar_r = []; all_residuals_r = []; all_radius_r = [];
all_galaxy_mass_r = []; all_galaxy_idx_r = [];

for g = 1:numGalaxies
    if ~isfield(RAR_Analysis(g), 'valid') || ~RAR_Analysis(g).valid
        continue;
    end
    if ~isfield(RAR_Analysis(g), g_bar_field) || ~isfield(RAR_Analysis(g), g_dyn_field)
        continue;
    end
    if ~isfield(Results(g), radius_field_to_use) || ~isfield(Results(g), 'M_enc_obs')
        continue;
    end
    
    % Get original data from Results (before filtering)
    r_orig = Results(g).(radius_field_to_use);
    M_orig = Results(g).M_enc_obs;
    if isempty(M_orig), M_orig = Results(g).M_enc_true; end
    v_orig = Results(g).v_obs;
    
    r_orig = r_orig(:); M_orig = M_orig(:); v_orig = v_orig(:);
    
    % Reconstruct the validity mask (matching main code logic)
    valid_reconstructed = (r_orig > 1e3) & (M_orig > 0) & ...
        isfinite(r_orig) & isfinite(M_orig) & isfinite(v_orig);
    
    if sum(valid_reconstructed) < 5
        continue;
    end
    
    % Compute accelerations from original data
    g_bar_orig = G * M_orig(valid_reconstructed) ./ r_orig(valid_reconstructed).^2;
    g_dyn_orig = v_orig(valid_reconstructed).^2 ./ r_orig(valid_reconstructed);
    
    % Apply acceleration range filter (matching RAR_Analysis filtering)
    valid_accel = (g_bar_orig > 1e-15) & (g_bar_orig < 1e-5) & ...
        (g_dyn_orig > 1e-15) & (g_dyn_orig < 1e-5) & ...
        isfinite(g_bar_orig) & isfinite(g_dyn_orig);
    
    if sum(valid_accel) < 5
        continue;
    end
    
    % Extract matched data
    r_matched = r_orig(valid_reconstructed);
    r_matched = r_matched(valid_accel);
    g_bar_matched = g_bar_orig(valid_accel);
    g_dyn_matched = g_dyn_orig(valid_accel);
    
    % Get galaxy total mass
    if ~isempty(M_orig(valid_reconstructed))
        gal_mass = max(M_orig(valid_reconstructed));
    else
        gal_mass = NaN;
    end
    
    % Compute residuals
    g_pred = ref_func(g_bar_matched);
    residuals = log10(g_dyn_matched) - log10(g_pred);
    
    % Validity check on residuals
    valid_res = isfinite(residuals) & isfinite(r_matched) & (r_matched > 1e3);
    
    if sum(valid_res) < 5
        continue;
    end
    
    % Store matched data
    all_g_bar_r = [all_g_bar_r; g_bar_matched(valid_res)];
    all_residuals_r = [all_residuals_r; residuals(valid_res)];
    all_radius_r = [all_radius_r; r_matched(valid_res)];
    all_galaxy_mass_r = [all_galaxy_mass_r; gal_mass * ones(sum(valid_res), 1)];
    all_galaxy_idx_r = [all_galaxy_idx_r; g * ones(sum(valid_res), 1)];
end

N_radius = length(all_residuals_r);
if N_radius < 50
    fprintf('[WARN] Insufficient radius-matched data (%d points); skipping radius plot\n', N_radius);
else
  
    fig = figure('Color','white','Position',[100 100 900 600],...
        'Name','Fig_4_4_Res_vs_Radius_C2_Empirical');
    
    % Scatter colored by galaxy mass - RESIDUALS ON X-AXIS
    valid_mass_r = isfinite(all_galaxy_mass_r) & (all_galaxy_mass_r > 0);
    if sum(valid_mass_r) > 0
        % Get global mass list (prefer galaxySummary; fallback to plotted galaxies)
        if exist('galaxySummary', 'var') && isfield(galaxySummary, 'M_total_obs_Msun')
            global_masses = [galaxySummary.M_total_obs_Msun];
            global_log_masses = log10(global_masses);
        else
            % Fallback: unique masses from plotted points (per galaxy)
            [unique_gal, ~, ic] = unique(all_galaxy_idx_r(valid_mass_r));
            unique_masses = accumarray(ic, all_galaxy_mass_r(valid_mass_r), [], @max);
            global_log_masses = log10(unique_masses);
        end
        
        [sorted_log, sort_order] = sort(global_log_masses);
        num_global = length(global_log_masses);
        
        cmap = parula(num_global);
        
        % Map each galaxy index to a color based on its mass rank
        color_per_galaxy = zeros(num_global, 3);
        for i = 1:num_global
            gal_idx = sort_order(i);          % galaxy with i-th smallest mass
            color_per_galaxy(gal_idx, :) = cmap(i, :);
        end
        
        % Assign colors to each point (only valid points)
        point_colors = color_per_galaxy(all_galaxy_idx_r(valid_mass_r), :);
        
        % Scatter plot using RGB colors
        sc = scatter(all_residuals_r(valid_mass_r), all_radius_r(valid_mass_r)/1e3, 20, ...
            point_colors, 'filled', ...
            'MarkerFaceAlpha', 0.6, 'MarkerEdgeAlpha', 0.3);
        
        colormap(parula);   
        
        cb = colorbar;
        ylabel(cb, 'log_{10}(M_{obs}/M_{☉})', 'FontSize', 14);
        global_min_log = min(global_log_masses);
        global_max_log = max(global_log_masses);
        
        % Ticks in normalized units (0 to 1) but labels in actual log-mass
        cb.Ticks = linspace(0, 1, 5);
        cb.TickLabels = arrayfun(@(x) sprintf('%.1f', x), ...
            linspace(global_min_log, global_max_log, 5), ...
            'UniformOutput', false);
    else
        % Fallback: single color if no mass data
        scatter(all_residuals_r, all_radius_r/1e3, 20, case_color, 'filled', ...
            'MarkerFaceAlpha', 0.5);
        fprintf('[WARN] No valid mass data for coloring; using single color\n');
    end
    
    hold on;
    set(gca, 'YScale', 'log');
    xlabel('Residual: log_{10}(g_{obs} / g_{RAR}) = log_{10}(g_{dyn}) - log_{10}[g_{bar}/(1-e^{-\surd(g_{bar}/a_0)})] [dex]', ...
        'FontSize', 12, 'FontWeight', 'bold', 'Interpreter', 'tex');
    ylabel('r_{obs} [km]', 'FontSize', 14, 'FontWeight', 'bold');
    
    % Reference lines
    xline(0, 'r--', 'LineWidth', 1.5);
    xline(median_res, 'b:', 'LineWidth', 1, 'Label', sprintf('Median = %.2f', median_res));
    
    % Statistical info box
    stats_text = {
        sprintf('N = %d points', N_radius),
        sprintf('σ = %.3f dex', sigma_tot),
        sprintf('Median = %+.2f dex', median_res)
    };
    annotation('textbox', [0.64 0.52 0.28 0.38], 'String', stats_text, ...
        'FontSize', 11, 'FontName', 'Consolas', 'BackgroundColor', 'white', ...
        'EdgeColor', case_color, 'LineWidth', 1.5, 'FitBoxToText', 'on');
    
    grid on; box on;
    
    exportgraphics(fig, fullfile(rarDir, 'Fig_4_4_Res_vs_Radius_C2_Empirical.png'), ...
        'Resolution', 300);
    close(fig);
end

%% --- 4.5 Empirical Mapping Acceleration Analysis ---
 
% analysis directory
ratioDir = fullfile(baseDir, 'AccelerationRatioAnalysis');
if ~exist(ratioDir, 'dir')
    mkdir(ratioDir);
end

% Initialize results structure
RatioAnalysis = struct();

%% STEP 4.5.1: PER-GALAXY RATIO CALCULATION
for g = 1:numGalaxies
    if ~isfield(RAR_Analysis(g), 'valid') || ~RAR_Analysis(g).valid
        RatioAnalysis(g).valid = false;
        continue;
    end
    
    % Extract Case 2 (Empirical Mapping) accelerations
    g_bar_C = RAR_Analysis(g).g_bar_C2;   % True baryonic acceleration
    g_dyn_C = RAR_Analysis(g).g_dyn_C2;   % Observed dynamical acceleration
    r_norm = RAR_Analysis(g).r_norm;      % Normalized radius
    
    % Column vectors
    g_bar_C = g_bar_C(:); g_dyn_C = g_dyn_C(:); r_norm = r_norm(:);
    
    % Compute the ratio R = g_dyn² / g_bar
    R_ratio = (g_dyn_C.^2) ./ g_bar_C;
    
    % Valid data mask
    valid_ratio = (g_bar_C > 1e-15) & (R_ratio > 1e-15) & (R_ratio < 1e-5) & ...
        isfinite(R_ratio) & isfinite(g_bar_C) & isfinite(g_dyn_C);
    
    if sum(valid_ratio) < 5
        RatioAnalysis(g).valid = false;
        continue;
    end
    
    % Apply mask
    g_bar_valid = g_bar_C(valid_ratio);
    g_dyn_valid = g_dyn_C(valid_ratio);
    R_valid = R_ratio(valid_ratio);
    r_norm_valid = r_norm(valid_ratio);
    
    %% Characterize ratio behavior
    [g_bar_sorted, sort_idx] = sort(g_bar_valid);
    R_sorted = R_valid(sort_idx);
    r_norm_sorted = r_norm_valid(sort_idx);
    
    % Regime masks
    high_g_mask = g_bar_sorted > a0_RAR;
    low_g_mask = g_bar_sorted <= a0_RAR;
    
    % Initialize statistics
    R_high_mean = NaN; R_low_mean = NaN;
    deviation_high = NaN; deviation_low = NaN;
    
    if sum(high_g_mask) > 3
        R_high_mean = mean(R_sorted(high_g_mask));
        g_bar_high_mean = mean(g_bar_sorted(high_g_mask));
        deviation_high = abs(R_high_mean - g_bar_high_mean) / g_bar_high_mean;
    end
    
    if sum(low_g_mask) > 3
        R_low_mean = mean(R_sorted(low_g_mask));
        deviation_low = abs(R_low_mean - a0_RAR) / a0_RAR;
    end
    
    % Overall statistics
    R_mean = mean(R_valid);
    R_median = median(R_valid);
    R_over_a0 = R_valid / a0_RAR;
    R_over_a0_mean = mean(R_over_a0);
    
    % Power-law slope
    slope_R_vs_gbar = NaN; intercept_R_vs_gbar = NaN; R2_powerlaw = NaN;
    if length(g_bar_valid) > 10
        log_g_bar = log10(g_bar_valid);
        log_R = log10(R_valid);
        p = polyfit(log_g_bar, log_R, 1);
        slope_R_vs_gbar = p(1);
        intercept_R_vs_gbar = p(2);
        R_fit = 10.^polyval(p, log_g_bar);
        SS_tot = sum((log_R - mean(log_R)).^2);
        R2_powerlaw = 1 - sum((log_R - log10(R_fit)).^2) / SS_tot;
    end
    
    %% Store results
    RatioAnalysis(g).valid = true;
    RatioAnalysis(g).galaxy = galaxy_names{g};
    RatioAnalysis(g).N_points = length(R_valid);
    RatioAnalysis(g).g_bar_C = g_bar_valid;
    RatioAnalysis(g).g_dyn_C = g_dyn_valid;
    RatioAnalysis(g).R_ratio = R_valid;
    RatioAnalysis(g).R_over_a0 = R_over_a0;
    RatioAnalysis(g).r_norm = r_norm_valid;
    RatioAnalysis(g).R_mean = R_mean;
    RatioAnalysis(g).R_median = R_median;
    RatioAnalysis(g).R_over_a0_mean = R_over_a0_mean;
    RatioAnalysis(g).R_high_g_mean = R_high_mean;
    RatioAnalysis(g).R_low_g_mean = R_low_mean;
    RatioAnalysis(g).slope_R_vs_gbar = slope_R_vs_gbar;
    RatioAnalysis(g).intercept_R_vs_gbar = intercept_R_vs_gbar;
    RatioAnalysis(g).R2_powerlaw_fit = R2_powerlaw;
    
    % Classification
    if isnan(slope_R_vs_gbar)
        regime_classification = 'Insufficient data';
    elseif slope_R_vs_gbar > 0.7
        regime_classification = 'Newtonian-dominated';
    elseif slope_R_vs_gbar < 0.3
        regime_classification = 'MOND-asymptotic';
    else
        regime_classification = 'Transition';
    end
    RatioAnalysis(g).regime_classification = regime_classification;
end

%% STEP 4.5.2: ENSEMBLE ANALYSIS — COLLECT ALL DATA BEFORE FIGURES
% Collect all valid data across galaxies — MUST RUN BEFORE FIGURES
all_g_bar_C = []; all_g_dyn_C = []; all_R_ratio = [];
all_R_over_a0 = []; all_r_norm = []; all_galaxy_idx = [];

for g = 1:numGalaxies
    if isfield(RatioAnalysis(g), 'valid') && RatioAnalysis(g).valid
        all_g_bar_C = [all_g_bar_C; RatioAnalysis(g).g_bar_C];
        all_g_dyn_C = [all_g_dyn_C; RatioAnalysis(g).g_dyn_C];
        all_R_ratio = [all_R_ratio; RatioAnalysis(g).R_ratio];
        all_R_over_a0 = [all_R_over_a0; RatioAnalysis(g).R_over_a0];
        all_r_norm = [all_r_norm; RatioAnalysis(g).r_norm];
        all_galaxy_idx = [all_galaxy_idx; g*ones(length(RatioAnalysis(g).R_ratio), 1)];
    end
end

N_total_ratio = length(all_R_ratio);
fprintf('[RESULT] Total valid measurements: %d\n', N_total_ratio);

if N_total_ratio == 0
    fprintf('[ERROR] No valid data for ensemble analysis!\n');
    return;
end

% Global statistics
global_R_median = median(all_R_ratio);
global_R_over_a0_median = median(all_R_over_a0);
global_R_over_a0_mean = mean(all_R_over_a0);

% Reference regime fractions
high_g_global = all_g_bar_C > a0_RAR;
low_g_global = all_g_bar_C <= a0_RAR;
fprintf('[RESULT] Low-g points (g_bar ≤ a₀): %d (%.1f%%)\n', sum(low_g_global), 100*sum(low_g_global)/N_total_ratio);

if sum(low_g_global) > 10
    R_low_global = all_R_ratio(low_g_global);
end

% Global power-law fit — compute locally for figure use
slope_global = NaN; intercept_global = NaN;
if length(all_g_bar_C) > 10
    log_g_all = log10(all_g_bar_C);
    log_R_all = log10(all_R_ratio);
    p_global = polyfit(log_g_all, log_R_all, 1);
    slope_global = p_global(1);
    intercept_global = p_global(2);
end

%% PREPARE MASS-BASED COLORING
% Extract galaxy masses for ALL figures (compute once)

global_log_masses = log10(galaxy_masses);
global_min_log = min(global_log_masses);
global_max_log = max(global_log_masses);

cmap = parula(numGalaxies);

% Map each galaxy index to a color based on its mass rank
[~, sort_order] = sort(global_log_masses, 'ascend');
color_per_galaxy = zeros(numGalaxies, 3);
for i = 1:numGalaxies
    gal_idx = sort_order(i);          
    color_per_galaxy(gal_idx, :) = cmap(i, :);
end

%% Figure 4.5.1: R/a₀ vs Radius — Radial Behavior of Mezzi Scale Factor ζ(r)
fig = figure('Color', 'white', 'Position', [100 100 1200 600], ...
    'Name', 'Fig_4_5_1_R_over_a0_RadialProfiles');
ax_main = axes('Position', [0.10 0.15 0.60 0.75]);
hold(ax_main, 'on');

% Collect radial data
all_r = []; all_R = [];
inner_R = []; outer_R = [];
mass_trend_low = []; mass_trend_high = [];
n_valid_galaxies = 0;

for g = 1:numGalaxies
    if isfield(RatioAnalysis(g), 'valid') && RatioAnalysis(g).valid
        r_n = RatioAnalysis(g).r_norm;
        R_norm = RatioAnalysis(g).R_over_a0;
        n_valid_galaxies = n_valid_galaxies + 1;
        all_r = [all_r; r_n(:)]; all_R = [all_R; R_norm(:)];
        inner_mask = r_n < 0.3; outer_mask = r_n > 0.3;
        inner_R = [inner_R; R_norm(inner_mask)];
        outer_R = [outer_R; R_norm(outer_mask)];
        if global_log_masses(g) < median(global_log_masses)
            mass_trend_low = [mass_trend_low; R_norm(:)];
        else
            mass_trend_high = [mass_trend_high; R_norm(:)];
        end
        % Plot with mass-based coloring 
        [r_sorted, idx] = sort(r_n);
        R_sorted = R_norm(idx);
        this_color = color_per_galaxy(g, :);
        plot(ax_main, r_sorted, R_sorted, 'LineWidth', 1.5, ...
            'Color', this_color, 'HandleVisibility', 'off');
    end
end

if n_valid_galaxies > 0 && ~isempty(all_R)
    % Calculate statistics
    median_inner_R = median(inner_R, 'omitnan');
    median_outer_R = median(outer_R, 'omitnan');
    transition_ratio = median_inner_R / median_outer_R;
    
    % Reference lines
    yline(ax_main, 1, 'r--', 'LineWidth', 2, 'DisplayName', 'MOND: R/a_0 = 1');
    yline(ax_main, 10, 'k:', 'LineWidth', 2, 'DisplayName', 'Newtonian: R/a_0 = 10');
    
    % COLORBAR
    cbar = colorbar(ax_main, 'Location', 'eastoutside');
    colormap(ax_main, parula);
    cbar.Label.String = 'log_{10}(M_{obs}/M_{☉})';
    cbar.Label.FontSize = 14;
    cbar.Label.FontWeight = 'bold';
    
    % Ticks in normalized units (0 to 1) but labels in actual log-mass
    cbar.Ticks = linspace(0, 1, 5);
    cbar.TickLabels = arrayfun(@(x) sprintf('%.1f', x), ...
        linspace(global_min_log, global_max_log, 5), 'UniformOutput', false);
    cbar.TickLength = 0;  
    
    xlabel(ax_main, 'Normalized Radius r/r_{max}', 'FontSize', 16, 'FontWeight', 'bold');
    ylabel(ax_main, 'R/a_0 = (g_{dyn}^2/g_{bar})/a_0', 'FontSize', 16, 'FontWeight', 'bold');
    legend(ax_main, {'MOND: R/a_0 = 1',  'Newtonian: R/a_0 = 10'}, ...
        'Location', 'best', 'FontSize', 11);
    
    set(ax_main, 'YScale', 'log', 'FontSize', 12, 'LineWidth', 1.2);
    grid(ax_main, 'on'); box(ax_main, 'on');
    
    exportgraphics(fig, fullfile(ratioDir, 'Fig_4_5_1_R_over_a0_RadialProfiles.png'), 'Resolution', 300);
    close(fig);
end

%% Figure 4.5.2: R/a₀ vs Radius — Binned Statistical Analysis
% Initialize collectors
all_r = []; all_R = []; all_gal_idx = []; all_gal_logmass = [];

for g = 1:numGalaxies
    if isfield(RatioAnalysis(g), 'valid') && RatioAnalysis(g).valid
        r_n = RatioAnalysis(g).r_norm;
        R_norm = RatioAnalysis(g).R_over_a0;
        
        % Validate data
        valid = isfinite(r_n) & isfinite(R_norm) & (R_norm > 0);
        r_n = r_n(valid); R_norm = R_norm(valid);
        
        if length(r_n) < 5
            continue;
        end
        
        % Store with galaxy tracking
        all_r = [all_r; r_n(:)];
        all_R = [all_R; R_norm(:)];
        all_gal_idx = [all_gal_idx; g * ones(sum(valid), 1)];
        all_gal_logmass = [all_gal_logmass; global_log_masses(g) * ones(sum(valid), 1)];
    end
end

N_total = length(all_R);
fprintf('[RESULT] Total valid R/a₀ data points: %d from %d galaxies\n', N_total, length(unique(all_gal_idx)));

%% Save results
save('EmpiricalMapping_ScalingFactor_Analysis.mat', ...
    'RatioAnalysis', 'all_g_bar_C', 'all_g_dyn_C', 'all_R_ratio', ...
    'slope_global', 'intercept_global', 'a0_RAR', ...
    'N_total_ratio', 'global_R_over_a0_mean');

 
%% ===== SECTION 5: BTFR ANALYSIS =====
%
% THEORETICAL FRAMEWORK:
% True physics: Newtonian, M ∝ v² (baryonic mass follows Keplerian scaling)
% Observed physics: Empirical BTFR with characteristic slope M ∝ v⁴
%
% Coordinate Scaling Factor ζ(r) where r_obs = r_true*ζ
%
% CASE DEFINITIONS (mathematical):
%   CASE T:      (M = M_true(r_true),    v = v_Newton(r_true))     [True Frame]
%   CASE E:      (M = M_true(r_true),    v = v_obs(r_obs))         [Empirical Mapping]
%
%% =========================================================================

fprintf(' BTFR Analysis in the True and the Observed Frames.\n');

% Structure to store results
BTFR = struct();
for g = 1:numGalaxies
   
    % Extract data
    r_obs_full = Results(g).r_obs(:);
    r_true_full = Results(g).r_true(:);
    M_obs_full = Results(g).M_enc_obs(:);
    M_true_full = Results(g).M_enc_true(:);
    v_obs_full = Results(g).v_obs(:);
    v_Newton_true_full = Results(g).v_Newton_true(:);
    Mezzi_full = Results(g).Mezzi_th(:);
    
    % Harmonize array sizes
    N = min([length(r_obs_full), length(r_true_full), length(M_obs_full), ...
        length(M_true_full), length(v_obs_full), length(v_Newton_true_full), length(Mezzi_full)]);
    r_obs = r_obs_full(1:N); r_true = r_true_full(1:N);
    M_obs = M_obs_full(1:N); M_true = M_true_full(1:N);
    v_obs = v_obs_full(1:N); v_Newton_true = v_Newton_true_full(1:N);
    Mezzi_local = Mezzi_full(1:N);
    
    % Valid data mask
    valid = (r_obs > 1e3*pc) & (M_obs > 0) & (r_true > 0) & (M_true > 0) & ...
        (v_obs > 0) & (v_Newton_true > 0) & isfinite(Mezzi_local) & ...
        (Mezzi_local > 0) & (Mezzi_local <= 1);
    
    if sum(valid) < 10
        BTFR(g).valid = false;
        fprintf('[SKIP] Galaxy %d: insufficient valid points\n', g);
        continue;
    end
    
    % Apply mask
    r_o = r_obs(valid); r_t = r_true(valid);
    M_o = M_obs(valid); M_t = M_true(valid);
    v_o = v_obs(valid); v_t_newton = v_Newton_true(valid);
    mezzi = Mezzi_local(valid);
    
    % Define outer region threshold (70% of max radius)
    outer_fraction = 0.3;
    r_threshold_o = outer_fraction * max(r_o);
    r_threshold_t = outer_fraction * max(r_t);
    
    % Select points where radius is larger than threshold
    outer_mask_o = (r_o >= r_threshold_o);
    outer_mask_t = (r_t >= r_threshold_t);
    N_outer_o = sum(outer_mask_o);
    N_outer_t = sum(outer_mask_t);
    
    if N_outer_o < 5 || N_outer_t < 5
        BTFR(g).valid = false;
        fprintf('[SKIP] Galaxy %d: insufficient outer points\n',g);
        continue;
    end
    
    %% ================================================================
    % Extract TOTAL mass at outermost point 
    M_tot_obs = M_o(end) / Msun;
    M_tot_true = M_t(end) / Msun;
    
    % Extract FLAT rotation velocity (mean of outer region)
    v_flat_obs = mean(v_o(outer_mask_o)) / 1e3;
    v_flat_true_T = mean(v_t_newton(outer_mask_t)) / 1e3;
    
    % CASE E: Empirical Mapping  
    % M_tot_true (global) vs v_flat_obs ( observed radii)
    v_flat_true_E = v_flat_obs;
    
    %% Convert to log-space
    % Method 2: Total mass  
    logM_obs_tot = log10(M_tot_obs);
    logv_obs_tot = log10(v_flat_obs);
    logM_T_tot = log10(M_tot_true);
    logv_T_tot = log10(v_flat_true_T);
    logM_E_tot = log10(M_tot_true);
    logv_E_tot = log10(v_flat_true_E);
    
    %% Store results
    BTFR(g).valid = true;
    BTFR(g).galaxy = galaxy_names{g};
    BTFR(g).N_outer = sum(valid);
    
    % Total mass method (single values)
    BTFR(g).tot.logM_obs = logM_obs_tot;
    BTFR(g).tot.logv_obs = logv_obs_tot;
    BTFR(g).tot.logM_T = logM_T_tot;
    BTFR(g).tot.logv_T = logv_T_tot;
    BTFR(g).tot.logM_E = logM_E_tot;
    BTFR(g).tot.logv_E = logv_E_tot;
    
    % compression metrics
    BTFR(g).Mezzi = mean(log10(mezzi(outer_mask_t)));
    BTFR(g).M_ratio = M_tot_true / M_tot_obs;
    
    %% Print summary
    fprintf('[BTFR] Galaxy %d/%d: %s\n', g, numGalaxies, galaxy_names{g});
%     fprintf('[RESULT]   v_obs=%.1f v_T=%.1f v_E=%.1f km/s\n', ...
%         v_flat_obs, v_flat_true_T, v_flat_true_E);
%     fprintf('[RESULT]   Total: M_obs=%.2e M_true=%.2e (ratio=%.2f)\n', ...
%         M_tot_obs, M_tot_true, BTFR(g).M_ratio);
end

%% 5.2 STEP B: GLOBAL ANALYSIS - CURVE SHAPE CHARACTERIZATION
%% =========================================================================
%% Aggregate data for both methods
% Total mass  
tot_logM_T = []; tot_logv_T = [];
tot_logM_E = []; tot_logv_E = [];

for g = 1:numGalaxies
    if ~isfield(BTFR(g), 'valid') || ~BTFR(g).valid
        continue;
    end
    
    % Total: append single values
    tot_logM_T = [tot_logM_T; BTFR(g).tot.logM_T];
    tot_logv_T = [tot_logv_T; BTFR(g).tot.logv_T];
    tot_logM_E = [tot_logM_E; BTFR(g).tot.logM_E];
    tot_logv_E = [tot_logv_E; BTFR(g).tot.logv_E];
end

N_gal = length(tot_logM_T);

%% Compute fits for BOTH methods (2 cases: T, E)
calc_fit = @(logv, logM) polyfit(logv, logM, 1);

% TOTAL MASS fits (2 cases)
p_tot_T = calc_fit(tot_logv_T, tot_logM_T);
p_tot_E = calc_fit(tot_logv_E, tot_logM_E);

%% Bootstrap confidence intervals (resample galaxies)  
n_boot = 1000;
boot_tot_T = zeros(n_boot, 1);
boot_tot_E = zeros(n_boot, 1);

for b = 1:n_boot
    gal_idx = randsample(N_gal, N_gal, true);
    
    % Total: collect single values
    b_tot_v_T = []; b_tot_M_T = [];
    b_tot_v_E = []; b_tot_M_E = [];
    
    for g = gal_idx'
        if isfield(BTFR(g), 'valid') && BTFR(g).valid
            % Total
            b_tot_v_T = [b_tot_v_T; BTFR(g).tot.logv_T];
            b_tot_M_T = [b_tot_M_T; BTFR(g).tot.logM_T];
            b_tot_v_E = [b_tot_v_E; BTFR(g).tot.logv_E];
            b_tot_M_E = [b_tot_M_E; BTFR(g).tot.logM_E];
        end
    end
    
    % Fit and store ONLY the slope (first element of polyfit output)
    if length(b_tot_v_T) > 3
        p_fit = polyfit(b_tot_v_T, b_tot_M_T, 1);
        boot_tot_T(b) = p_fit(1);
        p_fit = polyfit(b_tot_v_E, b_tot_M_E, 1);
        boot_tot_E(b) = p_fit(1);
    end
end

% Confidence intervals (2 cases)
ci_tot_T = prctile(boot_tot_T, [2.5 97.5]);
ci_tot_E = prctile(boot_tot_E, [2.5 97.5]);

%% Compute scatter around reference slopes  
calc_scatter_ref = @(logM, logv, slope) std(logM - (slope * logv + mean(logM - slope*logv)));

% Total (2 cases)
s_tot_T_s2 = calc_scatter_ref(tot_logM_T, tot_logv_T, 2.0);
s_tot_T_s4 = calc_scatter_ref(tot_logM_T, tot_logv_T, 4.0);
s_tot_E_s2 = calc_scatter_ref(tot_logM_E, tot_logv_E, 2.0);
s_tot_E_s4 = calc_scatter_ref(tot_logM_E, tot_logv_E, 4.0);

%% Store results
Global_BTFR_Results.tot.slope_T = p_tot_T(1);
Global_BTFR_Results.tot.slope_E = p_tot_E(1);
Global_BTFR_Results.tot.ci_T = ci_tot_T;
Global_BTFR_Results.tot.ci_E = ci_tot_E;
Global_BTFR_Results.tot.scatter_s2 = [s_tot_T_s2, s_tot_E_s2];
Global_BTFR_Results.tot.scatter_s4 = [s_tot_T_s4, s_tot_E_s4];
Global_BTFR_Results.tot.N_gal = N_gal;
Global_BTFR_Results.N_galaxies = N_gal;

save('BTFR_4Case_Analysis.mat', 'BTFR', 'Global_BTFR_Results', '-v7.3');

%% 5.3 STEP C: GENERATE FIGURES WITH SIDE PANEL TEXT
%% =========================================================================

btfrDir = fullfile(baseDir, 'BTFR_Figures');
if ~exist(btfrDir, 'dir')
    mkdir(btfrDir);
end

%% 5.3.2 PART 1: SLOPE EXTRACTION AND BASIC STATISTICS
%% ========================================================================
% Helper function for confidence intervals
calc_ci = @(boot_samples) prctile(boot_samples(~isnan(boot_samples)), [2.5 97.5]);

% Extract slopes and CIs for all cases (2 cases: T, E)
slopes_tot = [p_tot_T(1), p_tot_E(1)];
ci_tot = [ci_tot_T; ci_tot_E];

case_labels = {'True Frame', 'Empirical Mapping'};

% Reference slopes
slope_newtonian = 2.0;
slope_btfr = 4.0;

%% 5.3.8 FIGURE 5.3.8: BTFR SLOPES AND SCATTER ANALYSIS
%% =========================================================================

fig = figure('Color', 'white', 'Position', [100 100 1000 700], ...
    'Name', 'Fig_5_3_8_BTFR_Slopes_Total_Mass_Only');

% --- 1. Define Data Order: T -> E ---
s_tot_values = [s_tot_T_s4, s_tot_E_s4];
slopes_vals  = slopes_tot([1, 2]); 
ci_vals      = ci_tot([1, 2], :);  

% Labels & Colors
case_names_long = {'True Frame', 'Empirical Mapping'};
case_colors     = [0.2 0.4 0.8;     % 1: True Frame
                   0.9 0.5 0.1];    % 2: Empirical Mapping

% Scenario descriptions
scenario_desc = {'M_{True} & v_{Newton} at r_{True}', ...   
                 'M_{True}(r_{True}) vs v_{obs}(r_{obs})'};  

% --- 2. Main Plotting ---
ax_main = axes('Position', [0.15 0.18 0.70 0.72]);
hold on;
x_pos = 1:2;

% Store handles for legend
h = gobjects(3,1);

% Plot total mass with slope values next to scatter point
for i = 1:2
    h(i) = errorbar(x_pos(i), slopes_vals(i), ...
        slopes_vals(i) - ci_vals(i,1), ...
        ci_vals(i,2) - slopes_vals(i), ...
        'o', 'LineWidth', 2.5, 'MarkerSize', 14, ...
        'MarkerFaceColor', case_colors(i,:), ...
        'Color', case_colors(i,:)*0.7, 'CapSize', 10, ...
        'DisplayName', sprintf('%s (scatter=%.3f)', case_names_long{i}, s_tot_values(i)));
    
    slope_text = sprintf('%.2f', slopes_vals(i));
    text(x_pos(i) + 0.1, slopes_vals(i), ...
        slope_text, 'HorizontalAlignment', 'left', ...
        'FontSize', 10, 'FontWeight', 'bold', 'Color', case_colors(i,:), ...
        'VerticalAlignment', 'middle');
end

% --- 3. Reference Lines ---
yline(2.0, '--', 'Color', [0.2 0.7 0.3], 'LineWidth', 2.5, 'Alpha', 0.6);
yline(4.0, '--', 'Color', [0.8 0.2 0.2], 'LineWidth', 2.5, 'Alpha', 0.6);
text(0.5, 2.15, 'Newtonian (M ∝ v²)', 'Color', [0.2 0.7 0.3], ...
    'FontSize', 10, 'FontWeight', 'bold');
text(0.5, 4.15, 'BTFR (M ∝ v⁴)', 'Color', [0.8 0.2 0.2], ...
    'FontSize', 10, 'FontWeight', 'bold');

% --- 4. Formatting ---
set(gca, 'XTick', x_pos, 'XTickLabel', case_names_long, ...
    'FontSize', 11, 'XTickLabelRotation', 0);
ylabel('Fitted Slope (log M vs log v)', 'FontSize', 13, 'FontWeight', 'bold');

for i = 1:2
    text(x_pos(i), 0.2, scenario_desc{i}, ...
        'HorizontalAlignment', 'center', ...
        'FontSize', 8, 'Color', [0.4 0.4 0.4], ...
        'Interpreter', 'tex');
end

ylim([0 6.5]);
xlim([0.5 2.5]);
grid on; box on;

% --- 5. Legend ---
h(3) = plot(NaN, NaN, 'k-', 'LineWidth', 2, ...
    'DisplayName', 'Error bars = 95% Confidence Interval');

lgd = legend(h, 'Location', 'north', 'FontSize', 10, ...
    'EdgeColor', [0.5 0.5 0.5]);
hold off;

filePath = fullfile(btfrDir, 'Fig_5_3_8_BTFR_Slopes_Total_Mass.png');
exportgraphics(fig, filePath, 'Resolution', 300);
close(fig);

%% 5.3.8 FIGURE 5.3.8: LOG V VS LOG M SCATTER WITH SLOPE FIT
%% =========================================================================
fig = figure('Color', 'white', 'Position', [100 100 1200 550], ...
    'Name', 'Fig_5_3_9_LogV_vs_LogM');

v_range = linspace(min([tot_logv_T; tot_logv_E]), max([tot_logv_T; tot_logv_E]), 100);
scenario_desc_plots = {'M_{True} & v_{Newton} at r_{True}', ...   
                       'M_{True}(r_{True}) vs v_{obs}(r_{obs})'};  

% --- Subplot 1: True Frame ---
ax_T = subplot(1,2,1);
hold(ax_T, 'on');

% Scatter
h_scatter_T = scatter(ax_T, tot_logv_T, tot_logM_T, 50, [0.2 0.4 0.8], 'filled', 'MarkerFaceAlpha', 0.7, 'MarkerEdgeColor', 'none', 'DisplayName', sprintf('Scatter vs v^{2} = %.2f dex', s_tot_T_s4));

% Fit Line (Black)
h_fit_T = plot(ax_T, v_range, polyval(p_tot_T, v_range), '-', 'Color', 'k', 'LineWidth', 2.5, ...
    'DisplayName', sprintf('Fit (slope=%.2f)', p_tot_T(1)));

% Reference Lines
v_range_ref = linspace(min(tot_logv_T), max(tot_logv_T), 100);
h_s2_T = plot(ax_T, v_range_ref, 2.0*v_range_ref + mean(tot_logM_T - 2.0*tot_logv_T), '--', 'Color', [0.2 0.7 0.3], 'LineWidth', 2.0, 'DisplayName', 'Newtonian (M \propto v^2)');
h_s4_T = plot(ax_T, v_range_ref, 4.0*v_range_ref + mean(tot_logM_T - 4.0*tot_logv_T), '--', 'Color', [0.8 0.2 0.2], 'LineWidth', 2.0, 'DisplayName', 'BTFR (M \propto v^4)');

% Text Descriptions
text(0.05, 0.85, scenario_desc_plots{1}, 'Units', 'normalized', 'Parent', ax_T, ...
    'BackgroundColor', [1 1 1 0.8], 'EdgeColor', 'none', 'FontSize', 11, 'Interpreter', 'tex', ...
    'FontWeight', 'bold', 'Color', [0.2 0.4 0.8]);

% Legend
legend(ax_T, [h_scatter_T, h_fit_T, h_s2_T, h_s4_T], 'Location', 'southeast', 'FontSize', 11);

grid(ax_T, 'on'); box(ax_T, 'on');
xlabel(ax_T, 'log_{10}(v_{flat}) [km/s]', 'FontSize', 16);
ylabel(ax_T, 'log_{10}(M_{tot}) [M_{☉}]', 'FontSize', 16);
% title(ax_T, 'True Frame', 'FontSize', 14, 'FontWeight', 'bold');

% --- Subplot 2: Empirical Mapping ---
ax_E = subplot(1,2,2);
hold(ax_E, 'on');

% Scatter
h_scatter_E = scatter(ax_E, tot_logv_E, tot_logM_E, 50, [0.9 0.5 0.1], 'filled', 'MarkerFaceAlpha', 0.7, 'MarkerEdgeColor', 'none', 'DisplayName', sprintf('Scatter vs v^{4}= %.2f dex', s_tot_E_s4));

% Fit Line (Black)
h_fit_E = plot(ax_E, v_range, polyval(p_tot_E, v_range), '-', 'Color', 'k', 'LineWidth', 2.5, ...
    'DisplayName', sprintf('Fit (slope=%.2f)', p_tot_E(1)));

% Reference Lines
v_range_ref = linspace(min(tot_logv_E), max(tot_logv_E), 100);
h_s2_E = plot(ax_E, v_range_ref, 2.0*v_range_ref + mean(tot_logM_E - 2.0*tot_logv_E), '--', 'Color', [0.2 0.7 0.3], 'LineWidth', 2.0, 'DisplayName', 'Newtonian (M \propto v^2)');
h_s4_E = plot(ax_E, v_range_ref, 4.0*v_range_ref + mean(tot_logM_E - 4.0*tot_logv_E), '--', 'Color', [0.8 0.2 0.2], 'LineWidth', 2.0, 'DisplayName', 'BTFR (M \propto v^4)');

% Text Descriptions
text(0.05, 0.85, scenario_desc_plots{2}, 'Units', 'normalized', 'Parent', ax_E, ...
    'BackgroundColor', [1 1 1 0.8], 'EdgeColor', 'none', 'FontSize', 11, 'Interpreter', 'tex', ...
    'FontWeight', 'bold', 'Color', [0.9 0.5 0.1]);

% Legend
legend(ax_E, [h_scatter_E, h_fit_E, h_s2_E, h_s4_E], 'Location', 'southeast', 'FontSize', 11);

grid(ax_E, 'on'); box(ax_E, 'on');
xlabel(ax_E, 'log_{10}(v_{flat}) [km/s]', 'FontSize', 16);
ylabel(ax_E, 'log_{10}(M_{tot}) [M_{☉}]', 'FontSize', 16);
% title(ax_E, 'Empirical Mapping', 'FontSize', 14, 'FontWeight', 'bold');

filePath_log = fullfile(btfrDir, 'Fig_5_3_9_LogV_LogM_Scatter.png');
exportgraphics(fig, filePath_log, 'Resolution', 300);
close(fig);

%% 5.3.10 FIGURE 5.3.10: BTFR β-OFFSET VS MEZZI SCALE FACTOR (ζ)
%% 5.3.11 DATA EXTRACTION: Compute β-offsets for Observed and Empirical Mapping
%% ------------------------------------------------------------------------
residuals_obs = [];   % β-offset: log M_o - 4 log v_o 
residuals_E   = [];   % β-offset: log M_t - 4 log v_o 
Mezzi_values  = [];   % ζ = 10^Mezzi
galaxy_masses = [];   % Stellar mass for categorization

for g = 1:numGalaxies
    if ~isfield(BTFR(g),'valid') || ~BTFR(g).valid
        continue;
    end
    
    % Observed frame β-offset
    if isfield(BTFR(g).tot,'logM_obs') && isfield(BTFR(g).tot,'logv_obs')
        residuals_obs = [residuals_obs; BTFR(g).tot.logM_obs - 4*BTFR(g).tot.logv_obs];
    else
        residuals_obs = [residuals_obs; NaN];
    end
    
    % Empirical Mapping β-offset
    if isfield(BTFR(g).tot,'logM_E') && isfield(BTFR(g).tot,'logv_E')
        residuals_E = [residuals_E; BTFR(g).tot.logM_E - 4*BTFR(g).tot.logv_E];
    else
        residuals_E = [residuals_E; NaN];
    end
    
    % Compression parameter ζ
    if isfield(BTFR(g),'Mezzi')
        Mezzi_values = [Mezzi_values; 10^BTFR(g).Mezzi];
    else
        Mezzi_values = [Mezzi_values; NaN];
    end
    
    % Galaxy mass
    if isfield(Results(g),'Mtot_obs')
        galaxy_masses = [galaxy_masses; Results(g).Mtot_obs/Msun];
    else
        galaxy_masses = [galaxy_masses; NaN];
    end
end

%% 5.3.12 MASS CATEGORIZATION: Three mass bins with distinct colors
%% ------------------------------------------------------------------------
log_masses = log10(galaxy_masses);
mass_bin_edges = [8.5, 9.5, 10.5, 12.0];
mass_bin_labels = {'Low Mass', 'Intermediate', 'High Mass'};
mass_bin_colors = [0.8 0.3 0.2;   
                   0.6 0.4 0.8;    
                   0.2 0.5 0.8];   

mass_category = zeros(size(log_masses));
for i = 1:length(log_masses)
    if ~isnan(log_masses(i))
        if log_masses(i) < mass_bin_edges(2)
            mass_category(i) = 1;
        elseif log_masses(i) < mass_bin_edges(3)
            mass_category(i) = 2;
        else
            mass_category(i) = 3;
        end
    end
end

%% 5.3.15 STATISTICAL ANALYSIS: Compute comprehensive metrics
%% ------------------------------------------------------------------------
compute_fit = @(alpha, res) struct( ...
    'slope', polyfit(log10(alpha), res, 1), ...
    'r', corr(log10(alpha), res), ...
    'mean_res', mean(res, 'omitnan'), ...  % Average β-offset
    'std_res', std(res, 'omitnan'), ...
    'N', sum(~isnan(alpha) & ~isnan(res)) );

valid_obs = ~isnan(Mezzi_values) & ~isnan(residuals_obs);
stats_obs = compute_fit(Mezzi_values(valid_obs), residuals_obs(valid_obs));
[~, pval_obs] = corr(log10(Mezzi_values(valid_obs)), residuals_obs(valid_obs));

valid_E = ~isnan(Mezzi_values) & ~isnan(residuals_E);
stats_E = compute_fit(Mezzi_values(valid_E), residuals_E(valid_E));
[~, pval_E] = corr(log10(Mezzi_values(valid_E)), residuals_E(valid_E));

%% 5.3.13 COMPUTE SCATTER AND SCATTER WITHIN EACH MASS CATEGORY
%% ------------------------------------------------------------------------
N_bins = length(mass_bin_labels);
bin_stats_obs = struct('sigma', nan(1,N_bins), 'scatter', nan(1,N_bins), ...
    'N', zeros(1,N_bins), 'centroid_zeta', nan(1,N_bins), 'centroid_res', nan(1,N_bins));
bin_stats_E = struct('sigma', nan(1,N_bins), 'scatter', nan(1,N_bins), ...
    'N', zeros(1,N_bins), 'centroid_zeta', nan(1,N_bins), 'centroid_res', nan(1,N_bins));

for b = 1:N_bins
    cat_mask = (mass_category == b);
    
    % --- OBSERVED FRAME ---
    valid_obs_bin = ~isnan(Mezzi_values) & ~isnan(residuals_obs) & cat_mask;
    if sum(valid_obs_bin) >= 3
        % Scatter around slope-4 reference (for comparison)
        bin_stats_obs(b).sigma = std(residuals_obs(valid_obs_bin), 'omitnan');
        bin_stats_obs(b).N = sum(valid_obs_bin);
        
        % Centroid for tag positioning
        bin_stats_obs(b).centroid_zeta = median(Mezzi_values(valid_obs_bin));
        bin_stats_obs(b).centroid_res = mean(residuals_obs(valid_obs_bin), 'omitnan');
        
        % SCATTER: scatter around the FIT LINE β(ζ) = slope·log₁₀(ζ) + β_fit
        fit_vals = polyval(stats_obs.slope, log10(Mezzi_values(valid_obs_bin)));
        bin_stats_obs(b).scatter = std(residuals_obs(valid_obs_bin) - fit_vals, 'omitnan');
    end
    
    % --- EMPIRICAL MAPPING FRAME ---
    valid_E_bin = ~isnan(Mezzi_values) & ~isnan(residuals_E) & cat_mask;
    if sum(valid_E_bin) >= 3
        bin_stats_E(b).sigma = std(residuals_E(valid_E_bin), 'omitnan');
        bin_stats_E(b).N = sum(valid_E_bin);
        bin_stats_E(b).centroid_zeta = median(Mezzi_values(valid_E_bin));
        bin_stats_E(b).centroid_res = mean(residuals_E(valid_E_bin), 'omitnan');
        
        % Scatter around fit line for Case E
        fit_vals = polyval(stats_E.slope, log10(Mezzi_values(valid_E_bin)));
        bin_stats_E(b).scatter = std(residuals_E(valid_E_bin) - fit_vals, 'omitnan');
    end
end


%% 5.3.16 FIGURE 5.3.16A: OBSERVED FRAME - β-OFFSET VS ζ
%% =========================================================================
marker_size_uniform = 80;

fig = figure('Color','white','Position',[150 150 950 550], ...
    'Name','Fig_5_3_16A_BetaOffset_Mezzi_Observed');
ax_obs = axes('Parent',fig,'Position',[0.12 0.15 0.80 0.75]);
hold(ax_obs,'on');

% Plot each mass category with distinct colors
h_bin_obs = gobjects(N_bins, 1);
for b = 1:N_bins
    cat_mask = (mass_category == b) & valid_obs;
    if sum(cat_mask) > 0
        h_bin_obs(b) = scatter(Mezzi_values(cat_mask), residuals_obs(cat_mask), ...
            marker_size_uniform, mass_bin_colors(b,:), ...
            'filled','MarkerFaceAlpha',0.8,'MarkerEdgeColor','white', ...
            'MarkerFaceColor', mass_bin_colors(b,:), 'Parent', ax_obs);
    end
end

% Fit line: β(ζ) = slope·log₁₀(ζ) + β_fit
xfit_obs = linspace(min(log10(Mezzi_values(valid_obs))), ...
    max(log10(Mezzi_values(valid_obs))), 100);
yfit_obs = polyval(stats_obs.slope, xfit_obs);
h_fit_obs = plot(10.^xfit_obs, yfit_obs, '--k', 'LineWidth', 2.0, 'Parent', ax_obs);

% Labels and title with BTFR equation
ylabel(ax_obs, 'Offset β = log_{10}M_{Obs} - 4log_{10}v_{Obs}  [dex]', 'FontSize', 16);
xlabel(ax_obs, 'Mezzi Scale factor \zeta  ', 'FontSize', 16);
grid(ax_obs,'on'); box(ax_obs,'on');

% === TAGS: Position near cluster centroid with scatter values ===
y_all = [residuals_obs(valid_obs); residuals_E(valid_E)];  
y_range = max(y_all) - min(y_all);
if y_range == 0, y_range = 1; end
base_offset = 0.3 * y_range;    
min_gap = 0.5 * base_offset;     

% Store tag positions for overlap check
tag_positions = [];   % [x, y] for each bin that gets a tag
for b = 1:N_bins
    if ~isnan(bin_stats_obs(b).scatter)   
        % Candidate position
        tag_x = bin_stats_obs(b).centroid_zeta;
        
        sign_offset = 1 - 2*mod(b,2);   
        tag_y = bin_stats_obs(b).centroid_res + sign_offset * base_offset;
        
        if ~isempty(tag_positions)
            dist = abs(tag_y - tag_positions(:,2));
            if any(dist < min_gap)
                tag_y = tag_y + sign_offset * base_offset;
            end
        end
        
        tag_positions = [tag_positions; tag_x, tag_y];
        
        tag_text = sprintf('scatter=%.2f dex', bin_stats_obs(b).scatter);
        text(ax_obs, tag_x, tag_y, tag_text, ...
            'HorizontalAlignment','center', 'FontSize',9, ...
            'FontWeight','bold', 'Color', mass_bin_colors(b,:), ...
            'BackgroundColor',[1 1 1 0.9], 'EdgeColor', mass_bin_colors(b,:), ...
            'LineWidth',1.0, 'Margin',3);
    end
end

%% === LEGEND  ===
beta_fit_obs = stats_obs.slope(2);   
beta_avg_obs = stats_obs.mean_res;   

h_dummy_avg_obs = plot(ax_obs, NaN, NaN, 'Color','k', 'LineStyle','none');

legend_entries_obs = { ...
    sprintf('Fit components: slope = %.2f, int = %.2f  |  r = %.2f', ...
        stats_obs.slope(1), beta_fit_obs, stats_obs.r), ...
    sprintf('Offset β = %.2f dex  (all galaxies)', beta_avg_obs), ...
    sprintf('Low Mass  (N=%2d):  scatter = %.2f dex', ...
        bin_stats_obs(1).N, bin_stats_obs(1).scatter), ...
    sprintf('Intermediate (N=%2d):  scatter = %.2f dex', ...
        bin_stats_obs(2).N, bin_stats_obs(2).scatter), ...
    sprintf('High Mass   (N=%2d):  scatter = %.2f dex', ...
        bin_stats_obs(3).N, bin_stats_obs(3).scatter)};

legend_handles_obs = [h_fit_obs, h_dummy_avg_obs, ...
    h_bin_obs(1), h_bin_obs(2), h_bin_obs(3)];

legend(legend_handles_obs, legend_entries_obs, ...
    'Location','best', 'FontSize',11, 'Box','on', ...
    'NumColumns',1, 'Interpreter','tex');

filePath = fullfile(btfrDir, 'Fig_5_3_16A_BetaOffset_Observed.png');
exportgraphics(fig, filePath, 'Resolution', 300);
close(fig);

%% 5.3.17 FIGURE 5.3.17E: EMPIRICAL MAPPING - β-OFFSET VS ζ
%% =========================================================================
fig = figure('Color','white','Position',[150 150 950 550], ...
    'Name','Fig_5_3_17E_BetaOffset_Mezzi_Empirical');
ax_E = axes('Parent',fig,'Position',[0.12 0.15 0.80 0.75]);
hold(ax_E,'on');

% Plot each mass category
h_bin_E = gobjects(N_bins, 1);
for b = 1:N_bins
    cat_mask = (mass_category == b) & valid_E;
    if sum(cat_mask) > 0
        h_bin_E(b) = scatter(Mezzi_values(cat_mask), residuals_E(cat_mask), ...
            marker_size_uniform, mass_bin_colors(b,:), ...
            'filled','MarkerFaceAlpha',0.8,'MarkerEdgeColor','white', ...
            'MarkerFaceColor', mass_bin_colors(b,:), 'Parent', ax_E);
    end
end

% Fit line
xfit_E = linspace(min(log10(Mezzi_values(valid_E))), ...
    max(log10(Mezzi_values(valid_E))), 100);
yfit_E = polyval(stats_E.slope, xfit_E);
h_fit_E = plot(10.^xfit_E, yfit_E, '--k', 'LineWidth', 2.0, 'Parent', ax_E);

% Labels and title
ylabel(ax_E, 'Offset β = log_{10}M_{True} - 4log_{10}v_{Obs}  [dex]', 'FontSize', 16);
xlabel(ax_E, 'Mezzi Scale Factor \zeta ', 'FontSize', 16);
grid(ax_E,'on'); box(ax_E,'on');

% === TAGS: Position near cluster centroid with scatter values ===
y_all = [residuals_obs(valid_obs); residuals_E(valid_E)]; 
y_range = max(y_all) - min(y_all);
if y_range == 0, y_range = 1; end
base_offset = 0.2 * y_range;   
min_gap = 0.4 * base_offset;   

tag_positions = [];   
for b = 1:N_bins
    if ~isnan(bin_stats_E(b).scatter)
        tag_x = bin_stats_E(b).centroid_zeta;
        
        sign_offset = 1 - 2*mod(b,2);  
        tag_y = bin_stats_E(b).centroid_res + sign_offset * base_offset;
        
        if ~isempty(tag_positions)
            dist = abs(tag_y - tag_positions(:,2));
            if any(dist < min_gap)
                tag_y = tag_y + sign_offset * base_offset;
            end
        end

        tag_positions = [tag_positions; tag_x, tag_y];
        
        tag_text = sprintf('scatter=%.2f dex', bin_stats_E(b).scatter);
        text(ax_E, tag_x, tag_y, tag_text, ...
            'HorizontalAlignment','center', 'FontSize',9, ...
            'FontWeight','bold', 'Color', mass_bin_colors(b,:), ...
            'BackgroundColor',[1 1 1 0.9], 'EdgeColor', mass_bin_colors(b,:), ...
            'LineWidth',1.0, 'Margin',3);
    end
end

%% === LEGEND ===
beta_fit_E = stats_E.slope(2);
beta_avg_E = stats_E.mean_res;

h_dummy_avg_E = plot(ax_E, NaN, NaN, 'Color','k', 'LineStyle','none');

legend_entries_E = { ...
    sprintf('Fit components: slope = %.2f, int = %.2f  |  r = %.2f', ...
        stats_E.slope(1), beta_fit_E, stats_E.r), ...
    sprintf('Offset β = %.2f dex  (all galaxies)', beta_avg_E), ...
    sprintf('Low Mass  (N=%2d):  scatter = %.2f dex', ...
        bin_stats_E(1).N, bin_stats_E(1).scatter), ...
    sprintf('Intermediate (N=%2d):  scatter = %.2f dex', ...
        bin_stats_E(2).N, bin_stats_E(2).scatter), ...
    sprintf('High Mass   (N=%2d):  scatter = %.2f dex', ...
        bin_stats_E(3).N, bin_stats_E(3).scatter)};

legend_handles_E = [h_fit_E, h_dummy_avg_E, ...
    h_bin_E(1), h_bin_E(2), h_bin_E(3)];

legend(legend_handles_E, legend_entries_E, ...
    'Location','best', 'FontSize',11, 'Box','on', ...
    'NumColumns',1, 'Interpreter','tex');

filePath = fullfile(btfrDir, 'Fig_5_3_17E_BetaOffset_True.png');
exportgraphics(fig, filePath, 'Resolution', 300);
close(fig);

save('BTFR_4Case_Analysis.mat', 'Global_BTFR_Results', '-append', '-v7.3');












% ===========================================================================
% FUNCTIONS  
% ===========================================================================

function [M_true_cum, M_obs_cum, mass_ratio] = MassProfileFromMass(r_obs, M_enc, ScaleFactor)

%MASSPROFILEFROMMASS Reconstruct true mass profile via geometric projection
%   [M_TRUE_CUM, M_OBS_CUM, MASS_RATIO] = MASSPROFILEFROMMASS(R_OBS, M_ENC, SCALEFACTOR)
%   reconstructs the true mass profile from observed coordinates by applying
%   a geometric projection based on the Mezzi scale factor.
%
%   PHYSICAL MODEL
%   ==============
%   For an axisymmetric system, the annular area element dA transforms 
%   according to the Jacobian of the coordinate transformation. Observed 
%   coordinates r_obs are related to true coordinates r_true by the Mezzi 
%   scale factor ζ(r) (represented by SCALEFACTOR):
%
%       r_true = r_obs / ζ(r_obs)
%
%   Consider a discrete annular shell bounded by radii r_j and r_(j+1). 
%   The observed and true areas of this shell are, respectively:
%
%       Delta_A_obs  = pi * (r_(j+1)^2 - r_j^2)                           (15)
%       Delta_A_true = pi * (r_true_(j+1)^2 - r_true_j^2)                 (16)
%
%   The area ratio R_j = Delta_A_true / Delta_A_obs quantifies the 
%   geometric distortion of the j-th shell. For a Mezzi scale factor 
%   with ζ < 1, the shells expand in true coordinates (R_j > 1).
%
%   DENSITY INVARIANCE APPROXIMATION
%   =================================
%   The reconstruction utilizes a local approximation of density invariance. 
%   We posit that the observed surface density Sigma_obs provides the 
%   correct mass density for the local shell, but the coordinate grid is 
%   compressed, mapping the true extent onto smaller observed radii.
%
%   By projecting this invariant density onto the true, expanded metric, 
%   the true surface density is approximated by the observed value:
%
%       Sigma_true approx Sigma_obs                                        (17)
%
%   Consequently, the mass element for the expanded true shell is 
%   approximated by applying the observed surface density to the larger 
%   true geometry:
%
%       Delta_M_true,j = Sigma_obs,j * Delta_A_true,j                     (19)
%
%   Substituting the area ratio R_j, this can be expressed as:
%
%       Delta_M_true,j = Delta_M_obs,j * R_j                              (19b)
%
%   Since R_j > 1 for the Mezzi effect, the true shell mass exceeds the 
%   observed shell mass.
%
%   ALGORITHM
%   =========
%   Practical implementation employs a discrete shell algorithm:
%     1. Discretize the observed profile into N-1 annular shells.
%     2. Compute observed surface density: Sigma_obs,j = Delta_M_obs,j / Delta_A_obs,j
%     3. Transform shell boundaries: r_true,j = r_obs,j / ζ(r_obs,j)
%     4. Compute true shell areas and area ratios R_j.
%     5. Redistribute mass: Delta_M_true,j = Sigma_obs,j * Delta_A_true,j
%     6. Reconstruct cumulative profile via Eq. (20).
%
%   INPUTS
%   ======
%   R_OBS        [N×1 double] Observed radii r_obs in strictly increasing order.
%   M_ENC        [N×1 double] Enclosed mass M(<r_obs) at observed radii.
%   SCALEFACTOR  [N×1 double] Mezzi scale factor ζ(r_obs). Values < 1 indicate
%                shells expand in true coordinates (R_j > 1).
%
%   OUTPUTS
%   =======
%   M_TRUE_CUM   [N×1 double] Cumulative mass M_true(<r_true) in true coordinates.
%   M_OBS_CUM    [N×1 double] Original observed cumulative mass.
%   MASS_RATIO   [scalar double] Mass ratio η quantifying the net geometric bias:
%
%                η = M_total,true / M_total,obs                            (21)
%
%   VALIDATION AND BOUNDARY CONDITIONS
%   ===================================
%   The innermost enclosed mass M_1 = M_obs(<r_obs,1) is treated as a central
%   boundary condition. The assumption that the central region is unaffected 
%   by geometric distortion ensures M_true(<r_true,1) = M_obs(<r_obs,1).
%
%   For the Mezzi effect with ζ < 1, the expansion of outer shells yields 
%   R_j > 1 throughout the profile, resulting in η > 1. The true total mass 
%   exceeds the observed estimate because the observed aperture encompasses 
%   a smaller physical region than the true galactic extent.
%
%   NOTES AND CAVEATS
%   =================
%   * Surface density is assumed uniform within each shell (valid for fine binning).
%   * Negative shell masses (numerical noise) are clipped to zero.
%   * Non-monotonic M_ENC generates warning.

    %% 0. Basic Input Validation and Dimension Handling
    r_obs = r_obs(:);
    M_enc = M_enc(:);
    ScaleFactor = ScaleFactor(:);
    N = length(r_obs);
    
    %% Check for decreasing M_enc ( Non Monotonicity )
    dM_check = diff(M_enc);
    if any(dM_check < 0)
        warning('MassProfileFromMass:DecreasingMass', ...
            'M_enc has %d decreasing segment(s). Max decrease: %.4e', ...
            sum(dM_check < 0), max(abs(dM_check(dM_check < 0))));
    end

    %% 1. Convert enclosed mass to surface density in observed shells
    dM_obs_shell = diff(M_enc);  
    dA_obs_shell = pi * (r_obs(2:N).^2 - r_obs(1:N-1).^2);
    
    Sigma_shell = dM_obs_shell ./ dA_obs_shell;
    r_obs_mid = (r_obs(1:N-1) + r_obs(2:N)) / 2;
    
    %% 2. Coordinate transformation
    r_true = r_obs ./ ScaleFactor;
    
    %% 3. Calculate true shell masses
    dA_true_shell = pi * (r_true(2:N).^2 - r_true(1:N-1).^2);
    dM_true_shell = Sigma_shell .* dA_true_shell;
    dM_true_shell = max(dM_true_shell, 0);
    r_true_mid = (r_true(1:N-1) + r_true(2:N)) / 2;
    
    %% 4. EXACT CUMULATIVE RECONSTRUCTION 
    M_true_cum = zeros(N, 1);  
    M_true_cum(1) = M_enc(1);
    
    % Dimensions match explicitly
    cumsum_result = cumsum(dM_true_shell);
    M_true_cum(2:N) = M_enc(1) + cumsum_result;
    
    %% 5. Observed mass
    M_obs_cum = M_enc;
    
    %% 6. Calculate mass ratio
    M_total_true = M_true_cum(end);
    M_total_obs = M_enc(end);
    mass_ratio = M_total_true / M_total_obs;
end


 
   


function [galaxy_names, observed_radius, observed_speed, e_Vobs, ...
         observed_Vgas, e_Vgas] = read_rotation_data(filename)
    fileID = fopen(filename, 'r');
    if fileID == -1
        error('Cannot open file: %s', filename);
    end

    galaxy_names     = {};
    observed_radius  = {};
    observed_speed   = {};
    e_Vobs           = {};
    observed_Vgas    = {};
    e_Vgas           = {};

    while ~feof(fileID)
        line = fgetl(fileID);
        if isempty(line) || startsWith(strtrim(line), '#') || ...
           startsWith(strtrim(line), 'ID') || startsWith(strtrim(line), '=')
            continue;
        end

        parts = strsplit(strtrim(line));
        if numel(parts) < 6
            continue;
        end

        galaxy_name = parts{1};
        radius      = str2double(parts{2});
        Vobs        = str2double(parts{3});
        eVobs       = str2double(parts{4});
        Vgas        = str2double(parts{5});
        eVgas       = 0;

        galaxy_index = find(strcmp(galaxy_names, galaxy_name), 1);
        if isempty(galaxy_index)
            galaxy_names{end+1}     = galaxy_name;
            observed_radius{end+1}  = [];
            observed_speed{end+1}   = [];
            e_Vobs{end+1}           = [];
            observed_Vgas{end+1}    = [];
            e_Vgas{end+1}           = [];
            galaxy_index = numel(galaxy_names);
        end

        observed_radius{galaxy_index}(end+1) = radius;
        observed_speed{galaxy_index}(end+1)  = Vobs;
        e_Vobs{galaxy_index}(end+1)          = eVobs;
        observed_Vgas{galaxy_index}(end+1)   = Vgas;
        e_Vgas{galaxy_index}(end+1)          = eVgas;
    end
    fclose(fileID);
end

function [gs_galaxy_names, M2L, e_M2L, M2Lb, L36, e_L36, Lbulge, ...
         Rdisk, Reff, Distance, e_D, Mbar, e_Mbar] = read_mass_data(filename)
    % Extracts Mbar  ) and e_Mbar  TOTAL baryonic mass (stars + gas)
    fileID = fopen(filename, 'r');
    if fileID == -1
        error('Cannot open file: %s', filename);
    end

    gs_galaxy_names = {};
    M2L      = [];
    e_M2L    = [];
    M2Lb     = [];
    L36      = [];
    e_L36    = [];
    Lbulge   = [];
    Rdisk    = [];
    Reff     = [];
    Distance = [];
    e_D      = [];
    Mbar     = [];   
    e_Mbar   = [];   

    while ~feof(fileID)
        line = fgetl(fileID);
        if isempty(line) || startsWith(strtrim(line), '#') || ...
           startsWith(strtrim(line), 'ID') || startsWith(strtrim(line), '=')
            continue;
        end

        parts = strsplit(strtrim(line));
        if numel(parts) < 21
            continue;
        end

        gs_galaxy_names{end+1} = parts{1};
        M2L(end+1)    = str2double(parts{2});
        e_M2L(end+1)  = str2double(parts{3});
        M2Lb(end+1)   = str2double(parts{4});
        Distance(end+1) = str2double(parts{6});
        e_D(end+1)      = str2double(parts{7});
        L36(end+1)    = str2double(parts{10});
        e_L36(end+1)  = str2double(parts{11});
        Lbulge(end+1) = str2double(parts{12});
        Reff(end+1)   = str2double(parts{13});
        Rdisk(end+1)  = str2double(parts{17});
        Mbar(end+1)   = str2double(parts{20});   
        e_Mbar(end+1) = str2double(parts{21});  
    end
    fclose(fileID);
end

function [enc_disk, e_enc_disk] = compute_enclosed_disk_mass(r, M_disk, e_M_disk, R_disk)
    if R_disk > 0 && M_disk > 0
        x = r / R_disk;
        profile_factor = 1 - (1 + x) * exp(-x);
        enc_disk = M_disk * profile_factor;
        e_enc_disk = e_M_disk * profile_factor;
    else
        enc_disk = 0;
        e_enc_disk = 0;
    end
end

function [enc_bulge, e_enc_bulge] = compute_enclosed_bulge_mass(r, M_bulge, e_M_bulge, R_eff)
    if R_eff > 0 && M_bulge > 0
        a = R_eff / 1.8153;
        profile_factor = (r / (r + a))^2;
        enc_bulge = M_bulge * profile_factor;
        e_enc_bulge = e_M_bulge * profile_factor;
    else
        enc_bulge = 0;
        e_enc_bulge = 0;
    end
end
 