Category Archives: Advanced EM

Advanced Electromagentism (Phys 341)

Precession of Mercurian Planets: Project Plan


For my computational physics project I would like to investigate the relationship between the precession and the eccentricity of a planet’s orbit due to general relativity as detailed in chapter 4, section 3 of Computational Physics, by Nicholas J. Giordano and Hisao Nakanishi. The rate of precession examined in this project will only be that which is caused by the general relativistic effect of the model planet’s host star.

Week 1: Preliminary Results

As a check to make sure my code will run correctly, I will first model the precession of Mercury’s orbit as caused by the Sun. I must first write a code that will model the movement of a Mercury through it’s orbit according to Kepler’s laws. The force law predicted by General Relativity is given by:


where G is the Gravitational constant of the universe, Ms is the mass of the Sun, Mm is the mass of Mercury, and a Constant*α gives the rate of precession for the orbit. For the case of Mercury,  \alpha=1.1*10^{-8} AU^{2} .

Applying the Euler-Chromer method to the above equation to approximate the x and y velocities and positions of the planet as it travels throughout its orbit yields:

 V_{x,i+1}= V_{x,i}-(4*(\pi^{2})*x_{i})/(r_{i})^{3}+(\alpha*4*(\pi^2)*x_{i})/(r{i})^{5})*\Delta(t)
 x_{i+1} = x_{i}+(V_{x,i+1})
 V_{y,i+1}= V_{y,i}-(4*(\pi^2)*y_{i})/(r_{i})^{3}+(\alpha*4*(\pi^{2})*y_{I})/(r_{I})^{5})*\Delta(t)
 y_{i+1} = y_{i}+(V_{y,i+1})

