% ============================================================
%  RFA_komplett.m  -  GNU Octave conversion of Mathcad 6 document
%  Daniel's cosmological model (Neue Kosmologie)
%  Converted from LaTeX (Mathpix OCR of Mathcad 6 PDF)
%  Date: 2026-06-13 (refined)  |  Location: Neue_Kosmologie/RFA_komplett.m
% ============================================================
%
%  Physics summary:
%    The model treats the universe as a self-gravitating sphere of
%    baryonic matter (density rho_b). The dark-energy-like term is
%    encoded in the dimensionless ratio q = Omega_m / Omega_b.
%    The time-scale-factor relation uses an asinh form (cf. LambdaCDM
%    but with a model-specific x(a) instead of Omega_Lambda).
%
%  Run from Neue_Kosmologie/:
%    octave-cli.exe --no-gui --no-window-system --path "." --eval "RFA_komplett()"
% ============================================================

function RFA_komplett()

clear all
close all
format long

% Output directory: same folder as this script
out_dir = fileparts(which('RFA_komplett'));
if isempty(out_dir), out_dir = pwd; end

% Capture all console output to text file (for PDF printing)
diary(fullfile(out_dir, 'RFA_komplett_output.txt'));
diary on;

% ============================================================
%  SECTION 1: Fundamental Constants (Grundwerte)
% ============================================================

Mpc     = 3.085677581e22;       % [m]  1 Megaparsec
G       = 6.67430e-11;          % [m^3/(kg*s^2)]
c       = 299792458;            % [m/s]
H0_si   = 67.4e3 / Mpc;        % [1/s]  (67.4 km/s/Mpc)
H0_kms  = 67.4;                 % [km/s/Mpc]  for display

Omega_b   = 0.04960;
Omega_m0S = 0.31882;            % starting estimate
Omega_CDM = 0.315;              % ΛCDM matter density (for comparison)

rho_c = 3 * H0_si^2 / (8 * pi * G);   % [kg/m^3] critical density
rho_b = Omega_b * rho_c;               % [kg/m^3] baryon density

% Preliminary
y_prelim = (4/3) * pi * rho_b * c^2;
q_prelim = Omega_m0S / Omega_b;
x_prelim = (1 - q_prelim * Omega_b) / Omega_b;

fprintf('\n=== Section 1: Fundamental Constants ===\n');
fprintf('H0        = %.4f km/s/Mpc  = %.6e 1/s\n', H0_kms, H0_si);
fprintf('rho_c     = %.6e kg/m^3\n', rho_c);
fprintf('rho_b     = %.6e kg/m^3\n', rho_b);

% ============================================================
%  SECTION 2: Step 1 output  (Massenfaktor zur Gegenwart)
% ============================================================

q0       = 6.427606904079293;       % Schritt 1 result (2026-06-13)
Omega_m0 = q0 * Omega_b;           % = 0.318809302442333

fprintf('\n=== Section 2: Step 1 Output ===\n');
fprintf('q0        = %.15f\n', q0);
fprintf('Omega_m0  = %.15f\n', Omega_m0);
fprintf('  (target:  0.318809302442333)\n');

% ============================================================
%  SECTION 3 & 4: Sphere parameters and R0 self-consistency
% ============================================================
%
%  M = (4/3)*pi*R0^3 * rho_b   (constant baryonic mass)
%
%  Self-consistency condition:
%    R0 = c * sqrt(1 - 1/q0^2) * t(1)
%
%  We iterate: start from the document value R0 = 1.285e26 m,
%  compute t(1), update R0, repeat until converged.

R0 = 1.285e26;   % [m]  initial estimate from document

fprintf('\n=== Section 3 & 4: R0 self-consistency ===\n');

