2D Field Lines and Electric Potential

This notebook plots the field lines and potential for a random/defined distribution of charges using NDSolve.  Adapted from the book Mathematica Graphics by Tom Wickham-Jones -  Section 19.4 Field Lines by Edan Scriven.  Generalised to random charge distribution by Tyrone Curtis.  January/February 2006.

Initial Setup

ClearAll["Global`*"]

Off[General::spell1];
Off[NDSolve::ndsz];
Off[NDSolve::mxst];
Off[Power::infy];
Off[General::indet];

To use the command LabelContourLines need to load the LabelContour package.  Change the next line to AppendTo[$Path,"path containing the folder Graphics containing the file LabelContour.m"] (unless the file is under the usual Mathematica Graphics directory containing the other Graphics packages).  This only needs to be done once each Mathematica session.

AppendTo[ $Path, "/Users/au/Ty/GoodFiles"];

Needs["Graphics`LabelContour`"]

Charge Setup and Potential calculation

Firstly define the plane in which the charges are to be located.

xmin = 0;
xmax = 6;
ymin = 0;
ymax = 6;

Currently the list of charges is generated randomly, however you can input your own charges in the form {{q1{x1,y1}},{q2,{x2,y2}},.....,{qn,{xn,yn}}}.

n = 6;

charges = {{1, {2.99, 3}}, {-1, {3.01, 
            3}}}

charges = 
        Table[{(-1)^Random[Integer], {Random[Real, {xmin, xmax}], 
              Random[Real, {ymin, ymax}]}}, {n}];

n = Len gth[charges]   ;

Separate the positive,negative and neutral charges in order to calculate field lines later.

positiveCharges = 
        Select[charges, #[[1]] > 0 &] ;
negativeCharges = 
        Select[charges, #[[1]] < 0 &];
neutralCharges = 
        Select[charges, #[[1]] == 0 &];

Potential calculation

potential[{q1,{x1,y1}}] is the individual contribution to the total electric potential due to a charge of strength q1 located at {x1,y1}.

potential[{q1_, {x1_, y1_}}]  :=  
      (q1/(Sqrt[(((x1 - x))^2) + (((y1 - y))^2)]))

Total electric potential is the sum of the contributions from the individual charges.  That is:

\(ϕ \((x, 
          y)\) = ∑\+\(i = q\)\%n q\_i\/\@\(\((x - x\_i)\)\^2 + \((y - y\_i)\)\^2\)\)

where ((Subscript[x, i]), (Subscript[y, i])) is the position of the (i^th)charge of strength (Subscript[q, i]).  The electric potential generated this way will be in units of Volts, assuming x and y are in Å and (Subscript[q, i])'s in units of electron charge.  Potential is in units of 4  π(Subscript[ϵ, 0]).

elecPotential = Total[potential /@ charges];

Contour Plot and Gradient Plot

chargePlot is a list of Graphics objects used to plot the positions of the charges - positive charges are plotted in blue, negative charges in red and neutral charges in green.

chargePlot = {RGBColor[0, 0, 1], Disk[#, 0.07] & /@ 
            positiveCharges[[All, 2, {1, 2}]], 
          RGBColor[1, 0, 0], Disk[#, 0.07] & /@ 
            negativeCharges[[All, 2, {1, 2}]], 
          RGBColor[0, 1, 0], Disk[#, 0.07] & /@ 
            neutralCharges[[All, 2, {1, 2}]]};

Contour Plot

Shows the variation of the potential throughout the plotting region.

c1 = 
        ContourPlot[elecPotential, {x, xmin, xmax}, {y, ymin, ymax}, 
          PlotPoints -> 256, ContourShading -> False, 
          Contours -> 10, Epilog -> chargePlot, 
          ImageSize -> {600, 600}, 
          DisplayFunction -> Identity];

c1 = 
        LabelContourLines[c1, LabelPlacement -> Automatic, 
          PointFactor -> 10];

c2 = 
        ContourPlot[elecPotential, {x, xmin, xmax}, {y, ymin, ymax}, 
          PlotPoints -> 256, ContourShading -> True, 
          Contours -> 10, Epilog -> chargePlot, 
          ImageSize -> {600, 600}, ColorFunction -> GrayLevel, 
          DisplayFunction -> Identity];

c2 = 
        LabelContourLegend[c2, HeightName -> "ϕ (Volts)", 
          LabelFont -> 11];

Gradient Plot

This plots the gradient field for the electric potential, so shows the electric field vectors at grid points in the plane.  The electric field vectors are shown in much greater detail further below when the field lines are plotted.

Needs["Graphics`PlotField`"]

PlotGradientField[-elecPotential, {x, xmin, xmax}, {y, ymin, 
          ymax}, ScaleFactor -> 1, MaxArrowLength -> 10, 
        PlotPoints -> 15, ColorFunction -> GrayLevel, 
        Epilog -> chargePlot, 
        ImageSize -> {500, 500}];

Electric Field lines

Electric field vectors

Next calculate the electric field vectors.  The electric field vectors are given by the following relationships:

(Subscript[E, x]) = -∂(ϕ/∂)x
(Subscript[E, y]) = -∂(ϕ/∂)y

{ex, 
          ey}  =  -{D[elecPotential, x], 
            D[elecPotential, y]};

If the field line is thought of as the trajectory of a massless particle under the influence of the electric field, then the differential equations to be solved to plot this trajectory are:

(dx/dt) = (Subscript[E, x ])
(dy/dt) = (Subscript[E, y])

FieldLine[{x,ex,x0},{y,ey,y0},{t,t1}] calculates the field lines given the differential equations ex and ey for variables x and y, with starting values x0 and y0, for t up to t1.  It does this using NDSolve, the output is in the form Line[list of line coordinates for a single line].

FieldLine[{x_, ex_, x0_}, {y_, ey_, y0_}, {t_, t1_}, opts___]:=Module[{exf, eyf, sol, t2},
{exf, 
                eyf} = {ex, ey} /. {x -> x[t], y -> y[t]};
sol = 
              NDSolve[{x'[t] == exf, y'[t] == eyf, 
                  x[0] == x0, y[0] == y0}, {x[t], y[t]}, {t, 0, 
                  t1}, MaxStepFraction -> 1/1000, opts];
sol = {x[t], y[t]} /. First[sol];  
If[VectorQ[sol /. t -> 0, NumberQ]  ,
t2 = Part[sol, 1, 0, 1, 1, 2]   ;
sol = 
                    ParametricPlot[Evaluate[sol], {t, 0, t2}, 
                      DisplayFunction -> Identity];
First[
                    Cases[sol, Line[_], Infinity]]   ,
Line[{{x0, y0}}]   
]
]

Create a stopping test for use in NDSolve so that once the field lines are near a charge the calculation stops.

sTest[cs_, 
        pn_] := (Sqrt[
            Min[# . # & /@  ((cs - #) & /@ 
                    pn)]] < Δ)

Δ = 0.07;

Set m to be the number of points around the charge at which to start calculating field lines from - this does not mean that there will be only m field lines emanating from each charge, just that the number of field lines whose starting point is around an individual charge is m.

m = 11;

linesPos creates a list of field lines whose initial points are around the positive charges.

linesPos = 
        Flatten[Table[
            FieldLine[{x, ex, 
                positiveCharges[[i, 2, 1]] + 0.07  Sin[tj]}, {y, ey, 
                positiveCharges[[i, 2, 2]] + 0.07  Cos[tj]}, {t, 200}, 
              StoppingTest :> 
                sTest[{x[t], y[t]}, negativeCharges[[All, 2]]]], {tj, 0, 
              2  Pi - ((2 Pi)/m), ((2 Pi)/m)}, {i, 1, 
              Length[positiveCharges]}]];

linesNeg creates a list of field lines whose initial points are around the negative charges.  Field lines generated from these points move up field, towards the positive charges.

linesNeg = 
        Flatten[Table[
            FieldLine[{x, -ex, 
                negativeCharges[[i, 2, 1]] + 0.07  Sin[tj]}, {y, -ey, 
                negativeCharges[[i, 2, 2]] + 0.07  Cos[tj]}, {t, 200}, 
              StoppingTest :> 
                sTest[{x[t], y[t]}, positiveCharges[[All, 2]]]], {tj, 0, 
              2  Pi - ((2 Pi)/m), ((2 Pi)/m)}, {i, 1, 
              Length[negativeCharges]}]];

Adding  arrows to field lines

In order to show the direction of the field, it is useful to place arrows along the field lines.

<< "Graphics`Arrow`"

AddArrow[Line[pts],d,num] accepts input of the form Line[list of coordinates along single line] and places arrows along the line pointing in the direction given by the coordinates of the line - therefore when putting arrows on the field lines from the negative charges later it is necessary to reverse the order of the coordinates of the field lines.

AddArrow[Line[pts_], d_, 
        num_:  5] := 
Module[{arrow = {}, 
          n = 0}, 
Fold[
If[
              First[#1] ≥  d && 
                n < num, 
n++; 
AppendTo[arrow, 
                Arrow[Last[#1], #2, HeadScaling -> Absolute, 
                  HeadCenter -> 0.5, 
                  HeadLength -> 
                    4]]; 
{0, #2}, 
{First[#1] + 
                  Sqrt[Apply[Plus, (((Last[#1] - #2))^2)]], #2}] &, {0, 
            First[pts]}, Rest[pts]]; 
arrow]

arrowsPos = Map[AddArrow[#,  .5] &, linesPos]  (* 
          puts arrows on the field lines from the positive charges *) ;

In order for the arrows on the field lines running from the negative charges to the positive charges to be pointing in the right direction, must reverse the coordinates of the line.  The function rev[Line[pts]] requires input of the form Line[list of coordinates on the line] and reverses the list of coordinates.

rev[Line[pts_]] := 
        Module[{line1}, 
line1 = 
            Reverse[pts]; 
line1 = 
            Line[line1]];

arrowsNeg = Map[AddArrow[#, 1] &, rev /@ linesNeg]   ;

Plot field lines

Finally all these graphics objects are combined to show the field lines throughout the region, along with the charges.

c3 is a contour plot of the region, which is then labelled using the LabelContourLines function from the Extend Graphics book.  This looks a bit cluttered, may have to change.

c3 = 
        ContourPlot[elecPotential, {x, xmin, xmax}, {y, ymin, ymax}, 
          ContourShading -> False, ContourSmoothing -> Automatic, 
          PlotRange -> {-2, 2}, Contours -> 7, 
          PlotPoints -> 50, 
          ContourStyle -> {RGBColor[0, 0, 1], Dashing[{0.01, 0.01}]}, 
          DisplayFunction -> Identity];

c3 = 
        LabelContourLines[c3, LabelPlacement -> Automatic, 
          PointFactor -> 0.1, 
          DisplayFunction -> Identity];

Then show the field lines along with the contour plot.

pic = 
        Show[c3, 
          Graphics[{Thickness[0.001], arrowsPos, linesPos, arrowsNeg, 
              linesNeg}], AspectRatio -> Automatic, Frame -> True, 
          Epilog -> chargePlot, ImageSize -> {600, 600}, 
          PlotRange -> {{xmin, xmax}, {ymin, ymax}}, 
          DisplayFunction -> $DisplayFunction];

Cleaning Up

ClearAll["Global`*"]  (*removes everything of context Global*)