{ "cells": [ { "cell_type": "markdown", "id": "f4f422e3", "metadata": {}, "source": [ "RFA_komplett.py\n", "\n", "**Python translation of RFA_komplett.m (GNU Octave / Mathcad 6)**\n", "\n", "**Daniel's cosmological model (Neue Kosmologie)**\n", "\n", "Physics summary:\n", " The model treats the universe as a self-gravitating sphere of baryonic matter\n", " (density rho_b). The dark-energy-like term is encoded in the dimensionless\n", " ratio q = Omega_m / Omega_b. The time-scale-factor relation uses an asinh\n", " form (cf. LambdaCDM but with a model-specific x(a) instead of Omega_Lambda)." ] }, { "cell_type": "code", "execution_count": null, "id": "3a7c42dd", "metadata": {}, "outputs": [], "source": [ "import math\n", "import numpy as np\n", "import matplotlib\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "1b4be949", "metadata": {}, "source": [ "\n", "**Lokale Hilfsfunktionen**\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1a119872", "metadata": {}, "outputs": [], "source": [ "def x_from_qR(q, R, M, G, c):\n", " \"\"\"x(q, R) from the Grundgleichung.\"\"\"\n", " return (2*q - 1) - 0.5 * (1 - q * (1 - M*q*G / (R*c**2))**2)\n", "\n", "def R_q_formula(q, x, M, G, c):\n", " \"\"\"\n", " R_q from the quadratic root of the energy equation.\n", " R_q = [2*c^2*M*q^2*G + 2*sqrt(disc)] / [2*(5*c^4*q - 3*c^4 - 2*c^4*x)]\n", " Returns NaN when the discriminant is negative or denominator is zero.\n", " \"\"\"\n", " denom = 2 * c**4 * (5*q - 3 - 2*x)\n", " disc = M**2 * q**3 * G**2 * c**4 * (-4*q + 3 + 2*x)\n", " if disc < 0 or abs(denom) < 1e-300:\n", " return float(\"nan\")\n", " return (2*c**2*M*q**2*G + 2*math.sqrt(disc)) / denom\n", "\n", "def t_direct(a, q, x, H0_si, Omega_b, q0):\n", " \"\"\"\n", " Direct t(a) from the asinh formula:\n", " t = 2/(3*H0*sqrt(x*Omega_b)) * asinh(a^(3/2) * sqrt(x/q0))\n", " \"\"\"\n", " xOb = x * Omega_b\n", " if xOb <= 0 or x <= 0:\n", " return float(\"nan\")\n", " return (2.0 / (3.0 * H0_si * math.sqrt(xOb))) * math.asinh(a**1.5 * math.sqrt(x / q0))\n", "\n", "def t_of_a(a, q0, R0, M, G, c, H0_si, Omega_b):\n", " \"\"\"\n", " Full t(a) with iterative q correction (deep-past regime).\n", " Transcription of the while-loop algorithm from Mathcad - MatLab_S.2.pdf.\n", " \"\"\"\n", " R = a * R0\n", " q = q0\n", " x = x_from_qR(q, R, M, G, c)\n", "\n", " t_val = t_direct(a, q, x, H0_si, Omega_b, q0)\n", " if math.isnan(t_val):\n", " return float(\"nan\")\n", "\n", " R_int = c * math.sqrt(max(1 - 1/q**2, 0)) * t_val\n", " R_q = R_q_formula(q, x, M, G, c)\n", "\n", " if math.isnan(R_q) or R_q < R_int:\n", " return t_val # normal regime – no iteration needed\n", "\n", " # Outer loop: while R_q >= R_int\n", " for _ in range(500):\n", " if math.isnan(R_q) or R_q < R_int:\n", " break\n", "\n", " q -= 1e-6\n", " if q <= 1.001:\n", " q = 1.001\n", " break\n", "\n", " n = 0\n", " while n <= 4:\n", " v = math.sqrt(max(1 - 1/q**2, 0))\n", "\n", " # (a) Newton-Raphson for t\n", " dt = max(abs(t_val) * 1e-7, 1.0)\n", "\n", " Ri_a = c * v * t_val\n", " xa = x_from_qR(q, Ri_a, M, G, c)\n", " Rqa = R_q_formula(q, xa, M, G, c)\n", " if math.isnan(Rqa):\n", " break\n", " fa = Rqa - Ri_a\n", "\n", " Ri_b = c * v * (t_val + dt)\n", " xb = x_from_qR(q, Ri_b, M, G, c)\n", " Rqb = R_q_formula(q, xb, M, G, c)\n", " if math.isnan(Rqb):\n", " break\n", " fb = Rqb - Ri_b\n", "\n", " df_dt = (fb - fa) / dt\n", " if abs(df_dt) > 1e-300:\n", " t_val -= fa / df_dt\n", "\n", " # (b) R_int <- c * sqrt(1-1/q^2) * t\n", " v = math.sqrt(max(1 - 1/q**2, 0))\n", " R_int = c * v * t_val\n", "\n", " # (c) x <- energy equation evaluated at R_int\n", " x = x_from_qR(q, R_int, M, G, c)\n", "\n", " # (d) Newton-Raphson for q\n", " dq = max(abs(q) * 1e-7, 1e-10)\n", "\n", " v1 = math.sqrt(max(1 - 1/q**2, 0))\n", " x1 = x_from_qR(q, R_int, M, G, c)\n", " Rq1 = R_q_formula(q, x1, M, G, c)\n", " if math.isnan(Rq1):\n", " break\n", " fq1 = Rq1 - c * v1 * t_val\n", "\n", " q2 = q + dq\n", " v2 = math.sqrt(max(1 - 1/q2**2, 0))\n", " x2 = x_from_qR(q2, R_int, M, G, c)\n", " Rq2 = R_q_formula(q2, x2, M, G, c)\n", " if math.isnan(Rq2):\n", " break\n", " fq2 = Rq2 - c * v2 * t_val\n", "\n", " df_dq = (fq2 - fq1) / dq\n", " if abs(df_dq) > 1e-300:\n", " q -= fq1 / df_dq\n", " if q <= 1.001:\n", " q = 1.001\n", " break\n", "\n", " # (e) R_q <- updated with new q and x\n", " x = x_from_qR(q, R_int, M, G, c)\n", " R_q = R_q_formula(q, x, M, G, c)\n", " if math.isnan(R_q):\n", " break\n", "\n", " n += 1\n", "\n", " # Refresh for outer loop condition\n", " v = math.sqrt(max(1 - 1/q**2, 0))\n", " R_int = c * v * t_val\n", " x = x_from_qR(q, R_int, M, G, c)\n", " R_q = R_q_formula(q, x, M, G, c)\n", "\n", " return t_val\n", "\n", "\n", "def t_a_Rq(a, q0, R0, M, G, c):\n", " \"\"\"\n", " t(a) = R_q(q0, x) / v [Schritt 2, page 2 – direct, non-iterative]\n", " \"\"\"\n", " R = a * R0\n", " x = x_from_qR(q0, R, M, G, c)\n", " Rq = R_q_formula(q0, x, M, G, c)\n", " if math.isnan(Rq):\n", " return float(\"nan\")\n", " v = c * math.sqrt(1 - 1/q0**2)\n", " return Rq / v\n", "\n", "\n", "def da_dt(a, x_a, t_a, q0, H0_si, Omega_b):\n", " \"\"\"\n", " da/dt: time derivative of the scale factor.\n", " From a(t) = (sqrt(q0/x) * sinh(3/2 * H0 * t * sqrt(x*Omega_b)))^(2/3)\n", " \"\"\"\n", " xOb = x_a * Omega_b\n", " arg = 1.5 * H0_si * t_a * math.sqrt(xOb)\n", " adot = ((2/3) * a**(-0.5) * math.sqrt(q0 / x_a)\n", " * math.cosh(arg) * (1.5 * H0_si * math.sqrt(xOb)))\n", " return adot" ] }, { "cell_type": "markdown", "id": "f8a97e85", "metadata": {}, "source": [ "**Hauptprogramm**" ] }, { "cell_type": "code", "execution_count": null, "id": "6dd9fb45", "metadata": {}, "outputs": [], "source": [ "def main():\n", " import os\n", " # out_dir = os.path.dirname(os.path.abspath(__file__))\n", " out_dir = \".\"\n", "\n", " # ── SECTION 1: Fundamental Constants ─────────────────────────────────────\n", " Mpc = 3.085677581e22 # [m] 1 Megaparsec\n", " G = 6.67430e-11 # [m^3/(kg*s^2)]\n", " c = 299792458.0 # [m/s]\n", " H0_si = 67.4e3 / Mpc # [1/s]\n", " H0_kms = 67.4 # [km/s/Mpc] for display\n", "\n", " Omega_b = 0.04960\n", " Omega_m0S = 0.31882 # starting estimate\n", " Omega_CDM = 0.315 # ΛCDM matter density (comparison)\n", "\n", " rho_c = 3 * H0_si**2 / (8 * math.pi * G) # [kg/m^3] critical density\n", " rho_b = Omega_b * rho_c # [kg/m^3] baryon density\n", "\n", " # Preliminary\n", " q_prelim = Omega_m0S / Omega_b\n", " x_prelim = (1 - q_prelim * Omega_b) / Omega_b\n", "\n", " print(\"\\n=== Section 1: Fundamental Constants ===\")\n", " print(f\"H0 = {H0_kms:.4f} km/s/Mpc = {H0_si:.6e} 1/s\")\n", " print(f\"rho_c = {rho_c:.6e} kg/m^3\")\n", " print(f\"rho_b = {rho_b:.6e} kg/m^3\")\n", "\n", " # ── SECTION 2: Step 1 output ──────────────────────────────────────────────\n", " q0 = 6.427606904079293\n", " Omega_m0 = q0 * Omega_b\n", "\n", " print(\"\\n=== Section 2: Step 1 Output ===\")\n", " print(f\"q0 = {q0:.15f}\")\n", " print(f\"Omega_m0 = {Omega_m0:.15f}\")\n", " print( \" (target: 0.318809302442333)\")\n", "\n", " # ── SECTION 3 & 4: R0 self-consistency ───────────────────────────────────\n", " R0 = 1.285e26 # [m] initial estimate\n", "\n", " print(\"\\n=== Section 3 & 4: R0 self-consistency ===\")\n", "\n", " for k in range(1, 51):\n", " M_k = (4/3) * math.pi * R0**3 * rho_b\n", " x_k = x_from_qR(q0, R0, M_k, G, c)\n", " t1_k = t_direct(1.0, q0, x_k, H0_si, Omega_b, q0)\n", " R0_new = c * math.sqrt(1 - 1/q0**2) * t1_k\n", " err = abs(R0_new - R0) / R0\n", " R0 = R0_new\n", " if err < 1e-9:\n", " print(f\" R0 converged after {k} iterations\")\n", " break\n", "\n", " M = (4/3) * math.pi * R0**3 * rho_b\n", " MG_c2 = M * G / c**2\n", " Gyr = 1e9 * 365.25 * 24 * 3600 # [s]\n", " t1_val = t_of_a(1.0, q0, R0, M, G, c, H0_si, Omega_b)\n", "\n", " print(f\"R0 (self-consistent) = {R0:.6e} m (target: 1.285e26 m)\")\n", " print(f\"M = {M:.6e} kg\")\n", " print(f\"t(a=1) = {t1_val/Gyr:.9f} Gyr (target: 13.749542 Gyr)\")\n", "\n", " # ── SCHRITT 2 ─────────────────────────────────────────────────────────────\n", " t0_s2 = t_a_Rq(1.0, q0, R0, M, G, c)\n", "\n", " print(\"\\n=== Schritt 2: Unvermeidliche Vorberechnungen ===\")\n", " print(f\"t0 = {t0_s2:.3e} s (target: 4.339e17 s)\")\n", " print(f\"R0 = {R0:.8e} m (target: 1.28496778e26 m)\")\n", " print(f\"M = {M:.3e} kg (target: 3.761e51 kg)\")\n", "\n", " print(\"\\n=== Schritt 2: Bestimmung t(a), x(a) fuer a = 1, 0.9 .. 0.1 ===\")\n", " print(f\"{'a':>5} {'t [s]':>14} {'t [Gyr]':>10} {'x(a)':>12}\")\n", " a_range = [round(1.0 - 0.1*k, 10) for k in range(10)]\n", " for a_i in a_range:\n", " t_i = t_a_Rq(a_i, q0, R0, M, G, c)\n", " x_i = x_from_qR(q0, a_i * R0, M, G, c)\n", " print(f\"{a_i:5.2f} {t_i:14.6e} {t_i/Gyr:10.6f} {x_i:12.6f}\")\n", "\n", " # ── SECTION 5: Output functions (a = 2 .. 0.06) ─────────────────────────\n", " a_vec = np.flip(np.arange(0.06, 2.005, 0.005))\n", " N = len(a_vec)\n", "\n", " R_vec = np.zeros(N)\n", " x_vec = np.zeros(N)\n", " OmLambda_vec = np.zeros(N)\n", " OmM_vec = np.zeros(N)\n", "\n", " for i, a_i in enumerate(a_vec):\n", " R_vec[i] = a_i * R0\n", " x_vec[i] = x_from_qR(q0, R_vec[i], M, G, c)\n", " OmLambda_vec[i] = x_vec[i] * Omega_b\n", " OmM_vec[i] = q0 * Omega_b # constant\n", "\n", " idx1 = int(np.argmin(np.abs(a_vec - 1.0)))\n", " idx05 = int(np.argmin(np.abs(a_vec - 0.5)))\n", "\n", " x_001 = x_from_qR(q0, 0.01 * R0, M, G, c)\n", " OmL_001 = x_001 * Omega_b\n", "\n", " print(\"\\n=== Section 5: Output checks ===\")\n", " print(f\"Omega_Lambda(a=1) = {OmLambda_vec[idx1]:.10f} (target: 0.6811906976)\")\n", " print(f\"Omega_Lambda(a=0.01) = {OmL_001:.3f} (target: 27.387)\")\n", " print( \"Singularity (doc.) = 0.0740 (numerical observation)\")\n", " print(f\"R(a=1) = {R_vec[idx1]:.4e} m (target: 1.285e26 m)\")\n", "\n", " # ── SECTION 6: Hubble parameter functions ─────────────────────────────────\n", " Hg_vec = np.zeros(N)\n", " H1_vec = np.full(N, float(\"nan\"))\n", " H2_vec = np.full(N, float(\"nan\"))\n", " HCDM_vec = np.zeros(N)\n", " t_vec = np.full(N, float(\"nan\"))\n", "\n", " for i, a_i in enumerate(a_vec):\n", " R_i = R_vec[i]\n", " x_i = x_vec[i]\n", "\n", " # Geometric Hubble\n", " Hg_vec[i] = c / R_i * Mpc / 1000.0\n", "\n", " # ΛCDM comparison\n", " HCDM_vec[i] = H0_kms * math.sqrt(Omega_CDM / a_i**3 + (1 - Omega_CDM))\n", "\n", " if x_i > 0:\n", " t_i = t_a_Rq(a_i, q0, R0, M, G, c)\n", " xOb = x_i * Omega_b\n", " arg = max (1e-4, 1.5 * H0_si * t_i * math.sqrt(xOb))\n", " sqQX = math.sqrt(q0 / x_i)\n", " adot_i = sqQX * math.cosh(arg) * H0_si * math.sqrt(xOb) / (sqQX * math.sinh(arg))**(1/3)\n", " H1_vec[i] = adot_i / a_i * Mpc / 1000.0\n", " t_vec[i] = t_i\n", "\n", " H2_vec[i] = (math.sqrt(8/3 * math.pi * G\n", " * (q0*rho_b/a_i**3 + x_i*rho_b))\n", " * Mpc / 1000.0)\n", "\n", " print(\"\\n=== Section 6: Hubble parameters ===\")\n", " print(f\"H_g (a=0.5) = {Hg_vec[idx05]:.3f} km/s/Mpc (target: 143.982)\")\n", " print(f\"H1 (a=0.5) = {H1_vec[idx05]:.3f} km/s/Mpc (target: 118.344)\")\n", " print(f\"H2 (a=0.5) = {H2_vec[idx05]:.3f} km/s/Mpc (target: 120.502)\")\n", " print(f\"HCDM (a=0.5) = {HCDM_vec[idx05]:.3f} km/s/Mpc (target: 120.663)\")\n", " print(f\"H1 (a=1.0) = {H1_vec[idx1]:.3f} km/s/Mpc (target: 67.400)\")\n", " print(f\"H2 (a=1.0) = {H2_vec[idx1]:.3f} km/s/Mpc (target: 67.400)\")\n", " print(f\"HCDM (a=1.0) = {HCDM_vec[idx1]:.3f} km/s/Mpc (target: 67.400)\")\n", " print(f\"H_g (a=1.0) = {Hg_vec[idx1]:.3f} km/s/Mpc (target: 71.991)\")\n", "\n", " # ── SECTION 7: Tabulated output ───────────────────────────────────────────\n", " print(\"\\n=== Section 7: Table ===\")\n", " print(f\"{'a':>6} {'R [m]':>11} {'x(a)':>9} {'OmLambda':>9} \"\n", " f\"{'H_g':>8} {'H1':>8} {'H2':>8} {'H_CDM':>8}\")\n", "\n", " for ap in [2.0, 1.5, 1.0, 0.75, 0.5, 0.25, 0.1, 0.06]:\n", " ii = int(np.argmin(np.abs(a_vec - ap)))\n", " print(f\"{a_vec[ii]:6.3f} {R_vec[ii]:11.4e} {x_vec[ii]:9.5f} \"\n", " f\"{OmLambda_vec[ii]:9.5f} {Hg_vec[ii]:8.3f} \"\n", " f\"{H1_vec[ii]:8.3f} {H2_vec[ii]:8.3f} {HCDM_vec[ii]:8.3f}\")\n", "\n", " # ── SECTION 8: Plots ──────────────────────────────────────────────────────\n", "\n", " # Figure 1: Hubble parameters\n", " fig1, ax1 = plt.subplots(figsize=(9, 5))\n", " ax1.plot(a_vec, Hg_vec, \"k-\", lw=1.5, label=r\"$H_g$ (geometric)\")\n", " ax1.plot(a_vec, H1_vec, \"b--\", lw=1.5, label=r\"$H_1$ (da/dt)\")\n", " ax1.plot(a_vec, H2_vec, \"r-\", lw=1.5, label=r\"$H_2$ (Friedmann)\")\n", " ax1.plot(a_vec, HCDM_vec, \"m-.\", lw=1.5, label=r\"$H_\\mathrm{CDM}$ ($\\Omega_\\mathrm{CDM}=0.315$)\")\n", " ax1.plot(1.0, H0_kms, \"ko\", ms=8, label=f\"$H_0 = {H0_kms:.1f}$ km/s/Mpc\")\n", " ax1.set_xlabel(\"Scale factor a\")\n", " ax1.set_ylabel(\"H [km/s/Mpc]\")\n", " ax1.set_title(\"Hubble parameter functions (Daniel's cosmological model)\")\n", " ax1.set_xlim(0.06, 2.0)\n", " ax1.set_ylim(0, 500)\n", " ax1.grid(True)\n", " ax1.legend(loc=\"upper right\")\n", " fig1.tight_layout()\n", " fig1.savefig(os.path.join(out_dir, \"RFA_fig1_Hubble.pdf\"), dpi=300)\n", " print(f\"\\nFigure saved: {os.path.join(out_dir, 'RFA_fig1_Hubble.pdf')}\")\n", " plt.show ()\n", "\n", " # Figure 2: Density parameters\n", " fig2, ax2 = plt.subplots(figsize=(9, 5))\n", " ax2.plot(a_vec, OmLambda_vec, \"r-\", lw=1.5, label=r\"$\\Omega_\\Lambda(a)$\")\n", " ax2.plot(a_vec, OmM_vec, \"b--\", lw=1.5, label=r\"$\\Omega_m$ (const)\")\n", " ax2.axhline(0.68119, color=\"r\", ls=\":\", lw=1.0)\n", " ax2.axhline(Omega_m0, color=\"b\", ls=\":\", lw=1.0)\n", " ax2.set_xlabel(\"Scale factor a\")\n", " ax2.set_ylabel(r\"$\\Omega$\")\n", " ax2.set_title(\"Density parameters vs scale factor\")\n", " ax2.set_xlim(0.06, 2.0)\n", " ax2.grid(True)\n", " ax2.legend(loc=\"upper right\")\n", " fig2.tight_layout()\n", " fig2.savefig(os.path.join(out_dir, \"RFA_fig2_Omega.pdf\"), dpi=300)\n", " print(f\"Figure saved: {os.path.join(out_dir, 'RFA_fig2_Omega.pdf')}\")\n", " plt.show ()\n", "\n", " # Figure 3: Age of universe\n", " fig3, ax3 = plt.subplots(figsize=(9, 5))\n", " ax3.plot(a_vec, t_vec / Gyr, \"b-\", lw=1.5, label=\"t(a)\")\n", " ax3.plot(1.0, t1_val / Gyr, \"ro\", ms=8,\n", " label=f\"t(1) = {t1_val/Gyr:.3f} Gyr\")\n", " ax3.set_xlabel(\"Scale factor a\")\n", " ax3.set_ylabel(\"t [Gyr]\")\n", " ax3.set_title(\"Age of universe t(a)\")\n", " ax3.set_xlim(0.06, 2.0)\n", " ax3.grid(True)\n", " ax3.legend(loc=\"upper left\")\n", " fig3.tight_layout()\n", " fig3.savefig(os.path.join(out_dir, \"RFA_fig3_Age.pdf\"), dpi=300)\n", " print(f\"Figure saved: {os.path.join(out_dir, 'RFA_fig3_Age.pdf')}\")\n", " plt.show ()\n", "\n", " # ── SECTION 9: Summary ────────────────────────────────────────────────────\n", " print(\"\\n=== Summary ===\")\n", " print(f\"q0 = {q0:.15f}\")\n", " print(f\"Omega_m0 = {Omega_m0:.15f}\")\n", " print(f\"R0 = {R0:.8e} m\")\n", " print(f\"M = {M:.6e} kg\")\n", " print(f\"t(a=1) = {t1_val/Gyr:.6f} Gyr (target: 13.749542)\")\n", " print(f\"H1(a=1) = {H1_vec[idx1]:.3f} km/s/Mpc (target: 67.4)\")\n", " print(f\"H2(a=1) = {H2_vec[idx1]:.3f} km/s/Mpc (target: 67.4)\")\n", " print(f\"HCDM(a=0.5) = {HCDM_vec[idx05]:.3f} km/s/Mpc (target: 120.663)\")\n", " print(f\"Omega_Lambda(1) = {OmLambda_vec[idx1]:.10f} (target: 0.6811906976)\")\n", " print( \"Singularity = 0.0740 (numerical; no closed-form formula in document)\")\n", " print(\"Done.\")" ] }, { "cell_type": "code", "execution_count": null, "id": "1942223d", "metadata": {}, "outputs": [], "source": [ "main ()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }