function HeleShaw_Cosmological_Analysis()
% =========================================================================
% HELE-SHAW COSMOLOGICAL ANALOG: ANALYTICAL SUITE
% Investigates macroscopic Silk damping via Taylor-Aris dispersion and
% Bianchi Type-I anisotropic expansion in viscous squeeze flow.
% =========================================================================

    % Run both analytical modules
    analyze_isotropic_damping();
    analyze_anisotropic_expansion();

end

% =========================================================================
% MODULE 1: 1D Radial Isotropic Expansion & Silk Damping
% =========================================================================
function analyze_isotropic_damping()
    fprintf('Initializing 1D Gap-Averaged Kinematic Simulation...\n');

    % Physical & Cosmological Parameters
    R_max = 0.05;           % Cell radius (m)
    h_gap = 0.001;          % Gap thickness (m)
    kappa = 1e-6;           % Molecular diffusivity (m^2/s)
    lambda_0 = 0.005;       % Initial primordial wavelength (m)
    k0 = 2 * pi / lambda_0; % Comoving wavenumber
    sigma = 0.015;          % Gaussian envelope standard deviation

    % Spatiotemporal Mesh
    t_span = linspace(0, 3, 150);
    r_mesh = linspace(0, R_max, 500);

    % Execute PDE Solver (Cylindrical Coordinates, m=1)
    sol = pdepe(1, @pdefun, @icfun, @bcfun, r_mesh, t_span);
    C = sol;

    % Data Extraction: Comoving Trajectory & Amplitude Decay
    a_t = sqrt(1 + t_span);                     % Radiation epoch scale factor
    r_track = lambda_0 .* a_t;                  % Comoving radial trajectory
    A0 = exp(-(lambda_0^2) / (2 * sigma^2));    % Initial peak amplitude

    A_sim = zeros(1, length(t_span));
    A_theory = zeros(1, length(t_span));

    for i = 1:length(t_span)
        A_sim(i) = interp1(r_mesh, C(i, :), r_track(i), 'spline');
        A_theory(i) = A0 * (a_t(i)^(-2 * kappa * k0^2)); % Theoretical 2D decay
    end

    % ---------------------------------------------------------------------
    % QUANTITATIVE COMPARISON: ideal Silk law vs. simulation
    % Reports the number quoted in the text for the agreement between
    % A_sim and the ideal law A = A0 * a^(-2 kappa k0^2).
    % ---------------------------------------------------------------------
    dev_abs = abs(A_sim / A0 - A_theory / A0);   % absolute, in units of A0
    [max_dev, i_max] = max(dev_abs);

    fprintf('\n--- MACROSCOPIC SILK DAMPING: DIAGNOSTICS ---\n');
    fprintf('k0                      : %.1f m^-1\n', k0);
    fprintf('kappa*k0^2              : %.4f s^-1\n', kappa * k0^2);
    fprintf('predicted exponent 2kk0^2: %.4f\n', 2 * kappa * k0^2);
    fprintf('scale factor reached a(T): %.4f\n', a_t(end));
    fprintf('A_sim/A0    at t = T    : %.5f\n', A_sim(end) / A0);
    fprintf('A_theory/A0 at t = T    : %.5f\n', A_theory(end) / A0);
    fprintf('\n>> MAX ABSOLUTE DEVIATION: %.4f  (= %.2f%% of A0)\n', ...
            max_dev, 100 * max_dev);
    fprintf('   occurring at t = %.3f s (a = %.3f)\n', ...
            t_span(i_max), a_t(i_max));
    fprintf('   [use this percentage in place of the estimate in Sect. 6.3]\n');

    % Regime diagnostics: this configuration tests the 2D reduction, and is
    % NOT the physically realizable design point (see design-window table).
    Pi   = 0.5 * h_gap^2 / kappa;                % Hbar h^2 / kappa at t = 0
    r_TA = sqrt(210) * kappa / (0.5 * h_gap);    % Taylor-Aris radius at t = 0
    fprintf('\n--- REGIME CHECK (t = 0) ---\n');
    fprintf('mixing parameter  Pi = Hbar h^2/kappa : %.3f\n', Pi);
    fprintf('Taylor-Aris radius r_TA               : %.1f mm\n', r_TA * 1000);
    fprintf('tracked crest spans r = %.1f to %.1f mm\n', ...
            r_track(1) * 1000, r_track(end) * 1000);

    % ---------------------------------------------------------------------
    % FIGURE 1: Passive Scalar Concentration Heatmap (Comoving Stretching)
    % ---------------------------------------------------------------------
    figure('Name', 'Fig 1: Heatmap', 'Color', 'w', 'Position', [50, 100, 600, 500]);
    [R, T] = meshgrid(r_mesh, t_span);
    pcolor(R, T, C);
    shading interp;
    colormap parula;
    colorbar;
    hold on;
    % Overlay the comoving trajectory of the peak
    plot(r_track, t_span, 'w--', 'LineWidth', 2);
    title('Passive Scalar Concentration: Comoving Stretching', 'Interpreter', 'latex');
    xlabel('Radial Distance, $r$ (m)', 'Interpreter', 'latex');
    ylabel('Time, $t$ (s)', 'Interpreter', 'latex');

    % ---------------------------------------------------------------------
    % FIGURE 2: Macroscopic Silk Damping (Amplitude Decay)
    % ---------------------------------------------------------------------
    figure('Name', 'Fig 2: Silk Damping', 'Color', 'w', 'Position', [670, 100, 600, 500]);
    plot(t_span, A_theory / A0, 'k--', 'LineWidth', 2); hold on;
    plot(t_span, A_sim / A0, 'r-', 'LineWidth', 2);
    title('Diffusive Damping of the Comoving Mode', 'Interpreter', 'latex');
    xlabel('Time, $t$ (s)', 'Interpreter', 'latex');
    ylabel('Normalized Amplitude, $A(t) / A_0$', 'Interpreter', 'latex');
    legend('Ideal Silk Damping ($A \propto a^{-2\kappa k_0^2}$)', ...
           'Gap-Averaged Simulation ($\kappa_{eff}$ Breakaway)', ...
           'Interpreter', 'latex', 'Location', 'southwest');
    grid on;

    % PDE Definitions (Nested)
    function [c, f, s] = pdefun(r, t, u, DuDr)
        c = 1;
        H_fluid = 1 / (2 * (1 + t)); % Hubble parameter analogous to a(t) = t^(1/2)
        v_r = H_fluid * r;

        % Taylor-Aris Dispersion Injection
        Pe = (v_r * h_gap) / kappa;
        kappa_eff = kappa * (1 + (Pe^2 / 210));

        f = kappa_eff * DuDr;
        s = -v_r * DuDr;
    end

    function u0 = icfun(r)
        u0 = exp(-(r.^2) / (2 * sigma^2)) .* cos(k0 * r);
    end

    function [pl, ql, pr, qr] = bcfun(~, ~, ~, ~, ~)
        pl = 0; ql = 1; % Neumann symmetry at origin
        pr = 0; qr = 1; % Zero flux gradient at boundary
    end