for k = 1:50
    M_k    = (4/3) * pi * R0^3 * rho_b;
    t1_k   = t_direct(1.0, q0, x_from_qR(q0, R0, M_k, G, c), H0_si, Omega_b, q0);
    R0_new = c * sqrt(1 - 1/q0^2) * t1_k;
    err    = abs(R0_new - R0) / R0;
    R0     = R0_new;
    if err < 1e-9
        fprintf('  R0 converged after %d iterations\n', k);
        break
    end
end

M     = (4/3) * pi * R0^3 * rho_b;
MG_c2 = M * G / c^2;

t1_val = t_of_a(1.0, q0, R0, M, G, c, H0_si, Omega_b);
Gyr    = 1e9 * 365.25 * 24 * 3600;   % [s]

fprintf('R0 (self-consistent)  = %.6e m  (target: 1.285e26 m)\n', R0);
fprintf('M                     = %.6e kg\n', M);
fprintf('t(a=1)                = %.9f Gyr  (target: 13.749542 Gyr)\n', t1_val / Gyr);

% ============================================================
%  SCHRITT 2: Blick in die Vergangenheit
%  Unvermeidliche Vorberechnungen + Bestimmung von t(a), x(a)
% ============================================================
%
%  t(a) := R_q(q0, x(a)) / v          [direct, non-iterative — PDF page 2]
%    q    = q0  (fixed)
%    R    = a * R0
%    x    = (2*q0-1) - 0.5*(1 - q0*(1 - M*q0*G/(R*c^2))^2)
%    R_q  = [2*c^2*M*q0^2*G + 2*sqrt(disc)] / [2*c^4*(5*q0-3-2*x)]
%    v    = c * sqrt(1 - 1/q0^2)
%    disc = M^2*q0^3*G^2*c^4*(-4*q0+3+2*x)

t0_s2 = t_a_Rq(1.0, q0, R0, M, G, c);

fprintf('\n=== Schritt 2: Unvermeidliche Vorberechnungen ===\n');
fprintf('t0  = %.3e s   (target: 4.339e17 s)\n', t0_s2);
fprintf('R0  = %.8e m   (target: 1.28496778e26 m)\n', R0);
fprintf('M   = %.3e kg  (target: 3.761e51 kg)\n', M);

% Bestimmung von t(a), x(a): a := 1, 0.9 .. 0.1
fprintf('\n=== Schritt 2: Bestimmung t(a), x(a) fuer a = 1, 0.9 .. 0.1 ===\n');
fprintf('%5s  %14s  %10s  %12s\n', 'a', 't [s]', 't [Gyr]', 'x(a)');
for a_i = 1.0 : -0.1 : 0.1
    t_i = t_a_Rq(a_i, q0, R0, M, G, c);
    x_i = x_from_qR(q0, a_i * R0, M, G, c);
    fprintf('%5.2f  %14.6e  %10.6f  %12.6f\n', a_i, t_i, t_i / Gyr, x_i);
end

% ============================================================
%  SECTION 5: Output functions  (a = 2 .. 0.06)
% ============================================================

a_vec = fliplr(0.06 : 0.005 : 2.0);
N     = length(a_vec);

R_vec        = zeros(1, N);
x_vec        = zeros(1, N);
OmLambda_vec = zeros(1, N);
OmM_vec      = zeros(1, N);

for i = 1:N
    a_i          = a_vec(i);
    R_vec(i)     = a_i * R0;
    x_vec(i)     = x_from_qR(q0, R_vec(i), M, G, c);
    OmLambda_vec(i) = x_vec(i) * Omega_b;
    OmM_vec(i)   = q0 * Omega_b;   % constant
end

% Spot checks
[~, idx1]  = min(abs(a_vec - 1.0));
[~, idx05] = min(abs(a_vec - 0.5));

x_001   = x_from_qR(q0, 0.01 * R0, M, G, c);
OmL_001 = x_001 * Omega_b;

