In [2]:
import numpy as np
%pylab inline


Populating the interactive namespace from numpy and matplotlib

Monte Carlo Simulations

Jesse Livezey

Most of the integrals you have run into in your physics course have been exactly solveable. Often times in research, there will be a system that leads you to an integral that you cannot solve exactly. Generally one of two things happens at this point:

  • You can approximate the integral by another integral that is exactly solvable
  • You can numerically solve the integral using numerical methods (Monte Carlo being one example)

Millions of Small Simulations

At the heart of Monte Carlo methods is the idea that an analytical integral can be approximated by the average of many individual numerical simulations. In this case, we'll look at a research problem that I worked on as an undergrad.

Imagine shooting electrons at a fairly thin sheet of metal with cylindrical holes punched in it in a square grid patter. You want to know what percentage of electrons make it through the sheet as a function of angle (and later energy). This can be represented as an integral in an abstract way: $$\%_{\theta,E} = \int_{\text{initial conditions}}\text{acceptance(\theta,x,E)dEdx$$.

The lack of nice symmetry in the grid of holes makes this integral impossible to do by hand (although we'll exploit what symmetry it has.) We'll have to start lots of little particles and run them through a particle tracking code.

Let's say that our grid of holes is 1 cm from center to center and each hole is .25 cm in diameter. The sheet of metal the holes are punched in is .5 cm thick.


In [ ]:
# We'll create a "Simulation" class to store out variables and do the particle tracking
# Start by filling in the __init__ and run methods, then move on to the othersgr
class Simulation(object):
    """Class for creating and running particle tracking code"""
    def __init__(self):
        pass
    def update_position(self):
        pass
    def check_positions(self):
        pass
    def run(self):
        pass
    def results(self):
        pass

In [3]:
# You'll need a loop to initialize simulations with different angles and start the simulations
thetas = np.linspace(0., np.pi, 10)
results = np.zeros_like(thetas)
n_particles = 1000

for theta in thetas:
    pass

plt.plot(results)


Out[3]:
[<matplotlib.lines.Line2D at 0x10ad06410>]

Adding in Fields


In [ ]:

In the project I was working on, we could change the voltage on the metal grid to screen out electrons with low energies. At the same time, some deviced were in dipole magnetic fields. When these fields are added, it becomes neccesary to also simulate the particles as a function of their initial energy, not just angle.

Create a new simulation class that incorporates an electric or magnetic field (or both)! The electric field should point towards the grid so that electrons are repelled. You can set it up so that your simulation is is .5 cm tall before it hits the grid. This will allow you to convert voltages to fields. The magnetic field would also point normal to the grid surface. Try some different B field values to get a sense for what happens.


In [ ]: