{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "11.4.1 Lippmann-Schwinger equation\n",
    "=============================\n",
    "**keys**: volume integral equation, volume potential, penetrable scatterer, Mie series, Luneburg lens"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from ngsolve.bem import *\n",
    "from ngsolve.krylovspace import GMRESSolver\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.special import spherical_jn, spherical_yn\n",
    "\n",
    "dxi = dx(bonus_intorder=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "## Scattering by a penetrable sphere\n",
    "\n",
    "The total field $u = u_{in} + u_s$, with incoming wave $u_{in} = e^{ikz}$ and outgoing scattered field $u_s$, solves\n",
    "\n",
    "$$\n",
    "\\Delta u + k^2 \\eta(x)^2 \\, u = 0 \\qquad \\text{in} \\, \\mathbb R^3,\n",
    "$$\n",
    "\n",
    "where the refractive index $\\eta$, the ratio of the outside to the local wave speed, equals 1 outside the scatterer $D$. Moving the material to the right hand side and subtracting $(\\Delta + k^2) u_{in} = 0$ gives, with the contrast $m = 1-\\eta^2$,\n",
    "\n",
    "$$\n",
    "(\\Delta + k^2) \\, u_s = k^2 (1-\\eta^2) \\, u = k^2 m u.\n",
    "$$\n",
    "\n",
    "The material acts as a source that vanishes outside $D$. Green's representation formula for the outgoing field $u_s$ gives $u_s = -V_k\\big((\\Delta + k^2)\\,u_s\\big) = -k^2 V_k(m u)$, and adding $u_{in}$ gives the Lippmann-Schwinger equation\n",
    "\n",
    "$$\n",
    "u + k^2 \\, V_k(m u) = u_{in} \\qquad \\text{in} \\, D.\n",
    "$$\n",
    "\n",
    "We multiply by a test function $v$, integrate over $D$, and solve for $u_h$\n",
    "\n",
    "$$\n",
    "(u_h, v) + k^2 \\big( V_k(m u_h), v \\big) = (u_{in}, v) \\qquad \\forall \\, v\n",
    "$$\n",
    "\n",
    "where $(\\cdot, \\cdot)$ denotes the integral over $D$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveLS(mesh, m, k, order):\n",
    "    fes = L2(mesh, order=order, complex=True)\n",
    "    u,v = fes.TnT()\n",
    "    with TaskManager():\n",
    "        M = BilinearForm(u*v*dx).Assemble()\n",
    "        V = HelmholtzSL(m*u*dxi, k)*v*dxi\n",
    "        rhs = LinearForm(exp(1j*k*z)*v*dxi).Assemble()\n",
    "        inv = GMRESSolver(M.mat + k**2*V.mat, M.mat.Inverse(), tol=1e-8)\n",
    "        gfu = GridFunction(fes)\n",
    "        gfu.vec.data = inv * rhs.vec\n",
    "    return gfu, inv.iterations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "Once $u_h$ is known, $u_s = -k^2 V_k(m u_h)$ can be evaluated anywhere, also inside $D$. We use the square $[-3,3]^2$ in the $xz$-plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "plane = WorkPlane(Axes((0,0,0), n=Y, h=X)).RectangleC(6,6).Face()\n",
    "disk = WorkPlane(Axes((0,0,0), n=Y, h=X)).Circle(0,0,1).Face()\n",
    "pmesh = Mesh(OCCGeometry(Glue([plane-disk, disk])).GenerateMesh(maxh=0.2))\n",
    "screen = pmesh.Boundaries(\".*\")\n",
    "\n",
    "def ScatteredField(gfu, m, k):\n",
    "    u = gfu.space.TrialFunction()\n",
    "    return -k**2*HelmholtzSL(m*u*dxi, k)(gfu, screen)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "For the unit sphere with constant $\\eta$, the Mie series is the exact solution. With $k_1 = k \\eta$,\n",
    "\n",
    "$$\n",
    "u = \\sum_l i^l (2l+1) \\, a_l \\, j_l(k_1 r) \\, P_l(\\cos \\theta) \\quad \\text{in} \\, D, \\qquad\n",
    "u_s = \\sum_l i^l (2l+1) \\, b_l \\, h^{(1)}_l(k r) \\, P_l(\\cos \\theta) \\quad \\text{outside},\n",
    "$$\n",
    "\n",
    "where $a_l, b_l$ follow from the continuity of $u$ and $\\partial_r u$ at $r = 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Mie(k, eta):\n",
    "    k1 = k*eta\n",
    "    l = np.arange(int(k1)+16)\n",
    "    j, dj = spherical_jn(l, k), spherical_jn(l, k, True)\n",
    "    h, dh = j+1j*spherical_yn(l, k), dj+1j*spherical_yn(l, k, True)\n",
    "    j1, dj1 = spherical_jn(l, k1), spherical_jn(l, k1, True)\n",
    "    b = (k1*dj1*j-k*dj*j1) / (k*dh*j1-k1*dj1*h)\n",
    "    return l, (j+b*h)/j1, b, k1\n",
    "\n",
    "def MieOnAxis(k, eta, zs):\n",
    "    l, a, b, k1 = Mie(k, eta)\n",
    "    r, c = abs(zs)[:,None], 1j**l*(2*l+1)*np.where(zs < 0, -1, 1)[:,None]**l\n",
    "    i, o = abs(zs) < 1, abs(zs) >= 1\n",
    "    u = np.exp(1j*k*zs)\n",
    "    u[i] = (c[i]*a*spherical_jn(l, k1*r[i])).sum(axis=1)\n",
    "    u[o] += (c[o]*b*(spherical_jn(l, k*r[o])+1j*spherical_yn(l, k*r[o]))).sum(axis=1)\n",
    "    return u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "k, eta = 2, sqrt(2)\n",
    "m = 1-eta**2\n",
    "zs = np.linspace(-3, 3, 300)\n",
    "zaxis = pmesh(0*zs, 0*zs, zs, BND)\n",
    "inner = abs(zs) < 0.95\n",
    "zin = zs[inner]\n",
    "uz_mie = MieOnAxis(k, eta, zs)\n",
    "err = lambda u, uex: np.linalg.norm(u-uex) / np.linalg.norm(uex)\n",
    "\n",
    "cases = [(0.3, 0), (0.2, 0), (0.3, 1)]\n",
    "print (\" maxh  order  ndof  its  u_in+u_s   u_h\")\n",
    "uzs = []\n",
    "for h, order in cases:\n",
    "    mesh = Mesh(OCCGeometry(Sphere((0,0,0), 1)).GenerateMesh(maxh=h)).Curve(2)\n",
    "    gfu, its = SolveLS(mesh, m, k, order)\n",
    "    with TaskManager():\n",
    "        uz = np.exp(1j*k*zs) + ScatteredField(gfu, m, k)(zaxis)[:,0]\n",
    "    uzs.append(uz)\n",
    "    uh = gfu(mesh(0*zin, 0*zin, zin))[:,0]\n",
    "    print (f\" {h:4.2f}  {order:5d}  {gfu.space.ndof:4d}  {its:3d}  \"\n",
    "           f\"{err(uz[inner], uz_mie[inner]):.2e}  {err(uh, uz_mie[inner]):.2e}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(zs, uz_mie.real, \"k\", label=\"Mie\")\n",
    "for (h, order), uz in zip(cases, uzs):\n",
    "    plt.plot(zs, uz.real, \"--\", label=f\"maxh={h}, order={order}\")\n",
    "plt.axvspan(-1, 1, color=\"0.9\"); plt.xlabel(\"$z$\"); plt.ylabel(\"Re $u$\"); plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "The table compares, for $k = 2$ and $\\eta = \\sqrt 2$, the total field $u_{in} + u_s$ and $u_h$ itself with the Mie series on the $z$-axis inside $D$. Total field on the plane:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "uscat = GridFunction(H1(pmesh, order=3, complex=True))\n",
    "with TaskManager():\n",
    "    uscat.Set(ScatteredField(gfu, m, k), definedon=screen)\n",
    "Draw (exp(1j*k*z)+uscat, pmesh, min=-2, max=2, animate_complex=True, order=3, euler_angles=[90,0,0]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "## Variable index: Luneburg lens\n",
    "\n",
    "The Luneburg lens is a sphere with refractive index $\\eta = \\sqrt{2 - r^2}$: the index decreases from $\\sqrt 2$ at the center to 1 at the surface, so it matches the surrounding medium without a jump. In the ray-optics limit, the rays of an incoming plane wave are bent inside the sphere and all meet in one point on the opposite side. At $k = 5$ the wavelength $2\\pi/5 \\approx 1.26$ is comparable to the radius, and the focus is inside the sphere near $z \\approx 0.8$. For larger $k$ it moves toward $z = 1$. `SolveLS` is reused unchanged, now with $k = 5$ and the coefficient function $m = 1 - \\eta^2 = r^2 - 1$ as contrast:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 5\n",
    "m = x**2+y**2+z**2-1\n",
    "mesh = Mesh(OCCGeometry(Sphere((0,0,0), 1)).GenerateMesh(maxh=0.25)).Curve(2)\n",
    "gfu, its = SolveLS(mesh, m, k, order=1)\n",
    "print (\"ndof =\", gfu.space.ndof, \", GMRES iterations:\", its)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    uscat.Set(ScatteredField(gfu, m, k), definedon=screen)\n",
    "Draw(exp(1j*k*z)+uscat, pmesh, animate_complex=True, order=3, euler_angles=[90,0,0]);\n",
    "ulens = abs(np.exp(1j*k*zs) + uscat(zaxis)[:,0])\n",
    "print (\"max |u| =\", ulens.max(), \"at z =\", zs[ulens.argmax()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.13.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}