fprintf('\n=== Section 5: Output checks ===\n');
fprintf('Omega_Lambda(a=1)    = %.10f  (target: 0.6811906976)\n', OmLambda_vec(idx1));
fprintf('Omega_Lambda(a=0.01) = %.3f  (target: 27.387)\n',  OmL_001);
fprintf('Singularity (doc.)   = 0.0740  (numerical observation; see Mathcad - MatLab_S.3.pdf)\n');
fprintf('R(a=1)               = %.4e m  (target: 1.285e26 m)\n', R_vec(idx1));

% ============================================================
%  SECTION 6: Hubble parameter functions
% ============================================================
%
%  H_g(a)   = c/R(a) * Mpc/1000                              [km/s/Mpc]
%  H1(a)    = (da/dt) / a  analytical: arg=3/2*H0*t_a_Rq*sqrt(x*Ob), da/dt from a(t) formula
%  H2(a)    = sqrt(8/3*pi*G*(q0*rho_b/a^3 + x(a)*rho_b))    (Friedmann)
%  H_CDM(a) = H0 * sqrt(Omega_CDM/a^3 + (1-Omega_CDM))      (standard ΛCDM)
%
%  Note: H1 != H2 because t_a_Rq (R_q/v) differs from the asinh formula.

Hg_vec   = zeros(1, N);
H1_vec   = zeros(1, N);
H2_vec   = zeros(1, N);
HCDM_vec = zeros(1, N);
t_vec    = zeros(1, N);

for i = 1:N
    a_i = a_vec(i);
    R_i = R_vec(i);
    x_i = x_vec(i);

    % Geometric Hubble
    Hg_vec(i) = c / R_i * Mpc / 1000;

    % ΛCDM comparison (standard Friedmann, flat, matter+Lambda)
    HCDM_vec(i) = H0_kms * sqrt(Omega_CDM / a_i^3 + (1 - Omega_CDM));

    if x_i > 0
        % H1: da/dt from analytical formula of a(t)=(sqrt(q0/x)*sinh(arg))^(2/3)
        %   evaluated at t = t_a_Rq(a)  [page 2 direct formula, not t_direct]
        %   da/dt = sqrt(q0/x)*cosh(arg)*H0*sqrt(x*Omega_b) / (sqrt(q0/x)*sinh(arg))^(1/3)
        %   with arg = 3/2*H0*t_a_Rq(a)*sqrt(x*Omega_b)
        t_i    = t_a_Rq(a_i, q0, R0, M, G, c);
        xOb    = x_i * Omega_b;
        arg    = 1.5 * H0_si * t_i * sqrt(xOb);
        sqQX   = sqrt(q0 / x_i);
        adot_i = sqQX * cosh(arg) * H0_si * sqrt(xOb) / (sqQX * sinh(arg))^(1/3);
        H1_vec(i) = adot_i / a_i * Mpc / 1000;
        t_vec(i)  = t_i;

        % H2: Friedmann form
        H2_vec(i) = sqrt(8/3 * pi * G * (q0*rho_b/a_i^3 + x_i*rho_b)) * Mpc / 1000;
    else
        H1_vec(i) = NaN;
        H2_vec(i) = NaN;
        t_vec(i)  = NaN;
    end
end

fprintf('\n=== Section 6: Hubble parameters ===\n');
fprintf('H_g  (a=0.5) = %.3f km/s/Mpc  (target: 143.982)\n', Hg_vec(idx05));
fprintf('H1   (a=0.5) = %.3f km/s/Mpc  (target: 118.344)\n', H1_vec(idx05));
fprintf('H2   (a=0.5) = %.3f km/s/Mpc  (target: 120.502)\n', H2_vec(idx05));
fprintf('HCDM (a=0.5) = %.3f km/s/Mpc  (target: 120.663)\n', HCDM_vec(idx05));
fprintf('H1   (a=1.0) = %.3f km/s/Mpc  (target:  67.400)\n', H1_vec(idx1));
fprintf('H2   (a=1.0) = %.3f km/s/Mpc  (target:  67.400)\n', H2_vec(idx1));
fprintf('HCDM (a=1.0) = %.3f km/s/Mpc  (target:  67.400)\n', HCDM_vec(idx1));
fprintf('H_g  (a=1.0) = %.3f km/s/Mpc  (target:  71.991)\n',  Hg_vec(idx1));

