{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Schubart Master Formula P(t) — reproducibility notebook\n",
    "\n",
    "This notebook regenerates the two figures on the reproducibility page of\n",
    "https://www.pt-schubart.com/reproducibility and the illustrative Monte-Carlo\n",
    "estimate of the aggregated output.\n",
    "\n",
    "**Important:** the parameter magnitudes below are *illustrative* order-of-\n",
    "magnitude placeholders chosen only to exercise the machinery. They are not\n",
    "measured device specifications. The point of this notebook is that every\n",
    "figure and number on the page is produced by code you can run and inspect,\n",
    "not by hand-drawn artwork.\n",
    "\n",
    "Requires: `numpy`, `matplotlib`.  Run:  `pip install numpy matplotlib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "ACCENT = '#15448c'; GREY = '#aab0b8'; INK = '#6b727c'\n",
    "plt.rcParams.update({'font.family': 'serif', 'font.size': 11,\n",
    "                     'axes.edgecolor': INK, 'text.color': INK,\n",
    "                     'xtick.color': INK, 'ytick.color': INK,\n",
    "                     'axes.labelcolor': INK})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 1 — coherent enhancement of the effective cross-section\n",
    "\n",
    "When N constituents respond in phase, scattering amplitudes add, so the\n",
    "effective cross-section scales as N^2 rather than N. This is the mechanism\n",
    "behind the coherent enhancement discussed in the deep dives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = np.logspace(0, 6, 200)\n",
    "fig, ax = plt.subplots(figsize=(6.4, 4.0))\n",
    "ax.loglog(N, N**2, color=ACCENT, lw=2.2, label=r'coherent  $\\sigma_{\\rm eff}\\propto N^{2}$')\n",
    "ax.loglog(N, N,   color=GREY,  lw=2.0, ls='--', label=r'incoherent  $\\sigma_{\\rm eff}\\propto N$')\n",
    "ax.set_xlabel('Number of coherently coupled scatterers  N')\n",
    "ax.set_ylabel(r'Effective cross-section  $\\sigma_{\\rm eff}/\\sigma_0$')\n",
    "ax.grid(True, which='major', alpha=0.12)\n",
    "for s in ['top', 'right']: ax.spines[s].set_visible(False)\n",
    "ax.legend(frameon=False, loc='upper left')\n",
    "fig.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 2 — Monte-Carlo convergence of the aggregated output\n",
    "\n",
    "P(t) is a product of factors, each uncertain over an order of magnitude. We\n",
    "model each as log-normal and estimate the mean of their product by Monte-\n",
    "Carlo, showing the running estimate and its shrinking 95% confidence band\n",
    "converging to the analytic mean (available in closed form for a product of\n",
    "log-normals, which lets us check the estimator is unbiased)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(20260726)\n",
    "# ln-medians and ln-sigmas of five ILLUSTRATIVE multiplicative factors.\n",
    "mus = np.array([-6.0, -46.0, -24.0,  3.7, -0.7])\n",
    "sig = np.array([ 0.5,   0.4,   0.6,  0.2,  0.3])\n",
    "Nmax = 100_000\n",
    "s = np.exp(rng.normal(mus[:, None], sig[:, None], size=(5, Nmax))).prod(axis=0)\n",
    "true_mean = np.exp((mus + 0.5 * sig**2).sum())   # analytic E[P]\n",
    "n = np.arange(1, Nmax + 1)\n",
    "mean = np.cumsum(s) / n\n",
    "sd   = np.sqrt(np.maximum(np.cumsum(s * s) / n - mean**2, 0))\n",
    "nk = np.unique(np.geomspace(2, Nmax, 120).astype(int))\n",
    "band = 1.96 * sd[nk - 1] / np.sqrt(nk) / true_mean\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6.4, 4.0))\n",
    "ax.fill_between(nk, 1 - band, 1 + band, color=ACCENT, alpha=0.13, label='95% confidence band')\n",
    "ax.plot(nk, mean[nk - 1] / true_mean, color=ACCENT, lw=1.8, label='running estimate')\n",
    "ax.axhline(1.0, color=INK, lw=1.0, ls=':', label='analytic mean')\n",
    "ax.set_xscale('log'); ax.set_ylim(0.7, 1.3)\n",
    "ax.set_xlabel('Number of Monte-Carlo samples  n')\n",
    "ax.set_ylabel(r'Estimate of $E[P]$  (normalised)')\n",
    "ax.grid(True, which='major', alpha=0.12)\n",
    "for sp in ['top', 'right']: ax.spines[sp].set_visible(False)\n",
    "ax.legend(frameon=False, loc='upper right')\n",
    "fig.tight_layout(); plt.show()\n",
    "\n",
    "print(f'analytic mean   : {true_mean:.4e}')\n",
    "print(f'MC estimate     : {mean[-1]:.4e}')\n",
    "print(f'relative error  : {abs(mean[-1] - true_mean) / true_mean:.4%}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The running estimate converges on the analytic mean and the relative error\n",
    "falls below one percent well before 100000 samples, confirming the estimator\n",
    "is unbiased and the aggregation is numerically well behaved. Change the seed\n",
    "or the `mus`/`sig` arrays to explore sensitivity."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"},
  "language_info": {"name": "python", "version": "3.x"}
 },
 "nbformat": 4,
 "nbformat_minor": 5
}