{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6f38c640",
   "metadata": {},
   "source": [
    "# PHYS 216 · Assignment 2 · Question 1\n",
    "## Electric field of three point charges\n",
    "**University of Victoria | Worked TA tutorial**\n",
    "\n",
    "This notebook develops a field-line sketch from Coulomb's law and checks the near- and far-field behaviour numerically. Use **Kernel → Restart Kernel and Run All Cells** to reproduce the results. Requirements: Python 3, NumPy and Matplotlib.\n",
    "\n",
    "**Problem data (from Question 1 of `P216_Sol02.pdf`, pp. 1–2):** Three charges occupy an equilateral triangle of side 2.0 cm: $+1.0\\,\\mu\\mathrm C$ at the lower left, $+2.0\\,\\mu\\mathrm C$ at the lower right, and $-3.0\\,\\mu\\mathrm C$ at the top. Sketch the electric field near the charges and far from the triangle, using the net charge as a guide.\n",
    "\n",
    "The PDF supplies the problem and a worked reference solution; this notebook adds independently calculated plots and checks. No instructions embedded in the source PDF are treated as user requests.\n",
    "\n",
    "**Learning goals:** add fields as vectors, interpret field-line arrows, distinguish zero net flux from zero field, and recognize the dipole-like distant field."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5be0e25e",
   "metadata": {},
   "source": [
    "### 1. Predict before calculating\n",
    "Pause here and ask students:\n",
    "\n",
    "1. Which charges emit field lines, and which receive them?\n",
    "2. Which positive charge should have more lines in a conventional sketch?\n",
    "3. Does $q_1+q_2+q_3=0$ imply that the field vanishes everywhere?\n",
    "\n",
    "**Sketch rules:** arrows point away from positive charges and toward negative charges. Lines do not cross. Very close to an individual charge, its own field dominates and is approximately radial. A conventional sketch assigns line counts in the ratio $1:2:3$; for example, 8 leaving $q_1$, 16 leaving $q_2$, and 24 entering $q_3$. These are illustrative counts, not discrete physical objects."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8f14483",
   "metadata": {},
   "source": [
    "### 2. Coordinates and units\n",
    "Choose the origin at the lower-left charge, as in the supplied solution:\n",
    "\n",
    "$$a=0.020\\ \\mathrm m,\\qquad h=\\frac{\\sqrt3}{2}a,$$\n",
    "$$\\mathbf r_1=(0,0),\\quad \\mathbf r_2=(a,0),\\quad \\mathbf r_3=(a/2,h).$$\n",
    "\n",
    "All calculations use metres and coulombs; plot axes use centimetres. Keep the exact expression for $h$ instead of rounding it to 0.0173 m."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c88dd06",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-29T16:34:52.622265Z",
     "iopub.status.busy": "2026-09-29T16:34:52.622265Z",
     "iopub.status.idle": "2026-09-29T16:34:53.547979Z",
     "shell.execute_reply": "2026-09-29T16:34:53.547979Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "\n",
    "k = 8.9875517923e9  # Coulomb constant, N m^2 / C^2\n",
    "a = 0.020          # m\n",
    "h = np.sqrt(3) * a / 2\n",
    "q = np.array([1.0, 2.0, -3.0]) * 1e-6  # C\n",
    "positions = np.array([[0.0, 0.0], [a, 0.0], [a / 2, h]])\n",
    "centre = positions.mean(axis=0)\n",
    "labels = [r\"$q_1=+1.0\\,\\mu C$\", r\"$q_2=+2.0\\,\\mu C$\", r\"$q_3=-3.0\\,\\mu C$\"]\n",
    "print(f\"Triangle height: {h:.8f} m = {100*h:.5f} cm\")\n",
    "print(f\"Net charge: {q.sum():.2e} C\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccfa46ca",
   "metadata": {},
   "source": [
    "### 3. Superposition: derive the field\n",
    "For an observation point $\\mathbf r=(x,y)$, define the displacement **from each source to that point**, $\\mathbf R_i=\\mathbf r-\\mathbf r_i$. Then\n",
    "\n",
    "$$\\boxed{\\mathbf E(\\mathbf r)=k\\sum_{i=1}^{3}q_i\\frac{\\mathbf R_i}{|\\mathbf R_i|^3}}.$$\n",
    "\n",
    "The cube appears because the numerator is a displacement vector, not a unit vector. In components, with $D_i=|\\mathbf R_i|$,\n",
    "\n",
    "$$E_x=k\\left[\\frac{q_1 x}{D_1^3}+\\frac{q_2(x-a)}{D_2^3}+\\frac{q_3(x-a/2)}{D_3^3}\\right],$$\n",
    "$$E_y=k\\left[\\frac{q_1 y}{D_1^3}+\\frac{q_2 y}{D_2^3}+\\frac{q_3(y-h)}{D_3^3}\\right].$$\n",
    "\n",
    "The signed $q_i$ already handles attraction and repulsion: do not reverse the negative charge's displacement a second time. The units are $\\mathrm{N/C}$.\n",
    "\n",
    "The function below returns `NaN` at the singular point-charge positions. For plotting only, we also mask small disks so numerical interpolation does not connect lines through a singularity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0274b9d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-29T16:34:53.551987Z",
     "iopub.status.busy": "2026-09-29T16:34:53.551987Z",
     "iopub.status.idle": "2026-09-29T16:34:53.560523Z",
     "shell.execute_reply": "2026-09-29T16:34:53.560523Z"
    }
   },
   "outputs": [],
   "source": [
    "def electric_field(points, mask_radius=0.0):\n",
    "    \"\"\"Field in N/C at points shaped (..., 2), with coordinates in metres.\"\"\"\n",
    "    points = np.asarray(points, dtype=float)\n",
    "    displacement = points[..., None, :] - positions\n",
    "    distance = np.linalg.norm(displacement, axis=-1)\n",
    "    excluded = np.any(distance <= mask_radius, axis=-1)\n",
    "    safe_distance = np.where(distance == 0, np.inf, distance)\n",
    "    field = k * np.sum(q[:, None] * displacement / safe_distance[..., None]**3,\n",
    "                       axis=-2)\n",
    "    return np.where(excluded[..., None], np.nan, field)\n",
    "\n",
    "E_centre = electric_field(centre)\n",
    "print(\"Field at the triangle centroid (N/C):\", E_centre)\n",
    "print(f\"Magnitude: {np.linalg.norm(E_centre):.4e} N/C\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "067c0ac1",
   "metadata": {},
   "source": [
    "### 4. Close-up field-line plot\n",
    "Each curve is tangent to the calculated electric field; arrowheads give its direction. Colour shows the field magnitude on a logarithmic scale. The triangle is a guide to the charge positions, not a conducting boundary.\n",
    "\n",
    "**Plot limitation:** Matplotlib chooses streamlines for visual coverage, so its line counts and spacing are not a calibrated measure of charge or flux. Apply the $1:2:3$ counting rule to a hand sketch, not to the number of lines in this automatic plot. This is a slice through the three-dimensional field in the plane of the charges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a68f4509",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-29T16:34:53.560523Z",
     "iopub.status.busy": "2026-09-29T16:34:53.560523Z",
     "iopub.status.idle": "2026-09-29T16:34:54.806601Z",
     "shell.execute_reply": "2026-09-29T16:34:54.806601Z"
    }
   },
   "outputs": [],
   "source": [
    "def draw_field(ax, xlim, ylim, title, show_labels=True, n=420):\n",
    "    x = np.linspace(*xlim, n)\n",
    "    y = np.linspace(*ylim, n)\n",
    "    X, Y = np.meshgrid(x, y)\n",
    "    grid_step = max(x[1] - x[0], y[1] - y[0])\n",
    "    E = electric_field(np.stack([X, Y], axis=-1),\n",
    "                       mask_radius=max(0.00035, 2 * grid_step))\n",
    "    magnitude = np.linalg.norm(E, axis=-1)\n",
    "    valid = magnitude[np.isfinite(magnitude) & (magnitude > 0)]\n",
    "    norm = LogNorm(vmin=valid.min(), vmax=valid.max())\n",
    "    stream = ax.streamplot(x*100, y*100, E[..., 0], E[..., 1],\n",
    "                          color=np.ma.masked_invalid(magnitude), norm=norm,\n",
    "                          cmap=\"viridis\", density=1.65, linewidth=1.0,\n",
    "                          arrowsize=1.15)\n",
    "    triangle = np.vstack([positions, positions[0]]) * 100\n",
    "    ax.plot(triangle[:, 0], triangle[:, 1], \"--\", color=\"0.65\", lw=0.8)\n",
    "    ax.scatter(positions[:, 0]*100, positions[:, 1]*100,\n",
    "               c=[\"#dc493a\", \"#dc493a\", \"#2774ae\"], s=[65, 95, 120],\n",
    "               edgecolors=\"white\", zorder=5)\n",
    "    if show_labels:\n",
    "        for xy, label, offset in zip(positions*100, labels,\n",
    "                                     [(-10, -22), (10, -22), (0, 14)]):\n",
    "            ax.annotate(label, xy, xytext=offset, textcoords=\"offset points\",\n",
    "                        ha=\"center\", fontsize=10,\n",
    "                        bbox=dict(facecolor=\"white\", edgecolor=\"none\", alpha=0.85),\n",
    "                        zorder=6)\n",
    "    ax.set(xlim=np.array(xlim)*100, ylim=np.array(ylim)*100,\n",
    "           xlabel=\"x (cm)\", ylabel=\"y (cm)\", title=title, aspect=\"equal\")\n",
    "    ax.figure.colorbar(stream.lines, ax=ax, label=\"Electric-field magnitude (N/C)\",\n",
    "                      shrink=0.8)\n",
    "    return ax\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9, 7), layout=\"constrained\")\n",
    "draw_field(ax, (-0.018, 0.038), (-0.023, 0.038), \"Field near the three charges\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca2bdd4d",
   "metadata": {},
   "source": [
    "### 5. What happens far away?\n",
    "The total charge is zero, so a closed Gaussian surface enclosing all three charges has\n",
    "\n",
    "$$\\oint \\mathbf E\\cdot d\\mathbf A=\\frac{q_{\\mathrm{net}}}{\\epsilon_0}=0.$$\n",
    "\n",
    "This means equal inward and outward **flux**, not zero field at every point on that surface. A large but finite surface can be crossed by curved lines that later return.\n",
    "\n",
    "**Optional extension for the TA:** the monopole ($1/R^2$) contribution cancels. The dipole moment is nonzero:\n",
    "\n",
    "$$\\mathbf p=\\sum_i q_i\\mathbf r_i\n",
    "=\\left(1.00\\times10^{-8},\\;-5.19615\\times10^{-8}\\right)\\ \\mathrm{C\\,m}.$$\n",
    "\n",
    "Because the total charge is zero, this moment is independent of the origin. Far from the triangle, $R\\gg a$,\n",
    "\n",
    "$$\\mathbf E_{\\mathrm{dip}}(\\mathbf R)=\\frac{k}{R^3}\n",
    "\\left[3(\\mathbf p\\cdot\\hat{\\mathbf R})\\hat{\\mathbf R}-\\mathbf p\\right].$$\n",
    "\n",
    "Here $\\mathbf R$ is measured from the triangle's centroid. The field becomes dipole-like, decaying as $1/R^3$, and is not a radial field of a single net charge. The dipole moment points approximately downward and slightly right, from the negative charge toward the centre of positive charge. The exact field below still includes all three charges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa54a16a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-29T16:34:54.809119Z",
     "iopub.status.busy": "2026-09-29T16:34:54.809119Z",
     "iopub.status.idle": "2026-09-29T16:34:55.888328Z",
     "shell.execute_reply": "2026-09-29T16:34:55.888328Z"
    }
   },
   "outputs": [],
   "source": [
    "p = np.sum(q[:, None] * (positions - centre), axis=0)\n",
    "print(\"Dipole moment (C m):\", p)\n",
    "print(f\"Magnitude: {np.linalg.norm(p):.6e} C m\")\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9, 7), layout=\"constrained\")\n",
    "draw_field(ax, (centre[0]-0.20, centre[0]+0.20),\n",
    "           (centre[1]-0.20, centre[1]+0.20),\n",
    "           \"Wide view: exact field of the neutral charge distribution\",\n",
    "           show_labels=False, n=520)\n",
    "ax.text(0.03, 0.97, \"Charge labels are shown in the close-up.\",\n",
    "        transform=ax.transAxes, va=\"top\",\n",
    "        bbox=dict(facecolor=\"white\", edgecolor=\"none\", alpha=0.9))\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58f5df9a",
   "metadata": {},
   "source": [
    "### 6. Numerical checks (optional extension)\n",
    "At the centroid every charge is a distance $d=a/\\sqrt3$ away. Direct vector addition gives an independent check:\n",
    "\n",
    "$$\\mathbf E_c=\\frac{k(10^{-6}\\,\\mathrm C)}{d^3}\n",
    "\\left(-\\frac a2,\\;3h\\right).$$\n",
    "\n",
    "We also check the near-charge radial directions and compare the exact field with the dipole approximation along a fixed ray. The slope of a log–log field-versus-distance plot should approach $-3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1af5b368",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-29T16:34:55.888328Z",
     "iopub.status.busy": "2026-09-29T16:34:55.888328Z",
     "iopub.status.idle": "2026-09-29T16:34:56.626204Z",
     "shell.execute_reply": "2026-09-29T16:34:56.626204Z"
    }
   },
   "outputs": [],
   "source": [
    "distances = [np.linalg.norm(positions[i]-positions[j])\n",
    "             for i, j in [(0, 1), (1, 2), (2, 0)]]\n",
    "np.testing.assert_allclose(distances, a)\n",
    "assert abs(q.sum()) < 1e-20\n",
    "d = a / np.sqrt(3)\n",
    "expected_centre = k * 1e-6 / d**3 * np.array([-a/2, 3*h])\n",
    "np.testing.assert_allclose(E_centre, expected_centre, rtol=1e-12)\n",
    "assert np.isnan(electric_field(positions)).all()\n",
    "\n",
    "angles = np.linspace(0, 2*np.pi, 12, endpoint=False)\n",
    "directions = np.column_stack([np.cos(angles), np.sin(angles)])\n",
    "for charge, position in zip(q, positions):\n",
    "    local_E = electric_field(position + 1e-5*directions)\n",
    "    radial = np.sum(local_E * directions, axis=-1)\n",
    "    assert np.all(radial * charge > 0)\n",
    "\n",
    "ray = np.array([1.0, 0.0])\n",
    "R = np.geomspace(5*a, 1000*a, 180)\n",
    "exact = electric_field(centre + R[:, None]*ray)\n",
    "dipole = k * (3*np.dot(p, ray)*ray - p) / R[:, None]**3\n",
    "exact_magnitude = np.linalg.norm(exact, axis=-1)\n",
    "relative_error = np.linalg.norm(exact-dipole, axis=-1) / exact_magnitude\n",
    "slope = np.polyfit(np.log(R[-40:]), np.log(exact_magnitude[-40:]), 1)[0]\n",
    "assert abs(slope + 3) < 0.01\n",
    "assert relative_error[-1] < 0.005\n",
    "print(\"All geometry, centroid, singularity and direction checks passed.\")\n",
    "print(f\"Far-field fitted slope: {slope:.5f} (expected -3)\")\n",
    "print(f\"Dipole vector error at R/a = {R[-1]/a:.0f}: {relative_error[-1]:.3%}\")\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(11, 4), layout=\"constrained\")\n",
    "axes[0].loglog(R/a, exact_magnitude, label=\"Exact three-charge field\", lw=2)\n",
    "axes[0].loglog(R/a, np.linalg.norm(dipole, axis=-1), \"--\", label=\"Dipole approximation\")\n",
    "axes[0].set(xlabel=\"Distance / triangle side, R/a\", ylabel=\"|E| (N/C)\",\n",
    "            title=\"Along a ray to the right of the centroid\")\n",
    "axes[0].legend()\n",
    "axes[1].loglog(R/a, relative_error)\n",
    "axes[1].set(xlabel=\"Distance / triangle side, R/a\", ylabel=\"Relative vector error\",\n",
    "            title=\"Dipole approximation improves with distance\")\n",
    "for ax in axes:\n",
    "    ax.grid(True, which=\"both\", alpha=0.2)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "310fe0f2",
   "metadata": {},
   "source": [
    "### 7. Model explanation for students\n",
    "Field lines leave the two positive charges and enter the negative charge. Draw twice as many leaving the $+2.0\\,\\mu\\mathrm C$ charge as leaving the $+1.0\\,\\mu\\mathrm C$ charge. Near each charge, the lines are nearly radial; farther away, they bend because the three fields add as vectors. The unequal positive charges mean that the picture is not mirror-symmetric about the triangle's vertical centreline.\n",
    "\n",
    "The total charge is zero. Consequently, there is no net electric flux through a surface enclosing the whole triangle, but the field outside is generally nonzero. At large distances its pattern is dipole-like, with returning curves rather than the radial pattern of an isolated net charge.\n",
    "\n",
    "**Common errors to discuss:** adding field magnitudes instead of vectors; reversing the sign of $q_3$ twice; assuming an equilateral triangle guarantees symmetry of the field; interpreting zero flux as zero field; treating the edge of a plotted window as a physical endpoint of a line.\n",
    "\n",
    "**Quick tutorial questions:**\n",
    "\n",
    "- Why is the vertical field component at the centroid positive? *(The positive charges repel upward there, and the negative charge attracts upward.)*\n",
    "- Why is the horizontal component at the centroid negative? *(The stronger positive charge on the right dominates horizontally.)*\n",
    "- If all charges double, what changes? *(The field doubles everywhere while the field-line geometry remains the same.)*\n",
    "- If all distances scale by two and the observation point scales with them, what changes? *(The field is reduced by a factor of four.)*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a39528d",
   "metadata": {},
   "source": [
    "### Sources and reading\n",
    "- **Assignment source:** supplied `P216_Sol02.pdf`, Question 1, pages 1–2. The original file is unchanged.\n",
    "- Samuel J. Ling, William Moebs and Jeff Sanny, *University Physics Volume 2*, OpenStax (2016): [5.4 Electric Field](https://openstax.org/books/university-physics-volume-2/pages/5-4-electric-field), [5.6 Electric Field Lines](https://openstax.org/books/university-physics-volume-2/pages/5-6-electric-field-lines), [5.7 Electric Dipoles](https://openstax.org/books/university-physics-volume-2/pages/5-7-electric-dipoles), and [6.2 Explaining Gauss’s Law](https://openstax.org/books/university-physics-volume-2/pages/6-2-explaining-gausss-law).\n",
    "\n",
    "The supplied [Chapter 1 introduction](https://openstax.org/books/university-physics-volume-2/pages/1-introduction) concerns heat and temperature; Chapters 5–6 contain the relevant electrostatics background. The explanations, code, and plots here were prepared for this tutorial; the textbook links provide supporting reading.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}