% ============================================================
%  SECTION 7: Tabulated output
% ============================================================

fprintf('\n=== Section 7: Table ===\n');
fprintf('%6s  %11s  %9s  %9s  %8s  %8s  %8s  %8s\n', ...
    'a', 'R [m]', 'x(a)', 'OmLambda', 'H_g', 'H1', 'H2', 'H_CDM');

for ap = [2.0, 1.5, 1.0, 0.75, 0.5, 0.25, 0.1, 0.06]
    [~, ii] = min(abs(a_vec - ap));
    fprintf('%6.3f  %11.4e  %9.5f  %9.5f  %8.3f  %8.3f  %8.3f  %8.3f\n', ...
        a_vec(ii), R_vec(ii), x_vec(ii), OmLambda_vec(ii), ...
        Hg_vec(ii), H1_vec(ii), H2_vec(ii), HCDM_vec(ii));
end

% ============================================================
%  SECTION 8: Plots
% ============================================================

% --- Figure 1: Hubble parameters ---
figure(1);
clf;
hold on;
plot(a_vec, Hg_vec,   'k-',  'LineWidth', 1.5, 'DisplayName', 'H_g (geometric)');
plot(a_vec, H1_vec,   'b--', 'LineWidth', 1.5, 'DisplayName', 'H_1 (da/dt)');
plot(a_vec, H2_vec,   'r-',  'LineWidth', 1.5, 'DisplayName', 'H_2 (Friedmann)');
plot(a_vec, HCDM_vec, 'm-.', 'LineWidth', 1.5, 'DisplayName', 'H_{CDM} (\Omega_{CDM}=0.315)');
plot(1.0, H0_kms, 'ko', 'MarkerSize', 8, 'MarkerFaceColor', 'k', ...
     'DisplayName', sprintf('H_0 = %.1f km/s/Mpc', H0_kms));
hold off;
grid on;
xlabel('Scale factor a');
ylabel('H [km/s/Mpc]');
title('Hubble parameter functions (Daniel''s cosmological model)');
legend('Location', 'NorthEast');
xlim([0.06, 2.0]);
ylim([0, 500]);

% --- Figure 2: Density parameters ---
figure(2);
clf;
hold on;
plot(a_vec, OmLambda_vec, 'r-',  'LineWidth', 1.5, 'DisplayName', '\Omega_\Lambda(a)');
plot(a_vec, OmM_vec,      'b--', 'LineWidth', 1.5, 'DisplayName', '\Omega_m (const)');
line(xlim, [0.68119 0.68119], 'Color', 'r', 'LineStyle', ':', 'LineWidth', 1.0);
line(xlim, [Omega_m0 Omega_m0], 'Color', 'b', 'LineStyle', ':', 'LineWidth', 1.0);
hold off;
grid on;
xlabel('Scale factor a');
ylabel('\Omega');
title('Density parameters vs scale factor');
legend('Location', 'NorthEast');
xlim([0.06, 2.0]);

% --- Figure 3: Age of universe t(a) ---
figure(3);
clf;
hold on;
plot(a_vec, t_vec / Gyr, 'b-', 'LineWidth', 1.5, 'DisplayName', 't(a)');
plot(1.0, t1_val / Gyr, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r', ...
     'DisplayName', sprintf('t(1) = %.3f Gyr', t1_val / Gyr));
hold off;
grid on;
xlabel('Scale factor a');
ylabel('t [Gyr]');
title('Age of universe t(a)');
legend('Location', 'NorthWest');
xlim([0.06, 2.0]);

% ============================================================
%  Save figures as PDF
% ============================================================