end

% =========================================================================
% MODULE 2: 3D Bianchi-I Anisotropic Kinematics
% =========================================================================
function analyze_anisotropic_expansion()
    fprintf('\nAnalyzing Bianchi Type-I Kinematics from CFD Data...\n');

    % Verify presence of CFD data files
    if ~isfile('X_axis.csv') || ~isfile('Y_axis.csv')
        fprintf('WARNING: X_axis.csv or Y_axis.csv not found. Skipping Bianchi-I analysis.\n');
        return;
    end

    % Load Data (Compatible with all MATLAB versions)
    data_X = readtable('X_axis.csv');
    data_Y = readtable('Y_axis.csv');

    x_coords = data_X.Points_0; v_x = data_X.U_0;
    y_coords = data_Y.Points_1; v_y = data_Y.U_1;

    if max(x_coords) > 1
        x_coords = x_coords / 1000; y_coords = y_coords / 1000;
    end

    % Spatial Cropping (Isolating bulk flow from edge effects, r < 0.85 R)
    crop_X = 0.05 * 0.85;
    crop_Y = 0.025 * 0.85;

    valid_X = ~isnan(x_coords) & ~isnan(v_x) & (x_coords <= crop_X);
    valid_Y = ~isnan(y_coords) & ~isnan(v_y) & (y_coords <= crop_Y);

    % Linear Regression for Hubble Parameters
    p_x = polyfit(x_coords(valid_X), v_x(valid_X), 1); H_x = p_x(1);
    p_y = polyfit(y_coords(valid_Y), v_y(valid_Y), 1); H_y = p_y(1);

    % Cosmological Metric Comparison
    theoretical_ratio = (0.025 / 0.05)^2;
    simulated_ratio = H_x / H_y;
    error_margin = abs(theoretical_ratio - simulated_ratio) / theoretical_ratio * 100;

    % ---------------------------------------------------------------------
    % FIGURE 3: Anisotropic Kinematics
    % ---------------------------------------------------------------------
    figure('Name', 'Fig 3: Bianchi-I Analog', 'Color', 'w', 'Position', [1290, 100, 600, 500]);
    plot(x_coords, v_x, 'b-', 'LineWidth', 2); hold on;
    plot(y_coords, v_y, 'r-', 'LineWidth', 2);
    title('Bianchi Type-I Analog: Anisotropic Expansion Metrics', 'Interpreter', 'latex');
    xlabel('Radial Distance from Origin (m)', 'Interpreter', 'latex');
    ylabel('Fluid Velocity (m/s)', 'Interpreter', 'latex');
    grid on;
    legend({sprintf('Major Axis X-Flow ($H_x = %.4f$ s$^{-1}$)', H_x), ...
            sprintf('Minor Axis Y-Flow ($H_y = %.4f$ s$^{-1}$)', H_y)}, ...
            'Interpreter', 'latex', 'Location', 'northwest');

    % Output findings
    fprintf('--- BIANCHI-I EXPANSION RESULTS ---\n');
    fprintf('H_x (Major Axis Expansion): %.5f s^-1\n', H_x);
    fprintf('H_y (Minor Axis Expansion): %.5f s^-1\n', H_y);
    fprintf('Theoretical Anisotropy Ratio (Hx/Hy): %.4f\n', theoretical_ratio);
    fprintf('Simulated Anisotropy Ratio (Hx/Hy): %.4f\n', simulated_ratio);
    fprintf('Extraction Error (Post-Cropping): %.2f%%\n', error_margin);
end