[Tutor] Optimizing a simple radioactive decay simulation

[Tim Wilson]

Hi everyone,

My students and I did a fun little simulation of C-14 decay today using=20
pennies to represent the atoms. We flipped the coins repeatedly until all=
=20
the pennies had "decayed" (by turning up heads) and used the results to=20
make a graph of the decay curve.

I thought it would be fun to replicate the demo on the computer so I coul=
d=20
expand the simulation to many more pennies. My first approach, which I wi=
ll=20
include below, was to create Coin and DecaySim classes that contain the=20
needed methods. This worked perfectly, but was quite slow. Since I'm=20
interested in using simulations to teach science, I wonder if anyone has=20
suggestions for optimizing this code. I realize that there is probably a=20
much simpler approach, but I found this one appealing because of the 1:1=20
mapping between it and the real-life penny simulation.

-Tim

--snip--
"""
decay.py by Tim Wilson <wilson@visi.com>

This program simulates radioactive decay using coins.

"""

import random

class DecaySim:
    """Radioactive decay simulator class.
   =20
    Class uses 'coins' to simulate decayed and undecayed atoms.
   =20
    """
   =20
    def __init__(self, numcoins):
        """Create decay simulator instance."""       =20
        self.contents =3D []
        for i in range(numcoins):
            self.contents.append(Coin(1))
   =20
    def __len__(self):
        """Return number of coins left in simulator."""
        return len(self.contents)
      =20
    def removeTails(self):
        """Remove all coins from the simulator that are 'tails'."""
        self.contents =3D [coin for coin in self.contents if coin.value =3D=
=3D 1]
       =20
    def decay(self):
        """Flip each coin in the simulator and see which ones decay."""
        for coin in self.contents:
            coin.flip()

class Coin:
    """Coin class for use in radioactive decay simulation."""
   =20
    def __init__(self, value):
        """Create a coin instance. 1 =3D heads, 2 =3D tails"""       =20
        self.value =3D value
       =20
    def flip(self):
        """Flip the coin."""
        self.value =3D random.randint(0, 1)

def main():
    """Run the decay simulator."""
    numcoins =3D raw_input("How many coins? ")
    ds =3D DecaySim(int(numcoins))
    halflives =3D 0
    while len(ds) > 0:
        print "Undecayed (heads): %s" % len(ds)
        ds.decay()
        halflives +=3D 1
        ds.removeTails()
    print "\n%s half-lives required." % halflives
   =20
if __name__ =3D=3D '__main__':
    main()

--snip--

Here's the output on my 1-GHz Athlon with 512 MB RAM.

wilsont@galileo:~/Documents> time python decay.py
How many coins? 1000000
Undecayed (heads): 1000000
Undecayed (heads): 500350
Undecayed (heads): 250901
Undecayed (heads): 125447
Undecayed (heads): 62783
Undecayed (heads): 31182
Undecayed (heads): 15479
Undecayed (heads): 7796
Undecayed (heads): 3972
Undecayed (heads): 2001
Undecayed (heads): 1016
Undecayed (heads): 515
Undecayed (heads): 287
Undecayed (heads): 138
Undecayed (heads): 76
Undecayed (heads): 33
Undecayed (heads): 18
Undecayed (heads): 12
Undecayed (heads): 6
Undecayed (heads): 6
Undecayed (heads): 2

21 half-lives required.

real    12m6.267s
user    8m26.820s
sys     0m1.190s

--=20
Tim Wilson
Twin Cities, Minnesota, USA
Science teacher, Linux fan, Zope developer, Grad. student, Daddy
mailto:wilson@visi.com | http://qwerk.org/ | public key: 0x8C0F8813