Equipotentials in 3D

Getting started:

As usual we begin by firiing up the kernel and downloading the relevant graphics
     routines.  Execute the following statements:

2 + 2

\(Needs["Graphics`PlotField3D`"]; 

    Needs["Graphics`ContourPlot3D`"]; \)

The Monopole:
Equipotentials of a Point Charge in 3D.

One common way to portray potentials in 3D is to plot the value of the potential
	along some plane (or ``slice'') in 3D.  In order to do so, we first define a function 
	that is the electric potential for a point charge at the origin.  Please note the syntax:
	on the LHS the variables have underbars and on the RHS they don't. Also note 
	that we are using the `:='  to assign the meaning of the function.  This is called
	``delayed assignment'', and we use it because we don't want Mathematica to
	actually evaluate the RHS yet.  We want it to do that only when we supply the 
	values of x, y, and x.

V[x_, y_, z_] := (1/(Sqrt[(x^2) + (y^2) + (z^2)]))

We have defined above the electric potential for a point charge at the origin.  We have 
          have set the magnitude of the charge to be 4πϵ=1/K, in order to make things easier.
          
           You might wonder why one bothers with equipotentials.  First, they are very handy
         when trying to consider energy, and potential energy differences.   Second, they are a
         bit easier to interpret than field lines.   The function below will plot the gradient of
         -V[x,y,z], with  little arrows  representing the value of the vector field at each point.
         (It will have a few error messages since V[x,y,z] is infinite at the origin.  Don't 
         worry about them).

PlotGradientField3D[-V[x, y, z],  {x, -2.5, 2.5}, 
       {y, -2.5, 2.5},  {z, -2.5, 2.5},  VectorHeads -> True, 
      ScaleFactor -> True, MaxArrowLength -> .5]

Not very helpful, is it?  Even if we connected the arrows into continuous field
           lines, it wouldn't be much better.  Equipotentials are surfaces, and by drawing sets
           of them we can get a better feel for the system.
	
 	   You no doubt recall how to plot a surface in 3D, using the command Plot3D.  
	 Here we are plotting the value of the potential on the z=0 plane.   This should be
	 infinite at the origin.  Luckily, Mathematica is smart enough to chop off the peak
	 so that we don't try to plot the value at infinity.  Evaluate the next command to
	 get the potential.

Plot3D[V[x, y, 0], {x, -3, 3}, {y, -3, 3}]

As you can see, the value of the potential grows steeply as we approach the
	 origin.  Graphs such as these are the basis of the intuition of ``higher potential''
	 as being ``up-hill''.

(1)  Plot the surface potential for the z=1 and z=2 planes.  How do they differ 
		from the z=0 potential surface?  Be sure to note the actual numbers on
		vertical axis when making a comparison.

If we draw a path on the above surface that is at constant `height' or potential,
	   we draw out a line at constant potential.  However, we know that equipotentials
	   are really surfaces in 3D.
	   
	     We can plot such surfaces using the ContourPlot3D command, defined in the
	   package we just loaded.  You may wish to ask for help (?ContourPlot3D) on 
	   this command now that it is loaded.  Briefly, we must give it a function to 
	   evaluate, a range for the x, y and z values, and finally, a value of the potential
	   (the function) for which we are trying to find the equipotential.  Let's try it for
	   the point charge, and calculate the surface on which the potential takes on the
	   value V(x,y,z)=1..  While Mathematica is calculating, make sure you know what
	   you expect it should plot:  what is the equipotential surface for a point charge?
	   (This will take 2-3 minutes to evaluate.)

ContourPlot3D[
      V[x, y, z], {x, -2.1, 2.}, {y, -2.1, 2.}, {z, -2.1, 2}, 
      Contours -> {1.}]

This, while nice, should not be wholly unexpected.

The Dipole Case:

Let's look at a case with a bit more structure, that of a dipole.  Below we define the
	  potential of a dipole, consisting of a negative charge at y=-1, and a positive charge
	  at y=1.  Execute the command.  Note that we have named it Vd so as not to confuse
	  it with the potential of a point charge.

Vd[x_, y_, z_] := 
      -(1/(Sqrt[(x^2) + (((y - 1))^2) + (z^2)])) + 
        (1/(Sqrt[(x^2) + (((y + 1))^2) + (z^2)]))

(2)  Plot the surface potential for the z=0 surface, and the y=0 surface, and explain 
		their shapes.  Why are they so different?  Without plotting it, can you deduce 
		which one the x=0 surface will resemble?

We'd like plot several equipotential surfaces.  In this case ContourPlot3D would be
	  inefficient since it calculates the value of the potential at a set gridpoints throughout
	  the specified volume before plotting the surface each and  every time we
	  plot a surface.  It will be more efficient to store the value of the potential at every grid
	  point in a table, and just refer to it each time we need it.  The command below makes
	  a table on a 3D grid and names this array of numbers ``dpot''.  Note the trailing 
	  semi-colon.  It tells Mathematica to not pring the result to the screen.  Since the 
	  result is a big list of numbers, this is a good idea.  
	  
	     You may also note that we have wrapped the function Vd in an N[] function.  This
	  forces a numerical evaluation of the list, so that it isn't stored as symbolic square
	  roots.   The plots will go faster that way.

dpot = 
      Table[N[Vd[x, y, z]], {x, -3.1, 3, .5}, {y, -3.1, 3, .5}, {z, 
          -3.1, 3, .5}];

The command below will plot the equipotentials for V=0.5 and V=-0.5.  It will
	about 2 minutes to complete, execute the command and relax.  (The Timing 
	function evaluates its arguments, and reports how long it takes to evaluate the
	command.)

Timing[
      ListContourPlot3D[dpot, 
        MeshRange -> {{-3.1, 3}, {-3.1, 3}, {-3.1, 3}}, 
        PlotPoints -> {5, 7}, Contours -> {-.5, .5}]]

By symmetry you can see that for each positive equipotiential, we will ahve an
	 identically shaped surface reflected about the y=0 plane.  
	 
	    Below we have already calculated the equipotentials for V=-.5, -.1,  0, .1, and
	    0.5.  However, some surfaces are obscured by others.  We can persuade 
	    Mathematica to replot the surface by doing the following:   (1) click on the figure
	    below, so that a bounding square appears, (2) place your cursor over the lower
	    right corner until a small diagonal double-headed arrow appears, (3) click the
	    mouse button and resizf the box.  Mathematica will then redraw the box without
	    recalculating the figure.  You can watch as each surface is constructed, starting
	    at the back and working toward the front.

Timing[
      ListContourPlot3D[dpot, 
        MeshRange -> {{-3.1, 3}, {-3.1, 3}, {-3.1, 3}}, 
        PlotPoints -> {5, 7}, 
        Contours -> {-.1, -.5, 0., .1, .5}]]

(3)  Define a new potential function for two positive point charges located at y=1 
		and y=-1.  Construct a table of values over the same range as above.  Plot
		separate figures for the equipotentials, using ListContourPlot3D.
			(a) V=0.75
			(b) V=1.75
			(c) V=2.00
		(Each one will take about 2-3 minutes to plot.)   
	
	(4)  Based upon the above plots, answer the following:
		(a)  What will the equiotential surfaces look like far away from the origin? 
			Why?
		(b)  Can an equipotential be broken up into different closed surfaces, or 
			must all points at the same voltage lie on a single, connected surface?
			
	(5)  Print out this assignment, and label the above equipotentials with the value
		they represent.  (The positive y-axis sticks out from the page.)