CE3 Gauss's Law, Divergence, and Curl_Complete_student.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from shapegen import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computational Exercise 3: Gauss's Law\n",
    "\n",
    "In this exercise, we'll use computational methods to look at the shape independence of Gauss's Law. Please create markdown cells as needed to include your answers to the questions in your notebook. To to this, click the `+` icon near the top-left of the notebook to add a cell directly below the currently selected cell. Next, open the `Code` dropdown menu and change it to `Markdown`. Now when you run the cell, it should render as regular text (like this one).  \n",
    "\n",
    "## Overview\n",
    "\n",
    "The purpose of this exercise is to look at how Gauss's law is independent of surface shape, and to visualize the electric field over these shapes. There are three sections:\n",
    "* Calculating the flux from a given $Q_{enc}$ over a sphere\n",
    "* Calculating the flux from a given $Q_{enc}$ over a rectangular prism\n",
    "* Calculating the flux from a given $Q_{enc}$ over a cylinder\n",
    "\n",
    "__To calculate the flux, you'll use a numerical integration technique.__ In a normal flux calculation, you'd use $\\phi = \\oint_S  \\mathbf{E} ⋅ d \\mathbf{a}$, which is essentially a sum over infinitely small areas. To do a numerical integration, we look at the flux through a given face of the object, and sum all of those to get an approximation of the flux, so the equation would look like $\\phi \\approx \\sum_{i=1}^{n} \\mathbf{E}_i ⋅ \\mathbf{A}_i$.\n",
    "\n",
    "## Helpful Information\n",
    "\n",
    "__This exercises relies on a package called shapegen, which was written specifically to make this exercise much, much easier.__ It handles all of the computationally-annoying parts of this, like turning different shapes into a bunch of triangles, giving you the areas and locations of those triangles, and (the icing on the cake) plotting these triangularized shapes. For full documentation, there should be a file, shapgen.html, with all of the information included posted in your class documents. You're also welcome to open shapegen.py in a text editor and take a look!\n",
    "\n",
    "This package uses what's known as an object-oriented programming methodology, where you make an instance of an object that has certain attributes and functions attached to it. For example, you could make a Number object which has the attribute value and the functions add and subtract. To create the object, you would use\n",
    "<code>n = Number(1)</code> to declare a Number object with the value <code>1</code>. Then <code>n.value</code> would return 1, and <code>n.add(2)</code> would make <code>n = 3</code>. Note that to get an attribute there are no parentheses and for a function call there are, whether or not that function takes an argument. For more information on object-oriented programming, try taking CSC-207.\n",
    "\n",
    "### Declaring a Shape Object\n",
    "To generate a sphere, you can use:\n",
    "<pre><code>s = Sphere(r, res=n)</code></pre>\n",
    "Where `r` is the radius of the sphere and n is the resolution used to generate it. The resolution is optional (known as a keyword argument), and defaults to 16, so <code>c = Sphere(r)</code> works just as well.\n",
    "\n",
    "For a rectangular prism, use:\n",
    "<pre><code>p = Prism(corner_pt, res=n)</code></pre>\n",
    "Where `corner_pt` is the location of one of the corners of the prism. So to make a 2x2x2 cube, the declaration would be<code>p = Prism((1, 1, 1))</code>. Note that the prism will be centered at the origin, so make sure none of the coordinate entries are 0.\n",
    "\n",
    "For a cylinder, use:\n",
    "<pre><code>c = Cylinder(r, l, res=n)</code></pre>\n",
    "Where `r` is the radius and `l` is the length.\n",
    "\n",
    "## Getting Areas and Locations\n",
    "\n",
    "Since you'll be summing over the surface of each object, you'll need to get the location and area of each surface face (or Tile, in the package). For any object, __use the function calls <code>areas = shape_name.get_areas()</code> and <code>centroids = shape_name.get_centroids()</code>__ to get a list of the areas and centroids, respectively.\n",
    "\n",
    "Each area is given as a normal vector of the form $(x, y, z)$, and each centroid is given as a position vector of the same form. The area vector is a vector that points normal to the surface face and has a magnitude equal to the area of the face. By giving these to you as an vector, you can take the dot product of the area with a field and get the flux travelling normal to the face. This shows up often in the vector calculus that Griffiths uses.\n",
    "\n",
    "### Dealing with Vector Math with numpy\n",
    "\n",
    "The numpy (np) package in python has a number of useful tools to help you handle vectors. You can absolutely finish this assignment without doing vector math (you can unpack each element of any vector just like any array), but it'll be much, much easier if you use vectors. Trust me, I tried.\n",
    "\n",
    "Some useful vector operations:\n",
    "* [np.subtract](https://numpy.org/doc/stable/reference/generated/numpy.subtract.html): Find the difference (separation vector) between two vectors\n",
    "* [np.linalg.norm](https://numpy.org/doc/stable/reference/generated/numpy.linalg.norm.html): Find the magnitude of a vector\n",
    "* [np.dot](https://numpy.org/doc/stable/reference/generated/numpy.dot.html): Take the dot product of two vectors\n",
    "* [np.cross](https://numpy.org/doc/stable/reference/generated/numpy.cross.html): Take the cross product of two vectors\n",
    "\n",
    "Also, remember that you can find a unit vector by dividing a given vector by its magnitude.\n",
    "\n",
    "## Plotting the Shapes\n",
    "To plot a shape, you can use code such as:\n",
    "<pre><code>ax = shape_name.plot(color=e, cmap='plasma')\n",
    "ax.set_title(\"Sample plot\")\n",
    "plt.show()\n",
    "</code></pre>\n",
    "The `shape_name.plot()` function returns an Axis object, so that you can add your own titles, axis limits, and use other plotting functions.\n",
    "\n",
    "Note that color and cmap are both keyword arguments, so you can leave them out if you just want to plot the object. __The color argument takes an array of values of length <code>shape_name.len</code>.__ You should get this array by iterating through while calculating the fluxes (which will also have the same length). The cmap argumemt takes a [colormap](https://matplotlib.org/3.1.1/gallery/color/colormap_reference.html), or will just default to plasma."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# a) Calculating the Electric Field\n",
    "\n",
    "__First, use the template code below make a function that returns the electric field from a point charge at a given location.__ This function should take as inputs the magnitude of the charge `q`, the location of the charge `q_loc1` as a position vector of the form `(x, y, z)`, and the point of interest, also a position vector `r0`. Think about Griffith's equation for the electric field vector (eq. 1.3, excluding the summation, which we'll deal with separately):\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathbf{E}(\\mathbf{r}) = \\frac{1}{4 \\pi \\epsilon_{0}} \\frac{q}{\\mathscr{r}^2} \\hat{\\mathscr{r}}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Since we're taking all of our inputs as vectors, you'll have an easier time if you use numpy vector functions. The function should return $\\mathbf{E}$ as a vector of the form $(E_x, E_y, E_z)$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def get_E(q, q_loc, r0):\n",
    "    \"\"\"Return the electric field vector E=(Ex,Ey,Ez) due to charge \n",
    "        q at r0.\"\"\"\n",
    "    e_0 = 1 #Put the permitivitty of free space constant here\n",
    "    \n",
    "###YOUR CODE FOR FINDING THE E-Field HERE\n",
    "\n",
    "    return E\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# b) Predict Total Flux from a Charge\n",
    "\n",
    "In order to test your code, you'll want to know what the total flux through the Gaussian surface should be for a given charge. Use the cell below to calculate that prediction. Then in parts c, d, and e, you can sum the flux through the surface and compare with your prediction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Declare info about q\n",
    "q_enc = 1 # Charge of q\n",
    "q_loc = (0, 0, 0) #L ocation of q\n",
    "\n",
    "###YOUR CODE FOR FINDING THE TOTAL FLUX THROUGH A SURFACE FROM THIS CHARGE\n",
    "\n",
    "flux = 1 # Fix\n",
    "\n",
    "\n",
    "print(\"Total predicted flux: {:.3e} N m^2/C\".format(flux))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# c) Flux Through a Sphere\n",
    "Now that we have a function for the E-field, you can use it to find the flux of that electric field through a gaussian surface. \n",
    "\n",
    "__Use the next block to:__\n",
    "* __Create the sphere object and get the areas and centroids__\n",
    "* __Iterate through those to find the total flux and a list of the E-field on each surface face__\n",
    "* __Print the total flux__\n",
    "* __Plot the sphere using the E-fields to color the faces__.\n",
    "\n",
    "*Once you've done this, you can check your answer by calculating the flux over the surface by hand.*\n",
    "\n",
    "*I encourage you to try changing the size or dimensions of the sphere, the location of the enclosed charge, and the resolution of shapegen. You might even try putting the charge outside of the sphere! As you do so, notice how the flux at the surface does or doesn't change.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Generate sphere object\n",
    "r = 1 #Radius\n",
    "s = Sphere(r, res = 200)\n",
    "\n",
    "# Get areas and locations of each triangle\n",
    "areas = s.get_areas()\n",
    "centroids = s.get_centroids()\n",
    "\n",
    "# Dummy definitions of numpy arrays containing the fluxes and magnitude of the E-fields\n",
    "fluxes = np.zeros(s.len) #array of zeros with same length as number of triangles\n",
    "mag_efields = np.zeros(s.len) \n",
    "\n",
    "###YOUR CODE FOR FINDING E-Field AND FLUX HERE\n",
    "\n",
    "    \n",
    "# Calculate Total Flux\n",
    "flux = np.sum(fluxes)\n",
    "\n",
    "print(\"Total flux: {:.3e} N m^2/C\".format(flux))\n",
    "\n",
    "# Plot using magnitude of E-fields, since triangles have different areas\n",
    "ax = s.plot(color=np.array(mag_efields), cmap='plasma')\n",
    "ax.set_title(\"Electric Field over a Spherical Gaussian Surface\",\n",
    "            fontsize=25)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# d) Flux Through a Cube\n",
    "\n",
    "__Repeat part (c) but this time for a cube (prism).__ Note how the flux through the cube compares with the previously calculated flux through a sphere.\n",
    "\n",
    "*Once again, I encourage you to try changing the size or dimensions of the cube, the location of the enclosed charge, and the resolution of shapegen. As you do so, notice how the flux at the surface does or doesn't change.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Generate prism Object\n",
    "corner_pt = (1, 1, 2)\n",
    "p = Prism(corner_pt,res = 20)\n",
    "\n",
    "###YOUR CODE FOR FINDING AND PLOTTING E-Field AND FLUX HERE\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# e) Flux Through a Cylinder\n",
    "\n",
    "__Repeat the steps in part (c) and (d) but this time for a cylinder.__ Note how the flux through the cylinder compares with the previously calculated flux through a sphere and cube.\n",
    "\n",
    "*Once again, I encourage you to try changing the size or dimensions of the cylinder, the location of the enclosed charge, and the resolution of shapegen. As you do so, notice how the flux at the surface does or doesn't change.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "r = 1 #Radius\n",
    "l = 4 #Length\n",
    "\n",
    "c = Cylinder(l, r)\n",
    "\n",
    "###YOUR CODE FOR FINDING AND PLOTTING E-Field AND FLUX HERE\n",
    "\n"
   ]
  }
 ],
 "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}