Where the initial conditions are (assuming  r_{1} to be at the planet’s aphelion):

 v_{1} = \sqrt(\frac{G*M_{s}*(1-e)}{a*(1+e)}
 r_{1} = (1+e)*a

where a is the semi-major axis length of the planet’s orbit.

In Computational Physics, Giordano and Nakanishi find the precession of Mercury by first finding the precession rate for a much larger value of  \alpha , because it is difficult to measure the precession for a value this small (it would require computing the motion for hundreds of thousands of years). Using different values of α give different values for the precession rate, and graphing these two parameters against one another will provide the desired constant C. Once C is obtained, I can multiply this number by α to get the true precession of Mercury caused by the Sun.

Week 2: Results with Data Analysis

Once this code is working correctly, I can modify the initial conditions for the orbital motion portion of the code to represent new orbits. (Changing the parameters a, e, but maintaining the value of the perihelion as that of Mercury, so that I may compare my results to that for Mercury.)

If I have extra time I will add in other planets to the original Mercury-Sun system to observe the resulting effects on Mercury’s precession.

Week 3: Finalize project: make sure the code’s comments are informative and helpful, and review the program/results for any mistakes.

Week 4: Provide constructive criticism to other students on their projects/Prepare for project presentation.

Week 5: Give my project presentation to the class.

Computational Physics by Nicholas J. Giordano and Hisao Nakanishi



Experimental vs. Theoretical Values

It was immediately apparent when I began collecting data for this experiment that the equations I derived to anticipate voltage and frequency values were not entirely accurate. Here I will calculate expected curve shapes and expected values of frequency and decay, and and compare them to experimentally measured values. The majority of the difference in theoretical vs measured data is the fact that there is some resistance in the circuit. It may be possible to discover this resistance from the future tests, if not from the current data.

Shape Of Potential Curve

1465 nF Capacitor

Figure 1. 1465 nF Capacitor and 996 mH Inductor in series.

In theory, the circuit I constructed was a LC circuit. This means there was no resistor in series with the capacitor and inductor components, and that the wires have negligible resistance. Theoretically the oscillations of the voltage after the capacitor discharged through the inductor should continue indefinitely. This is analogous to how an object on a spring would oscillate forever without drag/friction forces acting upon it. That being said, it is readily apparent that the oscillations did not continue indefinitely. In fact they took around a 10th of a second to dissipate (in the case of the 1465 nF capacitor and 996 mH inductor circuit above). This means that there is damping occurring on the oscillations. So while theoretically I built an LC circuit, experimentally it behaved more like a damped LC circuit, or an RLC circuit, as if there was a resistor in series with the capacitor and inductor. This means that the relevant governing equation was not V\left( t\right)=V_{0}sin \left( \omega _{0}t-\delta \right)\left , but rather V\left( t\right) =V_{0}e^{-\beta t}\sin \left( \omega _{1}t-\delta \right) . Below (Figure 2) is the demonstration of expected shape (red) and measured shape (blue) of the voltage curve vs time as seen in my “Preliminary Data” post:

1465nF, 996mH, 128 Hz

Figure 2


To review,

 f = \frac{\omega}{2\pi}         (1)

In Simple Harmonic Motion

 f = \frac{\omega_{0}}{2\pi} ; \omega_{0}=\frac{1}{\sqrt{LC}}.         (2)

But as we just discussed, because the circuit is not ideal and has an equivalent resistance, these oscillations are actually described by Damped SHM, the equations being as follows:

V\left( t\right) =V_{0}e^{-\beta t}\sin \left( \omega _{1}t-\delta \right)        (3)

\omega _{1}\equiv \sqrt {\omega _{0}^{2}-\beta ^{2}}        \beta \equiv \frac{R}{2L}       (4)


 f= \frac{\sqrt{\frac{1}{LC} - (\frac{R}{2L})^2}}{2\pi}        (5)

Where, solving for R, we can estimate the equivalent resistance of the circuit:

 R=\sqrt{\frac{4L}{C}-(4 \pi L f)^2}          (6)

 Below is a table of values (Figure 3). On the left are all of the the experimentally measured values from this project. The expected theoretical frequency was calculated using Equation (2). The resistance of the circuit was estimated using Equation (6).


Figure 3

While the measured frequencies are not always very close to the predicted value, they are all well within an order of magnitude of each other, meaning the relationship defined by Equation (2) is clearly evident.

The estimated values of resistance, on the other hand, are extremely improbable. The resistivity of small connection wires is on the order of 10^-2 Ohms or lower. This indicates either a mistake in my derivation of resistance as a function of capacitance, inductance, and frequency, or a relationship that I don’t quite understand. It is possible that the equivalent resistance of the circuit changes as a function of voltage, or the change in current, in a way I do not know about.

Overall, the graphs of frequency as a function of L or C (Figure 4) were the correct shape.

Figure 5

Figure 4

They reflect the predicted behavior of a  \frac{1}{\sqrt{x}} relationship, whose plot looks very similar (Figure 5). In our case the actual relationship is  f = \frac{1}{2 \pi \sqrt{x}}.

Screen Shot 2014-05-20 at 9.01.56 PM

Figure 5 (Mathematica workbook: )

Exponential Decay and Circuit Resistance

When the experimentally measured decay coefficient β is plotted as a function of inductance, the data points are sparse, and one should be cautious of drawing too many conclusions. That being said, the points I collected (Figure 6) do seem to match the expected values (Figure 7):


Figure 13

Figure 6


Figure 7

Finding Resistance Using β

The decay coefficient β is related to L and R. Thus, since we know β and L experimentally, it should be fairly straightforward to calculate R empirically. In theory the resistance of the circuit should be the same for each measurement of β and L. Below is a table of measured values of β and L (Figure 8):

Screen Shot 2014-05-20 at 9.23.11 PM

Figure 8

A value of 0.0381 Ω is a reasonable order of magnitude for a circuit with theoretically no inherent resistance. This order of magnitude is also large enough to be largely responsible for the decay of the voltage function across the capacitor.

Can Significant Conclusions Be Drawn From This Data?

With scatter plots with only three data points, the answer is “no.” To confidently confirm that my data was following anticipated trends I would need many, many more points. The most conclusive and illuminating data that I did generate wast the Voltage vs. Time graphs. These graphs were supposed to be undamped sinusoidal waves, but instead exhibited very clear damped oscillation. The collection of these curves was easy to repeat and allowed me to generate very reliable curves. These graphs confirmed that there must indeed be a significant equivalent resistance in the circuit to damp the oscillations.

Applications Of Findings

This experiment could potentially be an effective way to anticipate the inherent resistance of a circuit. My first method of using L, C, and f to find the resistance, clearly did not work, but finding the decay coefficient and from there calculating R, could very well be effective. From my limited data I got a reasonable solution, but would have to generate far more data points to confirm this resistance. It would be particularly interesting to use this method in an instructional lab if there was some way to actually measure the resistance of the circuit (maybe known values of wire resistivity, or a very sensitive meter). It would provide students with data analysis experience, such as fitting curves and using Origin Pro, as well as having hands on experience with damped harmonic systems.

Changes To Experiment

The biggest change I would make to this experiment would be to have a far greater array of capacitors and inductors at my disposal. I did not have enough to produce satisfying data. One problem I kept running into when collecting data was that one of my inductors was several orders of magnitude bigger than the others. It could have been handy to have just made my own inductors out of coils of wire. This way I would be able to change the inductance at will. Ideally I would also exchange many capacitors for one variable capacitor. It would streamline the data collection and allow for more methodical, precise data collection. I would also use components and measurement devices that could handle more than 15-20V due to the fact that such low voltage has a greater window for interference.

Sources Of Error

I have no idea how much my measuring devices and power source effect the overall resistance of the circuit. I also do not know if there are other sources of induction (albeit minor) within my circuit setup, either within components (such as the switch) or from the wires of my measurement tools coiling around themselves. There were many times when the discharge curve across the capacitor would take a strange shape that I can only attribute to some form of saturation, either the probe or capacitor. I do not know how this my have impacted my data.


Final Data

In this post I will present my actual data findings with minimal interpretation. In the next blog post, “Conclusion,” I will interpret the data I collected and compare it to  expected values.

Preliminary Data Recap:

In my previous post I had measured the discharge of several different capacitor values in series with an inductor.

Below were the resulting plots of Voltage vs Time (Figures 1-3):

1465 nF Capacitor

Figure 1

1007 nF Capacitor

Figure 2

47 nF Capacitor

Figure 3

For the 1465, 1007, and 47 nF capacitors the measured frequency of oscillation was 128, 155, and 730 Hz, respectively.

Voltage As a Function Of Time  and Frequency of Oscillation

It is helpful to layer these three plots on top of each other (Figure 4). In doing so we can readily compare the curves to each other.

Figure 4

Figure 4

From Figure 4 it is evident that the higher the capacitor value, the lower the frequency of oscillation. This is not surprising, given the equation that the frequency and inductor*capacitor product are inversely related: \omega_{0} = \frac{1}{\sqrt{LC}}.

How are frequency and capacitance, or frequency and inductance related? I recorded the frequency of oscillation in several different LC circuits to try and illustrate this relationship. For three different inductor values I discharged five different capacitor values. The resulting plot (Figure 5) is below:

Figure 5

Figure 5

Because it is harder to see the shape of the blue line in Figure 5, I plotted it separately (Figure 6):

Figure 6

Figure 6

Figures 5 and 6 experimentally illustrate the inverse square root relationship between frequency and capacitance/inductance.

I then plotted the frequency of oscillation with respect to inductor value (Figure 7):

Hz vs uH

Figure 7

 Exponential Rate Of Decay

If you take a look at Figure 4 again, notice that despite changing capacitor values, the exponential decay of each curve is relatively the same. The decay for each curve is the same because the decay coefficient β is a function of resistance and inductance, not capacitance(\beta\equiv \frac{R}{2L}). In Equation (1) you can see the decay term of the voltage, e^(-βt):

V(t)=V_0e^{-\beta t}cos(\omega_1t-\delta)           (1)

In order to measure the decay constant I used Origin Pro to fit an exponential curve to the peaks of the voltage plots. The formula used for this function fit was as follows:


Where R0 is the decay coefficient, -β

Figure 8 below is an example of this function fit:

996 mH exp fit

Figure 8

I fit an exponential function to the 996 mH and 1007 nF, 47 nF oscillations as well, resulting in the table below (Figure 9):

 beta values

Figure 9

Theoretically the β value should be the same for all of them.

In order to explore how the decay coefficient β changes I needed to change the inductor values. Because the voltage curves were more manageable at greater LC values, I used the largest capacitor I had and varied the inductor value.  The results were as follows (Figures 10-11):

37 uH exp fit

Figure 10

171 uH exp fit

Figure 11

The resulting decay coefficients are as follows (Figure 12):

Screen Shot 2014-05-20 at 3.59.39 PM

Figure 12

Which can also be represented in as a scatter plot, despite having so few points(Figure 13):

Figure 13

Figure 13








Critique: Modeling Electric and Magnetic Fields

At the end of the preliminary data section, the Electric Field of a bar electret is stated to be approximately that of two opposing charged point charges at either end of the bar. However, little explanation is given as to why the bar can be modeled that way. If there was a justification given for this model, such as modeling many charges arranged in a line to show that the net fields cancel aside from the two charges at either end, the data would significantly benefit from this. Finally, it is unclear in the final data whether or not you are modeling a bar electret or an electric dipole.

Additionally, it is unclear why you chose the specifications for the dipole bar magnet, as well as the surrounding constants and values for the other scenarios. Although it is clear from a knowledgeable individual’s standpoint why the charges were assigned the equations as they were, a layperson may not be able to interpret the equations correctly.  Additionally, including the full derivations for each scenario on the preliminary data in mathematica would streamline the analysis of the data significantly, as the constants themselves would become more clear with this. Additionally, the equations that yield the specific graphs in the final data you are showing should be wrote, so that I we can associate the data and figures easily.

Finally, I suggest making the equations and derivations slightly more robust. By that I mean carrying out the derivations for these geometries in matter, as to better replicate the situations which a student may come across. Currently, this issue is approached without giving thought to the effects of polarization and the various materials involved. This may prove to both improve the depth and breadth of the project. Another concern that I have, although it may be a fault with mathematica, is that the field of the electret, does not show the field’s shape clearly. I suggest using a different graphing method, perhaps a contour graph, which would allow the interpretation of field lines much more intuitively. As such we would then expect the graph to take the shape of an imperfect dipole, as stated in your data.


Preliminary Data

Measuring and Modeling Physical RCL Circuits

Extended Derivations of LC and RLC circuit behavior.

EM circuit DC

To derive the equations governing the voltage across the capacitor I start with Kirchoff’s Loop Law:

\Delta V_{L}+\Delta V_{C}=0           (1)

Substituting in for the values of voltage across a capacitor and inductor we get:

 L\dfrac {dI\left( t\right) } {dt}+\dfrac {Q\left( t\right) } {C} =0         (2)

The current can be expressed as the rate of change of charge on the capacitor:

 I\left( t\right) =\dfrac {dQ\left( t\right) } {dt}         (3)

 L\dfrac {d^{2}Q\left( t\right) } {dt^{2}}+\dfrac {Q\left( t\right) } {c}=0         (4)

This expression is analogous to the equations governing Simple Harmonic Motion (SHM):

 F_{spring}=-kx           (5)

F_{net}=ma= m\dfrac {d^{2}x\left( t\right) } {dt^{2}}         (6)

m\dfrac {d^{2}x\left( t\right) } {dt^{2}}+kx =0           (7)

The solution to this second order differential equation can be written as.

x\left( t\right) =x_{0}\sin \left( \omega _{0}t-\delta \right)           (8)

Similarly the general solution for charge across the capacitor can be given by:

Q\left( t\right) =Q_{0}\sin \left( \omega _{0}t-\delta \right)           (9)

Thus the Voltage across the capacitor, \Delta V=\dfrac {Q\left( t\right) } {C}, can be written as:

\Delta V\left( t\right)=\dfrac {Q_{0}} {C}sin \left( \omega _{0}t-\delta \right)\left                   (10)

\Delta V\left( t\right)=V_{0}sin \left( \omega _{0}t-\delta \right)\left                   (11)


RLC circuit (damped oscillation)

EM circuit DC R


The RCL circuit can be similarly derived by adding a third term for the voltage across the resistor to Kirchoff’s Loop Law.

\Delta V_{L}+\Delta V_{R}+\Delta V_{C}=0       (12)

 L\dfrac {d^{2}Q\left( t\right) } {dt^{2}}+R\dfrac {dQ\left( t\right) } {dt}+\dfrac {Q\left( t\right) } {c}=0      (13)

\dfrac {d^{2}Q\left( t\right) } {dt^{2}}+\dfrac {R} {L}\dfrac {dQ\left( t\right) } {dt}+\dfrac {1} {LC}Q\left( t\right) =0      (14)

Here we define

\dfrac {1} {LC}\equiv \omega _{0}^{2}    and    \dfrac {R} {L}\equiv 2\beta       (15)

To simplify our second order differential equation to:

\dfrac {d^{2}Q\left( t\right) } {dt^{2}}+2\beta \dfrac {dQ\left( t\right) } {dt}+\omega _{0}^{2}Q\left( t\right) =0      (16)

Again, this equation is analogous to Damped Simple Harmonic Motion, described below:

 F_{spring}=-kx -b\dfrac {dx\left( t\right) } {dt}         (17)

F_{net}=ma= m\dfrac {d^{2}x\left( t\right) } {dt^{2}}       (18)

m\dfrac {d^{2}x\left( t\right) } {dt^{2}}+kx -b\dfrac {dx\left( t\right) } {dt}  =0       (19)

The general solution to this differential equation, describing an underdamped system (where \beta  < \omega _{0}), is as follows:

\omega _{1}\equiv \sqrt {\omega _{0}^{2}-\beta ^{2}}       (20)

V\left( t\right) =V_{0}e^{-\beta t}\sin \left( \omega _{1}t-\delta \right)       (21)

This solution describes a sinusoidal function whose amplitude decays as a function of e^{-\beta t}.

While in theory my experimental setup is an LC circuit, and should be described by SHM, in reality there should be some form of resistance in the circuit, causing the oscillations to decay over time as described in Damped SHM.

Experimental Setup:


The above diagram represent the physical circuit I built to observe the voltage oscillations of an LC circuit. This circuit has two functions. The first, when the switch connects the loop to the left, charges the capacitor with the voltage source. The second, when the switch completes the right loop, discharges the capacitor through the inductor. The right and left loops are independent of each other because of the switch. Connected in parallel across the capacitor is the PicoScope. This device records the voltage across the capacitor over time. Below are photographs of what experimental setup actually looked like:

DC power source, signal generator, oscilloscope, circuit (and multimeter)

DC power source, signal generator, oscilloscope, circuit (and multimeter).

Circuit Setup

Circuit setup. LC circuit with DC power source, switch, capacitor, inductor, and oscilloscope.


Voltage difference and reference signal

Oscilloscope with Channel A: Voltage difference, Channel B: reference signal.

First Data:

This project held a few challenges of its own before I even got to collecting data. Completely unfamiliar with the Picoscope (digital oscilloscope) probe and software, it took me a few days to figure out how to efficiently collect data with the setup.

My first attempt at taking data did not go well. With the oscilloscope I could see the voltage difference across the capacitor that I wanted, but when I switched the circuit to discharge the capacitor, the voltage dropped to zero immediately. I was looking for some sort of sinusoidal curve, possibly decaying very quickly (in ideal conditions the curve would never decay).

Searching for answers as to why I wasn’t seeing anything, attempted to estimate how long any capacitor would take to discharge.

The charge on a capacitor as a function of time (in an RC circuit) is given by Q=Q_{0}e^{-tk}, where k=\dfrac {1} {RC}. At the time I was using components of unknown value. I visited Larry Doe’s workshop in Blodget and measured the resistance, capacitance, and inductance values of all my components. The capacitor I had picked for my initial measurements had a value of 0.1 μF. What this means is the charge on the capacitor had discharged e^{-1} of it’s initial value (about 37% of it’s initial value) at t=RC . In this scenario the R is the resistance of the circuit. Using an estimation of 0.0001 Ω for 5 cm copper wire as the equivalent resistance of the circuit, and 10^-7 F as the capacitance, the charge on the capacitor should reach one third it’s initial value in 10 ps. To me this illustrated the incredibly short time period at which that particular LC combination would oscillate, and potentially decay.

Second Attempt:

To maximize the time it took the system to discharge I used the highest valued capacitor and inductor I had. The resulting curve was far more manageable. Below is the data I collected with the oscilloscope. The blue curve you see is voltage across the capacitor as a function of time. The total time shown for each image is 100 ms. The red curve is an artificially generated curve for reference. It represents the expected shape of the blue curve. In an ideal LC circuit there is no resistor, and thus no decay over time. Theoretically an ideal LC circuit would oscillate indefinitely.

Figure 1. Capacitor: 1465 nF   Inductor: 996 mH

1465nF, 996mH, 128 Hz

1465nF, 996mH, Oscillation: 128 Hz

Figure 1.1 With reference wave

1465nF, 996mH, 128 Hz

Oscillation: 128 Hz

I then used successively lower and lower capacitor values to approach my initial test where the capacitor discharge was far too quick to measure. The duration of measurement, initial voltage across the capacitor, and inductor value were held constant.

Figure 2. Capacitor: 1007 nF   Inductor: 996 mH

1007nF, 996mH, 155Hz

1007nF, 996mH, Oscillation: 155 Hz

Figure 2.1 With reference wave

1007nF, 996mH, 155Hz

Oscillation: 155 Hz

Figure 3. Capacitor: 47 nF   Inductor: 996 mH

47nF, 996mH, 730 hZ

47nF, 996mH, Oscillation: 730 hZ

Figure 3.1 With reference wave

47nF, 996mH, 730 hZ

Oscillation: 730 hZ

Figure 4. Capacitor: 1.9 nF   Inductor: 996 mH

1.9nF, 996mH, 4.23 kHz

1.9nF, 996mH, Oscillation: 4.23 kHz

Figure 4.1 With reference wave

1.9nF, 996mH,4.23 kHz

Oscillation: 4.23 kHz



Conclusion: Electric Fields of Spherical Objects

Gauss’s Law is a powerful tool for understanding electric fields for all types of configurations. The most general case of Gauss’s Law is given by $ \oint \! \textbf{E} \cdot \mathrm{d} \textbf{a} = \frac{1}{\epsilon_0} Q_{enc} $. In this equation, $ \textbf{E} $ represents the electric field vector, $ \mathrm{d} \textbf{a} $ represents the differential area vector of the Gaussian surface through which the electric field is pointed, $ \epsilon_0 $ represents the permittivity of free space, and $Q_{enc}$ represents the charge enclosed by a Gaussian surface.

As stated in my project plan, the general case can be used to solve for the electric field of a point charge, which is given by $\textbf{E} = \frac{1}{4 \pi \epsilon_0} \frac{q}{r^2} \mathbf{\hat{r}} $. The unit vector indicates that the electric field for a point charge points in the radial direction. In my project, I managed to use this equation to model the electric fields for both positive and negative point charges.

I also worked to model the electric field for hollow spherical conductors. The above equation still holds, but with one caveat:

    \begin{displaymath} \textbf{E} = \left\{ \begin{array}{lr} 0 & : r \leq R\\ \frac{1}{4 \pi \epsilon_0} \frac{q}{r^2} \mathbf{\hat{r}} & : r > R \end{array} \right. \end{displaymath}

In the above, $R$ represents the radius of the sphere. This piecewise function describes the fact represented in Gauss’s Law that a electric fields only exist in regions in which a charge is enclosed. In conductors, all of the charge exists at the outer surface, so there is no electric field on the inside of a hollow spherical conductor. When the distance away from the conductor is greater than its radius, it obeys the second piece of the piecewise function.

This image shows the electric field produced by a positive point charge in free space.

This image shows the electric field produced by a positive point charge in free space.

My project was initially intended to be a project that produced visualizations of the above equations. From that point, I intended to develop more complex systems, work out the math for them, and then visually model them as well. However, due to my inexperience with Mathematica 9, I encountered many issues. The RegionFunction proved to be my saving grace, as described in my results post, but it happened too late for me to be able to model more complex systems.


This image shows the electric field produced by a negative point charge. Compare to Figure 1.

This image shows the electric field produced by a negative point charge in free space. Compare to the figure for a positive point charge.

My results do show good visual models for electric fields of single positive and negative point charges, and a positively charged spherical conductor. These visual models are incredibly useful tools for understand electric fields. Looking at an equation (Gauss’s Law) can only do so much for one’s understanding; a visual strongly aids this.

This shows the electric field produced by a spherical conductor with positive charge on the surface.

This shows the electric field produced in free space by a spherical conductor with positive charge on the surface.

I wish that Mathematica had even more tools for visually modeling these electric fields. I tried to create an animation that showed the graphing of the vector field arrows coming out of the point charge, but this did not work. I don’t believe Mathematica’s VectorFunction3D  works with its Animate function, but it may just require different inputs. I also tried to create a Manipulate object to demonstrate changes in the electric field as charge changed, but this was also not successful. Part of me believe this is an issue of scaling the image, and part of my believes that this function is also not compatible with VectorFunction3D.

In the future, this project could be expanded by modeling more complex systems. One way to complicate this would be to introduced additional spherical conducting shells with charge on them. Another complicating factor would be to put a point charge inside of a spherical shell. These systems could be modeled with differing signs and values of charges. Additionally, an infinite charged plane could be added to show other ways the electric field might change.

Additional investigations could also include modeling the electric fields of dielectrics instead of conductors. This analysis would be especially complex because it would include polarizations and bound charges. Additionally, this could include the modeling of the electric displacement in addition to the electric field.

Sources and Resources:

  1. Introduction to Electrodynamics, David Griffiths, 4th Edition
  2. Mathematica 9: Student Edition, created by Wolfram
  3. Collaborators: Brian Deer and Tewa Kpulun

Van Allen Belts Modeling – Conclusion and Final Thoughts

At the conclusion of this project I have to say I came out of this project having learned a lot about organizing a project, the Van Allen Belts and how to work with Mathematica. But first I would like to discuss the conclusions that can be drawn from my final data.

Unfortunately I can’t draw any major conclusions about the Van Allen Belts and their location depending on core size as I was not able to incorporate that particular information into my models. I was however successfully able to model the vector fields of a magnetic dipole using Mathematica and an initial equation given by our textbook, a task that proved to be difficult enough on its own to figure out. In the end I did prove that the equation was valid and created a vector field that resembled the expected field of the dipole. Even with this however, I question the validity of the models that were obtained using Mathematica. As I mentioned in my “Final Data” post, the vectors do not seem to vary with respect to the radius away from the origin which doesn’t match the actual system I was attempting to model. It should be decreasing with respect to distance. In the end the project was never really able to combine the Van Allen Belt mapping and the modeling of Earth’s magnetic field into the single cohesive unit.

Furthermore the model itself had a lot of assumptions that needed to be put into place in order to make it easier to develop or simply because I could not figure out how to account for it within Mathematica. The assumptions were listed throughout my final data but they included using a constant magnetic field isolated from the sun’s solar wind to model out the Van Allen belts and assuming the Earth to be a sphere. In the long term however, these were all reasonable assumptions in astronomical scales.

If there was one thing I was somewhat disappointed about, it was the learning curve to reacquaint myself with Mathematica. It had been awhile since I had used Mathematica and it took some time to recall the best functions to use for particular tasks. Furthermore when trying to learn how to do a new task; like converting between Spherical and Cartesian coordinates with Mathematica, it was a long process of trial and error to achieve a result you were looking for. More often than not you would get something more like this:

No Mathematica...this isn't rotating my graph...good try though.

Occasionally, the program has a cruel sense of humor

So a great deal of the project was spent playing around with Mathematica more that playing around with the physics. I do certainly wish I had a bit more time to play around with the physics and even more time to play with the Mathematica program. After a few breakthroughs, I was quickly understanding how to make the program much more cooperative and finding relevant information in the help menus faster. With a little more time, I fixed and accounted for a few of the issues that had been plaguing my models.

Regardless, I was a bit overzealous with the scope of my project. It had a lot of components that were much more difficult than I initially expected them to be. What I ended with were two separate projects, a research summary of Van Allen Belts and a 3D model of Earth’s magnetic field with vague notions of where the Van Allen Belts would be located. Still I am proud of what was accomplished. Modeling the vector field of a dipole was a difficult task and a properly organized model was the end. It was also plotted with proper constants though it needed a little tweaking to confirm that the structure was as expected. Finally from this point it could easily be picked up and completed at a later date with those willing to work on it needing only to manipulate the dynamo equation into an easier form.



Final Remarks

In this post, I will discuss the results of the final data, what I learned, and suggestions for further work.

I began by using one of the most fundamental and well known laws in electromagnetism: Gauss’s Law. It allows for the relatively simple modeling of electric fields of various geometries, where the direct formula for the electric field would have been difficult, if not impossible, to use.

I used this law to derive some the electric fields of some more common physical geometries dealt with in undergraduate physics, which the exception of one configuration: the electric dipole. Through this process I not only re-affirmed my knowledge of the principals involved but learned how to communicate them in an effective way. I then used Mathematica to model the derived electric fields to provide a visual representation.

This was particularly challenging as I do not had a significant background with computer programing at all, let alone with Mathematica. From this work I significantly improved my working knowledge and understanding of this software, not only through the programing but through explaining what was being done throughout the code.

The final results displayed how the electric field of the particular configurations of a sphere, cylinder and ideal dipole work and these results agreed with known models. The basic principles showed that electric fields point away from positively charged objects and towards negatively charged ones, which is most clearly seen in the model of the dipole.

Further work with this topic ideally would include the modeling of more complex systems, such as multiple spheres or cylinders of varying charge or charge densities. It would also be interesting to compare the magnetic field configurations of similar objects, which unfortunately was not able to be done as originally planned due to limitations of Mathematica and available working material.


Van Allen Belt Modeling – Final Data

Background Information:

Earth’s Van Allen Belts are a naturally occurring phenomena due to the presence of the Earth’s magnetic fields. Specifically, they are two zones of space where large concentrations of high energy particles are trapped by the magnetic field lines of our planet. The zones are defined by our field lines; toroidal in shape, and encompass the planet at particular radii from the surface. Following the field lines, the regions begin from Earth’s Southern Pole and loop to its Northern Pole  This is illustrated in the diagram below.

Cross-Section of the Van Allen Belts

Cross-Section of the Van Allen Belts

The particles of these regions include fast moving electrons,and ions from the solar wind along with protons and other remnants of cosmic ray interactions with our atmosphere. All of these particles are electrically charged and thus their motion becomes highly constrained as they pass through our magnetic field; or magnetosphere. As they interact with our magnetic field lines, a force is exerted and influences their direction of motion. It attempts to make the particles rotate about the lines rather than allowing them to continue drifting in a straight line. Thus “sliding” begins to occur; the particles gain a spiral trajectory during which they orbit the field lines while vertically moving along them. Since the magnetic field is stronger as you get closer to Earth however the particles will stop and eventually reverse their motion as they approach the poles. This bouncing between the hemispheres establishes the permanent belts.

The constituent particles and stability of the belts however depends on the location.and thus the strength of the magnetic field. The inner belt is compact; ranging only from 1,700 km to 13,000 km above the Earth’s surface. It is made up primarily of high energy protons created during cosmic ray interactions with the atmosphere. Of the two belts, it is the most stable due to the higher magnitude of the magnetic field lines at that proximity to the Earth.  While it gains high energy particles at a comparatively slow rate, it holds onto the particles for an extended period of time and reaches a much higher intensity than the other belt.

The outer belt ranges from 20,000 km to 40,000 km from the Earth’s surface and consists primarily of the ions and electrons that make up the solar wind. Due to the weaker magnetic field at that distance from Earth, the zone is tenuous and loses its particles very quickly. The constant supply of particles from the solar wind however means that the belt has plenty of material to be created from, though its dependence on solar wind particles means that its intensity fluctuates proportionally with solar wind activity. With this background information established, the information needed to created valid assumptions was gained.

2D Vector Fields:

Unfortunately, this also marks the point where my project strayed away from modeling the Van Allen Belts to a more general modeling of the Earth’s magnetic field. The initial task of modeling the magnetic field of the Earth was necessary to define the locations of the Belts with specific field lines that would be seen in the vector plot. Now with some research I found that the Earth’s field could be generally modeled as a magnetic dipole oriented with the south end pointed in the positive z direction. This model is a very basic assumption and ignores the effects of the solar wind on Earth’s magnetosphere, however the area where the Van Allen Belts are present does not get severely effected by these factors except in extreme circumstances. This means a model of a dipole will be generally accurate for our purposes. Thus the first task was to create a model of the magnetic field of the dipole.

The equation for the magnetic field of a dipole with respect to the radius was found to be

    \[ \vec{B_{dip}}(\vec{r})=\frac{(\mu_0)(m)}{(4\pi*r^3)}(2cos(\theta)\hat{r}+sin(\theta)\hat{\theta})\]

This equation was obtained from Griffiths’ Introduction to Electrodynamics. To insure the equation worked I began by applying a constant free version, evaluating only the variables. However, the equation was given in spherical coordinates which Mathematica can’t plot efficiently. To compound that issue, the command CoordinateTransform does not properly evaluate the points which led to several attempts to obtain the proper Cartesian conversion of our equation. Finally, I was able to find the command TransformedField which took into account the variables being used and allowed for a direct conversion. The code and direct conversion can be seen below.


From these equations, I created a vector field of the xz plane with VectorPlot which can be seen below.

xz plane of dipole magnetic field

Simple vector model of the xz plane of the magnetic field of a dipole

With this successful, I then applied the proper constants to the vector field components and attempted to recreated the field. These constants were the permeability of free space ( $\mu_0$ = $4\pi$*$10^{-7}$ $N$/$A^2$) and the magnetic dipole moment of the Earth ($m_{Earth}$ = $7.79$*$10^{22}$ $A*m^{2}$).The circle in the center of the graph represents the Earth, under the assumption that the Earth has a constant radius of $6.37*10^6$ m.


Model of the xz plane of the dipole’s magnetic field using Earth defined constants

Note that in this plot, the vectors seem to lack the fine structure of the previous field however the vectors are still properly oriented. They continue out from the southern hemisphere and seem to eventually loop back around and end in the northern hemisphere. I attempted a variety of scales and vector size values but no matter what, when the constants were applied to the equation they never had the same shape as the previous vector field. The other concerning aspect of the field was the fact that the vectors don’t seem to weaken as the field continues away from the Earth. It could be that the scale was just off so it was not evident, but even with the shifting of the scales, the vectors still seemed to remain the same size. The general direction and shape was obtained with these equations so I continued on to the next step.

3D Vector Fields:

Applying the same methods as before, I used the vector field components from the obtained equation and using the command VectorPlot3D, I created a vector field of the general equation without the constants. This field can be seen below.


Simple vector model of the magnetic field of a dipole

Without constants the shape of the vector field directly resembles the magnetic field of a dipole in 3 dimensions. Unfortunately, we also have the vectors not varying as expected as the distance from the origin increases. Regardless, after the success of obtaining the proper vector field shape and direction I decided to apply the proper constants to see if the shape would be altered as it was in the 2D graph. The constant applied vector field can be seen below. Once again the constants remain the same and the sphere in the center of the plot represents the Earth, assuming it is a sphere with a constant radius.


Vector model of the magnetic field of a dipole using Earth defined constants

As we can see, the same lack of detail is present when the constants are added to the vector field. However the shape and direction of the vector field are still generally correct though the fine structure is again lost. The vectors still don’t seem to vary with distance as well, which is troubling for actually plotting the field but the model does seem to accurately plot the field itself.

Further Plans: 

Unfortunately it is at this point where I ran out of time and needed to call it quits on the project, with only the magnetic field of the dipole modeled. While I will discuss what I would have liked to change and go into more detail on complications within my “Conclusion”, I would like to take a moment to discuss where I wanted to continue on with this project. Specifically, next on the list would have been adjusting the vector field  23 degrees from the z-axis to account for the Earth’s natural tilt. I experimented with a few methods to do this but none of them were truly successful.

Next I would have liked to incorporate a model of the dynamo theory to account for the Earth’s magnetic field. One of the biggest assumptions made within this project was that the Earth’s magnetic field simply existed as a giant bar magnet in the surface. The magnetic field is actually generated by the motion of liquid iron present in the Earth’s outer core. Using this knowledge, I was hoping to create an animation that would allow the observer to alter the size of the planet and then the percent of the outer core took up of the interior. This would allow the observer to have a general model about how the magnetic fields change with different core sizes and thus the location of the Van Allen Belts. Knowing this, scientists could prep future scientific equipment for radiation exposure even if they needed to orbit within the Van Allen belts.

Here are my Mathematica Notebooks in case there were any questions on my code: 2D Vector Fields, 3D Vector Fields, and General Experimenting



Griffiths, David, Introduction to Electrodynamics

Stern, David, Radiation Belts (

Dynamo Theory (

Magnets and Electromagnetism (




Relativistic E&M: Final Conclusions

While I did succeed in creating a model that allowed the user to view the field at different speeds, it only managed to show the direction and relative magnitude compared to the rest of the current fields. If I were to continue with this project, I would attempt to improve it in several ways.

-First, I would try to find a way to make Mathematica account for the absolute magnitude of the fields so that the uniform fields would appear to change as the speed was varied. This might be accomplished by using equipotential surfaces rather than vector fields.

-Second, I would write code such that the user could simply input (or maybe copy and paste) a particular E-field equation for a rest frame, and the program would interpret the equation and produce the vector plots with sliders as shown in my model.

-Third, I would model more discrete situations rather than infinite distributions of charge or current so that the changes in field would be more visible.

-Fourth, I might consider a different type of plot, since Mathematica seems to deal with this kind of plot poorly. For instance, a static plot of E versus v and B versus v at a single point in space might better show the changing relative magnitudes of the fields.

That being said, this project was successful in several regards. It emphasizes the interplay of electric and magnetic behaviors when considered in moving reference frames, and at least in the case of the point charge, it allows for an examination of the changing shape of the fields at different speeds. It also sets up the framework for being able to model more complex systems by providing a proof-of-concept using simpler cases. Therefore, this project successfully modeled the fundamentals of relativistic electromagnetism, and provides a good foundation for a more in-depth exploration.