print(1, fullfile(out_dir, 'RFA_fig1_Hubble.pdf'),  '-dpdf', '-r300');
print(2, fullfile(out_dir, 'RFA_fig2_Omega.pdf'),   '-dpdf', '-r300');
print(3, fullfile(out_dir, 'RFA_fig3_Age.pdf'),     '-dpdf', '-r300');
fprintf('\nFigures saved:\n');
fprintf('  %s\n', fullfile(out_dir, 'RFA_fig1_Hubble.pdf'));
fprintf('  %s\n', fullfile(out_dir, 'RFA_fig2_Omega.pdf'));
fprintf('  %s\n', fullfile(out_dir, 'RFA_fig3_Age.pdf'));

% ============================================================
%  SECTION 9: Summary
% ============================================================

fprintf('\n=== Summary ===\n');
fprintf('q0              = %.15f\n',   q0);
fprintf('Omega_m0        = %.15f\n',   Omega_m0);
fprintf('R0              = %.8e m\n',  R0);
fprintf('M               = %.6e kg\n', M);
fprintf('t(a=1)          = %.6f Gyr  (target: 13.749542)\n', t1_val / Gyr);
fprintf('H1(a=1)         = %.3f km/s/Mpc  (target: 67.4)\n',   H1_vec(idx1));
fprintf('H2(a=1)         = %.3f km/s/Mpc  (target: 67.4)\n',   H2_vec(idx1));
fprintf('HCDM(a=0.5)     = %.3f km/s/Mpc  (target: 120.663)\n', HCDM_vec(idx05));
fprintf('Omega_Lambda(1) = %.10f  (target: 0.6811906976)\n',  OmLambda_vec(idx1));
fprintf('Singularity     = 0.0740  (numerical; no closed-form formula in document)\n');
fprintf('Done.\n');
fprintf('Console output saved: %s\n', fullfile(out_dir, 'RFA_komplett_output.txt'));
diary off;


end  % function RFA_komplett

% ============================================================
%  LOCAL FUNCTIONS
% ============================================================

% --- x(q, R) from the Grundgleichung ---
function val = x_from_qR(q, R, M, G, c)
    val = (2*q - 1) - 0.5 * (1 - q * (1 - M*q*G / (R*c^2))^2);
end

% --- R_q from the quadratic root of the energy equation ---
%  R_q = [2*c^2*M*q^2*G + 2*sqrt(disc)] / [2*(5*c^4*q - 3*c^4 - 2*c^4*x)]
%  Verified against Mathcad - MatLab_S.2.pdf (2026-06-13)
function Rq = R_q_formula(q, x, M, G, c)
    denom = 2 * c^4 * (5*q - 3 - 2*x);
    disc  = M^2 * q^3 * G^2 * c^4 * (-4*q + 3 + 2*x);
    if disc < 0 || abs(denom) < 1e-300
        Rq = NaN;
        return
    end
    Rq = (2*c^2*M*q^2*G + 2*sqrt(disc)) / denom;
end

% --- Direct t(a) from the asinh formula ---
%   t = (2 / (3*H0*sqrt(x*Omega_b))) * asinh(a^(3/2) * sqrt(x/q0))
function t = t_direct(a, q, x, H0_si, Omega_b, q0)
    xOb = x * Omega_b;
    if xOb <= 0 || x <= 0
        t = NaN;
        return
    end
    t = (2 / (3 * H0_si * sqrt(xOb))) * asinh(a^(3/2) * sqrt(x / q0));
end

