An electric field line simulator straight in your browser!

Try the simulator here

Note that there are some limitations to this simulation, leaving it a bit unphysical when it comes to lines that might end up in open space, wrong field line density and other quirks. All in all it should be correct to a decent approximation, though.

For the programming enthusiasts out there, the simulation is created using Javascript and HTML5. You can check out the source code here.

Enjoy!

[TerminalIPythonApp] Unrecognized flag: '-wthread'

Without this flag, the Mayavi plots lock up the UI and hangs. If you want to get the possibility back to rotate and play around with the plots, just start IPython the following way from now on:

ipython --pylab=qt

This will launch IPython with the Qt backend and threading. Using only –pylab does not include threading. For easy and quick access, add the following to a file named .bashrc in your home folder:

alias pylab='ipython --pylab=qt'

From now on you can launch IPython just by typing

pylab

in a terminal.

ImportError: No module named enthought.mayavi.mlab

Now, to avoid switching the import statement every time you switch systems, you can make Python check if one of the versions is installed during import. If it is not, we’ll tell it to try the other. This is done in this simple command:

try:
    from enthought.mayavi.mlab import *
except ImportError:
    from mayavi.mlab import *

Just replace any other similar import statements the same way and your code should once again be working across all your installations.

So here you go. A pendulum on a cart.

I won’t go into details about deriving and solving the Lagrange equation for this system now, but I might get back to that at a later time. And of course there will be some added controls to this so you could set the initial positions, velocities and stuff. Sometime later, I hope.

The game is called Nanoparticles and was written by me about half a year ago. It is based on Coulomb’s law with a few modifications and Newtonian physics. But instead of having you calculate anything, the game does it all for you in the background. You can download the game for Ubuntu, Linux, Windows and Nokia phones here (Mac and Android versions will be available sometime in the future):

Download Nanoparticles

Your mission is to make sure the positive particles you control does not crash into the negative ones. That would cause it to annihilate. To do this, you need to place out other positive and negative particles to keep your particle away for as long as possible.

The game is also open source, so if you’re interested in programming and C++, you can have a look at the source code too. More information about the game and how to obtain the source code is found here.

But now, let’s have a look at matplotlib’s animation capabilities. The script below shows a very easy approach to animation in matplotlib. This results in an animation of the standing wave shown here:

The script is as follows:

from matplotlib.pylab import *  # pylab is the easiest approach to any plotting
import time                     # we'll do this rendering i real time

ion()                           # interaction mode needs to be turned off

x = arange(0,2*pi,0.01)         # we'll create an x-axis from 0 to 2 pi
line, = plot(x,x)               # this is our initial plot, and does nothing
line.axes.set_ylim(-3,3)        # set the range for our plot

starttime = time.time()         # this is our start time
t = 0                           # this is our relative start time

while(t < 5.0):                 # we'll limit ourselves to 5 seconds.
                                # set this to while(True) if you want to loop forever
    t = time.time() - starttime # find out how long the script has been running
    y = -2*sin(x)*sin(t)        # just a function for a standing wave
                                # replace this with any function you want to animate
                                # for instance, y = sin(x-t)

    line.set_ydata(y)           # update the plot data
    draw()                      # redraw the canvas

I’ve commented each line, since the script is fairly small, but let’s outline a few important things.

First of all, you need to disable interactive mode with

ion()

This makes it possible to animate, but disables all the button controls in matplotlib’s figure window.

Furthermore, we need to assign our plot to a variable (or pointer if you like), named line. This is done by plotting some dummy data:

line, = plot(x,x)

Note the comma after “line”. This is placed here because plot returns a list of lines that are drawn. Since we draw only one line we unpack only this by placing a comma and nothing else after “line”. If you plot multiple lines at once you may unpack them together by issuing a command like this:

line_a,line_b,line_c = plot(...

This plotting makes our axes aligned to the dummy data (in this case y = x). To make sure we have enough room, we manually set the limits of our y axis:

line.axes.set_ylim(-3,3)

Finally, we do our calculations and eventually change the y-data for each time step. After this, we draw everything back onto the canvas. This is done through the following two lines:

    line.set_ydata(y)           # update the plot data
    draw()                      # redraw the canvas

That’s it. Your animation should now look like the one above.

Saving the animation to a  file

Saving everything to file is fairly simple with the savefig command. You can save all the frames as images and then convert them to video using your favorite video editor. I recommend using ffmpeg from command line or Blender (which also does 3D). There are surely easier tools out there, but I find the work flow using ffmpeg and Blender quite quick, and they are also useful tools for many other tasks.

I’ve made the script animation ready below, and there’s really just one important thing to note: Earlier we showed the animation in real time. Saving all the images in real time lags badly and many frames are dropped. Because of this, we now set a delta time between each frame and iterates over each and every frame to make sure we have enough data for a smooth video.

Feel free to use the script below as it suits you:

from matplotlib.pylab import *  # pylab is the easiest approach to any plotting

ion()                           # interaction mode needs to be turned off

fig = figure(figsize=(16,9),dpi=80)     # create a figure in 1280 x 720 pixels
                                # this is not as simple as it could be, but
                                # it is like this because matplotlib is focused
                                # on print, not screen

x = arange(0,2*pi,0.01)         # we'll create an x-axis from 0 to 2 pi
line, = plot(x,x)               # this is our initial plot, and does nothing
line.axes.set_ylim(-3,3)        # set the range for our plot

t = 0                           # this is our relative start time
dt = 0.04
i = 0
while(t < 5.0):                 # we'll limit ourselves to 5 seconds.
                                # set this to while(True) if you want to loop forever
    y = -2*sin(x)*sin(t)        # just a function for a standing wave
                                # replace this with any function you want to animate
                                # for instance, y = sin(x-t)

    line.set_ydata(y)           # update the plot data
    draw()                      # redraw the canvas

    t = t + dt                  # increase the time
    i = i + 1                   # increase our counter

    # save the figure with a 4-digit number
    savefig("outdata/blah" + '%04d' % i + ".png")

This post is also arriving a bit late to be useful for anyone this semester, but what it shows might be useful for others. Or maybe for someone taking the FYS1120 course at a later time. So I’ll post it here still.

In mandatory exercise 2 we needed to load a huge file containing data from an AM radio signal. This was available a Matlab file and as a text file with all the values on each row. Those of us using Python realized quite quickly that most of the time spent on this exercise went to load the data before performing any computations on the data.

However, in proper hindsight, I thought it might have been a better idea to save the file as binary data instead of as a text file. After all, binary data saves space and usually is a bit quicker to load.

Doing this in NumPy is extremely simple. Just load the text file data and save it back in NumPy’s binary format:

from numpy import *
data = loadtxt("input_data.txt")
save("binary_data.npy", data)

Now, loading the file in your application is just as simple:

from numpy import *
data = load("binary_data.npy")

The time you’ll save on doing this for large data sets is extreme. Loading the data set we recieved in mandatory exercise 2 as a text file took me about 2 minutes and 21 seconds, while loading the data from the binary file format took only 0.06 seconds! Yes, that is 60 milliseconds.

On top of this, loading the data as text using “loadtxt” practically wasted 2 GB of memory on my computer, while loading it as a binary file used only about one hundred megabytes; just about the size of the binary data file.

Now I just wish I thought of this a couple of weeks earlier. But I guess that is what hindsight is for.

Four charges with different magnitude plotted in 3D using Mayavi

In one of my earlier posts about Mayavi, I wrote about how you could visualize 2D field line plots using the flow function. At the end of that post I added that Mayavi is actually best at 3D plotting, and to follow up on that I’ll show you some of these plots with a few example Python scripts you might try out on your own.

First of all, you might want to know how to install Mayavi. For those lucky ones of you who have freed yourself and jumped on the Linux bandwagon, installing Mayavi should be quite easy. If you are using Ubuntu in particular, you may just install the package mayavi2 using either Synaptic or apt-get. If you are on Windows or Mac, you may either install Enthought’s own Python distribution (EPD) or give a shot at compiling on your own. Just note that EPD is quite expensive, even though all its components are open source, but if you are a student or academic user you could go ahead and download the academic version for free. It is basically the same as the commercial one, but with an academic license. (Kudos to Enthought for both making Mayavi open source, building an business model around it and still providing a great solution for students!)

Now, Enter 3D!

The way you do your plots in Mayavi depends on what you want to express. Most likely, you would prefer to show some simple plots giving just the necessary amount of information to tell you how the electric field behaves around your charges. A simple example of this is shown below:

The simple plots often give you a great perspective about what happens in the electric field

On the other hand, you might want to give a strong visualization to show off the density and beauty of an electric field. In such a case, increasing the resolution of the flow/streamline seeds gives you a greater amount of field lines, which could result in plots like this:

An highger seed resolution gives a more dense plot.

Note that the plot shows the same four charges as above, but from a different angle and, of course, with more field lines.

The greatest part of using Mayavi to visualize these plots, however, comes from the fact that you may rotate the plot in real time in the Mayavi scene view. This gives you great control and insight of what you are plotting as you may rotate it as if it was a physical object in front of you. To show how interesting this may be, I’ve created a short video rotating around the same plot as above, recorded using Mayavi’s animation features:

You may even animate the charges individually, showing you what happens when a charge moves through space.

Note that the video below shows some artifacts around the charges. I believe I could have tweaked the settings a bit more to avoid these, but I decided it was good enough for the purpose of showing Mayavi’s capabilities.

If you would like to test out these plots on your own, you can download the source code here:

• Field lines around a dipole
• Field lines around multiple charges
• Field lines around multiple charges, animated

And if you couldn’t get enough of those field line plots, here are all the above and a some more, stacked in one set:

3D plot for exercise 2d

There are only two days left before the first mandatory exercise in Fys1120 is due. Last time I posted about Mayavi I wrote in a comment that I would give you some hints on how to use it with the oblig. Actually, there are not so many useful applications for Mayavi in this particular oblig, since most plots are 2D and are better created using scitools.

In fact, the 3D plots in this oblig might also be plotted using scitools, but I decided to use Mayavi to test it out for this particular purpose.

In my code I stored the position vector for the particle in an array r, which is an array created using the NumPy function r = zeros((n,3),float). If you have stored the position in a similar array, you may use the code below to plot the position in exercise 1d, 2a and 2d. After plotting, use your mouse to find a camera angle you’re satisfied with and hit the save-button in the toolbar.

from enthought.mayavi.mlab import *
fig = figure(fgcolor=(0,0,0), bgcolor=(1,1,1)) # set the background and foreground of our figure
plot3d(r[:,0], r[:,1], r[:,2], tube_radius=None, color=(0.807, 0.10, 0.023))
axesm = axes()
axesm.title_text_property.bold = False
axesm.label_text_property.bold = False
axesm.title_text_property.italic = False
axesm.label_text_property.italic = False
axesm.title_text_property.font_family = "times"
axesm.label_text_property.font_family = "times"
outlinem = outline()

A decent tip for the LaTeX-users out there is to save your file as “.pdf”. This way you may import it to a LaTeX document without losing quality. Note that the tube_radius parameter of the plot3d function must be set to None to avoid humongous PDF files.

Below you see the output I got for the different exercises, but don’t panic if your results are different. There might be some errors in my own code for all I know.

Mayavi renders great field line plots.

While searching for a good Python module to visualize electric fields, I found Mayavi. Developed by Enthought, Mayavi is a very good module for visualizing a huge range of different scientific data sets. Everything from surfaces, flows and streamlines to bar charts, 3D plots and contour surfs are beautifully drawn on screen and exported to several file formats, such as PDF, PNG, EPS and more.

What I needed it for, however, was to visualize electric field lines in the course FYS1120 at the University of Oslo. We were told to use Matlab with the streamline and quiver functions, but even so, I wanted to use Python and decided to do a search and see if something similar was possible with Python. It took me some time to figure out how to use the scitools package to do streamline plots, but eventually I made it. However, these were a bit tedious to get working correctly and looked only about as good as the Matlab plots.

I continued my search and found Mayavi as an alternative. It seemed like it was a bit out of the scope of my quest, but I decided to give it a serious try in any case.

I’m glad I did. The quality of the plots in Mayavi are amazing compared to Matlab’s. Let’s take the following field line plot done in Matlab as an example:

The field line plot in Matlab.

This is done with the streamline function, and if you have been working with field lines in electromagnetism, you’ll notice that the plot doesn’t really follow the definition of field lines.

In the field line definition, the tangent of any point on the field line is the same as the electric field in that point. At the same time, the density of field lines should represent the strength of the field. That is, the electric field in a specified area is proportional to the number of field lines in the same area.

This is where this Matlab approach with the streamline function doesn’t work well. Actually, it misses completely. As you see, the density of the field lines is bigger in the edges, further away from the charge, while it should be bigger in the middle, close to the charge.

Methodically, the streamline function allows us to select from which points we want to draw field lines, and we could of course generate the streamlines in a circle around the charge instead of in a grid like now. However, the plot will still be a bit “rough” in the edges and doing this for every point requires some extra amount of logic for the programmer.

Giving Mayavi a try at the job shows us how much prettier and more correct the representation in Mayavi is:

The same field line plot in MayaVi.

Now that’s what I call a field line plot! Giving our lonely charge a few positive and negative friends to play with yields another pretty picture:

A field line plot using Python with MayaVi, showing four charges

You might notice that this plot is also not perfect in terms of the field line definition, but it is very close and a lot better than what we would have gotten from Matlab’s equivalent plot.

How to do it

Below you see the full source code for this script. I’ve included comments on each line instead of describing the script in detail here. However, I do want you to notice one important thing:

Mayavi is very good at visualizing things in 3D space, but it does not seem to be intended for 2D use. Therefore, we need to do some hacks to simulate a 2D space. It works flawlessly when you know how to do this, but it makes everything seem a bit harder than necessary. The payoff, on the other hand, is humongous, and it is one-time thing to learn how to do.

As well, it makes it easily possible to scale our plot up to a 3D-plot instead, which I will show you in the next post about Mayavi.

Enjoy!

from numpy import *
from enthought.mayavi.mlab import *

N = 100 # the bigger the number, the more calculations your script has to make, but you save MayaVi's Runge Kutta method a lot of trouble and the total plot time is lowered
a = 0. # the lowest coordinate
b = 1. # the highest coordinate
dt = b / N; # the step we need to make to get to each position on our grid

q = [1., -1., 1., -1.] # the values of our four charges
qpos = [[0.56, 0.56, 0.50], # and their positions
        [0.26, 0.76, 0.50],
        [0.66, 0.16, 0.50],
        [0.66, 0.86, 0.50]]

x,y,z = mgrid[a:b:dt, a:b:dt, 0.:1.:0.5] # here is the trick - we want a 2d plot, but MayaVi works best in 3D space. Let's just keep the z-axis to a minimum (two steps)
Ex, Ey, Ez = mgrid[a:b:dt, a:b:dt, 0.:1.:0.5] # create a few arrays to store the vector field in, using meshgrid

for i in range(N): # iterate over all rows
    for j in range(N): # and all columns
        Ex[i,j] = 0.0 # set the value of each point to 0, initially
        Ey[i,j] = 0.0
        for num in range(len(q)): # for each charge, calculate the electric field it provides
            rs = ((x[i,j] - qpos[num][0])**2 + (y[i,j] - qpos[num][1])**2) # the distance from the point to the charge, squared
            r = sqrt(rs)
            q1x = q[num] * (x[i,j] - qpos[num][0]) / (r * rs) # the x-component of the field
            q1y = q[num] * (y[i,j] - qpos[num][1]) / (r * rs) # and the y-component
            # this is $\frac{q}{r^2} \hat{\mathbf r}$ on component form
            Ex[i,j] = q1x + Ex[i,j] # now, add this to the electric field in this point, together with the contribution from the other charges
            Ey[i,j] = q1y + Ey[i,j]

fig = figure(fgcolor=(0,0,0), bgcolor=(1,1,1)) # set the background and foreground of our figure
#obj = quiver3d(x, y, z, Ex, Ey, Ez, line_width=1) # uncomment this if you want a quiver plot
streams = list() # create a list to hold all our streamlines (or flows if you speak MayaVi)

for s in range(len(q)): # for each charge, create a streamline seed
    stream = flow(x,y,z,Ex, Ey, Ez, seed_scale=0.5, seed_resolution=1, seedtype='sphere') # the seed resolution is set to a minimum initially to avoid extra calculations
    stream.stream_tracer.initial_integration_step = 0.01 # the integration step for the runge kutta method
    stream.stream_tracer.maximum_propagation = 20.0 # the maximum length each step should reach - lowered to avoid messy output
    stream.stream_tracer.integration_direction = 'both' # integrate in both directions
    stream.seed.widget.center = qpos[s] # set the stream widget to the same position as the charge
    stream.seed.widget.radius = dt * 2 # and its radius a bit bigger than the grid size
    stream.seed.widget.theta_resolution = 30 # make the resolution high enough to give a fair number of lines
    stream.seed.widget.phi_resolution = 1 # but we are looking at a plane for now, so let's not have any resolution in the z-direction
    stream.seed.widget.enabled = False # hide the widget itself
    streams.append(stream) # and eventually, add the stream to our list for convenience

xlab = xlabel("x") # set the labels
ylab = ylabel("y")
showm = show() # show everything
axesm = axes() # add some axes
axesm.axes.y_axis_visibility = False # remove the z-axis (named y for some MayaVi reason)
fig.scene.z_plus_view() # let's look at it from the top!
fig.scene.parallel_projection = True # and we don't need any projection when looking at it in 2D
print "Done."