{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c77bb8c5",
   "metadata": {},
   "source": [
    "# PHYS 216 — Assignment 3: all worked solutions\n",
    "**Gauss's law, conductors, and electric potential\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0a4ae8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.patches import Circle, PathPatch\n",
    "from matplotlib.path import Path\n",
    "\n",
    "EPS0 = 8.8541878128e-12            # Vacuum permittivity, F/m.\n",
    "K = 1 / (4*np.pi*EPS0)            # Coulomb constant, N m^2/C^2.\n",
    "E_CHARGE = 1.602176634e-19         # Positive magnitude of electron charge, C.\n",
    "plt.rcParams.update({\"font.size\": 11})\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be95fd46",
   "metadata": {},
   "source": [
    "## Problem 1 — Parallel conducting plates\n",
    "**Given:** area $A=400\\,\\mathrm{cm^2}=0.0400\\,\\mathrm{m^2}$, separation $d=0.020\\,\\mathrm m$, and $N=1.0\\times10^{12}$ electrons moved from left to right. Neglect edge effects, as requested.\n",
    "\n",
    "### (a) Surface charge densities\n",
    "The left plate loses negative charge and becomes positive; the right plate becomes negative:\n",
    "$$Q_L=+Ne,\\quad Q_R=?,\\qquad \\sigma_L=\\frac{Q_L}{A},\\quad\\sigma_R=?.$$\n",
    "Thus $\\boxed{\\sigma_L\\approx?,\\ \\sigma_R\\approx?\\ \\mathrm{C/m^2}}$.\n",
    "\n",
    "### (b) Field between the plates\n",
    "Each ideal sheet contributes $|\\sigma|/(2\\epsilon_0)$. Using Gauss's law find the electric field between two plates:\n",
    "$$\\boxed{\\mathbf E= ? \\hat{\\mathbf x}\\approx ? \\hat{\\mathbf x}\\ \\mathrm{N/C}}.$$\n",
    "Here $+x$ points right. The field points from the positive left plate toward the negative right plate. In this idealization the exterior fields cancel and the fields inside the metals are zero.\n",
    "\n",
    "### (c) Potential difference\n",
    "$$V_R-V_L=-\\int_L^R\\mathbf E\\cdot d\\mathbf l=-Ed.$$\n",
    "Therefore $\\boxed{V_R-V_L\\approx-9.0\\times10^3\\ \\mathrm V}$. The magnitude is about $9.0\\,\\mathrm{kV}$, and the **left plate is at higher potential**. Choosing $V_R=0$ is a convenient reference, not an additional physical assumption."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef48c6a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem 1: convert units before calculating.\n",
    "N = 1.0e12                       # Number of transferred electrons.\n",
    "A = 400 * 1e-4                   # 1 cm^2 = 10^-4 m^2.\n",
    "d = 2.0e-2                       # Separation in metres.\n",
    "Q_left = N * E_CHARGE             # Losing electrons leaves positive charge.\n",
    "Q_right = -Q_left\n",
    "sigma_left = Q_left / A\n",
    "sigma_right = Q_right / A\n",
    "E_plate = sigma_left / EPS0       # Both sheets contribute to the gap field.\n",
    "V_right_minus_left = -E_plate*d   # Potential decreases in the field direction.\n",
    "print(f\"Charges: left {Q_left:+.5e} C; right {Q_right:+.5e} C\")\n",
    "print(f\"Surface densities: {sigma_left:+.5e}, {sigma_right:+.5e} C/m^2\")\n",
    "print(f\"Gap field: {E_plate:.5e} N/C, to the right\")\n",
    "print(f\"V_right - V_left = {V_right_minus_left:.5e} V\")\n",
    "\n",
    "x_gap = np.linspace(0, d, 200)\n",
    "V_gap = E_plate*(d-x_gap)          # Set the right plate's potential to zero.\n",
    "fig, axes = plt.subplots(1, 2, figsize=(10, 3.8), layout=\"constrained\")\n",
    "axes[0].plot(x_gap*100, np.full_like(x_gap, E_plate/1e5))\n",
    "axes[0].set(xlabel=\"Position from left plate (cm)\", ylabel=\"E_x (10^5 N/C)\",\n",
    "            title=\"Uniform field in the gap\", ylim=(0, 5.5))\n",
    "axes[1].plot(x_gap*100, V_gap/1000)\n",
    "axes[1].set(xlabel=\"Position from left plate (cm)\", ylabel=\"V (kV)\",\n",
    "            title=\"Potential, with V_right = 0\")\n",
    "for ax in axes:\n",
    "    ax.grid(alpha=0.25)\n",
    "plt.show()\n",
    "assert np.isclose(Q_left + Q_right, 0)\n",
    "assert V_right_minus_left < 0\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4199026e",
   "metadata": {},
   "source": [
    "## Problem 2 — Field from a potential\n",
    "Given $V(x,y,z)=-xy^2z+4xy$, with coordinates in metres and potential in volts. The polynomial coefficients carry the implied units.\n",
    "\n",
    "### (a) Differentiate the potential\n",
    "$$\\mathbf E=-\\nabla V,\\qquad\n",
    "\\boxed{E_x=y^2z-4y,\\quad E_y=2xyz-4x,\\quad E_z=xy^2}.$$\n",
    "\n",
    "### (b) Field at $(1,1,1)$\n",
    "Plug in the numbers for the given point in each electric field component and fill in the below question marks with numbers:\n",
    "$$\\boxed{\\mathbf E=(? \\hat{\\mathbf i}+? \\hat{\\mathbf j}+?\\hat{\\mathbf k})\\ \\mathrm{N/C}}.$$\n",
    "$$|\\mathbf E|=\\sqrt{?}\\approx?\\ \\mathrm{N/C},\\qquad\n",
    "\\boxed{\\hat{\\mathbf E}=\\frac{?\\hat{\\mathbf i}+?\\hat{\\mathbf j}+?\\hat{\\mathbf k}}{\\sqrt{?}}}.$$\n",
    "The unit vector gives the direction separately from the magnitude.\n",
    "\n",
    "### (c) External work\n",
    "For a slow move with negligible change in kinetic energy,\n",
    "$$W_{\\rm ext}=\\Delta U=q[V_f-V_i].$$\n",
    "Calculate the potential values and find the work done by the agent. In the section below, include the name of each quantity in the corresponding print(\"...\") statement so you can easily identify it in the output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5f214e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def potential_xyz(points):\n",
    "    \"\"\"Potential in volts for one point or an array with final axis [x,y,z].\"\"\"\n",
    "    points = np.asarray(points, dtype=float)\n",
    "    x, y, z = np.moveaxis(points, -1, 0)\n",
    "    return -x*y**2*z + 4*x*y\n",
    "\n",
    "def field_xyz(points):\n",
    "    \"\"\"Electric field -grad(V), in N/C.\"\"\"\n",
    "    points = np.asarray(points, dtype=float)\n",
    "    x, y, z = np.moveaxis(points, -1, 0)\n",
    "    return np.stack([y**2*z-4*y, 2*x*y*z-4*x, x*y**2], axis=-1)\n",
    "\n",
    "start = np.zeros(3)\n",
    "finish = np.ones(3)\n",
    "q_test = 2.0e-6                   # Positive test charge, C.\n",
    "E_point = field_xyz(finish)\n",
    "E_magnitude = np.linalg.norm(E_point)\n",
    "E_direction = E_point/E_magnitude\n",
    "W_external = q_test*(potential_xyz(finish)-potential_xyz(start))\n",
    "# All of the field components are being calculated in this section, like E_point, E_magnitude, W_external.\n",
    "# Replace each ___ with the variable for the requested quantity.\n",
    "print(\"Field components (N/C):\", ___)\n",
    "print(\"Field magnitude (N/C):\", ___)\n",
    "print(\"Unit direction:\", ___)\n",
    "print(\"External work (J):\", ___)\n",
    "\n",
    "# Gaussian quadrature integrates these polynomial line integrals accurately.\n",
    "nodes, weights = np.polynomial.legendre.leggauss(8)\n",
    "t = (nodes+1)/2                   # Map integration nodes from [-1,1] to [0,1].\n",
    "weights = weights/2\n",
    "\n",
    "def field_integral_segment(p0, p1):\n",
    "    delta = np.asarray(p1, float)-np.asarray(p0, float)\n",
    "    path_points = np.asarray(p0, float) + t[:, None]*delta\n",
    "    return np.sum(weights * (field_xyz(path_points) @ delta))\n",
    "\n",
    "straight = field_integral_segment(start, finish)\n",
    "vertices = [start, [1,0,0], [1,1,0], finish]\n",
    "broken = sum(field_integral_segment(vertices[i], vertices[i+1]) for i in range(3))\n",
    "print(f\"Integral E dot dl: straight = {straight:.6f} V; segmented = {broken:.6f} V\")\n",
    "np.testing.assert_allclose([straight, broken], [-3., -3.], atol=1e-12)\n",
    "np.testing.assert_allclose(-q_test*np.array([straight, broken]), W_external)\n",
    "\n",
    "# Check the analytical gradient independently with central finite differences.\n",
    "step = 1e-5\n",
    "numeric_gradient = np.array([\n",
    "    (potential_xyz(finish+step*u)-potential_xyz(finish-step*u))/(2*step)\n",
    "    for u in np.eye(3)\n",
    "])\n",
    "np.testing.assert_allclose(E_point, -numeric_gradient, rtol=1e-8)\n",
    "np.testing.assert_allclose(E_magnitude, np.sqrt(14))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efdd6441",
   "metadata": {},
   "source": [
    "## Problem 3 — Conducting sphere and concentric shell\n",
    "Given $a=0.020\\,\\mathrm m$, $b=0.050\\,\\mathrm m$, $c=0.060\\,\\mathrm m$, sphere charge $Q_s=+5.0\\,\\mu\\mathrm C$, and shell net charge $Q_{\\rm sh}=-3.0\\,\\mu\\mathrm C$.\n",
    "\n",
    "### (a) Shell surface charges\n",
    "A Gaussian surface in the shell metal has $E=0$, so its enclosed charge must vanish:\n",
    "$$Q_s+Q_b=0\\quad\\Rightarrow\\quad\\boxed{Q_b=?\\,\\mu\\mathrm C}.$$\n",
    "Charge conservation on the shell gives\n",
    "$$Q_b+Q_c=Q_{\\rm sh}\\quad\\Rightarrow\\quad\\boxed{Q_c=?\\,\\mu\\mathrm C}.$$\n",
    "The inner solid sphere carries its $+5.0\\,\\mu\\mathrm C$ on its own surface.\n",
    "\n",
    "### (b) Field in all four regions\n",
    "For a concentric configuration, spherical symmetry makes the field radial and constant in magnitude on a spherical Gaussian surface. Thus $4\\pi r^2 E_r=Q_{\\rm enc}/\\epsilon_0$. Replace each question mark with corresponding charge value.\n",
    "$$\\boxed{E_r(r)=\\begin{cases}\n",
    "0,&0\\le r<a,\\\\\n",
    "k(?)/r^2 ,&a<r<b,\\\\\n",
    "0,&b<r<c,\\\\\n",
    "k(?)/r^2,&r>c.\n",
    "\\end{cases}}$$\n",
    "Both nonzero regions point radially outward. \n",
    "\n",
    "### (c) Potential, with $V(\\infty)=0$\n",
    "The electric potential can be found using $$V_r-V_\\infty=-\\int_\\infty^r\\mathbf E\\cdot d\\mathbf l.$$ \n",
    "Start from $\\infty$ and integrate inward, and matching boundary values will provide the answer. Consider $V(\\infty)=0$.\n",
    "For $r>c$: \n",
    "$$V_r-V_\\infty=-\\int_\\infty^r\\mathbf k\\frac{Q_{tot}}{r^2} d\\mathbf r.$$ \n",
    "so \n",
    "$$V(r)=k\\frac{Q_{\\mathrm{tot}}}{r},\\qquad r>c.$$\n",
    "Continue integrating inward through the different regions. Inside each conductor, the electric field is zero, so the potential remains constant throughout that conductor; it is not necessarily zero\n",
    "### (d) Graphs and boundaries\n",
    "$V$ is continuous everywhere, including $a,b,c$, but its radial derivative changes abruptly. $E_r$ jumps at **all three surfaces** because each carries surface charge. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a68830b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "a, b, c = 0.020, 0.050, 0.060  # Three radii, in metres.\n",
    "Q_sphere, Q_shell = 5.0e-6, -3.0e-6\n",
    "Q_inner = -Q_sphere              # Ensures zero field in the shell metal.\n",
    "Q_outer = Q_shell-Q_inner        # Preserve the shell's net charge.\n",
    "V_shell = K*Q_outer/c\n",
    "V_sphere = V_shell + K*Q_sphere*(1/a-1/b)\n",
    "\n",
    "# Use masks rather than evaluating 1/r at r=0.\n",
    "def radial_field(r):\n",
    "    r = np.asarray(r, dtype=float)\n",
    "    if np.any(r < 0):\n",
    "        raise ValueError(\"Radius must be nonnegative.\")\n",
    "    E = np.zeros_like(r)\n",
    "    gap = (r > a) & (r < b)\n",
    "    outside = r > c\n",
    "    E[gap] = K*Q_sphere/r[gap]**2\n",
    "    E[outside] = K*Q_outer/r[outside]**2\n",
    "    E[(r == a) | (r == b) | (r == c)] = np.nan  # Surface limits differ.\n",
    "    return E\n",
    "\n",
    "def radial_potential(r):\n",
    "    r = np.asarray(r, dtype=float)\n",
    "    if np.any(r < 0):\n",
    "        raise ValueError(\"Radius must be nonnegative.\")\n",
    "    V = np.full_like(r, V_shell)\n",
    "    V[r <= a] = V_sphere\n",
    "    gap = (r > a) & (r < b)\n",
    "    outside = r > c\n",
    "    V[gap] = V_shell + K*Q_sphere*(1/r[gap]-1/b)\n",
    "    V[outside] = K*Q_outer/r[outside]\n",
    "    return V\n",
    "\n",
    "print(f\"Shell inner charge = {Q_inner*1e6:+.1f} microC\")\n",
    "print(f\"Shell outer charge = {Q_outer*1e6:+.1f} microC\")\n",
    "print (\"*Replace each ___ with the variable for the requested quantity below. V_shell and V_sphere are calculated above.*\")\n",
    "print(f\"Shell potential = { ___} V\")\n",
    "print(f\"Sphere potential = { ___} V\")\n",
    "\n",
    "# Draw each field region separately to avoid false connecting lines at jumps.\n",
    "fig, axes = plt.subplots(2, 1, figsize=(9, 7), sharex=True, layout=\"constrained\")\n",
    "for low, high in [(0,a), (a,b), (b,c), (c,0.10)]:\n",
    "    r = np.linspace(low+1e-9, high-1e-9, 350)\n",
    "    axes[0].plot(r*100, radial_field(r)/1e7, color=\"tab:red\", lw=2)\n",
    "r = np.unique(np.r_[np.linspace(0, 0.10, 1200), a, b, c])\n",
    "axes[1].plot(r*100, radial_potential(r)/1e6, color=\"tab:blue\", lw=2)\n",
    "for ax in axes:\n",
    "    ax.axvspan(0, a*100, color=\"0.9\")       # Shade the solid conductor.\n",
    "    ax.axvspan(b*100, c*100, color=\"0.9\")   # Shade the shell metal.\n",
    "    for radius in [a,b,c]:\n",
    "        ax.axvline(radius*100, color=\"0.5\", ls=\":\", lw=1)\n",
    "    ax.set_xlim(0,10)\n",
    "    ax.grid(alpha=0.2)\n",
    "axes[0].set(ylabel=\"E_r (10^7 N/C)\", title=\"Concentric conductors: field and potential\")\n",
    "axes[1].set(xlabel=\"Radius r (cm)\", ylabel=\"Potential (MV)\")\n",
    "axes[0].text(1, 9, \"Sphere\", ha=\"center\")\n",
    "axes[0].text(5.5, 9, \"Shell\", ha=\"center\")\n",
    "plt.show()\n",
    "\n",
    "# Verify charge conservation and potential matching at all interfaces.\n",
    "np.testing.assert_allclose(Q_inner+Q_outer, Q_shell)\n",
    "np.testing.assert_allclose(V_shell + K*Q_sphere*(1/a-1/b), V_sphere)\n",
    "np.testing.assert_allclose(V_shell + K*Q_sphere*(1/b-1/b), V_shell)\n",
    "np.testing.assert_allclose(K*Q_outer/c, V_shell)\n",
    "probe = np.array([0.01,0.03,0.055,0.08])\n",
    "step = 1e-7\n",
    "minus_dVdr = -(radial_potential(probe+step)-radial_potential(probe-step))/(2*step)\n",
    "np.testing.assert_allclose(radial_field(probe), minus_dVdr, rtol=1e-8, atol=1e-3)\n",
    "surface_charge_densities = np.array([Q_sphere,Q_inner,Q_outer])/(4*np.pi*np.array([a,b,c])**2)\n",
    "field_jumps = np.array([K*Q_sphere/a**2, -K*Q_sphere/b**2, K*Q_outer/c**2])\n",
    "np.testing.assert_allclose(field_jumps, surface_charge_densities/EPS0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3dc19cc",
   "metadata": {},
   "source": [
    "### (e) Move the inner sphere off-centre\n",
    "The charges and outer spherical geometry stay fixed, and the sphere does not touch the shell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "646591e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Qualitative cross-section for part (e); lengths here are in cm.\n",
    "# Create a figure and one set of axes.\n",
    "# figsize=(8,8) makes the plot square, and layout=\"constrained\"\n",
    "# automatically adjusts spacing around labels and the title.\n",
    "fig, ax = plt.subplots(figsize=(8, 8), layout=\"constrained\")\n",
    "\n",
    "\n",
    "# ------------------------------------------------------------\n",
    "# GEOMETRY\n",
    "# ------------------------------------------------------------\n",
    "\n",
    "# Position of the centre of the small charged sphere.\n",
    "# The conducting shell is centred at (0,0), but the charged sphere\n",
    "# is shifted upward by 1.5 cm.\n",
    "offset = np.array([0., 1.5])\n",
    "\n",
    "\n",
    "# Convert the radii a, b, and c to centimetres.\n",
    "# ra = radius of the inner charged sphere\n",
    "# rb = inner radius of the conducting shell\n",
    "# rc = outer radius of the conducting shell\n",
    "ra, rb, rc = a*100, b*100, c*100\n",
    "\n",
    "\n",
    "# Draw the entire outer conducting shell as a grey circle\n",
    "# with radius rc centred at the origin.\n",
    "ax.add_patch(\n",
    "    Circle(\n",
    "        (0, 0),             # centre of circle\n",
    "        rc,                 # radius\n",
    "        facecolor=\"0.85\",   # grey interior\n",
    "        edgecolor=\"black\",  # black boundary\n",
    "        zorder=1            # drawing layer\n",
    "    )\n",
    ")\n",
    "\n",
    "\n",
    "# Draw a white circle of radius rb on top of the grey circle.\n",
    "# This creates the empty cavity inside the conductor.\n",
    "ax.add_patch(\n",
    "    Circle(\n",
    "        (0, 0),\n",
    "        rb,\n",
    "        facecolor=\"white\",\n",
    "        edgecolor=\"black\",\n",
    "        zorder=2\n",
    "    )\n",
    ")\n",
    "\n",
    "\n",
    "# Draw the small positively charged sphere inside the cavity.\n",
    "# Its centre is offset from the centre of the conducting shell.\n",
    "ax.add_patch(\n",
    "    Circle(\n",
    "        offset,             # centre at (0, 1.5 cm)\n",
    "        ra,                 # radius of charged sphere\n",
    "        facecolor=\"#ffd9d5\",\n",
    "        edgecolor=\"tab:red\",\n",
    "        zorder=4\n",
    "    )\n",
    ")\n",
    "\n",
    "\n",
    "# ------------------------------------------------------------\n",
    "# ELECTRIC FIELD INSIDE THE CAVITY\n",
    "# ------------------------------------------------------------\n",
    "\n",
    "# Construct 18 field lines around the charged sphere.\n",
    "# theta is the angular position around the small sphere.\n",
    "for theta in np.linspace(0, 2*np.pi, 18, endpoint=False):\n",
    "\n",
    "    # Unit vector pointing in the direction theta.\n",
    "    u = np.array([np.cos(theta), np.sin(theta)])\n",
    "\n",
    "\n",
    "    # Starting point of the field line.\n",
    "    # It lies on the surface of the small charged sphere.\n",
    "    start = offset + ra*u\n",
    "\n",
    "\n",
    "    # Projection of the displacement of the small sphere's centre\n",
    "    # onto the direction u.\n",
    "    projection = np.dot(offset, u)\n",
    "\n",
    "\n",
    "    # Calculate how far we must move from the displaced sphere's\n",
    "    # centre along direction u before reaching the cavity wall.\n",
    "    #\n",
    "    # This comes from solving the intersection condition\n",
    "    #\n",
    "    # |offset + distance*u| = rb\n",
    "    #\n",
    "    # because the cavity wall is a circle of radius rb centred at (0,0).\n",
    "    distance = (\n",
    "        -projection\n",
    "        + np.sqrt(\n",
    "            projection**2\n",
    "            + rb**2\n",
    "            - np.dot(offset, offset)\n",
    "        )\n",
    "    )\n",
    "\n",
    "\n",
    "    # Position where the field line reaches the inner wall\n",
    "    # of the conducting shell.\n",
    "    end = offset + distance*u\n",
    "\n",
    "\n",
    "    # Straight-line distance between the start and end points.\n",
    "    # This is used only to choose reasonable Bezier control points.\n",
    "    length = np.linalg.norm(end - start)\n",
    "\n",
    "\n",
    "    # First Bezier control point.\n",
    "    # It forces the field line initially to leave the charged sphere\n",
    "    # approximately perpendicular to its surface.\n",
    "    cp1 = start + 0.32*length*u\n",
    "\n",
    "\n",
    "    # Second Bezier control point.\n",
    "    # end/rb is a radial unit vector normal to the cavity surface.\n",
    "    # Therefore this makes the field line approach the conductor\n",
    "    # approximately perpendicular to the inner wall.\n",
    "    cp2 = end - 0.32*length*(end/rb)\n",
    "\n",
    "\n",
    "    # Four points defining a cubic Bezier curve:\n",
    "    # start point, control point 1, control point 2, end point.\n",
    "    vertices = np.array([start, cp1, cp2, end])\n",
    "\n",
    "\n",
    "    # Create the cubic Bezier path.\n",
    "    #\n",
    "    # MOVETO = move to the first point.\n",
    "    # CURVE4 = cubic Bezier curve using the remaining three points.\n",
    "    curve = Path(\n",
    "        vertices,\n",
    "        [\n",
    "            Path.MOVETO,\n",
    "            Path.CURVE4,\n",
    "            Path.CURVE4,\n",
    "            Path.CURVE4\n",
    "        ]\n",
    "    )\n",
    "\n",
    "\n",
    "    # Draw the electric-field line.\n",
    "    ax.add_patch(\n",
    "        PathPatch(\n",
    "            curve,\n",
    "            fill=False,         # don't fill the curve\n",
    "            color=\"tab:blue\",   # blue field lines\n",
    "            lw=1.2,             # line width\n",
    "            zorder=3\n",
    "        )\n",
    "    )\n",
    "\n",
    "\n",
    "    # Function returning the position on the cubic Bezier curve\n",
    "    # for a parameter t between 0 and 1.\n",
    "    #\n",
    "    # t = 0 gives the start point.\n",
    "    # t = 1 gives the end point.\n",
    "    def bezier(t):\n",
    "\n",
    "        return (\n",
    "            (1-t)**3 * start\n",
    "            + 3*(1-t)**2*t * cp1\n",
    "            + 3*(1-t)*t**2 * cp2\n",
    "            + t**3 * end\n",
    "        )\n",
    "\n",
    "\n",
    "    # Add a small arrowhead to the field line.\n",
    "    # xytext is slightly earlier on the curve than xy,\n",
    "    # so the arrow points from the charged sphere toward the conductor.\n",
    "    ax.annotate(\n",
    "        \"\",\n",
    "        xy=bezier(0.58),\n",
    "        xytext=bezier(0.48),\n",
    "        arrowprops=dict(\n",
    "            arrowstyle=\"->\",\n",
    "            color=\"tab:blue\"\n",
    "        ),\n",
    "        zorder=3\n",
    "    )\n",
    "\n",
    "\n",
    "# ------------------------------------------------------------\n",
    "# ELECTRIC FIELD OUTSIDE THE CONDUCTING SHELL\n",
    "# ------------------------------------------------------------\n",
    "\n",
    "# Outside a spherical conductor, the electric field is radial.\n",
    "# The off-centre position of the internal charge does not change\n",
    "# the spherical symmetry of the EXTERNAL field.\n",
    "for theta in np.linspace(0, 2*np.pi, 16, endpoint=False):\n",
    "\n",
    "    # Radial unit vector in the direction theta.\n",
    "    u = np.array([np.cos(theta), np.sin(theta)])\n",
    "\n",
    "\n",
    "    # Draw an outward-pointing arrow beginning at the outer\n",
    "    # surface of the conductor and extending outward.\n",
    "    ax.annotate(\n",
    "        \"\",\n",
    "        xy=7.4*u,        # arrow end\n",
    "        xytext=rc*u,     # arrow starts at conductor surface\n",
    "        arrowprops=dict(\n",
    "            arrowstyle=\"->\",\n",
    "            color=\"tab:orange\",\n",
    "            lw=1.3\n",
    "        )\n",
    "    )\n",
    "\n",
    "\n",
    "# ------------------------------------------------------------\n",
    "# TEXT LABELS\n",
    "# ------------------------------------------------------------\n",
    "\n",
    "# Label the inner positively charged sphere.\n",
    "# The text \"+5 microC\" indicates its charge.\n",
    "ax.text(\n",
    "    0,\n",
    "    1,\n",
    "    \"+5 microC\\nE = 0\",\n",
    "    ha=\"center\",        # horizontal alignment\n",
    "    va=\"center\",        # vertical alignment\n",
    "    zorder=5\n",
    ")\n",
    "\n",
    "\n",
    "# Label the conducting material.\n",
    "# Inside the metal of a conductor in electrostatic equilibrium,\n",
    "# the electric field is zero.\n",
    "ax.text(\n",
    "    0,\n",
    "    -5.5,\n",
    "    \"Shell metal: E = 0\",\n",
    "    ha=\"center\",\n",
    "    fontsize=9,\n",
    "    zorder=5\n",
    ")\n",
    "\n",
    "\n",
    "# ------------------------------------------------------------\n",
    "# AXIS SETTINGS\n",
    "# ------------------------------------------------------------\n",
    "\n",
    "# Set:\n",
    "# x-axis limits,\n",
    "# y-axis limits,\n",
    "# equal x/y scaling,\n",
    "# x-axis label,\n",
    "# y-axis label,\n",
    "# and plot title.\n",
    "ax.set(\n",
    "    xlim=(-8, 8),\n",
    "    ylim=(-8, 8),\n",
    "    aspect=\"equal\",\n",
    "    xlabel=\"x (cm)\",\n",
    "    ylabel=\"y (cm)\",\n",
    "    title=\"Off-centre sphere: qualitative field-line sketch\"\n",
    ")\n",
    "\n",
    "\n",
    "# Display the completed figure.\n",
    "plt.show()\n",
    "\n",
    "\n",
    "# Print a message after the plot is generated.\n",
    "# It also emphasizes that the curved field lines inside the cavity\n",
    "# are qualitative/schematic rather than an exact numerical\n",
    "# electrostatic solution.\n",
    "print(\n",
    "    \"All numerical checks passed. \"\n",
    "    \"Off-centre gap curves are schematic.\"\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cd0d68b",
   "metadata": {},
   "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
}