From 3b598a9fe0ce6132188e4df4913fccf45c45b7ae Mon Sep 17 00:00:00 2001 From: Remco de Boer <29308176+redeboer@users.noreply.github.com> Date: Thu, 14 Sep 2023 16:58:16 +0200 Subject: [PATCH] DOC: visualize Riemann sheets for T matrix (#18) * DOC: add issue button * ENH: improve cutting mechanism for Riemann sheet plots --- docs/_config.yml | 1 + docs/lecture17.ipynb | 314 +++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 303 insertions(+), 12 deletions(-) diff --git a/docs/_config.yml b/docs/_config.yml index 11c11b9..e6598e6 100644 --- a/docs/_config.yml +++ b/docs/_config.yml @@ -13,6 +13,7 @@ execute: timeout: 600 html: + use_issues_button: true use_repository_button: true launch_buttons: diff --git a/docs/lecture17.ipynb b/docs/lecture17.ipynb index f7ca0e2..25ec75a 100755 --- a/docs/lecture17.ipynb +++ b/docs/lecture17.ipynb @@ -63,7 +63,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "tags": [] + }, "source": [ "## Riemann sheets" ] @@ -103,7 +105,6 @@ " funcs: list[Callable],\n", " func_unicode: str,\n", " boundaries: tuple[complex, float] | tuple[complex, complex] = (0, 1),\n", - " T_cutoff: float | None = None,\n", " resolution: int | tuple[int, int] = 50,\n", " colorize: bool = True,\n", " mask: Callable[[np.ndarray, np.ndarray], bool] | None = None,\n", @@ -111,13 +112,17 @@ " X, Y = create_meshgrid(boundaries, resolution)\n", " Z = X + Y * 1j\n", " T = [f(Z) for f in funcs]\n", - " if T_cutoff is not None:\n", - " for t in T:\n", - " t[np.abs(t) > T_cutoff] = np.nan\n", " if mask is not None:\n", - " the_mask = mask(X, Y)\n", + " the_mask = np.full(Z.shape, False)\n", + " for t in T:\n", + " the_mask |= mask(Z, t)\n", + " if np.all(the_mask):\n", + " raise ValueError(\"All points were masked away\")\n", " X[the_mask] = np.nan\n", " Y[the_mask] = np.nan\n", + " Z[the_mask] = np.nan\n", + " for t in T:\n", + " t[the_mask] = np.nan\n", "\n", " vmax = max(max(t.imag.max(), t.real.max()) for t in T)\n", " colorscale = [[i, c] for i, c in enumerate(DEFAULT_PLOTLY_COLORS)]\n", @@ -176,7 +181,17 @@ " return np.meshgrid(\n", " np.linspace(x1, x2, num=x_res),\n", " np.linspace(y1, y2, num=y_res),\n", - " )" + " )\n", + "\n", + "\n", + "def cut_t(\n", + " cutoff: float | tuple[float, float]\n", + ") -> Callable[[np.ndarray, np.ndarray], bool]:\n", + " if isinstance(cutoff, tuple):\n", + " re_cut, im_cut = cutoff\n", + " else:\n", + " re_cut, im_cut = cutoff, cutoff\n", + " return lambda z, t: (np.abs(t.real) > re_cut) | (np.abs(t.imag) > im_cut)" ] }, { @@ -213,7 +228,7 @@ " lambda z: +1 / np.sqrt(z),\n", " ],\n", " func_unicode=\"1/±√z\",\n", - " T_cutoff=9,\n", + " mask=cut_t(10),\n", ")" ] }, @@ -246,7 +261,7 @@ " ],\n", " func_unicode=\"log z\",\n", " boundaries=(0, np.e**2),\n", - " T_cutoff=10,\n", + " mask=cut_t((np.e, np.nan)),\n", ")" ] }, @@ -334,6 +349,9 @@ "cell_type": "code", "execution_count": null, "metadata": { + "jupyter": { + "source_hidden": true + }, "mystnb": { "code_prompt_show": "Show Riemann surface of ⁺√z" }, @@ -348,7 +366,7 @@ "plot_riemann_surfaces(\n", " funcs=[sp.lambdify(z, SignedSqrt(z).doit())],\n", " func_unicode=\"⁺√z\",\n", - " mask=lambda x, y: (np.abs(y) < 1e-5) & (x > 0),\n", + " mask=lambda z, t: (np.abs(z.imag) < 1e-5) & (z.real > 0),\n", " resolution=(30, 301),\n", ")" ] @@ -468,6 +486,9 @@ "cell_type": "code", "execution_count": null, "metadata": { + "jupyter": { + "source_hidden": true + }, "tags": [ "hide-input", "full-width" @@ -484,9 +505,278 @@ " boundaries=(240 - 40j, 320 + 40j),\n", " colorize=False,\n", " resolution=(50, 401),\n", - " mask=lambda x, y: np.abs(y) < 1e-6,\n", + " mask=lambda z, t: np.abs(z.imag) == 0,\n", ")" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## T-matrix definition" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "mystnb": { + "code_prompt_show": "Define expression classes" + }, + "tags": [ + "hide-input", + "scroll-input" + ] + }, + "outputs": [], + "source": [ + "@implement_doit_method\n", + "class S(UnevaluatedExpression):\n", + " is_commutative = True\n", + " is_real = False\n", + "\n", + " def __new__(cls, s, m, Mρ, GV, fπ, g0, sign=+1, **hints) -> Sigma:\n", + " return create_expression(cls, s, m, Mρ, GV, fπ, g0, sign, **hints)\n", + "\n", + " def evaluate(self) -> sp.Expr:\n", + " s, m, Mρ, GV, fπ, g0, sign = self.args\n", + " return 1 - 2 * sp.I * Sigma(s, m) / (16 * sp.pi) * T(*self.args)\n", + "\n", + " def _latex(self, printer, *args) -> str:\n", + " s = printer._print(self.args[0])\n", + " sign = self.args[-1]\n", + " number = \"I\" if sign < 0 else \"II\"\n", + " return f\"S_{{{number}}}({s})\"\n", + "\n", + "\n", + "@implement_doit_method\n", + "class T(UnevaluatedExpression):\n", + " is_commutative = True\n", + " is_real = False\n", + "\n", + " def __new__(cls, s, m, Mρ, GV, fπ, g0, sign=+1, **hints) -> Sigma:\n", + " return create_expression(cls, s, m, Mρ, GV, fπ, g0, sign, **hints)\n", + "\n", + " def evaluate(self) -> sp.Expr:\n", + " s, m, Mρ, GV, fπ, g0, sign = self.args\n", + " return 1 / (1 / V1(s, m, Mρ, GV, fπ) - G(s, m, g0, sign))\n", + "\n", + " def _latex(self, printer, *args) -> str:\n", + " s = printer._print(self.args[0])\n", + " sign = self.args[-1]\n", + " number = \"I\" if sign < 0 else \"II\"\n", + " return f\"T_{{{number}}}({s})\"\n", + "\n", + "\n", + "@implement_doit_method\n", + "class V1(UnevaluatedExpression):\n", + " is_commutative = True\n", + " is_real = False\n", + "\n", + " def __new__(cls, s, m, Mρ, GV, fπ, **hints) -> Sigma:\n", + " return create_expression(cls, s, m, Mρ, GV, fπ, **hints)\n", + "\n", + " def evaluate(self) -> sp.Expr:\n", + " s, m, Mρ, GV, fπ = self.args\n", + " return -(2 * p2(s, m)) / (3 * fπ**2) * (\n", + " 1 - GV**2 / fπ**2 * 2 * s / (s - Mρ**2)\n", + " ) - GV**2 / fπ**4 * p2(s, m) * h(Mρ**2 / (2 * p2(s, m)))\n", + "\n", + " def _latex(self, printer, *args) -> str:\n", + " s = printer._print(self.args[0])\n", + " return Rf\"V_1\\left({s}\\right)\"\n", + "\n", + "\n", + "@implement_doit_method\n", + "class h(UnevaluatedExpression):\n", + " is_commutative = True\n", + "\n", + " def __new__(cls, a, **hints) -> Sigma:\n", + " return create_expression(cls, a, **hints)\n", + "\n", + " def evaluate(self) -> sp.Expr:\n", + " a = self.args[0]\n", + " return -sp.Mul(\n", + " sp.Rational(2, 3),\n", + " (1 + 6 * a + 3 * a**2),\n", + " evaluate=False,\n", + " ) + a * (2 + 3 * a + a**2) * sp.log(1 + 2 / a)\n", + "\n", + " def _latex(self, printer, *args) -> str:\n", + " a = printer._print(self.args[0])\n", + " return Rf\"h\\left({a}\\right)\"\n", + "\n", + "\n", + "@implement_doit_method\n", + "class p2(UnevaluatedExpression):\n", + " is_commutative = True\n", + "\n", + " def __new__(cls, s, m, **hints) -> Sigma:\n", + " return create_expression(cls, s, m, **hints)\n", + "\n", + " def evaluate(self) -> sp.Expr:\n", + " s, m = self.args\n", + " return s / 4 - m**2\n", + "\n", + " def _latex(self, printer, *args) -> str:\n", + " s = printer._print(self.args[0])\n", + " return Rf\"p^2\\left({s}\\right)\"\n", + "\n", + "\n", + "a, Mρ, GV, fπ = sp.symbols(\"a m_rho, G_V f_pi\")\n", + "_exprs = [\n", + " S(s, m, Mρ, GV, fπ, g0, sign=-1),\n", + " T(s, m, Mρ, GV, fπ, g0, sign=-1),\n", + " T(s, m, Mρ, GV, fπ, g0, sign=+1),\n", + " V1(s, m, Mρ, GV, fπ),\n", + " h(a),\n", + " p2(s, m),\n", + "]\n", + "Math(aslatex({e: e.doit(deep=False) for e in _exprs}))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "gv = sp.Symbol(\"g_v\")\n", + "substitutions = {\n", + " fπ: 87.3,\n", + " GV: sp.sqrt(gv**2 * fπ**2) / 2,\n", + " gv: 1,\n", + " m: 139,\n", + " Mρ: 770,\n", + " g0: -3,\n", + "}\n", + "Math(aslatex(substitutions))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "T_exprs = [\n", + " T(s, m, Mρ, GV, fπ, g0, sign).doit().xreplace(substitutions).xreplace(substitutions)\n", + " for sign in [-1, +1]\n", + "]\n", + "T_funcs = [sp.lambdify(s, expr) for expr in T_exprs]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "hide-input", + "scroll-input" + ] + }, + "outputs": [], + "source": [ + "x = np.linspace(500, 1_100, num=200)\n", + "y = np.linspace(1e-5, 150, num=100)\n", + "X, Yn = np.meshgrid(x, -y)\n", + "X, Yp = np.meshgrid(x, +y)\n", + "Zn = X + Yn * 1j\n", + "Zp = X + Yp * 1j\n", + "Tn = T_funcs[1](Zn**2)\n", + "Tp = T_funcs[0](Zp**2)\n", + "\n", + "vmax = 100\n", + "colorscale = [[i, c] for i, c in enumerate(DEFAULT_PLOTLY_COLORS)]\n", + "sty = lambda t: dict(\n", + " cmin=-vmax,\n", + " cmax=+vmax,\n", + " colorscale=\"RdBu_r\",\n", + " surfacecolor=t.imag,\n", + ")\n", + "Sn = go.Surface(x=X, y=Yn, z=Tn.real, **sty(Tn), name=\"Unphysical\")\n", + "Sp = go.Surface(x=X, y=Yp, z=Tp.real, **sty(Tp), name=\"Physical\", colorbar_title=\"Re T\")\n", + "y = Yp[0]\n", + "z = x + y * 1j\n", + "line = go.Scatter3d(\n", + " x=x,\n", + " y=y,\n", + " z=T_funcs[0](z**2).real,\n", + " marker=dict(size=0),\n", + " line=dict(color=\"darkgreen\", width=1),\n", + ")\n", + "fig = go.Figure(data=[Sn, Sp, line])\n", + "fig.update_layout(height=550, width=600)\n", + "fig.update_scenes(\n", + " xaxis_title_text=\"Re s\",\n", + " yaxis_title_text=\"Im s\",\n", + " zaxis_title_text=\"Im T\",\n", + " zaxis_range=[-vmax, +vmax],\n", + ")\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "hide-input", + "scroll-input" + ] + }, + "outputs": [], + "source": [ + "sty = lambda t: dict(\n", + " cmin=-vmax,\n", + " cmax=+vmax,\n", + " colorscale=\"RdBu_r\",\n", + " surfacecolor=t.real,\n", + ")\n", + "Sn = go.Surface(x=X, y=Yn, z=Tn.imag, **sty(Tn), name=\"Unphysical\")\n", + "Sp = go.Surface(x=X, y=Yp, z=Tp.imag, **sty(Tp), name=\"Physical\", colorbar_title=\"Im T\")\n", + "y = Yp[0]\n", + "z = x + y * 1j\n", + "line = go.Scatter3d(\n", + " x=x,\n", + " y=y,\n", + " z=T_funcs[0](z**2).imag,\n", + " marker=dict(size=0),\n", + " line=dict(color=\"darkgreen\", width=1),\n", + ")\n", + "fig = go.Figure(data=[Sn, Sp, line])\n", + "fig.update_layout(height=550, width=600)\n", + "fig.update_scenes(\n", + " xaxis_title_text=\"Re s\",\n", + " yaxis_title_text=\"Im s\",\n", + " zaxis_title_text=\"Re T\",\n", + " zaxis_range=[-vmax, +vmax],\n", + ")\n", + "fig.show()" + ] } ], "metadata": { @@ -505,7 +795,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.10.13" } }, "nbformat": 4,