% --- Full t(a) with iterative q correction (deep-past regime) ---
%  Transcription of the while-loop algorithm from Mathcad - MatLab_S.2.pdf
function t_val = t_of_a(a, q0, R0, M, G, c, H0_si, Omega_b)
    R = a * R0;
    q = q0;
    x = x_from_qR(q, R, M, G, c);

    t_val = t_direct(a, q, x, H0_si, Omega_b, q0);
    if isnan(t_val)
        return
    end

    R_int = c * sqrt(max(1 - 1/q^2, 0)) * t_val;
    R_q   = R_q_formula(q, x, M, G, c);

    if isnan(R_q) || R_q < R_int
        return   % normal regime — no iteration needed
    end

    % Outer while R_q >= R_int
    for outer = 1:500
        if isnan(R_q) || R_q < R_int, break; end

        q = q - 1e-6;
        if q <= 1.001, q = 1.001; break; end

        n = 0;
        while n <= 4

            v = sqrt(max(1 - 1/q^2, 0));

            % (a) Newton-Raphson for t
            dt   = max(abs(t_val) * 1e-7, 1.0);

            Ri_a = c * v * t_val;
            xa   = x_from_qR(q, Ri_a, M, G, c);
            Rqa  = R_q_formula(q, xa, M, G, c);
            if isnan(Rqa), break; end
            fa   = Rqa - Ri_a;

            Ri_b = c * v * (t_val + dt);
            xb   = x_from_qR(q, Ri_b, M, G, c);
            Rqb  = R_q_formula(q, xb, M, G, c);
            if isnan(Rqb), break; end
            fb   = Rqb - Ri_b;

            df_dt = (fb - fa) / dt;
            if abs(df_dt) > 1e-300
                t_val = t_val - fa / df_dt;
            end

            % (b) R_int <- c * sqrt(1-1/q^2) * t
            v     = sqrt(max(1 - 1/q^2, 0));
            R_int = c * v * t_val;

            % (c) x <- energy eq evaluated at R_int
            x = x_from_qR(q, R_int, M, G, c);

            % (d) Newton-Raphson for q
            dq  = max(abs(q) * 1e-7, 1e-10);

            v1  = sqrt(max(1 - 1/q^2, 0));
            x1  = x_from_qR(q, R_int, M, G, c);
            Rq1 = R_q_formula(q, x1, M, G, c);
            if isnan(Rq1), break; end
            fq1 = Rq1 - c * v1 * t_val;

            q2  = q + dq;
            v2  = sqrt(max(1 - 1/q2^2, 0));
            x2  = x_from_qR(q2, R_int, M, G, c);
            Rq2 = R_q_formula(q2, x2, M, G, c);
            if isnan(Rq2), break; end
            fq2 = Rq2 - c * v2 * t_val;

            df_dq = (fq2 - fq1) / dq;
            if abs(df_dq) > 1e-300
                q = q - fq1 / df_dq;
            end
            if q <= 1.001, q = 1.001; break; end

            % (e) R_q <- updated with new q and x
            x   = x_from_qR(q, R_int, M, G, c);
            R_q = R_q_formula(q, x, M, G, c);
            if isnan(R_q), break; end

            n = n + 1;
        end

        % Refresh for outer while condition check
        v     = sqrt(max(1 - 1/q^2, 0));
        R_int = c * v * t_val;
        x     = x_from_qR(q, R_int, M, G, c);
        R_q   = R_q_formula(q, x, M, G, c);
    end
end

% --- t(a) = R_q(q0, x) / v  [Schritt 2, page 2 — direct, non-iterative] ---
function t_val = t_a_Rq(a, q0, R0, M, G, c)
    R     = a * R0;
    x     = x_from_qR(q0, R, M, G, c);
    Rq    = R_q_formula(q0, x, M, G, c);
    if isnan(Rq)
        t_val = NaN;
        return;
    end
    v     = c * sqrt(1 - 1/q0^2);
    t_val = Rq / v;
end

% --- da/dt: time derivative of scale factor ---
%   From a(t) = (sqrt(q0/x) * sinh(3/2*H0*t*sqrt(x*Omega_b)))^(2/3)
function adot = da_dt(a, x_a, t_a, q0, H0_si, Omega_b)
    xOb  = x_a * Omega_b;
    arg  = 1.5 * H0_si * t_a * sqrt(xOb);
    adot = (2/3) * a^(-0.5) * sqrt(q0 / x_a) ...
           * cosh(arg) * (1.5 * H0_si * sqrt(xOb));
end