Physical Thought

Adobe Reader on Linux

2012-01-05T16:02:00.002Z

This afternoon I tried to print a chapter of my thesis to read through and figure out where I needed to make changes. I opened it in Adobe Reader on Linux, as I had done hundreds of times before, and tried to print the relevant page range. I was then informed that the document could not be printed!

Now, this is a little worrying for a thesis—something which is inherently meant for printing! The weird part was that I'd printed bits of it before, to test layout, colour, etc.

I tried to print other PDFs through Adobe, which all worked fine. I tried to print the thesis from Okular, which also worked fine. I then gave up and Googled the error message (document could not be printed). It appears that there was some kind of bug in Adobe Reader 9.4.2 that caused some documents to become unprintable. I downloaded (and, obviously, installed) Adobe Reader 9.4.6 and the problem was gone.

The moral of the story is, I suppose, always keep up to date on your Adobe Reader version. Obviously this is difficult on Linux since the versions are a way behind, and apparently they're going to stop making them altogether [warning: unsubstantiated rumour]. Unfortunately, none of the other PDF viewers I've tried on Linux actually work very well (both Okular and Evince have various rendering problems, though not the same problems as each other, which is odd because last time I checked they both used the same rendering engine library...)

Still, I'd rather work with PDF than any of the much less portable formats people regularly send around for no good reason (MS Word documents, OpenOffice documents, etc.)

Assassin's Creed

2012-01-04T16:01:00.004Z

I recently spent quite a while playing Assassin's Creed and Assassin's Creed II (I haven't finished II yet...)

The games are fantastic! The gameplay in ACII is rather better than the original, but the storyline continues sensibly and it's seriously addictive.

The addictive nature is helped (?) by the lack of obvious breakpoints; you tend to just keep playing, rather than thinking "this is a good place to stop".

I still have quite a bit of ACII to play through, then Brotherhood and Revelations, so I suspect my evenings for the next week or two are looking quite full!

Musings of a Ph.D. student

2011-09-21T14:06:00.000+01:00

As September 2011 draws to a close, I am coming up on the final six months of my Ph.D. and felt that now would be a good time to make a post about some of the things I've managed to do, and some of the things I have left to finish.

One of the "goals" of most Ph.D. students is to go somewhere nice for a conference or "school" (essentially a week-long course / conference where the goal is to educate or inform about some aspect of theory, experiment or technology). Over the course of my Ph.D. I've managed to go to Italy twice, Oxford once and Glasgow once...

The first trip happened as the first year of my Ph.D. was almost over; this should be a fairly familiar story for any UK experimental particle physicist, since I am referring to the STFC RAL Summer School in Oxford. This is a two-week long smorgasbord of quantum field theory, and is simultaneously thoroughly enjoyable and complete hell.

Let me explain that one; if you're interested in particle physics and you're doing an experimental Ph.D. the summer school is the closest you get to a decent understanding of the theory behind it all. The problems for the first week are interesting and solvable, while the problems for the second week are interesting and impossible unless you are a theorist! The hell aspect comes from having two weeks of lectures, problems, solutions, more lectures, too much food, and far too much alcohol.

Shortly after the summer school, I had a week-long trip to a place called Bertinoro, in the Forlì-Cesena province of the Emilia-Romagna region of Italy (near Bologna). This was for the 'European School of Computing' (ESC 09) which was a week of lectures on high-performance computing for particle physics. Some of the things discussed here were quite interesting, and my understanding of things like virtual memory improved dramatically over the course of the week. One of the interesting aspects was the food; each night we got to try a different local restaurant!

My third trip was to Siena, Tuscany (Italy). This was for the IPRD10 conference - a week-long conference on particle detectors. Here, the food was spectacular, but more expensive than in Bertinoro. The weather in Siena was wonderful, and it is a beautiful place to walk around. The same cannot be said for Pisa, though. Staying in Pisa was necessary since I flew into the airport there, and had a flight back from there. Once you've seen the tower (and it really does have quite a lean to it!) there's not an awful lot left.

One of the highlights of Siena was finding what I thought was a little restaurant on a side street. On walking inside, I discovered that it went quite a long way back, opening out onto a balcony looking out over the Tuscan countryside. This proved to be an ideal place for a lunch of lobster linguine!

The final trip was to Glasgow, for the IoP NPPD 2011 (Institute of Physics Nuclear & Particle Physics Division) Conference. Since several of us from Warwick were going, we rented a house in Glasgow for the duration of the conference. The house was amazing, and I think it helped to keep me sane, having a "house" to go home to rather than a hotel room and another restaurant dinner. We managed to find a Waitrose on the way between the house and Glasgow University, where the conference was held, so all of our evening meal and breakfast needs were taken care of. The cooking facilities in the house were second to none, and the entire experience was a huge success.

I managed to give talks at both the IPRD10 and NPPD11 conferences, on both occasions to a small but interested audience, so that is one of the Ph.D. checkboxes ticked.

I've spent three years now working on software for track reconstruction in liquid Argon TPCs (time-projection chambers). This software has been in a mixture of C++ and Python. As I write, I'm currently finishing some work on one such algorithm and heading towards a sensible analysis to round off my thesis... which I started writing last week.

So, in terms of things left to do, I need to formulate some coherent idea of what my final analysis will actually involve, get all the pieces of software together and working, and then run them on large data samples. All this while simultaneously writing a thesis...

Regarding the thesis itself, I spent a little while working on stylistic issues such as changing the chapter title pages from the default to something a little more colourful and modern. In doing so, I purchased the Futura Std Medium font from adobe.com to use as the chapter title font in my thesis (and then had to figure out how to make XeTeX pick up this font and use it).

I guess now I just have six months of being very busy, working simultaneously on analysing physics events, writing chapters of my thesis, writing any papers that seem appropriate along the way, and trying to remain sane.

Gears of War 3

2011-09-21T11:58:00.001+01:00

So I'm not what you'd call a "big gamer" (I'm neither obese, nor do I play a lot of games). I do, however, own an Xbox 360, and I was therefore eagerly awaiting Gears 3, which was released yesterday.

Now, my introduction to the Gears saga began with Horde on Gears of War 2, while visiting friends with Xboxes. This experience was sufficient to make me buy a 360 and a copy of Gears 2. Once I'd completed that, I bought Gears 1. I haven't actually *finished* Gears 1, but that is mostly due to a lack of time, rather than, say, lack of interest.

Anyway, yesterday was eagerly anticipated; I even went so far as to pre-order Gears 3 from Amazon (though they didn't manage to get it to me by the release day, which I consider to be a pretty major problem considering the concept of pre-ordering, and the concept of Amazon Prime... oh well!) Instead, I bought a copy from Sainsburys (and will sell the other copy to my housemates) and I spent much of the evening playing.

The first point I want to make is that, as I write this, I haven't finished Gears 3. In fact, I've only played through Act 1 on Normal difficulty, and played waves 1 through 33 of Horde on Normal difficulty. Then, extreme tiredness kicked in and I went to bed.

As far as gameplay is concerned, it feels much the same as Gears 2. The graphics have been polished a little, but it's basically the familiar experience; today's Penny Arcade pretty much sums it up.

I won't go into details about the storyline, but the first act has an interesting split to it, and the Horde mode has been enhanced somewhat, which makes it a little more interesting than the "wait, shoot, wait" cycle of Gears 2 Horde.

I haven't tried out the Beast mode yet, but that sounds like it should be fun; maybe I'll write more about that another day!

In all, Gears 3 delivered exactly what I expected... the third chapter in the Gears story, the familiar gameplay in familiar environments... having played the previous versions, I'm not having too much trouble playing through on Normal difficulty. I guess Hardcore and Insane are going to be more tricky; and probably even more fun.

As for achievements, the usual complement return (Waves 1-10 on Horde, Waves 1-50 on Horde, story progression in the campaign, etc.) though I haven't spotted any as ridiculously hard to achieve as Seriously 2.0 on Gears 2... yet!

More Android Development

2011-03-12T18:27:00.001Z

So far, I've managed to set up Eclipse to use the Android dev. tools and completed the 'Hello World' tutorial. This proved relatively simple, except for the matter of getting the entire emulator visible on a Macbook display. My resolution is 1280x800 and the Android emulator wants to display something 800 pixels high. Luckily, I discovered that you can pass the emulator a command line option to scale the device on screen:

emulator -scale 0.75

The full list of options is available at http://developer.android.com/guide/developing/tools/emulator.html

There's a 'Notepad Tutorial' to try next, so I'll probably do enough of that to get a decent idea of the components available, then start digging around the reference documentation for the components I'll actually need for one of the ideas I have (more on this later!)

Android Development

2011-03-12T02:04:00.001Z

It's been a while since I posted here, and since I've decided to spend a little time on a new project, I thought I might use the opportunity to write some new posts.

The project is to get acquainted with the Android SDK and write one or two apps for my Android phone (HTC Desire, since you didn't ask).

So far, it's going well. I'm taking the easy route, as outlined at http://developer.android.com/sdk/installing.html so I have downloaded the SDK starter, and downloaded Eclipse. I'm now installing the 'ADT' (Android Developer Tools?) into Eclipse.

I've never really bothered with Eclipse before, since I prefer to do most of my software development with gvim and a terminal window, but it is apparently the easiest way to get started with Android development. We'll see...

Comment Moderation

2010-06-30T23:27:00.000+01:00

In the last few days I've noticed a number of spam comments appearing on here. I've turned on comment moderation in order to prevent these receiving any airtime, but please accept my apologies if you are a genuine commenter and your comment is not published immediately; I'll approve the comments just as soon as I can.

Wine Reviews

2010-06-26T22:14:00.000+01:00

Some semi-regular readers may have noticed that I occasionally post wine reviews on this blog. Today, I created a new blog, Distilled Taste, where I, along with two friends, will review and discuss beers, ciders, wines and spirits from around the world. If you're interested in that kind of thing, why not take a look?

The Moon

2010-05-23T01:41:00.000+01:00

This is the best image from fifteen photographs of the Moon that I took this evening. My flat has a skylight thing in the roof, which opens and tilts in various ways (it would be great for using a telescope, actually...) so I decided to set up my camera and tripod and take some photos of the moon, since it was a clear night.

Wine: Shelter Bay Sauvignon Blanc 2009

2010-04-30T21:02:00.000+01:00

This sauvignon blanc is from Marlborough, New Zealand, and possesses strong gooseberry aromas with hints of tropical fruit. The fruit theme continues through to the taste, which is long and smooth. This is a very clean, dry white, with refreshing characteristics ideal for Spring and Summer drinking. I had it with tagliatelle carbonara, but it would go well with fish, chicken or light cheeses.

Must look like I'm running away, to you at your faster pace

2010-03-19T22:48:00.001Z

(Post title inspired by Tim's blog post, over at http://timfightsbears.blogspot.com/)

So I thought I'd make a more traditional "blog-style" post; talking about things I've done, and such...

UK Neutrino Network Meeting
I attended the March 17th UK Neutrino Network Meeting (link here) in Sussex, this week. As an experimental physicist, most of the theory talks (so, most of the session) went right over my head, but as someone with an interest in theory and phenomenology it's always nice to hear about what's going on in the field; I had no idea so many people were working so actively on leptogenesis. Seems that the Harrison-Perkins-Scott tri-bi-maximal mixing stuff is pretty popular, too, which is nice.

The more experimentally-oriented talks proved quite interesting; it's good to see there's a lot of work going on, despite the complete lack of funding!

Paper
I've been working towards writing a paper on the stuff I've been doing for my Ph.D. since summer 2009. It's looking like it might be finished in the more imminent future, for the first time since then; the basic software framework is in place and all the algorithms are working to varying degrees, so unless another big surprise is lurking just around the corner, I should be able to actually make some progress. I can't really say any more about this at present, but stay tuned for details when the thing is finally published!

Food
I haven't really had time to cook elaborate meals, lately, but tonight I made a sort of variant on Chicken Basque. I had in mind to do some kind of Coq au vin kind of thing, and I was looking through cookbooks to find some ideas. What I ended up doing was to fry chicken drumsticks and thighs in a little groundnut oil until the skin went golden, then put them aside on a plate, and fry onions, chorizo, garlic and red pepper in the same oil. I then added some sun-dried tomato, and about 150g basmati rice, before pouring in about 250ml of red wine (Cono Sur Pinot Noir, Chile), 200ml chicken stock, black pepper, paprika and herbs. Once this was all simmering again, the chicken pieces went back in, and the whole thing went in the oven at about 170C (180C if you don't have a fan oven) for 45 minutes. The result was tasty.

Sleep
I haven't been sleeping much. This is bad. I need more exercise, and to stop working earlier in the evening!

The Road To Quantum Field Theory: Updated Post List

2010-03-04T20:14:00.000Z

Once again, I'm posting an index, essentially, of the posts in the `Road to QFT' series, so far. The previous list was here, but I'll reproduce all of the links in this post too.

Pre-Series Posts
These were posts on related topics (special relativity, quantum mechanics and particle physics)

The `Road To Quantum Field Theory' Series
Posts in the `Road to QFT' series, in chronological order.

Out-of-Series Posts
Posts on related topics which were not intended as part of the `Road to QFT' series, but may be of interest anyway.

Basic Relativity

Roadmap
As before, I'll post a rough list of things still to be covered. By comparing the previous roadmap to the most recent two posts, you'll notice I've done things slightly out of the order I mentioned there, so this is really a very approximate guide.

Relativistic QM: The Klein-Gordon Equation & Spin-Zero Particles
Relativistic QM 2: The Dirac Equation & Spin-Half Particles
Quantum Electrodynamics (QED)
Quantum Chromodynamics (QCD)
SU(2) and Electroweak Unification
Quark Flavour Mixing
Spontaneous Symmetry Breaking and the Higgs Mechanism
Beyond the Standard Model: Neutrino Mass Terms & Neutrino Mixing

This is still a very ambitious set of things to cover. In particular, the posts so far have only brushed the surface, providing the background needed to understand relativistic quantum mechanics and gauge theories. The real work is yet to come! Note also that I am not a theorist, so my coverage of some of the more advanced topics may be less than completely thorough! I've added sections on quark flavour mixing and neutrino masses & mixing to the roadmap since they are of particular interest to me.

Gauge Invariance in Quantum Mechanics

2010-03-02T21:37:00.070Z

This post continues the 'Road to Quantum Field Theory' series.

The Lorentz force law for a non-relativistic particle of charge q, moving with velocity v in both electric (E) and magnetic (B) fields, is given by:

\mathbf{F} = q\mathbf{E} + q\mathbf{v}\times\mathbf{B}

This can be derived from Hamilton's equations using the classical Hamiltonian given below.

H = \frac{1}{2m}\left(\mathbf{p}-q\mathbf{A}\right)^2 + qV

The Schrödinger equation for a charged particle in an EM field is,

\left[\frac{1}{2m}\left(-i\nabla - q\mathbf{A}\right)^2 + qV \right] \psi(\mathbf{x},t) = i\frac{\partial \psi(\mathbf{x},t)}{\partial t}

obtained from the Hamiltonian through the substitution

\mathbf{p} \rightarrow -i\nabla

as is usual for the quantum mechanical momentum operator. We can identify operator combinations:

\mathbf{D} = \nabla - iq\mathbf{A}

D^0 = \frac{\partial}{\partial t} + iqV

which replace the operators

\nabla, \frac{\partial}{\partial t}

when we move from the free-particle Schrödinger equation to the electromagnetic field case.

The solutions for the wavefunctions of the Schrödinger equation describe completely the behaviour of a particle under the influence of the potentials A and V, but these potentials are not unique, as I showed in the previous post in this series. Instead, they can be changed by a gauge transformation:

\mathbf{A}\rightarrow\mathbf{A}^\prime=\mathbf{A}+\nabla\chi

V\rightarrow V^\prime = V - \frac{\partial\chi}{\partial t}

Maxwell's equations for E and B will remain invariant under these transformations (this was the topic of the previous post), but will the physics described by the Schrödinger equation be the same if we make these changes to the potentials there?

The answer is no! The Schrödinger equation is not gauge invariant, since the same wavefunction cannot satisfy both the original and transformed versions! All is not lost, however. The wavefunction itself is not directly observable, whereas E and B are. If we do not require the wavefunction to remain invariant when the potentials undergo a transformation, then we are free to change the wavefunction such that we retain invariance of physical observables. Write the Schrödinger equation in terms of some transformed wavefunction,

\psi\rightarrow\psi^\prime

i.e.

\left[\frac{1}{2m}\left(-i\nabla - q\mathbf{A}^\prime\right)^2 + qV^\prime\right] \psi^\prime(\mathbf{x},t) = i\frac{\partial \psi^\prime (\mathbf{x},t)}{\partial t}

Now, the form in which we have written the equation above is identical to that of the original Schrödinger equation, except that we are writing it in terms of primed quantities (psi', A', V') instead of unprimed (psi, A, V). Both equations describe the same physics, so if we can find such a psi', then the Schrödinger equation is gauge covariant; that is, it maintains the same form under a gauge transformation.

We know the relationship between A, V and A', V', so we can write down the transformation for psi':

\psi^\prime (\mathbf{x},t) = e^{iq\chi(\mathbf{x},t)}\psi(\mathbf{x},t)

where

\chi

is the same space- and time- dependent function appearing in the transformations of A and V. We can verify that this results in the gauge covariance of the Schrödinger equation:

(-i\nabla - q\mathbf{A}^\prime)\psi^\prime=\left[-i\nabla-q\mathbf{A}-q(\nabla\chi)\right]e^{iq\chi}\psi

=q(\nabla\chi)e^{iq\chi}\psi + e^{iq\chi}(-i\nabla\psi) + e^{iq\chi}(-q\mathbf{A}\psi) - q(\nabla\chi)e^{iq\chi}\psi

The first and last terms cancel, leaving:

(-i\nabla - q\mathbf{A}^\prime)\psi^\prime = e^{iq\chi}(-i\nabla - qA)\psi

which can be written as:

(-i\mathbf{D}^\prime \psi^\prime) = e^{iq\chi}(-i\mathbf{D}\psi)

The space-time dependent phase factor feels the action of

\nabla

, but passes through the combined D' operator, converting it to D, so it becomes clear that D'psi' is related to Dpsi in the same way psi' is related to psi! Similarly,

(iD^{0\prime}\psi^\prime = e^{iq\chi}(iD^0\psi)

We can now erite:

\frac{1}{2m}(-i\mathbf{D}^\prime)^2 \psi\prime = e^{iq\chi}\frac{1}{2m}(-i\mathbf{D})^2 \psi = e^{iq\chi}iD^0\psi = iD^{0\prime}\psi^\prime

which demonstrates the correct relationship between psi and psi'. The question remains whether the same physics is described by both psi and psi'. We can see that it is in a number of ways:

1. The probability density is given by

|\psi|^2 = \psi^\dagger \psi

, which is equivalent to

|\psi\prime|^2 = \psi^{\prime\dagger} \psi^\prime = \psi^\dagger e^{-iq\chi}e^{+iq\chi}\psi = \psi^\dagger \psi

2. The probability current,

\psi^\dagger(\nabla\psi) - (\nabla\psi)\psi^\dagger

is not invariant under a gauge transformation, but when we replace

\nabla \rightarrow \mathbf{D} ~~,~~ \frac{\partial}{\partial t} \rightarrow D^0

then,

\psi^\dagger\prime(\mathbf{D}\psi^\prime = \psi^\dagger e^{-iq\chi} e^{iq\chi}(\mathbf{D}\psi) = \psi^\dagger\mathbf{D}\psi

(and similarly for the other term).

Hence, identical physics is described by both the wavefunction psi, and the gauge transformed wavefunction psi'; the gauge invariance of Maxwell's equations presents as a gauge covariance in quantum mechanics, provided that we transform not only the potentials, but the wavefunction also:

\mathbf{A}\rightarrow \mathbf{A}^\prime = \mathbf{A} + \nabla\chi

V\rightarrow V^\prime = V - \frac{\partial \chi}{\partial x}

\psi\rightarrow \psi^\prime = i^{iq\chi}\psi

Finally, we note that the new differential operators,

\mathbf{D} = \nabla - iq\mathbf{A} ~~,~~D^0 = \frac{\partial}{\partial t} + iqV

can be written in a Lorentz covariant form:

D^\mu = \partial^\mu + iqA^\mu

This allows us to write:

-iD^{\prime\mu}\psi^\prime = e^{iq\chi}(-iD^\mu \psi)

and it follows that an equation can involving the operator

\partial^\mu

can be made to be gauge invariant by under the combined transformations:

A^\mu \rightarrow A^{\prime\mu} = A^\mu - \partial^\mu \chi

\psi \rightarrow \psi^\prime = e^{iq\chi}\psi

provided the minimal substitution is also made:

\partial^\mu \rightarrow D^\mu = \partial^\mu + iqA^\mu

This provides a simple mechanism to get the wave equation for a particle in an electromagnetic field from the equation for a free particle - by making the above substitution. This forms the basis of the gauge principle, i.e. the form of an interaction is determined by an insistence on local gauge invariance.

That's enough for now; next time we'll move on to study the Klein-Gordon equation for relativistic spin-zero particles.

Farscape Intro Randomiser

2010-02-08T15:01:00.000Z

While watching Farscape Season 2, Episode 4 (Crackers Don't Matter) last night, I wrote a Python script to randomise elements of the Farscape intro. Here's an example:

My name is John Crichton, a wormhole hit and I got shot through a military commander. Now I'm lost in some distant part of the universe on an insane radiation wave, full of escaped prisoners. I'm being hunted by a living ship. Doing everything I can. I'm just looking for a way home.

import random

adjectives = [
 ('insane', 'an'),
 ('living', 'a'),
 ('escaped', 'an'),
 ]

nouns = [
 ('ship', 'a'),
 ('radiation wave', 'a'),
 ('military commander', 'a'),
 ('prisoner', 'a'),
 ('wormhole', 'a')
 ]


parts = []
parts.append(random.choice(nouns))
nouns.remove(parts[0])
parts.append(random.choice(nouns))
nouns.remove(parts[1])
parts.append(random.choice(adjectives))
adjectives.remove(parts[2])
parts.append(random.choice(nouns))
nouns.remove(parts[3])
parts.append(random.choice(adjectives))
adjectives.remove(parts[4])
parts.append(random.choice(nouns))
nouns.remove(parts[5])
parts.append(random.choice(adjectives))
adjectives.remove(parts[6])
parts.append(random.choice(nouns))
nouns.remove(parts[7])

result = 'My name is John Crichton, '
result += parts[0][1] + ' ' + parts[0][0]
result += " hit and I got shot through "
result += parts[1][1] + ' ' + parts[1][0]
result += ". Now I'm lost in some distant part of the universe on "
result += parts[2][1] + ' ' + parts[2][0] + ' ' 
result += parts[3][0]
result += ", full of "
result += parts[4][0] + ' ' + parts[5][0] + 's'
result += ". I'm being hunted by "
result += parts[6][1] + ' ' + parts[6][0] + ' ' + parts[7][0]
result += ". Doing everything I can. I'm just looking for a way home."

print result

Orange Tiramisu

2010-01-26T17:05:00.001Z

At the weekend I made (yet another) tiramisu. This time, though, I used orange juice & Cointreau instead of coffee & amaretto, and I grated the zest of one orange into the egg & mascarpone mix. The result was very tasty - the orange and the dark chocolate grated on top worked well together. I recommend it.

Tiramisu

2009-12-23T22:19:00.001Z

It seems I make Tiramisu every Christmas. My recipe is pretty much the same as the one in Delia Smith's Winter collection, but where she uses three egg yolks and two egg whites, I also use the remaining egg white. Oh, and amaretto instead of dark rum. Delia's recipe is available at http://www.deliaonline.com/recipes/cuisine/european/italian/tiramisu.html

I either make this in a single large flat dish, or in cocktail glasses. This year I made it in cocktail glasses - to be served individually. They're topped with chopped dark chocolate and a generous dusting of cocoa powder.

Lorentz-Covariant and Gauge-Invariant Electrodynamics

2009-12-08T21:38:00.000Z

Over a year ago I started a series of blog posts aimed at describing quantum field theory, beginning with Lagrangian mechanics and working all the way up. It's been a while since I made any posts in that series. Indeed, the last post was a roadmap of sorts; a list of the topics still to be covered.

The roadmap hinted that the next step is a relativistic description of electromagnetism. Actually, what we really want to do is describe electromagnetism as a gauge theory. That is, a theory which is gauge invariant. By the end of this post, what we mean by gauge invariance will become obvious.

So, here goes.

Maxwell's Equations
Any treatment of electromagnetism must start with Maxwell's equations, presented below.

(1)\,\,\,\,\,\nabla\cdot\mathbf{E} = \rho_{em}

(2)\,\,\,\,\,\nabla\times\mathbf{E} = -\frac{\partial\mathbf{B}}{\partial t}

(3)\,\,\,\,\,\nabla\cdot\mathbf{B} = 0

(4)\,\,\,\,\,\nabla\times\mathbf{B} = \mathbf{j}_{em}

Here,

\rho_{em}

and

\mathbf{j}_{em}

are the electric charge density and electric current density, respectively, and we are working in Heaviside-Lorentz units. The densities act as the sources of the electric and magnetic fields.

Taking the divergence of (4) leads to a problem; the continuity equation for electric charge states that

(5)\,\,\,\,\,\frac{\partial \rho}{\partial t} + \nabla\cdot\mathbf{j} = 0

(note that I will omit the "em" subscript on the densities from now on).

Since

\nabla\cdot\nabla\times\mathbf{B} = 0\Rightarrow\nabla\cdot\mathbf{j}=0

This is only true when the charge density does not vary with time. In general, Ampere's law must be modified to be

(6)\,\,\,\,\,\nabla\times\mathbf{B}=\mathbf{j}+\frac{\partial\mathbf{E}}{\partial t}

This is now consistent with (5), the continuity equation.

The continuity equation states that the rate of change of charge in some volume is due entirely to the flux of current through its surface. That is to say, electric charge is conserved. Since this volume can be made arbitrarily small, we can say that electric charge is locally conserved. This means that we can't balance things out by having a negative charge appear out of nowhere, as long as a positive charge is made the other side of the universe (or a negative charge disappears on the other side of the universe); this conservation law has to apply on an arbitrarily small scale.

We can introduce the vector potential

A_\mu

in place of the fields E and B:

(7)\,\,\,\,\,\mathbf{B}=\nabla\times\mathbf{A}

(8)\,\,\,\,\,\mathbf{E}=-\nabla V-\frac{\partial \mathbf{A}}{\partial t}

This defines the three-vector potential A, and the scalar potential V (sometimes called the electrostatic potential, and sometimes denoted

\phi

. Equations (2) and (3) are then satisified automatically.

The gauge invariance of electromagnetism is based on the fact that, for a given E and B, the equations above do not uniquely determine A and V. In fact, one can apply transformations to A and V which leave E and B completely unchanged; and therefore leave all of the physics unchanged. These are known as gauge transformations, and the invariance of Maxwell's equations under these transformations is known as gauge invariance

We are free to change A as follows;

(9)\,\,\,\,\,\mathbf{A}\rightarrow \mathbf{A}^\prime = \mathbf{A}+\nabla\chi

where

\chi

is an arbitrary scalar function. Since

\nabla\times(\nabla\chi) = 0

, this does not change B.

To preserve E when making the above transformation, we must simultaneously transform the scalar potential V:

(10)\,\,\,\,\,V\rightarrow V^\prime=V-\frac{\partial\chi}{\partial t}

The transformations (9) and (10) can be written in a single expression using the four-potential,

(11)\,\,\,\,\,A^\mu = (V, \mathbf{A})

Note now that the differential operators,

(12)\,\,\,\,\,\left(\frac{\partial}{\partial t}, -\nabla\right)

form a four-vector operator

\partial^\mu

. The gauge transformation can then be specified by:

(13)\,\,\,\,\,A^\mu \rightarrow A^{\prime\mu} = A^\mu - d^\mu\chi

Maxwell's equations can be written in a manifestly Lorentz covariant form (meaning, they have the same form under Lorentz transformations) using the four-current:

(14)\,\,\,\,\,\mathbf{j}^\mu = (\rho,\mathbf{j})

The continuity equation can then be written:

(15)\,\,\,\,\,\partial_\mu j^\mu = 0

which in turn allows us to write equations (1) and (6) as:

(16)\,\,\,\,\,\partial_\mu F^{\mu\nu} = j^\nu

The equation above introduced the Field Strength tensor (or the Faraday tensor),

(17)\,\,\,\,\,F^{\mu\nu}=\partial^\mu A^\nu - \partial^\nu A^\mu

Under a gauge transformation (13), the Faraday tensor (17) remains unchanged; it is gauge invariant, and so, therefore, are Maxwell's equations in the form (16). As I already said above, Maxwell's equations are Lorentz covariant in the form (16), leading to a Lorentz covariant and Gauge invariant theory.

Now, one can write

(18)\,\,\,\,\,\Box A^\nu - \partial^\nu(\partial_\mu A^\mu) = j^\nu

Equation (6) can be inferred (and indeed, was inferred by Maxwell) from the local charge conservation requirement expressed by the continuity equation (15). The field equations (16) automatically include the continuity equation. Mathematically speaking, the Faraday tensor (17) is a sort of "four-dimensional curl". This is unchanged by the gauge transformation (13), suggesting that gauge invariance in electromagnetism is related to charge conservation. The reality is a little more complicated.

It was shown in 1949 by Wigner that the principle that no physical quantity can depend on the absolute value of the electrostatic potential V, combined with conservation of energy, implies that charge conservation holds. This relates charge and energy conservation to an invariance under the transformation of the electrostatic potential by a constant. Charge conservation alone does not require the more general space-time dependent transformation of the gauge invariance we've discussed above.

Changing the electrostatic potential by a constant amount is an example of a global transformation. Invariance under this global transformation is related to the conservation of electric charge, but it is not sufficient to obtain all of electromagnetism. Instead, we must impose a local change in the electrostatic potential V (the time-derivative term in (10), which is compensated by a simultaneous change in the vector potential A, leaving Maxwell's equations ultimately unchanged. By including these magnetic effects, the global invariance under a change in V, related to conservation of electric charge, can be extended to a local invariance.

The concept of local gauge invariance is important for the development of quantum field theory in subsequent posts. Next time, a tour of gauge invariance in quantum mechanics, and the gauge principle.

Much of this post was based on Chapter 3 of Aitchison & Hey (Gauge Theories in Particle Physics, 3rd Edition, Vol. I - From Relativistic Quantum Mechanics to QED). It is a book I strongly recommend purchasing if you have an interest in quantum field theory. The second volume discusses QCD and Electroweak theory, and is another good purchase.

Mince Pies

2009-12-06T16:42:00.002Z

Here's my recipe for mince pies. I don't tend to make the mincemeat myself; get a decent jar from a supermarket. You'll need about 400g for 12 mince pies, and you can always pour a little Cognac into it to add a bit of a kick!

175g (6oz) Plain flour
40g (1.5oz) White fat (I use White flora in place of lard)
40g (1.5oz) Unsalted butter
Pinch of salt
Cold water

The key to making pastry is to keep everything as cold as possible. Use the fats straight from the fridge. Add the fats to the flour in a large bowl, then run your hands under the cold tap for a while and dry them. Next, using only your fingertips to avoid transferring too much heat into the pastry, rub the fats into the flour.

Once this is done, add the very cold water, a small amount at a time, and use a fork to mix it into the flour. What you want is for the pastry to just come together and leave the bowl clean, without being too wet or sticky. Add the water slowly and stop as soon as the pastry holds together.

Put the lump of pastry into a plastic bag, and leave in the fridge for at least half an hour. This helps to cool the pastry down again after working with it.

When you're ready to make the pies, preheat the oven to 200C (190 for a fan oven) and grease the tins. You can either use separate small tins, or one of those tins with twelve hollowed out bits, more commonly used for muffins and yorkshire puddings!

Flour the work surface, your hands, and the rolling pin, then roll out the pastry, being sure to roll in all directions (turn the pastry round as you roll to do this). Then, use large pastry cutters (or the top of a pint glass) to cut out 12 bases. Put these into the tins and gently push them into shape. Fill with mincemeat to just below the level of the pastry.

If necessary, re-roll your pastry and use a smaller cutter (or a smaller glass!) to cut the lids. Place the lids on top of the pies, glaze with a little milk, and make a small hole in the centre using a pair of scissors or a barbecue skewer.

Bake at 200C for 25 to 30 minutes. Allow to cool in the tins for 10 to 15 mins before transferring to a wire cooling rack. To serve, dust with icing sugar.

The mince pies will keep in an airtight container for about a week, or can be frozen (but don't dust with icing sugar until you're ready to serve them!)

makefilter: Colour for make output

2009-11-09T17:05:00.001Z

I wrote a Python program today which runs make (and passes command-line arguments to it) and filters the output, adding ANSI/VT100 colour codes to highlight certain parts.

The goal was twofold; to distinguish different parts of the gcc / g++ command, and to distinguish commands from error messages and other output lines.

The highlighting is performed by directly emitting the VT100 control sequences. Eventually I might change this to use ncurses / terminfo / something more portable, but for now this will have to do.

There's a webpage relating to the program at http://www.warwick.ac.uk/~phrebd/makefilter.html and a screenshot below (the same screenshot as on the webpage!)

I'm aware of several other projects to achieve similar effects. The closest is pretty_make.py, available at http://phil.freehackers.org/pretty-make/index.html. There is also colorgcc, http://schlueters.de/colorgcc.html, which does a similar thing for gcc error messages (the down side is it's written in Perl).

Update: Version 1.3 of makefilter has been released today (Saturday November 14, 2009). This includes a few minor improvements and bugfixes. The screenshot above has also been updated accordingly.

Google C++ Testing Framework

2009-09-29T21:06:00.001+01:00

I was looking around for C++ unit testing frameworks, earlier. I'd already disqualified CppUnit and Boost.Test for various reasons, so I started out with Wikipedia's List of Unit Testing Frameworks, from which I quickly gravitated towards the Google C++ Testing Framework.

Google Test (as it is apparently known, for short) is a decent stab at a unit testing framework. Admittedly, I use Gmail, Blogger, and Google Maps, and I'm generally quite a fan of all things Google, but still, I had already thrown out two unit test frameworks based on their unsuitability (if that is a word). Google Test looks to be well designed, based on the xUnit concept, and provides a lot of nice stuff, all of which seems to be well documented. What sold it, for me, was the multitude of build systems. Google Test builds under XCode, Visual Studio, GNU make, autotools, Scons, and probably some other things I forgot already.

The build system was my primary reason for discarding CppUnit; I wanted something I could integrate neatly into an existing project's build system, and not have to introduce an extra system dependency on a unit test framework. So, since Google Test appeared to have the most flexible build system, I started to incorporate it into my project.

At that point, I went home, so I have yet to see how well it works for actual testing, but if it's as good as the build integration, it will be a pleasure to work with!

Bad Examples

2009-09-23T12:46:00.004+01:00

I realised today (while browsing examples in the CppUnit Cookbook, documentation for the CppUnit Unit Testing framework) that there exists, on the Internet, a plague of incredibly bad examples.

Allow me to explain. I'm trying to figure out how to use CppUnit, and I'm trying to do it quickly. Like most programmers, I go straight for the examples, because I can understand code very quickly, without having to read a lot of text explaining the history of, and relationship between, various unit testing frameworks.

The code examples look good at first glance - they tell me everything I need to know about how to use the framework; I just derive from a certain class, implement a couple of methods, and stick it all into another class to run the tests.

Ah, but wait. The first example, the one with the class I'm supposed to derive from, doesn't tell me which header to include. I don't know where the CppUnit::TestCase class is declared, and I'm going to have to hunt through the documentation (or grep the source) to find out, because at no point does it give a complete example of this stuff, with the headers!

Now, I understand that headers often get in the way of demonstrating a single point, and that writers often leave out common headers to save space and to make it clearer what is going on in the example of interest. However, when the header is non-standard, and is absolutely critical to what follows, it's important to keep the include statement visible.

Ok, rant over. Now I'm going to go back to hunting for the declaration of CppUnit::UnitTest.

Edit: Crisis over. I dug deeper into the documentation, eventually finding a class hierarchy with information about the headers in which things are defined, but still, having complete examples would've saved me a lot of time.

Wine: Domaine de la Vougeraie 'Les Petits Noizons' 2006 (Pommard)

2009-09-21T23:52:00.003+01:00

A strong offering from Bourgogne, this 100% Pinot Noir has a dark red-brown colour. Initially it does not offer much by way of fragrance, but after leaving to the air for a little, strong notes of plum, cherry, melon and blackcurrant practically leap out of the glass.

The first taste similarly reveals little but a smooth, velvet character. Subsequent tastings open up the robust, fruity character of this wine, exhibiting many rich berry flavours; blackberry and raspberry in particular. There's a very smooth finish with a hint of very dark chocolate.

The entire experience is a velvety smooth, rich and full wine. There appears to be a very slight nutty character both to the taste and fragrance, though this may be a mental inference from the name, refering to a walnut plantation.

In all, an incredible wine, with a very long finish. It would go well with rare steak, lamb or other red meats. Avoid strong flavours such as cheese, which would detract from the fullness of this wine.

Reversing Lines In A File

2009-07-15T11:50:00.002+01:00

As is often the case, an interesting little text manipulation task came up in the office today; given a file containing lines of text, reverse it (i.e. put the last line first, the penultimate line second, etc.)

I have a feeling there's a neat little command-line utility that already does this (but can't remember what it is...), but my mind soon went to a short Python script. Here's my first stab at this:


#!/usr/bin/env python
import sys

lines = []
for line in open(sys.argv[1]):
 lines.append(line.strip())

lines.reverse()
for line in lines:
 print line

After this, I tried to shorten it a bit. Putting the lines into a list could be done nicely with a list comprehension,


#!/usr/bin/env python
import sys

lines = [line.strip() for line in open(sys.argv[1])]

lines.reverse()
for line in lines:
 print line

Finally, I discovered the reversed() function, which allows you to create reverse iterators for any Python sequence. Here's where the list comprehension and Python's iteration stuff really comes into its own, leading to a wonderful two-liner (not counting the module import and the #!)


#!/usr/bin/env python
import sys
for line in reversed([line.strip() for line in open(sys.argv[1])]):
 print line

Wine: Santa Rita Merlot Reserva 2007

2009-06-20T19:54:00.002+01:00

This merlot from Chile's Maipo valley has all the hallmarks of a decent Bordeaux. It has a wonderful ruby-red colour and rich plum aromas with a sweet, decadent note. The robust taste carries strong blackcurrant flavours with undertones of plum and blackberry. A smooth finish rounds off this classic Chilean offering. Strong and alcoholic, this wine is warming and has a good kick to it.

Drink alone (that is, without food, not without other people) or with food that will stand up to the strength of this wine; steak, blue cheese, mildly spiced or high-garlic foods.

Wine: Vistasur Cabernet Sauvignon 2007

2009-06-17T17:47:00.003+01:00

This Chilean Cabernet has a rich cherry aroma and a deep red colour. It tastes light and smooth, with notes of ripe summer berries.The finish is powerful yet fleeting. Goes well with blue cheese!

Wine: Vistasur Sauvignon Blanc 2008

2009-06-12T22:51:00.003+01:00

I had this wonderful Chilean white with salmon (and again with tuna). Chilled in the fridge, it was refreshing after warm June days.

The colour is delicate and clear with a hint of gold, and the fresh fruity aromas rise from the glass. There's a definite peach scent, with something that reminds me of citrus blossom. The taste is also fruity, again with prominent peach flavours. I also got a definite citrus note from the taste, which I initially pinned as lemon, then changed my mind to lime, then back to lemon again. It is a crisp, dry wine with a light finish, leaving the mouth refreshed.

Wine: Cono Sur Carmenere 2008

2009-06-06T19:13:00.004+01:00

(My first stab at a wine review. In my efforts to remain neutral, I didn't read the label before writing. However, I do have something of a predilection for Chilean reds).

The first thing that struck me about this wine is the vibrant purple-red colour, indicating a young wine (of course, the 2008 vintage was a big clue here!) After that, the aromas; plummy richness with a hint of autumn berries. The wine has a very smooth, full-bodied taste, with the characteristic Chilean warmth; hints of spice. The taste lingers, giving the wine a wonderfully long finish, which can best be described only as sublime.

STLSearch

2009-06-01T16:48:00.004+01:00

I almost don't want to post this, since it'll inevitably spawn a number of misinformed, misguided and generally miserable comments. Oh well.

Today I wrote a Python program called STLSearch. It allows you to quickly bring up the C++ standard library documentation (courtesy of Rogue Wave) for a standard library class or header. There is more information here and it can be downloaded here.

As I said on the detailed page, I know it's bad to call the C++ standard library the "STL", but I wanted a nice snappy name for the program. StandardLibrarySearch.py just doesn't sound so good. So if you don't like it, tough luck.

Also, I'm aware that this is probably completely irrelevant, reinventing the wheel, has been done a million times before, integrated with bigger and better things, is a part of everyone's favourite development environment, etc. etc. I don't care. I find this useful, as I do most of my programming in vim, and I am publishing it in the hope that someone else might find it useful too. If you don't have anything nice to say, don't say anything at all.

GCC vs. ICC: Error Messages

2009-06-01T11:53:00.007+01:00

Here's an error message I was recently confronted with from GCC (g++, actually):

g++ -I. -Wall -O2 -Wextra -Werror -std=c++98 -pedantic-errors -Wfatal-errors -Wwrite-strings -ftrapv -fno-merge-constants -fno-nonansi-builtins -fno-gnu-keywords -fstrict-aliasing -c -o HitMap.o HitMap.cpp
HitMap.cpp: In function ‘std::ostream& KDTree::operator<<(std::ostream&, const KDTree::HitMap&)’:
HitMap.cpp:12: error: no match for ‘operator<<’ in ‘std::operator<< [with _Traits = std::char_traits<char>](((std::basic_ostream<char, std::char_traits<char> >&)((std::ostream*)KDTree::operator<<(((std::ostream&)((std::ostream*)ost)), ((const KDTree::HitPoint&)(& i.std::_Rb_tree_const_iterator<_Tp>::operator-> [with _Tp = std::pair<const KDTree::HitPoint, KDTree::HitInfo>]()->std::pair<const KDTree::HitPoint, KDTree::HitInfo>::first))))), ((const char*)": ")) << i.std::_Rb_tree_const_iterator<_Tp>::operator-> [with _Tp = std::pair<const KDTree::HitPoint, KDTree::HitInfo>]()->std::pair<const KDTree::HitPoint, KDTree::HitInfo>::second’
compilation terminated due to -Wfatal-errors.

And here's ICC's (the Intel C++ compiler) take on the same error:

icc -I. -c -o HitMap.o HitMap.cpp
HitMap.cpp(12): error: no operator "<<" matches these operands
operand types are: std::basic_ostream<char, std::char_traits<char>> << const KDTree::HitInfo
ost << i->first << ": " << i->second << ", ";
^

compilation aborted for HitMap.cpp (code 2)

Now, the GCC error is detailed, sure. It also makes it nearly impossible to figure out, at a glance, where the actual error is (in this case, I had no operator<< defined for the HitInfo class).

The moral of the story? If you can't be bothered to decypher long, cryptic GCC errors, try ICC. Chances are, it'll be more direct.

I wonder why GCC can't do this? Sometimes, having all that extra information is useful. Mostly, it just gets in the way. Perhaps there should be a command-line option to make gcc display more / less verbose C++ errors. (And yes, I know there are tools to filter standard library and template instantiation related errors to reduce the verbosity, and that there are a whole host of issues surrounding this, and—as I said above—the verbosity is sometimes really useful)

Anyway, just thought I'd make a post about this. Mostly as a distraction from having to figure out the bigger (runtime!) error in this damn program.

Bash History Editing

2009-05-14T21:35:00.002+01:00

This post sets a challenge; the goal is to come up with a bash one-liner (which can include multiple statements, pipes, etc. to arbitrary level) which removes "incriminating" lines from the .bash_history file, without looking too incriminating itself.

For example, you sneak onto someone's system and type xhost + so you can run X apps on their display, but you want to hide how you worked your magic, so they can be in awe that bit longer. Of course, I don't condone such activities, but it's an interesting hypothetical. The problem is to remove the line of interest (xhost +) from the .bash_history file, without leaving something equally obvious; like vim .bash_history

One way, obviously, would be to simply run vim, then open the file within vim rather than on the command-line. I'm looking for sneakier approaches; things that will look really impressive, but almost incomprehensible, to the untrained eye. Here's my first stab at a one-liner to remove xhost + (amongst other things) from .bash_history... I'd like to see what any readers can come up with; feel free to leave your offerings in comments to this post.

VAR=".$(basename $SHELL)_h*" sed "/x.+\+/d" -i $VAR

Rules:
1. All offerings must be valid bash commands or sequences thereof
2. All offerings must work when typed into a single line on the bash prompt
3. Nothing destructive, please; don't do anything more than remove a few lines from .bash_history
4. While you can probably come up with a million variants on the "open it with vim" approach, these are trivial. I don't want to see those.
5. If you're replying to someone else's comment, try to remember to treat them with the respect and decency you would expect to receive. In other words, keep things polite, and clean.
6. The more obscure, obfuscated or oblique, the better.

On Traffic Lights and the BBC (Unrelated)

2009-04-17T09:02:00.003+01:00

Today, a rant. Tomorrow, the world?

A couple of weeks ago, I read James May's column in the Telegraph Motoring section, on Saturday. In it, he mentioned the apparent increase in red traffic lights. This morning I experienced the same phenomenon; a junction which does not usually incur much of a wait was suddenly the greatest bottleneck of my journey. The lights were on red when I approached, turned green just long enough to let two cars through, then went back to red for an inordinately long amount of time. When I eventually got through, the next set of lights were, of course, also red.

I wouldn't mind, but while those lights were red, I had a clear view of the road in front. Was it full of cars from other directions? No. Was there a queue as those cars filtered onto the next road? No. Was there any reason at all for the lights to be on red? No.

Meanwhile, sitting in the car, I was listening to BBC Radio 4. I'd already had the BBC News channel on the TV while I had breakfast, so I wasn't surprised to hear their story that MI5 are looking for a 'chief scientific advisor'. Fair enough.

Now, is anyone else getting annoyed by the BBC's tendency to sensationalise mediocre news by completely ignoring the facts (or, in this case, common sense)? I sure am. It happened with the LHC switch on (something I wrote about here), and it's happening again today. The BBC, on the TV, the radio and their website (article here) open with a comparison to 'Q', the guy who produces all the gadgets for James Bond. Is this what a scientific advisor really does? No. Is this even remotely within the realm of possibility? No. Is there any justification for bringing this fantasy into a news story? No. Does it get the attention of the public, who know no better? Yes. And of course, that's all the BBC want. To hell with reporting factual news, or trying to inform the public about what scientists really do. They just want a good story.

Granted, somewhere down the page of the online article, well beyond the point where most people stop reading, there's a quote from someone saying that it is unlikely the post will involve developing weapons systems for Aston Martins. Well, no shit. What a surprise!

It is this kind of 'journalism' I didn't used to expect from the BBC, but now it seems to be commonplace. Once upon a time, they reported the news. Now, I don't know what they do.

First Barbecue!

2009-04-14T23:11:00.003+01:00

Today marks my first barbecue of 2009; pork chops in a tomato, basil and oregano marinade, with barbecued potatoes (new potatoes, boiled first, then covered in garlic and olive oil, and finished on the barbecue, skewered) and vegetables [and of course, a nice bottle of wine—Vacqueyras—to go with it all].

Barbecued food tastes incredible, and I always enjoy this time of year, when it's warm enough and dry enough—not to mention light enough—to have a barbecue, and maybe even to eat outside!

So, I thought I'd make this blog post to announce that summer is on the way! That is all.

PDF and the Internet

2009-03-16T15:29:00.003Z

I have become increasingly annoyed of late with the number of people (including supposedly professional organisations) who publish documents on the Internet in formats not accessible to a large number of people. I refer, of course, to the prevalence of MS Word and PowerPoint documents online.

The solution is simple, the word `portable' is even in the name of the Portable Document Format. It is also easy to make PDF files these days; either directly from Microsoft Office using their Save as PDF extension, using Acrobat to convert files to PDF, or using PDF Printers and similar. There really is no excuse for providing read-only information in formats other than PDF, online.

The advantages to PDF are numerous; the primary advantage being convenience! I can read PDF files directly from my web browser with the Adobe Reader plugin. Furthermore, I can open and print PDF files on Windows, Linux and Mac OS X. The same cannot be said of MS Word documents, or PowerPoint presentations.

I could almost understand if the people putting such documents online were individuals with few resources, but when large, professional organisations provide only Word documents as a source of information, you have to wonder how professional they really are.

Update

2009-03-15T17:47:00.003Z

Well it has been almost a month since I wrote anything here. In fact, my last post was a day before I spent the weekend in Oxford and London, back in February.

Since then, I've mostly been working, going to the gym, and sleeping. I'm currently working on tracking and reconstruction, a part of which involves creating a piece of software in C++ to follow a particle track (parametrised by start coordinates, end coordinates and track energy or dE/dx) through a 3D volume made up of small cubic cells, depositing energy in each cell; the amount of which depends on the length of the track through that cell.

If you're thinking that sounds a lot like GEANT, you're right. Only it doesn't need to be as complicated as GEANT. Still, it doesn't work very well right now, and I spent the best part of last week chasing down bugs (actually, chasing bugs round in a circle!)

Other than that, I'm reading a lot about track reconstruction, pattern recognition, track finding and fitting approaches from the visual scans performed in bubble chambers to the Kalman Filter and the great number of approaches it has spawned. I need to look at the LHC track finding (well, ATLAS and CMS, where they have a very high density of tracks, less so LHCb I think) a little more than I already have, too.

Outside of the office, I've really only had time to cook, eat and sleep. Such is the life of a Ph.D. student!

Writing Scientific Documents Using LaTeX: Updated

2009-02-19T14:30:00.003Z

I submitted a fifth edition (5/E) of `Writing Scientific Documents Using LaTeX` to CTAN yesterday (and another version today which fixes a couple of minor mistakes in the examples given).

As usual, the document and associated files can be found at:
http://www.ctan.org/tex-archive/info/intro-scientific/

Programming Languages

2009-02-12T10:58:00.005Z

It seems that I "collect" programming languages! Not content with learning a single language, or one of each "kind", I play with a number of different languages and try to find the best (or most interesting) one for any given task.

If I'm creating something that has to be done reasonably quickly, I stick with the ones I know best: C, C++, Python and PHP. Otherwise, while Python dominates my recreational programming, I can often be found exploring other languages.

Here's a list of languages I know (to some extent) and could write useful software in:

C
C++
Python
PHP
Bash (if you count shell scripts)
Perl (to a lesser extent, and mostly write-only)
SQL (Well it is a 'language'...)
BASIC (Oh yeah!)
JavaScript
Mathematica (A language, according to Wikipedia's List of programming languages)
Maple (Also in the List of programming languages)
sed

And now the ones I have played with and could probably write something useful in, given a good source of documentation:

Java
C#
Visual Basic (Eugh!)
bc
LabVIEW
Tcl

And the languages I've looked at and written a few short trial programs in and intend to learn properly:

Common Lisp
IA32 ASM
Pike
Ruby
Smalltalk

And languages I've had some exposure to but have no intention of using (if I can help it):

Fortran
Pascal

And finally the languages I haven't touched, but want to learn:

Objective Caml
Haskell
Erlang
Lua
PostScript (Yes, really.)
Scheme

The lists of languages I've played with but want to learn more of, or haven't touched and want to learn from scratch are substantially larger than the list of languages I already know or could use with a little documentation, so I have a lot of "work" ahead of me!

EDIT: I've added a list of markup languages I know, in response to a comment on the original post.

Markup Languages:

HTML / XHTML
XML
CSS (Does this count?)
TeX / LaTeX
bbCode (Technically markup, right?)

Hevea

2009-02-11T21:46:00.002Z

Back in 2006 I converted my nmap tutorial to LaTeX, then realised that I wanted an HTML version of it too. This turned out to be rather a difficult goal to achieve, as the only software I could find to do this at the time was latex2html.

latex2html generates multi-page output, and doesn't do too well with some of the more complicated LaTeX features. Even though I tried to keep the nmap tutorial quite simple, latex2html struggled.

For a single-page version, I wrote a couple of Python scripts to convert LaTeX to bbCode (used on forum systems) and to convert bbCode to HTML. The output was pretty horrific, but it sufficed.

Now, three years on, I've lost my python scripts, and I wanted to regenerate the single-page HTML version; so I Googled "latex to html".

A few clicks later I was at the HeVeA homepage; a LaTeX to HTML translator written in Objective Caml.

Moments later, I'd downloaded a Slackware 12.0 package for ocaml from linuxpackages.net, installed it, built HeVeA, and modified my nmap-tutorial makefile to run HeVeA on the TeX file.

My first surprise was that it all just worked; no errors and only one warning produced by HeVeA. The second surprise was that the output looks so good!

So, if you want to generate HTML from LaTeX files, I can highly recommend taking a look at HeVeA. Furthermore, I've added yet another language (Ocaml) to my list of languages to learn (LISP is still top of that list...)

Nmap-Tutorial.com: Site Down

2009-02-09T14:43:00.002Z

At the end of last week, I took down the server hosting the nmap-tutorial.com website. At this time, I have no plans to put it back up, since I cannot justify the cost of running such a server.

I would like to find a home for the Nmap Tutorial, but I'm not prepared to pay for it any more. I'm looking into one or two possibilities, but for the time being the site will remain down. If you'd like to offer to host it, by all means get in touch!

In the mean time, the tutorial is still available as a thread on Security-Forums.com and a number of other mirrors.

Mathematica Home Edition

2009-02-06T11:38:00.002Z

It seems that Wolfram have created a version of Mathematica 7 targeted at home users. It cannot be used for commercial, nonprofit, government or academic purposes, and comes with a price tag of $295.

From their press release,

Mathematica Home Edition contains the full functionality of the recently released Mathematica 7, which integrates powerful technologies including image processing, parallel computation, and rich data on astronomy, geography, life sciences, and more. Mathematica Home Edition also includes Mathematica's innovative user interfaces that make it possible for new users to achieve high-impact results, like dynamic applications, instantly.

I thought I'd post this in case anyone wants to check it out, but my thoughts are that the license is somewhat limiting; it'd be useful for personal projects and not much else, and is still prohibitively expensive for such use!

Anyone who really needs the functionality of Mathematica is likely to have access to a site-license or a huge budget, but if you don't, and you really want to play with Mathematica (as long as it's not for academic, commercial, nonprofit or government use; which doesn't leave much!) then the Home Edition is a "mere" $295.

Merging PDF Files in Linux

2009-02-03T15:09:00.004Z

The following bash script uses gs (GhostScript) to merge several PDFs into one file. The basic command is

gs -dBATCH -dNOPAUSE -q -sDEVICE=pdfwrite -sOutputFile=out.pdf in1.pdf in2.pdf ...

To use the bash script, chmod +x it and run as
./mergepdf output.pdf in1.pdf in2.pdf ...
(In other words, it just saves you remembering all the gs options!)


#!/bin/bash

output=$1
echo "Output: $output"
shift
echo "Input: $@"
echo -n "Merging PDFs..."
gs -dBATCH -dNOPAUSE -q -sDEVICE=pdfwrite -sOutputFile=$output $@
echo " Done."

Note that this script does no error checking; in particular it doesn't make sure you passed it any options, and it doesn't give you any hints as to how to use it. These are fairly simple things to add, though, so I'll leave them as an exercise to the reader!

Constness

2009-02-01T17:08:00.005Z

Whenever I talk to someone about C++, the topic of const qualifiers almost inevitably arises. It seems many people don't understand what the const qualifier does, and with good reason; const can mean different things in different places. Fundamentally, though, it means that you're telling the compiler "this cannot be changed, it's a constant".

So, I'm now going to go through almost every 'type' of constness you're likely to encounter, explaining what it does, and perhaps why, in the context of several C++ declarations.

const int a = 4;
Here, the const means that the variable a (assigned the value 4) cannot have its value changed; that is, once 4, it's always 4. Subsequent attempts to assign a new value will be compile-time errors. This is perhaps the most basic form of constness, and is relatively uncommon as it has limited use beyond defining numerical constants, string literals, etc.

const int *b = &a;
This form of const declares a pointer to a const int. That means that the value pointed to is const, and you cannot change it using, say
*b = 5;
However, the pointer itself remains non-const, so you could do something like:
const int c;
b = &c;
which is valid, and makes b point to c instead of a. Furthermore, you can make a pointer to a const int point to an integer which isn't implicitly const. All this means is that you cannot change the value of the integer through the pointer (even though you can change it directly). Consider this code:
int a = 4;
const int* b = &a;
*b = 6; // Compile error!
a = 6; // Perfectly acceptable

int * const b = &a;
Given code such as,
int a = 4;
int c = 6;
int * const b = &a;
then b is a const pointer to a (non-const) int. You can change the value of the integer pointed to, e.g.
*b = 7; // sets a to 7
but you cannot change where b points:
b = &c; // Compile error!

This is an important distinction! a const before the asterisk makes the type pointed to const; a const after the asterisk makes the pointer itself const.

const int * const b = &a;
e.g.
int a = 4;
int c = 6;
const int * const b = &a;
Here, the value of the variable pointed to by b cannot be changed through b, and the location b points to also cannot be changed:
*b = 5; // Compile error!
b = &c; // Compile error!

const member function
Consider the following class:
class Employee {
public:
Employee(const std::string& name) : e_name(name) {};
const std::string& getName() const { return e_name; };
void setName(const std::string& name) { e_name = name; };
private:
std::string name;
};

Note that the const qualifier is used liberally here. Note also that references are used extensively. Although this post isn't really about references, I'll mention them briefly in what follows.

Now, we can construct an Employee object as follows:
Employee me("Andrew");
we could also construct it using something like
std::string name = "Andrew";
Employee me(name);

Notice that the constructor takes a reference to a const std::string object. The reference is so that we avoid copying the string in the function call, and the const ensures the caller that we're not going to change it. Here, all we want to do is construct our internal e_name string from the one passed in. If a reference (or pointer) was not used, we'd copy the string once when making the function call, then copy it again when constructing the internal e_name string. The const qualifier makes sure we can't forget what we're doing, and try to directly modify the string passed in by reference. That'd cause a compile-time error.

Next, we'll take a look at getName(). This returns a reference to a const std::string. Again, the reference is because we have the name stored internally, so we want to avoid making copies by letting other people use this directly. Of course, we still want to have control over what changes can be made to the e_name variable, so we return a const reference, preventing anyone "outside" the class from modifying the string! Finally, notice that there is an additional const qualifier after the function name and parameter list (but before the function body)? This tells the compiler that this member function does not (indeed, cannot) modify any of the member variables. This is useful for two reasons:

1. get... methods should be logically const, and shouldn't change the internal state of the class (the exception is reference-counted resource handles, but see the mutable keyword!)

2. It allows us to call such a method if we only have a const Employee, or a reference/pointer to a const Employee.

You see, if you have a const object, the only methods (member functions) you can call on it are those marked const! This is the const that comes after the function declaration but before the definition of the body of the function. Of course, usually you'd separate the declaration into a header file and the definition into a source file, so the distinction is clear. Here, to conserve space, I wrote the definitions inline.

And, last but not least, notice the setName() function. It takes a reference to a const std::string, for the same reasons as the constructor. It doesn't have to return anything, but we can't mark this member function as const, because it changes the internal state -- it changes the value of e_name!

Thus concludes our tour of constness. Note that you can combine return types, argument types and the const modifier to a member function, resulting in a very large number of occurrences of the word const in a single line, for example here is a (member) function which returns a const int, takes a reference to a const string, and a const pointer to a const char, and does not modify the state of the object upon which it is called:

class Stuff {
...
const int getIndexOf(const std::string& str1, const char * const str2) const;
...
};

My recommendation? Make use of the const qualifier as much as possible. It will help you in the long run and will make your code safer (more resilient against unexpected behaviour).

Multicast UDP Host Discovery in Python

2009-01-24T18:02:00.005Z

I'm in the process of writing a file transfer application in Python (don't ask; I don't like to run services such as FTP, SMB, etc. and a lot of people on Windows machines don't want to use SSH...) and I wanted the application to be able to automatically discover other instances of itself on the network.

I chose Python in part because I love it, and also because it's easy to get running on Windows, Linux and OS X, the three OS' I'm targeting here. With the help of the Twisted networking engine, UDP multicast is pretty straightforward. Here's a working Python/Twisted class for host discovery:

#!/usr/bin/env python

# Transpyre
# Network File Transfer Program

from twisted.internet.protocol import DatagramProtocol
from twisted.internet import reactor
from twisted.application.internet import MulticastServer
from twisted.internet import task
import config

class HostDiscovery(DatagramProtocol):
 def __init__(self):
  self.addresses = []
  self.ping_str = 'PING Transpyre:' + config.nodename
  self.pong_str = 'PONG Transpyre:' + config.nodename

 def startProtocol(self):
  print 'Host discovery: listening'
  self.transport.joinGroup(config.multicast_addr)
  # self.ping() # Not needed any more; the refresh loop calls self.ping() on initialisation!
 
 def refreshList(self):
  print 'Refreshing list...'
  self.addresses = []
  self.ping()
 
 def ping(self):
  self.transport.write(self.ping_str, (config.multicast_addr, config.udp_port))

 def pong(self):
  self.transport.write(self.pong_str, (config.multicast_addr, config.udp_port))

 def datagramReceived(self, datagram, address):
  # Check if we received a dgram from another Transpyre instance
  parts = str(datagram).split(':',1)
  if len(parts) == 2:
   print datagram
   if parts[0] == 'PING Transpyre':
    # Send a PONG
    self.pong()
   if parts[0] == 'PONG Transpyre':
    try:
     # Try to find address in list
     idx = self.addresses.index(address[0])
    except:
     # Add it if it isn't there!
     self.addresses.append(address[0])
  print self.addresses

discovery = HostDiscovery()
reactor.listenMulticast(config.udp_port, discovery)
refresh = task.LoopingCall(discovery.refreshList)
refresh.start(config.refresh_delay) # Refresh list at rate defined in config.py
reactor.run()

It sends a UDP packet to the multicast address on port 8192 (by default) with the text 'PING Transpyre:some host name' every 30 seconds and collects the responses ('PONG Transpyre:some other host name'). Each time, the current address list is cleared and regenerated from the packets returned. The packet payload (the text) is checked to determine if the multicast event was indeed another instance of this application, and the PING/PONG mechanism prevents excess packets flying around in some kind of chain reaction!

Here's config.py which is required by the above file

# config.py
# Edit this to reflect your network settings
# Possibly the GUI will edit this for you...

# Node name - reported to other hosts
nodename = 'Loki'

# Download directory
downloads = '/stuff/Transpyre'

# Refresh address list after every N seconds
refresh_delay = 30

# IPv4 Multicast Address - Don't change this if you want it to work with other standard installations
multicast_addr = '224.0.0.1'

# Transpyre UDP Port - Don't change this if you want it to work with other standard installations
udp_port = 8192

# Transpyre TCP Port - Don't change this if you want it to work with other standard installations
tcp_port = 16384

Some Lists

2009-01-24T17:41:00.003Z

After reading this post over at xkcd I decided to try to find some phrases of my own with no results on Google (at the time of posting).

Here goes...

"passenger-side toothbrush holder"
"novel fossilised cow dung"
"airbags are for decoration only"
"marketing a new approach to personal hygiene"
"inconclusive proof for the existence of dog"
"positively neutral in charge"
"blind man saves life of guide dog"
"She'll know when it's time to feed the gimp"
"Nine hundred and ninety nine point nine recurring red balloons"
"Golden Retriever retrieves gold"

And some things I'd hoped would have no hits, but surprisingly did...

"Smarter than your average raptor"
"unimaginably imaginary"
"Poodle in a revolving door"

Happy Googling, everyone!

Soldering LEDs!

2009-01-13T19:16:00.003Z

I've spent most of the last two days soldering (then re-soldering) LEDs and trying to get a power supply working.

What I wanted to do was to get 16 white LEDs all switched on at the same time. Easy, right?

For what I'm doing, I needed about 35cm between each LED, so after acquiring the LEDs in the first place (RS part number 654-2291) I soldered them all together into a single chain with 40cm wires between each one, taking care to connect the LEDs the correct way round!

I wanted to stick a single resistor onto the end of this chain and push 10 mA through the LEDs (the nominal current each LED requires). The threshold forward voltage per LED is 2.7 V, which gives a minimum power supply voltage of 16 * 2.7 = 43.2 V plus the voltage drop across the current-limiting resistor. Since many (commonly available) power supplies don't provide voltages that high, I had to split the chain into two strings of eight, and drive those two strings in parallel. Using a power-supply that provides a constant current (or provides current-limiting on its own), there should in theory be no need for a resistor in the circuit, since I could just set the power supply to go no higher than 10 mA. However, it was suggested (by the people in the electronics workshop) that I use a 100 Ohm resistor on each branch to protect the LEDs should anyone try to push more current through!

So, the circuit looks something like this:

and is currently powered by a 0 - 32V power supply with a current limit set at 0.01 A (10 mA). Strangely enough, the LEDs switch on at around 20 V on the power supply (and certainly less than the expected 21.6 V), which leads me to believe the stated threshold voltage is not quite accurate! Of course, the current drawn appears to also be much less than 10 mA, and for the purposes I require, that the LEDs are on (and that the brightness can be controlled by adjusting the supply voltage) is sufficient!

A success, so far. Tomorrow I have to try to paint some cardboard black.

sed Magic: Removing Lines

2009-01-12T17:29:00.007Z

Problem: Remove lines from a file f2 if they contain a word in file f1.

Solution:

for word in `cat f1`; do
sed -i "/$word/d" f2;
done

sed is a stream editor, the -i option performs the editing in-place (i.e. in the original file).

This came up in the office today and I thought the solution was rather neat (if equally straightforward). Alternatives in most scripting languages (except maybe Perl?) are longer and/or uglier.

Some Random Stuff

2009-01-11T23:36:00.003Z

I don't generally make posts about non-physics-or-programming-related topics, but I thought I'd just outline some of the things that I've done, bought or considered over the last few months.

I think the main thing has to be my new monitor. I used to run dual 19" Samsung TFTs at home. One of them died in November, so I decided to treat myself to the Samsung T240, a 24" widescreen. I'm now using that as my only monitor, until I can get enough desk space to add a second one, and a graphics card that will cope with both! The surviving 19" is now at my parents' house so I can use it when I go home for the weekend (I took it back when I went for Christmas), so I'll probably get either another widescreen 24" (the same as this one) or two smaller monitors to go either side of it and have a triple-monitor setup, next.

Around the same time, I got Amazon Prime. This amazing offer gives unlimited free next day delivery on anything sold by Amazon directly (not marketplace sellers) and costs only £49 a year. Since getting Amazon prime, I've more than made back that much in saved delivery costs (but of course, Amazon win too since I've bought a lot more stuff than I usually would have!) If you spend a lot of time and money at Amazon, I highly recommend Prime. The service is amazing, and the ability to order something and have it arrive the next day (Saturdays included) without paying for delivery is great. Amazon is now my first port-of-call when buying something online.

Since I stayed on at Warwick to do a Ph.D. a large number of people I know have moved on to other places, and those that remained at Warwick have become more broadly geographically distributed. As a result, I'm often quite bored these days. Especially during the evenings and weekends, although to be perfectly honest I'm not sure I have all that much time to do things anyway. Given that I've decided to (try to) go swimming every morning before work, and to the gym three times a week (Tuesdays, Thursdays and Saturdays) there really isn't all that much time for anything beyond eating and sleeping in the evenings, and doing cleaning and washing at the weekends. I guess this is working life!

Shortly before Christmas I bought a wired USB Apple Keyboard. The main reason was because I was taking my laptop home for Christmas, but not my deskop (my laptop is an Apple Macbook). When the keyboard arrived, however, I liked it so much that I'm now using it as my primary desktop keyboard at home. It's a little weird because I haven't (yet) remapped the keys to correspond to the Apple layout, so the labels don't necessarily correspond to the symbol that is generated when you press a key, but since I hardly ever look at the labels, it doesn't matter enough for me to have been sufficiently annoyed to invest the time playing with xmodmap.

I think I mentioned in a post a few days/weeks ago that I was considering getting rid of the hosting for Nmap-Tutorial.com. I still haven't decided what to do about that, but given that I'm considering a move to Kenilworth, a holiday in the U.S. and several weekends in European cities, I should probably start to look at reducing my unnecessary expenses.

On that note, I've been looking into savings accounts at various banks, especially tax free savings such as ISAs. They often have conditions about how much (or little) interest you get paid when you do or don't make a withdrawal in a given month, or about minimum payments, access and notice. To that end, if any reader of this blog has comments about savings accounts they have seen (or use personally) I'd love to hear them. Basically I'm looking for something into which I can stick whatever I think I can afford every month, then withdraw enough to pay for a week in America, or whatever.

And on that note I think I'll wrap up this post. I just thought it would be good to explain what's going on in the real world, for me, amidst the almost certainly confusing series of posts about quantum field theory, and the random posts on programming and playing with Mathematica.

The Road To Quantum Field Theory

2009-01-05T14:31:00.002Z

Dedicated readers will know that I've been presenting topics from theoretical physics with the goal of leading from Lagrangian mechanics all the way to quantum field theories. There were a few posts I made on special relativity and relativistic quantum mechanics before the 'QFT Series' of posts really began, and a number of other posts interspersed with those in the series. The most recent post introduced a mathematical notation useful in special relativity and relativistic quantum mechanics, notably four-vectors, tensors, etc.

Near the beginning I presented a list of topics I was going to cover. At this stage, I think it is appropriate to present a list of topics I've already covered, with links to the posts covering them, then to briefly outline the steps which still remain.

Pre-Series posts on Special Relativity, Quantum Mechanics and Particle Physics

The 'Road to QFT' Series

Out-of-Series Posts

Basic Relativity

Roadmap

Relativistic Electrodynamics (Covariant Formulation of Maxwell's Equations)
Relativistic Quantum Mechanics: From Schrödinger to Klein-Gordon
Relativistic Quantum Mechanics II: The Dirac Equation and Spin
Symmetries and the Gauge Principle
Quantum Electrodynamics (QED)
Aside on Quantum Chromodynamics (QCD)
SU(2) and Electroweak Unification
Aside on Superconductivity (maybe...)
Spontaneous Symmetry Breaking and the Higgs Mechanism

You may notice that the imminent parts of the roadmap are becoming more clearly defined, while the still-distant final topics are broad and loosely defined. For example, there's likely to be more than one post on QED, but until I get there I don't know how many posts it will take!

The Mathematics of Special Relativity

2009-01-03T16:49:00.013Z

As promised, here is a post outlining the mathematical objects used in special relativity (SR). I'm not going to present fully formal definitions, and I'll almost certainly miss out some things I'll later rely on to explain more advanced concepts, but I'll get to those when the time comes. As for the formal definitions, find a good text on special relativity.

Contravariant Four-Vectors
A contravariant four-vector is a four-component object which transforms according to the rule

x^{\prime\mu}=L^\mu_\nu x^\nu

where L is the Lorentz transformation (LT) matrix.

For a transformation from a frame S to a frame S' moving with respect to S with velocity v in the positive x direction, define

\begin{align*}
\beta &= \frac{v}{c}\\ \gamma &= \frac{1}{\sqrt{1-\beta^2}}\\ L &= \left(\begin{array}{cccc}\gamma&-\gamma\beta&0&0\\ -\gamma\beta&\gamma&0&0\\0&0&1&0\\0&0&0&1\end{array}\right)\end{align*}

For example, the space-time four-vector

x^\mu = \left(\begin{array}{c}ct\\x\\y\\z\end{array}\right)

transforms as

x^{\prime\mu} = \left(\begin{array}{c}\gamma ct-\gamma\beta x\\-\gamma\beta ct+\gamma x\\y\\z\end{array}\right)

= \left(\begin{array}{c}c\gamma\left(t-\frac{vx}{c^2}\right)\\\gamma\left(x-vt\right)\\y\\z\end{array}\right)

i.e. the usual form of the SR Lorentz transforms.

Contravariant four-vectors are represented with a greek-letter index in the superscript position. In particle physics and quantum field theory it is conventional for time to be the zeroth component, as in all of my posts. In some areas of relativistic quantum mechanics, one will see time as the fourth component. This alters the definition of the LT matrix and the metric tensor, but of course the physics remains the same. I will always use time as the zeroth component in these posts.

Covariant Four-Vectors and Rank-Two Tensors
In addition to contravariant four-vectors, there are mathematical objects called covariant four-vectors (sometimes also called one-forms). It is possible to convert between covariant and contravariant vectors using an object known as a metric tensor.

Covariant vectors are written with a single greek-letter index in the subscript position. Covariant and contravariant vectors are both rank-one tensors; objects where the rank is defined by the number of indices (equivalently the number of times one must apply the Lorentz transforms).

A metric tensor defines the structure of the space-time in which we are working. For special relativity, we use Minkowski space, and the metric tensor is defined as g:

g = \left(\begin{array}{cccc}1&0&0&0\\0&-1&0&0\\0&0&-1&0\\0&0&0&-1\end{array}\right)

g is a rank-two tensor representing the structure of space-time. In general relativity (GR) different metric tensors can be used for different space-time geometries.

The metric tensor allows us to convert contravariant vectors into covariant vectors, and vice versa:

\begin{align*}
x_\mu &= g_{\mu\nu}x^\nu\\
x^\mu &= g^{\mu\nu}x_\nu\end{align*}

Note that in these cases,

g_{\mu\nu}

has both indices either raised or lowered, whereas in the Lorentz transforms, L has one raised index and one lowered index. The positions of the indices determine the result (covariant or contravariant) and whether the result is physically meaningful. The Einstein summation convention states that one should sum over repeated indices where one is raised and the other is lowered.

In Minkowski space, the conversion between co- and contra-variant vectors is simple:

x^\mu = \left(\begin{array}{c}ct\\x\\y\\z\end{array}\right)

x_\mu = \left(\begin{array}{cccc}ct&-x&-y&-z\end{array}\right)

Note that this leads to the useful property that the contraction of a covariant and a contravariant four-vector is a Lorentz-invariant scalar quantity. For example, using the space-time four-vector:

x^\mu x_\mu = x^\mu g_{\mu\nu}x^\nu = c^2 t^2 - x^2 - y^2 - z^2

which should be recognisable as the Lorentz-invariant space-time interval of special relativity.

This property holds for all four-vectors and in all inertial reference frames.

Rapidity
Special Relativity is sometimes formulated in terms of rapidity

\theta

rather than velocity v.

Here,

\mathrm{tanh}\,\theta = \beta

\Rightarrow \mathrm{cosh}\,\theta = \gamma = \frac{1}{\sqrt{1-\beta^2}}

.

The Lorentz transformation matrix is then written (often called

\Lambda

rather than L):

\Lambda =\left(\begin{array}{cccc}\mathrm{cosh}\,\theta&-\mathrm{sinh}\,\theta&0&0\\-\mathrm{sinh}\,\theta&\mathrm{cosh}\,\theta&0&0\\0&0&1&0\\0&0&0&1\end{array}\right)

for the equivalent of

L = \left(\begin{array}{cccc}\gamma&-\gamma\beta&0&0\\-\gamma\beta&\gamma&0&0\\0&0&1&0\\0&0&0&1\end{array}\right)

The usefulness of this approach lies in the fact that while velocities in SR do not add linearly, i.e.

v_{12} \neq v_1 + v_2

Rapidities do add linearly:

\theta_{12} = \theta_1 + \theta_2

Useful Objects
I'm going to list a few of the vectors and tensors which appear often in particle physics, so that when I need to use them in future posts, they won't look completely new. I won't explain much about them here.

Energy-Momentum four-vector

p^\mu = \left(\begin{array}{c}\frac{E}{c}\\p_x\\p_y\\p_z\end{array}\right) = \left(\begin{array}{c}\frac{E}{c}\\\mathbf{p}\end{array}\right)

Four-Current

j^\mu = \left(\begin{array}{c}\rho\\\mathbf{j}\end{array}\right)

(Electromagnetic) Four-Potential

A^\mu = \left(\begin{array}{c}\frac{\phi}{c}\\\mathbf{A}\end{array}\right)

where A is the vector potential of electromagnetism.

Covariant Derivative

\partial_\mu = \left(\begin{array}{c}\frac{1}{c}\frac{\partial}{\partial t}\\ \nabla\end{array}\right)

Faraday Tensor

F^{\mu\nu} = \partial^\mu A^\nu - \partial^\nu A^\mu

F = \left(\begin{array}{cccc}0&-\frac{E_1}{c}&-\frac{E_2}{c}&-\frac{E_3}{c}\\
\frac{E_1}{c}&0&-B_3&B_2\\
\frac{E_2}{c}&B_3&0&-B_1\\
\frac{E_3}{c}&-B_2&B_1&0\end{array}\right)

Note that

F^{\mu\nu}F_{\mu\nu} = 2\left(B^2 - \frac{E^2}{c^2}\right)

is Lorentz invariant!

2008

2008-12-31T12:55:00.005Z

Welcome to this, my final post of 2008. This year has seen a number of important events in science, notably the switch-on of the LHC at CERN (and the subsequent switch-off due to problems with magnets and helium). As far as I'm aware, things should be back online soon (ish), and although I'm not directly involved in anything going on at the LHC, I'm watching it with interest.

In other news, I started my Ph.D. in particle detector research & development, looking specifically at the use of Liquid Argon based detectors. So far, I've written some C++ code to simulate point light sources and attempt to reconstruct their locations, got very annoyed with the ROOT software, figured out how to use ROOT a little better, only to get annoyed with the next slightly-more-complicated thing, and attended a lot of postgraduate lectures on particle physics.

While I initially thought the lectures weren't going to be much use, it turns out I've learned quite a lot from them, and many things I already knew have become much more concrete in my mind, so I have to say they were a success, even if they did take up a huge amount of time!

In addition to all of that, I set out in September and October to make posts on quantum field theory, beginning with classical Lagrangian mechanics, moving through Special Relativity (SR) into relativistic quantum mechanics (RQM) and finally into quantum field theory (QFT). I got about as far as the relativity part before I ran out of time. Hopefully in 2009 I'll make time to continue this (hopefully) thorough treatment of theoretical physics.

If you've been following the posts on those topics, I'd appreciate any feedback, thoughts, comments or suggestions you may have. Are the posts understandable, and at what level? For instance, if you're a physics undergraduate, do you have a decent chance of following the material presented? If you're not a physicist at all, do you get the idea of what's going on, even if you can't follow every step of the way? I have no intention to "dumb down" those posts; I want them to be as in-depth as possible, but I may be convinced to make a parallel series of posts explaining the same material in a less mathematical, more phenomenological way.

As some readers may know, I wrote a detailed (at the time) tutorial for the Nmap port scanner in the years 2003 - 2006. The 2006 version is kept at http://www.nmap-tutorial.com, a site which I maintain and host. Given that I haven't updated it since 2006, I'm debating whether to keep it there at all, given the cost of hosting the site. There are reasons for keeping it, not least of which is that the same server provides a useful dumping ground for files I want to publish on this blog (which aren't images), such as the Mathematica notebook from the Fourier Transform of a Cat post. There are also reasons for getting rid of it, notably the time and money required to maintain a server. I'm not entirely sure what will happen here, if I can find alternative hosting for the Nmap Tutorial, that would be the ideal solution, but that doesn't seem terribly likely at this time.

On the topic of websites, I created a site at http://www.warwick.ac.uk/~phrebd/ to provide information useful mostly to Warwick students and staff, particle physicists and programmers. I wanted to keep certain links and information together in one place and not in a blog format.

And now, on a more personal note, I'm thinking of moving to Kenilworth after June 2009. If any Warwick Ph.D. students are reading this, and wouldn't mind living in Kenilworth (my minimum requirements are somewhere with off-road parking, a decent sized kitchen, and a dishwasher) then get in touch.

Finally, I'd like to wish anyone who bothers to read this a happy new year, and I'll try to post more often in 2009!

Wolfram release Mathematica 7

2008-11-19T09:13:00.003Z

Mathematica 7 has been released today. The release continues the emphasis on bringing together aspects of Mathematica, improving and optimising the way it works, and putting more power and flexibility at the user's fingertips. Better image processing, support for parallel computation and high-performance computing are also mentioned in the Wolfram press release, a section of which reads,

November 18, 2008, Champaign, Illinois--Wolfram Research today
announced Mathematica 7, a major release that accelerates the
drive to integrate and automate functionality as core Mathematica
capabilities, adding image processing, parallel high-performance
computing (HPC), new on-demand curated data, and other recently
developed computational innovations--in total over 500 new
functions and 12 application areas.

I haven't had time to check out Mathematica 7 for myself yet, but I'll definitely take a look as soon as I can.

Basic Relativity

2008-11-10T22:57:00.012Z

I realise I haven't made any posts here recently; I've been far too busy with my Ph.D. and various other commitments. In lieu of the post I wanted to make on special relativity, defining everything carefully and introducing the tensor formulation, I'm going to present some basic concepts of Special Relativity, almost exactly as I taught them to a first year undergraduate examples class. These people had lectures on relativity already, so the stuff below is really a recap and a worked example.

Lorentz Transforms
Consider a coordinate frame S' moving with velocity v in the positive x-direction of a frame S. The Lorentz transforms between the two frames (x, y, z, t in S and x', y', z', t' in S') are given by:


\begin{align*}
x^\prime &= \gamma(x-vt) = \frac{1}{\sqrt{1-\frac{v^2}{c^2}}}(x - vt)\\
y^\prime &= y\\
z^\prime &= z\\
t^\prime &= \gamma\left(t - \frac{vx}{c^2}\right) = \frac{1}{\sqrt{1-\frac{v^2}{c^2}}}\left(t - \frac{vx}{c^2}\right)
\end{align*}

The reverse transforms, from S' back to S, are given by:


\begin{align*}
x &= \gamma(x^\prime+vt^\prime) = \frac{1}{\sqrt{1-\frac{v^2}{c^2}}}(x^\prime + vt^\prime)\\
y &= y^\prime\\
z &= z^\prime\\
t &= \gamma\left(t^\prime + \frac{vx^\prime}{c^2}\right) = \frac{1}{\sqrt{1-\frac{v^2}{c^2}}}\left(t^\prime + \frac{vx^\prime}{c^2}\right)
\end{align*}

Velocity Transformation
If an object is moving with a velocity u' in frame S' (as defined above), we need to determine a transform for the velocity in frame S. Start out with the fact that the velocity in S is

u = \frac{\Delta x}{\Delta t}

and substitute in the transforms above:


u = \frac{\Delta x}{\Delta t} &= \frac{\gamma(\Delta x^\prime + v\Delta t^\prime)}{\gamma\left(\Delta t^\prime + \frac{v\Delta x^\prime}{c^2}\right)}

Divide through by

\Delta t^\prime


u= \frac{\frac{\Delta x^\prime}{\Delta t^\prime} + v}{1+\frac{v\Delta x^\prime}{\Delta t^\prime c^2}}

Since

\frac{\Delta x^\prime}{\Delta t^\prime}=u^\prime


u = \frac{u^\prime + v}{1+\frac{u^\prime v}{c^2}}

which is the velocity transform between u and u'!

Time Dilation
Consider a particle with a lifetime

\tau

in its rest frame (i.e. in an inertial frame in which the particle has zero momentum). The lifetime t, in a frame moving with velocity v is given by:


t = \gamma\tau

Since

\gamma

is always greater than or equal to 1, t is always as large as, or larger than,

\tau

. This is the phenomenon of time dilation and the relationship can be used to determine the time dilation in any inertial frame.

Length Contraction
If the length measured in the lab frame is given by

L_0

, the length seen by a particle moving at speed v in the lab frame is given by:


L = \frac{L_0}{\gamma}

The length is therefore always smaller than (or exactly the same as) the length in the lab frame, and this is the phenomenon of length contraction.

Be careful never to use both time dilation and length contraction in the same (part of a) physics problem! If you're working in the lab frame, use time dilation to alter the particle lifetime, and leave the length as measured in the lab frame. If you're in the particle rest frame, use the particle lifetime as given (since they are defined as lifetime in a particle rest frame) and use length contraction to determine the length seen by the particle from the length as measured in the lab.

The classic atmospheric muon problem, below, illustrates how to use these concepts.

The Atmospheric Muon Problem
Muons are produced in cosmic-ray interactions high in the Earth's atmosphere (at, say 8000m). The muon lifetime is

2.2\times 10^{-6}~\mathrm{s}

.

Muons produced in this way are detected at the Earth's surface. For a muon travelling at 0.998 times the speed of light,
(a) Calculate how far the muon would travel using only classical (non-relativistic) physics.

Here, we use the classical

s = ut

relationship. The maximum distance the muon can travel is given by:


s = 0.998c \tau = 0.998\times 3\times 10^{8} \times 2.2\times 10^{-6} = 658.68~\mathrm{m}

In this situation, the muon cannot reach the Earth's surface.

(b) In the rest frame of the muon, use relativistic physics to show that the muon can indeed reach the surface.

In the muon rest frame, the lifetime is still

2.2\times 10^{-6}~\mathrm{s}

, but the length it has to travel is shortened by length contraction:


L_\mu = \frac{1}{\gamma}L_E

The Lorentz factor is given by:


\gamma = \frac{1}{\sqrt{1-\frac{v^2}{c^2}}} = 15.8

giving the length seen by the muon as:


L_\mu = \frac{8000}{15.8} = 506.32~\mathrm{m}

We already worked out that a muon travelling at 0.998c can travel over 600m before it decays using its lifetime as stated, in part a, so now the muons will make it to the surface of the Earth!

(c) In the Earth frame, use relativistic physics to show that the muon can make it to the surface.

Here, the length is as measured, 8000m, but the lifetime of the muon must be dilated:

t = \gamma \tau

We worked out the Lorentz factor above as 15.8, so now we have:


t = 15.8 \times 2.2\times 10^{-6} = 3.48 \times 10^{-5}

The maximum distance such a muon can travel is then:


s = u t = 0.998\times 3\times 10^8 \times 3.48\times 10^{-5} = 10419~\mathrm{m}

The muon can now travel over 10 km in the Earth frame, well over the 8 km it needs to hit the surface, so again we can detect such muons on the Earth.

This example demonstrates several interesting aspects of relativistic physics. Firstly, that it works! We can detect muons produced in this way, which we wouldn't be able to detect if those relativistic effects didn't occur! Secondly, you can analyse a problem using either length contraction or time dilation, but you need to choose your frame carefully. Don't contract or dilate quantities that were measured in the frame you're using; only alter those measured in other frames. Note that lengths are often (but not always) measured in a lab frame, and that particle lifetimes are always given in the rest frame of the particle concerned!

An Unhelpfulness of Optics Books

2008-10-28T22:29:00.003Z

The collective noun is probably something more like "library" or simply "collection", but I think an Unhelpfulness better describes their applicability to my current work.

As I mentioned in my previous post, I was looking for an optics book which would help me with figuring out how to map a set of point sources onto a close imaging plane, and to be able to reconstruct the positions of those sources.

Having looked through all the books I mentioned in that post, as well as Optics (4th Edition), E. Hecht (Pearson, 2003), I've reached the (in retrospect, fairly obvious) conclusion that none of them had the treatment I wanted.

They deal with wave optics and primarily in the far-field regime, for coherent light sources. In my case, the sources are not coherent, and are pretty damn close to the image plane. I resorted to the principle of superposition and to the inverse square law for light intensity, which appears to be producing adequate results, though I would still be interested in any texts on point source reconstruction from images on nearby planes.

The Quest for an Optics Book

2008-10-19T13:17:00.004+01:00

I'm searching for a book on optics; one which covers topics such as reconstructing point sources from an imaging plane a short distance away from them, where the imaging plane consists of an array of detectors (which do not cover the entire surface). It would also be nice if I could find something that covered things like dealing with edge effects and how best to reconstruct points at the edge of the detector area.

I've been to the University of Warwick Library, and acquired three books on optics. At first glance, these seem to cover some of the introductory topics as well as a few of the things I will need specifically. Since, currently, I know very little about optics (despite owning a digital SLR camera!) I need to find a gentle introduction followed by some very technical stuff.

As I read through these to assess their suitability, I may end up posting brief reviews of each one, describing their content and relevance for my purposes. In the mean time, here are the books I've found so far. If anyone can suggest an alternative, or something which contains exactly what I want, please feel free to leave a comment!

Optical Imaging and Aberrations - Part I: Ray Geometrical Optics, Virendra N. Mahajan (SPIE, 1998)
Astronomical Image and Data Analysis, Jean-Luc Starck and Fionn Murtagh (Springer, 2002)
Introduction to Optics and Optical Imaging, Craig Scott (IEEE Press, 1998)

Contraction of Four-Vectors: A C++ Exercise

2008-10-07T12:02:00.001+01:00

NOTE: The examples as typed have not been compiled or checked, however they are directly based on examples I talked through while teaching C++, so there should not be any particularly glaring errors or omissions. Note that I've avoided talking too much about object-oriented design, advanced concepts or why certain things are done the way they are. This article is intended mainly as an introduction to the abstraction mechanisms available with C++ class-based design, and is not intended to introduce any new concepts to the world of C++. In particular, I do not claim to be the first person to have thought of such concepts (as a number of comments have suggested in the past); I'm merely providing an example for those interested!

Without going into the physics behind it, if we have two objects A and B, each with four components, we can define their contraction as:

A_{\mu}B^{\mu} = A_0 B_0 - (A_1 B_1 + A_2 B_2 + A_3 B_3)

I intend now to use this concept to develop a C++ program from a `first draft' to a finished object-oriented program. This is indented primarily as an introductory-level tutorial for those new to C++ but familiar with many aspects of C or similar languages, so I will not go into a lot of detail about the basics, instead focussing on the C++ specific stuff.

The first approach could look something like this:


#include <iostream>

int main()
{
 // Create some variables
 double a0, a1, a2, a3;
 double b0, b1, b2, b3;
 
 // Prompt and read eight components
 std::cout << "Enter the four components of A, separated by spaces:" << std::endl;
 std::cin >> a0 >> a1 >> a2 >> a3;
 std::cout << "Enter the four components of B, separated by spaces:" << std::endl;
 std::cin >> b0 >> b1 >> b2 >> b3;
 
 // Calculate contraction and display result
 double result = a0*b0 - a1*b1 - a2*b2 - a3*b3;
 std::cout << "A.B = " << result << std::endl;
 
 return 0; 
}

This is a fairly simple program, achieving exactly what we set out to do. However, from the perspective of simplicity, there's a lot of overhead in the main program which detracts from the understanding of what is actually going on. We can make things slightly clearer by writing the contraction as a function; then we'll be able to call it repeatedly from our code, instead of defining what a contraction of four-vectors means at every point it is used:


double contract(double a0, double a1, double a2, double a3,
 double b0, double b1, double b2, double b3)
{
 return a0*b0 - a1*b1 - a2*b2 - a3*b3;
}

Now, our main() looks like this:


int main()
{
 // Create some variables
 double a0, a1, a2, a3;
 double b0, b1, b2, b3;
 
 // Prompt and read eight components
 std::cout << "Enter the four components of A, separated by spaces:" << std::endl;
 std::cin >> a0 >> a1 >> a2 >> a3;
 std::cout << "Enter the four components of B, separated by spaces:" << std::endl;
 std::cin >> b0 >> b1 >> b2 >> b3;
 
 // Calculate contraction and display result
 double result = contract(a0, a1, a2, a3, b0, b1, b2, b3);
 std::cout << "A.B = " << result << std::endl;
 
 return 0; 
}

No major improvement there, but we've moved the contraction code out into a function so that we can reuse it. Still, it's pretty ugly to pass eight doubles into the function. Much nicer would be to pass two four-vectors; after all, that is what the contraction works with, at a higher level:


#include <iostream>

class FourVector {
 public:
  double x0;
  double x1;
  double x2;
  double x3;
};

double contract(FourVector a, FourVector b)
{
 return a.x0*b.x0 - a.x1*b.x1 - a.x2*b.x2 - a.x3*b.x3;
}

int main()
{
 // Create some variables
 Fourvector a, b;
 
 // Prompt and read eight components
 std::cout << "Enter the four components of A, separated by spaces:" << std::endl;
 std::cin >> a.x0 >> a.x1 >> a.x2 >> a.x3;
 std::cout << "Enter the four components of B, separated by spaces:" << std::endl;
 std::cin >> b.x0 >> b.x1 >> b.x2 >> b.x3;
 
 // Calculate contraction and display result
 double result = contract(a, b);
 std::cout << "A.B = " << result << std::endl;
 
 return 0;
}

Things are a little neater now, but still not ideal. For a start, we have public member variables in our class. This is generally considered bad design, so we should make them private and instead provide alternative mechanisms for accessing them. In this case, the class lends itself perfectly towards the use of operator overloading. The idea here is to define what it means to "multiply" two four-vectors together. We can define operator*() as a member-function of the FourVector class:


#include <iostream>

class FourVector {
 public:
  double x0;
  double x1;
  double x2;
  double x3;
  double operator*(const FourVector& b);
};

double FourVector::operator*(const FourVector& b)
{
 return x0*b.x0 - x1*b.x1 - x2*b.x2 - x3*b.x3;
}

int main()
{
 // Create some variables
 Fourvector a, b;
 
 // Prompt and read eight components
 std::cout << "Enter the four components of A, separated by spaces:" << std::endl;
 std::cin >> a.x0 >> a.x1 >> a.x2 >> a.x3;
 std::cout << "Enter the four components of B, separated by spaces:" << std::endl;
 std::cin >> b.x0 >> b.x1 >> b.x2 >> b.x3;
 
 // Calculate contraction and display result
 double result = a * b;
 std::cout << "A.B = " << result << std::endl;
 
 return 0;
}

As you can see, we have now defined enough code that the statement result = a * b; works as expected, producing the contraction of the two. The member function, when called, is bound to object a (in this case), so x0, x1, x2, x3 all refer to a's variables. We have to refer to the variables in b (passed as an argument to the operator) explicitly as b.x0 etc.

We still haven't hidden the member variables. To do this, we need to add a couple more things: a constructor, and a stream operator.


#include <iostream>

class FourVector {
 public:
  FourVector() : x0(0), x1(0), x2(0), x3(0) {};
  FourVector(double i0, double i1, double i2, double i3)
   : x0(i0), x1(i1), x2(i2), x3(i3) {};
  double operator*(const FourVector& b);
  friend std::istream& operator>>(std::istream& in, FourVector& v);
 private:
  double x0;
  double x1;
  double x2;
  double x3;
};

double FourVector::operator*(const FourVector& b)
{
 return x0*b.x0 - x1*b.x1 - x2*b.x2 - x3*b.x3;
}

std::istream& operator>>(std::istream& in, FourVector& v)
{
 in >> v.x0 >> v.x1 >> v.x2 >> v.x3;
 return in;
}

int main()
{
 // Create some variables
 Fourvector a, b;
 
 // Prompt and read eight components
 std::cout << "Enter the four components of A, separated by spaces:" << std::endl;
 std::cin >> a;
 std::cout << "Enter the four components of B, separated by spaces:" << std::endl;
 std::cin >> b;
 
 // Calculate contraction and display result
 double result = a * b;
 std::cout << "A.B = " << result << std::endl;
 
 return 0;
}

There are a few things to note here. First, the default constructor:

FourVector() : x0(0), x1(0), x2(0), x3(0) {};

Constructors are named identically to the class they construct, and have no return type. This particular one takes no arguments, has an initialiser list which sets all the member variables to zero, and has an empty function body -- {};

The other constructor allows one to construct a FourVector object by passing four doubles into the constructor:


double E, px, py, pz;
// Some code which sets E and p goes here
FourVector EnergyMomentum(E, px, py, pz);

Finally, operator>> is a little different. Because the order is stream >> FourVector, we can't make the operator>> a member of FourVector (since it would expect the FourVector object to be passed first!). Furthermore, we need access to the private member variables, or we have nothing to read into. This is a bit of a dilemma, until we use the `friend' statement. This tells the compiler to let a function which is defined outside the class have access to internal class members; precisely what we want here!

You'll notice that the mechanism by which we got here is quite convoluted; we've added constructors and defined operators, and the class is far from complete -- for it to be generally useful we'd need an output operator, other mathematical operations and more. The amount of code has increased, as has the complexity of the overall solution. But look at the code inside main()... this is now considerably simpler. If we ignore comments and don't prompt the user, we can now write exactly the code which corresponds to the mathematical concept:


int main()
{
 FourVector a, b;
 std::cin >> a >> b;
 double result = a * b;
 std::cout >> result >> std::endl;
 return 0;
}

The best part is that the FourVector class is reusable, and sufficiently powerful to allow this abstraction to be carried forward to develop a whole system for tensor calculus, expressed in terms of the mathematical objects involved!

More Mathematica Visualisations: Non-Orientable Foo

2008-10-02T22:51:00.002+01:00

Here are some more random visualisations I created in Mathematica. I uploaded the videos to YouTube a week ago with the intention of posting links here, but never got round to it. The code I used to generate most of these is fairly straightforward. For the Möbius Strip Traversal:

R := 4
w := 1
Animate[Show[
ParametricPlot3D[{(R + s Cos[t/2]) Cos[t], (R + s Cos[t/2]) Sin[t],
s Sin[t/2]}, {s, -w, w}, {t, 0, 2 \[Pi]}, Boxed -> False,
Axes -> False, Mesh -> 3,
PlotStyle ->
Directive[Opacity[0.9], Orange, Specularity[White, 40]],
BoundaryStyle -> Directive[Black]],
Graphics3D[Sphere[{R Cos[v], R Sin[v], 0}, 0.5]]], {v, 0, 2 \[Pi],
0.05}, DisplayAllSteps -> True]

And for the rotating Klein Bottle,

Kx = 6 Cos[u] (1 + Sin[u]);

Ky = 16 Sin[u];

Kr = 4 (1 - Cos[u]/2);

X = If[\[Pi] <>
Kx + Kr Cos[v + \[Pi]],
Kx + Kr Cos[u] Cos[v]];

Y = If[\[Pi] <>
Ky,
Ky + Kr Sin[u] Cos[v]];

Z = Kr Sin[v];

ParametricPlot3D[{X, Y, Z}, {u, 0, 2 \[Pi]}, {v, -\[Pi], \[Pi]},
Boxed -> False, Axes -> False,
PlotStyle ->
Directive[Opacity[0.75], Orange, Specularity[White, 40]]]

Animate[ParametricPlot3D[{X, Y, Z}, {u, 0,
2 \[Pi]}, {v, -\[Pi], \[Pi]}, Boxed -> False, Axes -> False,
PlotStyle ->
Directive[Opacity[0.75], Orange, Specularity[White, 40]],
ViewPoint -> {5 Cos[\[Theta]], 5 Sin[\[Theta]], 5}], {\[Theta], 0,
2 \[Pi], 0.1}, DisplayAllSteps -> True]

Möbius Strip

Klein Bottle
(Technically, a Klein-Bottle immersed in 3D space)

Decadence Factors: Curve Fitting

2008-09-28T11:47:00.008+01:00

Some months ago I made a post about Decadence Factors. These I defined to be the total electricity consumption of a given country, divided by its population (i.e. electricity consumption per capita). The data was obtained from the CIA World Factbook.

Going back over my old posts, a couple of days ago, I noticed that the `top 20' graph in my earlier post had what looked like an exponential decay curve, so I set out to find a fit for the decadence vs. position in list, with the most decadent country sitting at position 1.

I decided to use Mathematica because I've been playing with it again recently, and it proved easy to import the tab-separated data I'd previously used in Gnuplot;

data = Import["decadence.csv", "Table"]

This imports the data into a nested list of the form {{Name,Factor},{Name,Factor},...}. In order to fit the data (or plot it as a list) I needed the decadence factors in a list without the country names. Easily done:

factors = Transpose[data][[2]]

With this achieved, I tried two kinds of fit. The first, a power series up to order

x^{120}

:

fit = Fit[factors, Table[x^n, {n, 0, 120}], x]

In the plot below, the fit curve is shown in red, and the original decadence data points are shown in blue. The power series is an extremely good fit, but almost useless for further analysis.

The second fit I tried was of the form

A\exp(-Bx)

. Mathematica fit the curve with the parameters:
A = 17032.5
B = 0.0272113

The fit command was:

fit2 = FindFit[factors, A Exp[-B x], {A, B}, x]

This wasn't quite such a good fit, but it's close enough for us to be able to do a little further work. Again, the red line is the fit curve, the blue points are the original data.

The nice thing about the exponential fit is that it's easy to invert. That is, given a decadence factor, we can work out an approximation of that country's decadence ranking; their position in the list! For a decadence factor F, and a rank R, the fit provided by Mathematica gives:


F(R) = 17032.5\exp(-0.0272113 R)

Inverting this, we get:


R(F) = \frac{-1}{0.0272113}\log\left(\frac{F}{17032.5}\right)

This approximation would place Sweden, with a decadence factor of 14825.2 in rank 5.10063. Sweden is actually ranked 5th, so the approximation agrees reasonably with the data.

Visualising Random Walks

2008-09-25T23:21:00.003+01:00

I've been playing with visualisation in Mathematica recently; particularly with the Animate[] function, allowing a representation (graphical or otherwise) to be viewed as some parameter, or list thereof, changes within some defined range. I used it to produce the earlier video I posted, showing how the Taylor series approximation of Exp[x] approaches the true Exp[x] function as the number of terms in the expansion increases.

Today I was trying out random walks in Mathematica. They turn out to be fairly easy to code; random number generation is built in, as is a mechanism of accumulating steps into a list. In order to generate a list of points corresponding to a 1000-step random walk with origin (0,0), the following command suffices:

Points = Accumulate[RandomChoice[{{-1, 0}, {1, 0}, {0, 1}, {0, -1}}, 1000]]

In order to visualise the progress of the walk step-by-step, in Mathematica, I used the following Animate command:

Animate[Graphics[{Red, Line[{Points[[1 ;; n]]}]},
PlotRange -> {{-Max[Points], Max[Points]}, {-Max[Points],
Max[Points]}}], {n, 1, 1000, 1}, DisplayAllSteps -> True]

I'll go through exactly what this is doing... The Graphics[...] part produces a 2D Graphics object. In this case, it consists of a line whose points are given by the list produced in the earlier command. The Points[[1;;n]] part allows Animate to work its magic...

By telling Animate that we want to see n run from 1 to 1000, in steps of 1, displaying each step as we go, we obtain a new Graphics object each step, with a line containing the first n points of the Points list we generated earlier. Adding a few options to keep the plot range constant across all of these new visualisations (otherwise, Mathematica would rescale the axes, and the resulting image would `flicker'), we get a video of each step of the walk being added on to the previous steps.

Unfortunately, for some reason I couldn't get this to directly export in that format. I managed to produce an AVI file which showed the Animate box, but didn't animate the behaviour correctly. So instead, I produced a Table[] of Graphics objects, and exported that using one image for each frame of the AVI, using the command:

Export["random.avi",
Table[Graphics[{Red, Line[{Points[[1 ;; n]]}]},
PlotRange -> {{-Max[Points], Max[Points]}, {-Max[Points],
Max[Points]}}], {n, 1, 1000, 1}]]

The result, after transcoding to MPEG-4 to dramatically reduce the filesize, and uploading to YouTube, is included below!

LISP: Taylor Expansion of Exponential Function

2008-09-24T14:23:00.008+01:00

The code below defines a LISP function to approximate

\exp(x)

using a Taylor expansion. I first had to define the functions factorial (x!) and pow (raise x to the power n) in order to define the Taylor expansion of exp(x),


\exp(x) \approx 1 + x + \frac{x^2}{2!} + \ldots + \frac{x^n}{n!} + \ldots

(defun factorial(x)
 (let ((acc 1))
(dotimes (i x)
  (setf acc (* acc (+ i 1))))
  acc))

(defun pow(x n)
 (let ((acc 1))
(dotimes (i n)
  (setf acc (* acc x)))
acc))

(defun taylor-exp(x)
 (let ((a 1) (b x) (n 1) (tolerance 0.00001))
(loop
  (setf a (+ a (/ (pow x n) (factorial n))))
  (setf n (+ n 1))
  (if (< (abs (- b a)) tolerance)   (return (* a 1.0))   (setf b a)))))

And here's a test run comparing the built-in (exp 2) with (taylor-exp 2):

(exp 2)
7.389056
(taylor-exp 2)
7.389055

With a little more fine-tuning of the tolerance parameter, the two could agree to arbitrary precision!

Here's a video produced using Mathematica, showing the convergence of the Taylor series to the true Exp[x]. It was produced using Animate[Plot[ Evaluate[{Exp[x], Normal[Series[Exp[x], {x, 0, n}]]}], {x, 0, 5}, PlotRange -> 200, PlotStyle -> Thick], {n, 1, 12, 2}]

Symmetries, Conserved Quantities and Noether's Theorem

2008-09-24T08:23:00.015+01:00

Following on from my posts on Lagrangian and Hamiltonian Mechanics [1, 2, 3, 4], I'd like to discuss one of the most amazing topics in physics. The consideration of the symmetries of physical laws and how those symmetries relate to conserved quantities, the fundamentally beautiful mathematics lying beneath, and the extent to which we can develop theories of the world around us from such simple concepts; these things continually inspire me and provoke my interest.

I'll start by showing that the formulation of Lagrangian and Hamiltonian Mechanics, thus far, allows us to determine several conservation laws. Consider, for example, the homogeneity of space. Space is homogeneous if the motion (or time-evolution) of a particle (or system thereof) is independent of absolute position. That is, the potential does not vary with absolute position (it can still vary with the vector distance between two particles, as an interaction potential, for example!)

If we make a transformation

\mathbf{r}\rightarrow\mathbf{r}+\delta\mathbf{r}

, then the Lagrangian will also transform as

L \rightarrow L+\delta L

. For a single particle, we can Taylor expand as follows:


L(\mathbf{r}+\delta\mathbf{r},\mathbf{v}) = L(\mathbf{r},\mathbf{v})+\frac{\partial L}{\partial x}\delta x+\frac{\partial L}{\partial y}\delta y+\frac{\partial L}{\partial z}\delta z

Which we can use to write

\delta L = \frac{\partial L}{\partial \mathbf{r}}\cdot\delta\mathbf{r}

\frac{\partial L}{\partial \mathbf{r}}

is a vector quantity; each component is the derivative of L with respect to the corresponding coordinate of r. For a single particle, then,

\frac{\partial L}{\partial \mathbf{r}}=\nabla L

.

Homogeneity of space requires that

\delta L = 0

. Since

\delta \mathbf{r}

is arbitrary (and therefore not necessarily zero), we have that


\frac{\partial L}{\partial q_i} = 0 ~~~~~~ (\star)

This holds only if L does not depend on absolute position, otherwise there would be a contribution

\delta L

from many of the possible choices of

\delta\mathbf{r}

. Spatial dependence of e.g. V(x) implies spatial variation of L, and momentum would not be conserved.

The Euler-Lagrange Equation applies for each coordinate in the vector r. The sum of these Euler-Lagrange Equations (ELEs) means that

(\star)

requires that:


\begin{multiline*}
\frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i} = 0 \\
\Rightarrow p_i = \frac{\partial L}{\partial \dot{q}_i} ~~\mathrm{remains~constant}
\end{multiline*}

We have, therefore, demonstrated the conservation of momentum as a result of requiring translational invariance. That is, any canonical momenta whose conjugate coordinates do not appear explicitly in the Lagrangian are conserved.

Turning once again to time symmetries, let us re-derive the conservation of energy. If the Lagrangian is homogeneous in time, i.e.

L(q,\dot{q})~~\mathrm{not}~~L(q,\dot{q},t)

, then:


\frac{dL}{dt} = \sum_i \frac{\partial L}{\partial q_i}\dot{q}_i + \sum_i\frac{\partial L}{\partial \dot{q}_i}\ddot{q}_i

As L does not depend explicitly on time, there is no term

\frac{\partial L}{\partial t}

on the RHS. Sunstituting

\frac{\partial L}{\partial q_i}

from the ELE,


\frac{dL}{dt} = \sum_i\dot{q}_i\frac{d}{dt}\frac{\partial L}{\partial\dot{q}_i}+\sum_i\frac{\partial L}{\partial\dot{q}_i}\ddot{q}_i = \sum_i\frac{d}{dt}\left(\dot{q}_i\frac{\partial L}{\partial\dot{q}_i}\right)


\begin{multiline*}
\Rightarrow \frac{d}{dt}\sum_i\left( \dot{q}_i\frac{\partial L}{\partial\dot{q}_i}-L \right) = 0 \\
\Rightarrow H = \sum_i \dot{q}_i\frac{\partial L}{\partial \dot{q}_i} - L~~~\mathrm{remains~constant}

The conservation of energy holds for any motion in a non-time-varying external field V(x)

We turn now to the isotropy of space, and show that angular momentum is conserved due to rotational invariance of the Lagrangian. Consider rotation by an angle

|\delta\theta|

(with a direction given by

\delta\theta

) about a vector. For small rotations,

\mathbf{r}\rightarrow \mathbf{r}+\delta\mathbf{r}

, with

\delta\mathbf{r} = \delta\theta\times r

. Each component of the velocity is also transformed by this rotation,

\delta\mathbf{v}=\delta\theta\times\mathbf{v}

.

For a single body, we now impose the requirement that the Lagrangian be unchanged under such a rotation (i.e. we require space to be isotropic).

\delta L = \sum_i\left( \frac{\partial L}{\partial q_i}\cdot\delta r_i + \frac{\partial L}{\partial\dot{q}_i}\cdot\delta v_i \right) = 0

We can replace

\frac{\partial L}{\partial v_i}

by the vector canonical momentum

p_i

, and

\frac{\partial L}{\partial q_i}

\dot{p}_i

, leaving:

\begin{multiline*}
\left(\dot{\mathbf{p}}\cdot\delta\mathbf{r} + \mathbf{p}\cdot\delta\mathbf{v}\right) = 0 \\
\Rightarrow \dot{\mathbf{p}}\cdot\left(\delta\theta\times \mathbf{r}\right) + \mathbf{p}\cdot\left(\delta\theta\times\mathbf{v}\right) = 0
\end{multiline*}

Since

\mathbf{a}\cdot\left(\mathbf{b}\times\mathbf{c}\right) = \mathbf{b}\cdot\left(\mathbf{c}\times\mathbf{a}\right)


\begin{multiline*}
\delta\theta\cdot\left(\left[\mathbf{r}\times\dot{\mathbf{p}}\right]+\left[\mathbf{v}\times\mathbf{p}\right]\right)=0\\
\Rightarrow\delta\theta\cdot\frac{d}{dt}\left(\mathbf{r}\times\mathbf{p}\right)=0\end{multiline*}

Since

\delta\theta

is arbitrary, this requires that

\mathbf{r}\times\mathbf{p}

does not change in time, hence angular momentum is a conserved quantity.

Hamilton's Equations
Using the ideas presented above, I'm going to take a moment to derive Hamilton's Equations, which will prove useful later on.

Consider changes in the Lagrangian L, according to

dL = \sum_i\frac{\partial L}{\partial\dot{q}_i}\,d\dot{q}_i + \sum_i\frac{\partial L}{\partial q_i}\,dq_i

This can be written:

dL=\sum_ip_i\,d\dot{q}_i+\sum_i\dot{p}_i\,dq_i

since

\frac{\partial L}{\partial q_i}=\dot{p}_i

and

\frac{\partial L}{\partial\dot{q}_i}=p_i

.
Using,

\sum_i p_i\,d\dot{q}_i = d\left(\sum_i p_i q_i\right) - \sum_i\dot{q}_i\,dp_i


d\left(\sum_i p_i\dot{q}_i - L\right) = -\sum_i\dot{p}_i\,dq_i+\sum_i\dot{q}_i\,dp_i

The argument of the differential on the left is the Hamiltonian, H,

H(q,p,t)=\sum_i p_i\dot{q}_i - L

, therefore:

dH = -\sum_i\dot{p}_i\,dq_i + \sum_i\dot{q}_i\,dp_i

From here, we can obtain Hamilton's Equations:

\begin{align*}
\dot{q}_i &= \frac{\partial H}{\partial p_i}\\
\dot{p_i} &= \frac{\partial H}{\partial q_i}
\end{align*}

For m coordinates (and m momenta), Hamilton's Equations form a system of 2m first-order differential equations, compared to the m second-order equations in the Lagrangian treatment.

The total time derivative,


\frac{dH}{dt}=\frac{\partial H}{\partial t}+\sum_i\frac{\partial H}{\partial q_i}\dot{q}_i+\sum_i\frac{\partial H}{\partial p_i}\dot{p}_i

Substituting Hamilton's equations for

\dot{q}_i, \dot{p}_i

, the last two terms cancel, so

\frac{dH}{dt} = \frac{\partial H}{\partial t}

and if H does not depend explicitly on time,

\frac{dH}{dt}=0

and energy is conserved!

Noether's Theorem
The three conserved quantities above were shown to be related to the invariance of the Lagrangian under some symmetry transformation:

Translational invariance (homogeneity of space) ==> Conservation of momentum
Rotational invariance (isotropy of space) ==> Conservation of angular momentum
Time invariance (homogeneity of time) ==> Conservation of energy

Noether's Theorem states that any differentiable symmetry of the Action (integral of the Lagrangian) of a physical system has a corresponding conservation law.

To every differentiable symmetry generated by local actions, there corresponds a conserved current.

`Symmetry' here, refers to the covariance of the form of a physical law with respect to a Lie group of transformations; the conserved quantity is known as a charge and the flow carrying it as a current (c.f. electrodynamics).

Noether's Theorem, which I will discuss in more detail at a later date, is another critical component used to build gauge theories. The key thing to remember right now is that a symmetry (invariance) of the Lagrangian corresponds to a conserved quantity; we can use this result to look for the underlying symmetry behind quantities we know to be conserved (for example, electric charge).

Feedback on `Writing Scientific Documents Using LaTeX'

2008-09-22T14:38:00.002+01:00

Since I submitted my article Writing Scientific Documents Using LaTeX to CTAN, I've had several emails offering comments, suggestions and improvements. I plan on taking all of these into account for the next edition, which should be available sometime soon. Meanwhile, I thought I'd just post a little about some of the suggestions here.

The single most common suggestion was that I replace the use of $$ ... $$ to delimit displayed mathematics with the \[ ... \] form. The dollar-variant is apparently deprecated, and the newer square bracket style preferred. The reason for this, as given in http://www.ctan.org/tex-archive/info/l2tabu/, is that the $$ form is a plain TeX command, and should be avoided in LaTeX due to inconsistencies. Check the l2tabu document for more details.

Another frequent comment was that the eqnarray environment is bad, and should always be replaced by the align or similar, from the amsmath package. I originally thought I'd leave a discussion of eqnarray in there for reference, but now it seems better to remove it entirely.

Other comments related to minor corrections and requests for additional explanation in some sections, expanded coverage of BibTeX and more details on table design.

As I said above, I plan to incorporate all of these suggestions into the 5th Edition of the article.

From Lagrangian to Hamiltonian Mechanics

2008-09-22T13:19:00.011+01:00

If the Lagrangian L does not depend explicitly on time, and varies with time only through the time-dependence of the coordinate q (and its time derivative), then we can define a constant of motion known as the Hamiltonian, H.

That is, if


\begin{align*}
L = L(q(t),\dot{q}(t)) \\
\mathrm{not}~~L = L(q(t), \dot{q}(t), t)
\end{align*}

then the Hamiltonian H can be defined as

H = \dot{q}\frac{\partial L}{\partial \dot{q}} - L

For the simplest example of a particle of mass m, moving through a potential V(x), with velocity v = dx/dt, we have

L = \frac{m\dot{x}^2}{2} - V(x) = T - V

. Using the above definition of the Hamiltonian,

H = \frac{m\dot{x}^2}{2} + V(x)

which can be identified as the total energy. Since H is a constant of motion, this result corresponds to conservation of energy!

More generally,

\frac{dH}{dt}=\frac{d}{dt}\left(\dot{x}\frac{\partial L}{\partial\dot{x}}\right)-\frac{\partial L}{\partial\dot{x}}\frac{d\dot{x}}{dt}-\frac{\partial L}{\partial x}\frac{dx}{dt}

= v\left( \frac{d}{dt}\frac{\partial L}{\partial \dot{x}} - \frac{\partial L}{\partial x}\right) = 0

since the last term inside the brackets vanishes (it's the Euler-Lagrange Equation!)

Note that if L depends explicitly on time, the above does not hold, and we lose conservation of energy; here, energy can be `used' in thermodynamically irreversible processes. Of course, if we expanded the Lagrangian to take into account a large enough system, we regain conservation of energy for the universe as a whole!

The Hamiltonian H is an integral of the motion, since

\frac{dH}{dt} = 0 \Rightarrow \int_{t_0}^{t_1}\frac{dH}{dt}\,dt = \left[H\right]_{t_0}^{t_1} = 0

It contains only first-order time derivatives of the coordinate q, whereas the Euler-Lagrange equation contains second-order time derivatives,

\frac{d}{dt}\frac{\partial L}{\partial \dot{q}}

.

As an example of determining the Hamiltonian from a Lagrangian, I'd like to look at the case of Special Relativity (SR). There are two reasons for this... first, the SR Lagrangian looks a bit different from the usual T - V form of classical mechanics. Secondly, Special Relativity will feature prominently in a number of future posts as I start to direct the methods of Lagrangian and Hamiltonian mechanics towards describing gauge theories.

The SR Lagrangian in a potential V(x) is given by

L = \frac{-mc^2}{\gamma} - V(x) = -mc^2\sqrt{1-\frac{v^2}{c^2}} - V(x)

where I have used

v = \dot{x} = \frac{dx}{dt}

and have referred to the rest mass as simply m, while many books use the notation m_0.

The Hamiltonian H can be found as

\begin{multiline*}
H = v\frac{\partial L}{\partial v} - L\\
= \frac{mc^2}{\sqrt{1-\frac{v^2}{c^2}}} + V(x)\\
= \gamma mc^2 + V(x)
\end{multiline*}

The last line is readily identified as the total relativistic energy of a particle of rest mass m in a potential V(x), so our interpretation of the Hamiltonian holds!

Canonical Momenta
I'd like to take a moment now to introduce some new nomenclature and notation, and to explain why it is useful here.

First, a reminder that the Euler-Lagrange equation for a coordinate q can be written as

\frac{d}{dt}\frac{\partial L}{\partial \dot{q}} = \frac{\partial L}{\partial q}

We can define a quantity p, known as the canonical momentum conjugate to the coordinate q, as follows

p=\frac{\partial L}{\partial \dot{q}}

and a quantity F, known as the canonical force conjugate to the momentum p, as

F = \frac{\partial L}{\partial q}

These follow the usual definitions such that

T = \frac{p^2}{2m}

and

F = -\nabla V

, since the Lagrangian L = T - V.

Written in this form, the Euler-Lagrange Equation directly represents Newton's 2nd Law:

F = \dot{\mathbf{p}}

For multiple coordinates

q_i

, the Hamiltonian H is given by


H = \sum_i p_i \dot{q}_i - L

where

p_i

is the canonical momentum conjugate to the coordinate

q_i

.

Introducing the concept of canonical momentum conjugate to a coordinate is essential for the future posts I'd like to make on this blog. The reason for this is that when I discuss quantum field theory, the concepts introduced above will play an essential role in the procedure. More on this later!

Meanwhile, although I provided a couple of classical mechanics examples for the use of the Lagrangian, in order to introduce the concept at a level that most people can understand from the world around them, I will refrain from posting similar examples making use of the Hamiltonian, since I wish to progress to more complicated topics.

In order to illuminate the path from here to a discussion of gauge theories, I'd like now to simply list some of the topics I must first cover. As such, the interested reader will notice any such posts! I anticipate that this list will be incomplete, and that I will need to cover certain topics in a lot more detail than others.

Symmetries and Noether's Theorem (some of this)
Special Relativity in four-vector formulation (lots of this)
Relativistic electrodynamics (not much of this)
Relativistic Quantum Mechanics (I've already discussed the Klein-Gordon equation, but I'll go over this and other aspects, again)
Spin & Relativity (Pauli matrices and the gamma matrices)
The Dirac Equation (in some detail!)
U(1) symmetries and the Gauge Principle
Quantum Electrodynamics
SU(2) and Electroweak Unification
Aside on Superconductivity (which may be omitted)
Spontaneous Symmetry Breaking and the Higgs Mechanism

I expect those posts to take some considerable time. I've spent at least two weeks discussing Lagrangian mechanics, and we have a long way to go before I can introduce the Higgs mechanism, which is the ultimate goal of this entire series of posts. The idea is that anyone with a basic grounding in physics and mathematics should be able to learn enough about Lagrangian & Hamiltonian Mechanics, Quantum Mechanics and Special Relativity to be able to appreciate the formulation of the gauge theories of the Standard Model of Particle Physics, and to go a little beyond that and get a look at the Higgs mechanism. Let's hope I make it to the end!

Lagrangian Mechanics: From the Euler-Lagrange Equation to Simple Harmonic Motion

2008-09-20T12:06:00.008+01:00

I already wrote about obtaining Newton's Laws from the Principle of Least Action. Now I'm going to analyse a simple mass-spring system; effectively just a case of substituting in a suitable potential for the spring.

Let us go through all of the steps, once more. We'll choose our coordinate x to be the displacement from the equilibrium length of the spring, as shown on the diagram. The kinetic energy of the moving mass is then just

T=\frac{m\dot{x}^2}{2}

. The potential of a spring stretched (or compressed) x metres from its equilibrium length is given by

V=\frac{kx^2}{2}

, where k is the spring constant (N m^-1).

So, for our mass-spring system, the Lagrangian is


L = \frac{1}{2}\left(m\dot{x}^2 - kx^2\right)}

Applying the Euler-Lagrange equation, we obtain


\begin{multiline*}
-kx - \frac{d}{dt}\left(m\dot{x}\right) = 0 \\
\Rightarrow \ddot{x} - \frac{k}{m}x = 0
\end{multiline*}

Comparing the last line above with the general form of Simple Harmonic Motion (SHM),

\ddot{q} - \omega^2 q = 0

we can see that the equation of motion rendered by applying the Euler-Lagrange equation to the Lagrangian of a mass-spring system provides Simple Harmonic Motion with a frequency


\omega = \sqrt{\frac{k}{m}}

This is as expected for the case of a mass-spring system!

Lagrangian Mechanics: From the Euler-Lagrange Equation to Newton's Laws

2008-09-19T14:08:00.006+01:00

Last week I wrote about the Euler-Lagrange Equation, and how it can be obtained from the Principle of Least Action. Today, I'd like to show that this equation is consistent with Newtonian Mechanics.

For reference, the Euler-Lagrange equation for some arbitrary coordinate q is:


\frac{\partial L}{\partial q} - \frac{d}{dt}\frac{\partial L}{\partial \dot{q}} = 0

L, the Lagrangian, is the kinetic energy of a system minus the potential.

Let us now consider the case of a free particle of mass m, moving in a potential V(x). The particle's instantaneous speed is given by

v = \frac{dx}{dt}=\dot{x}

. The kinetic energy is the familiar

T=\frac{1}{2}m\dot{x}^2

.

Applying the Euler-Lagrange Equation, we have:


\begin{align*}
\frac{\partial L}{\partial x} = -\frac{\partial V(x)}{\partial x} \\
\mathrm{and}~~\frac{\partial L}{\partial \dot{x}} = m\dot{x} \\
\Rightarrow \frac{d}{dt}\frac{\partial L}{\partial \dot{x}} = m\ddot{x}
\end{align*}

Putting these together, we have:


\begin{align*}
-\frac{\partial V(x)}{\partial x} - m\ddot{x} = 0 \\
\Rightarrow m\ddot{x} = -\frac{\partial V}{\partial x}
\end{align*}

The right-hand side is the gradient of a potential energy (in 1D). Force can be defined in terms of the gradient of a potential V:


F = -\nabla V

And since

\ddot{x} = \frac{d^2 x}{dt^2}

is acceleration, a, the result of applying the Euler-Lagrange equation to a classical-mechanical Lagrangian is the familiar form of Newton's Second Law:


F = ma

In other words, applying the Euler-Lagrange equation to a suitable Lagrangian provides an equation of motion!

On Science in Education and in the Media

2008-09-18T00:04:00.001+01:00

It is a rare, but not entirely unprecedented, occasion that I find myself writing not about scientific fact, programming or mathematics, but expressing an opinion. Indeed, this time it may prove even less favourable than my post calling for the elimination of PowerPoint presentations in physics lectures. The topic for today is a mixed bag of sensationalism, media scaremongering, poor reporting and poorer education.

My thoughts turned to these topics last week as the media fully embraced the coming apocalypse. I refer, of course, to the coverage of the first beam of protons to traverse the Large Hadron Collider (LHC) at CERN. Reports of the end of the world at the hands of a man-made black hole were clearly exaggerated (after all, we're still here). That the claims of black-hole catastrophe had little scientific merit, and no expert had voiced such concerns, did not appear to bother the mass media, as a number of newspapers and even the BBC led with reports of the possibility of a planet-swallowing black hole being produced at the LHC on September 10, 2008.

What did not occur to any members of the media to point out was that no collisions were planned for that day, so there was absolutely no chance of producing a black hole on September 10, whatsoever. If you're reading this after about mid-October, you can join me in laughing at the media for completely failing to warn the public bout the collisions of protons and anti-protons now taking place routinely, at centre-of-mass energies of up to 10 TeV (assuming everything went according to plan).

What annoyed me further still were the constant references to the LHC as a `big bang' machine. ALICE, the project that will use the LHC to study quark-gluon plasmas created by the collision of protons with Lead nuclei, recreating the conditions of the early universe, is only one of four current LHC experiments; proton-Lead collisions will account for at most a third of the LHC's operational time. Meanwhile, LHCb will look at the physics of B-hadrons, particles containing one or more `bottom' quarks, at the high energies available (up to 14 TeV at full power), while CMS and ATLAS will look for whatever new physics they can find at these unprecedented centre-of-mass energies. It is ATLAS and CMS which are geared up to detect the Higgs boson, a particle which would confirm the theoretical Higgs Mechanism of spontaneous symmetry breaking to provide massive particles to the otherwise massless gauge theories of the Standard Model of Particle Physics. I shall avoid getting into the topic of the separation of science and religion, broken by the continual reference to the Higgs boson as some kind of `God particle'.

I realise that, had the headlines been along the lines of "Scientists successfully guide protons around 27 km loop", there would have been much less public interest, but that is certainly no reason to sensationalise, to play on a fear of the unknown, or to exploit the natural curiosity many people express about the origins of our universe. As a physicist, I've spent much of the past week convincing people that the world is not going to end (at least, not as a result of the LHC producing black holes) and that the research being carried out is much broader than trying to make a Big Bang of our own!

Onward, then, to the issue of science education. I can't help thinking that the public understanding of the physics going on at the LHC should be much deeper than it currently is. I know that the vast majority of people on the planet will never again hear the terms `gauge theory' or `spontaneous symmetry breaking', but surely there are other issues which a better science education would open the public's mind to.

Two examples spring to mind, and it is a coincidence, precipitated in part by media coverage, that they are related: Energy, and Global Warming (or should I say `Climate Change'?) I'll write first about energy; specifically electrical power generation.

Power generation is covered in secondary science courses in a broadly Geographical context. The student learns about the renewables vs. non-renewables, Carbon emission, inefficiency of, say, wind power, and even the relative cost of different sources of energy. Nuclear power is given brief coverage too. The problem is, none of this is science! A lot of it is politics, or economics at best.

How would I change it? I'd use the issue of power generation to launch into the science behind it; thermodynamics, turbines, efficiency and electromagnetic induction. Solar power to introduce the photoelectric effect and basic quantum physics, nuclear power to introduce radioactivity, fission and nuclear binding energy. Even fusion, though there is little new science here, and magnetic confinement is beyond the scope of secondary education! Instead, how about inertial confinement (laser fusion) to illustrate the concepts of energy, intensity, power transmission and constructive interference in waves?

Given such an education, emphasising the actual science, the public should be better able to make informed decisions, carry out accurate and productive debate, and think for themselves, rather than relying on the media to misinform.

Science is taking a central role in modern society, and will only become more central in the future. It is negligent to soften the focus of secondary science education, even for those who have no intention of pursuing a career in a science-related discipline. The public has a right to understand where science is heading and how it affects them.

I shall avoid saying anything about global warming, since it is a topic which has been demonstrated to cause `heated' debate, and I'd rather avoid arguments in which neither side has sufficient facts to back up their claims, here. Is it a good thing to try to reduce our CO2 emissions, to become more energy efficient, and less reliant on limited reserves of fossil fuels? Yes, of course it is. But when a large number of scientists cannot even agree on the cause of climate change, it is a poorly informed media pressing the issue home to the general public.

Steering now away from global warming, I will close by reiterating my central points. First, I am decidedly unimpressed by the media coverage of the LHC last week. Not only was it mostly misinformed, it was an obvious case of sensational reporting, serving only to inflate readership (or viewers, for the televised media). In particular, I found the BBC coverage sadly disappointing.

Secondly, I fear that compulsory science education does not do enough to provide the public with the knowledge they are increasingly going to need in order to understand the issues of the world around them. Finally, I make the suggestion that topical issues should be used to lead into hard science, making it relevant and capturing the interest of those not of a scientific disposition.

Writing Scientific Documents Using LaTeX: Permanent home on CTAN

2008-09-16T23:53:00.003+01:00

Following the suggestion of someone who commented on my previous post about my article Writing Scientific Documents Using LaTeX, I submitted the PDF and TeX source to CTAN, the Comprehensive Tex Archive Network. After changing the license to the LPPL, it was accepted into the archive.

As a result, the PDF can now be found permanently at the following location:
http://www.ctan.org/get/info/intro-scientific/scidoc.pdf

The CTAN entry is at http://www.ctan.org/tex-archive/info/intro-scientific/

Lagrangian Mechanics: From the Principle of Least Action to the Euler-Lagrange Equation

2008-09-14T23:40:00.005+01:00

A quantity known as the Action, A is defined as


A = \int L \,dt

where L is known as the Lagrangian. The simplest Lagrangian is given as the kinetic energy T of a system minus the potential energy V of the system:


L(x,t) = T(x,t) - V(x,t)

By extremising the Action (minimising it, in this case) we can obtain the Euler-Lagrange Equation, a key equation for Lagrangian and Hamiltonian mechanics. The concept is also used extensively in quantum mechanics and particle physics, particularly when dealing with gauge theories.

Consider an extremised path x(t) from a point x(t0) to a point x(t1), and some excursion from this path, as shown in the figure below. The excursion is given by some small function a(t), and the velocity is changed accordingly:


\begin{align*}
x(t) \rightarrow x(t) + a(t) \\
v(t) \rightarrow v(t) + \dot{a}(t)
\end{align*}

If we take x(t) to be the extremal path from x(t0) to x(t1), with the end-points fixed, and a(t) to be some small but general excursion from that path which must pass through the end-points, we can assert that:


a(t_0) = a(t_1) = 0

The Lagrangian will be changed as a result of these excursions. To first-order in small a(t), the Lagrangian transforms as:


\begin{align*}
L(x,v) \rightarrow L(x+a, v+\dot{a}) \\
= L(x, v) + a(t)\frac{\partial L}{\partial x} + \dot{a}(t)\frac{\partial L}{\partial v}
\end{align*}

The Action therefore transforms according to

A \rightarrow A + \delta A

where:


\delta A = \int_{t_0}^{t_1}\,dt \left( a(t)\frac{\partial L}{\partial x} + \frac{da}{dt}\frac{\partial L}{\partial v} \right)

The second term in the brackets above can be integrated by parts:


\int_{t_0}^{t_1}\,dt \frac{da}{dt}\frac{\partial L}{\partial v} = \left[ a(t)\frac{\partial L}{\partial v} \right]_{t_0}^{t_1} - \int_{t_0}^{t_1}dt\,a(t)\frac{d}{dt}\frac{\partial L}{\partial v}

Since a(t0) = a(t1) = 0, the integrated part (in square brackets) vanishes, leaving the following form for delta A:


\delta A = \int dt\,a(t)\left(\frac{\partial L}{\partial x} - \frac{d}{dt}\frac{\partial L}{\partial v}\right)

For arbitrary a(t), to minimise the action (by having delta A zero), the part in brackets must be zero. This is the Euler-Lagrange Equation, rewritten for a generalised coordinate q (in place of x) and using the time-derivative of q in place of the velocity v, above:


\frac{\partial L}{\partial q} - \frac{d}{dt}\frac{\partial L}{\partial \dot{q}} = 0

As I mentioned above, the Euler-Lagrange Equation is critical for the understanding of several aspects of quantum mechanics and gauge theories. Since I plan on making posts about some of these topics in the future, I felt I should begin by explaining Lagrangian mechanics, rather than jumping in at the deep end.

Two-Column Layouts in LaTeX

2008-09-14T00:33:00.006+01:00

I was recently asked about using LaTeX to create the two-column layout often seen in scientific journals. I used this layout in some lab reports and project reports.

The first thing to do is to include the multicol package in your preamble:

\usepackage{multicol}

After that, to place text within a multi-column layout, start the multicols environment (note the s at the end!) and specify the number of columns you'd like:

\begin{multicols}{2}

Finally, end the multi-column layout with:

\end{multicols}

It's as simple as that! Any sections, subsections, equations etc. between the two will now be automatically formatted into two columns.

There are, however, a few caveats and points to note for formatting considerations.

The first thing to note is that the abstract is often formatted as a single column. To do this, simply place the \begin{multicols}{2} command below your abstract, but above the first \section command.

Secondly, I usually typeset references within the two-column layout. To do this, end the multicols environment after your references section.

Finally, and most importantly, graphics and tables do not always work well with multi-column environments. In order to achieve reasonable results, I used a few workarounds.

If the graphic is large, or has small text, set it normally and it will float outside of the column layout, spanning both columns. For smaller figures, you can create an adjusted figure or table environment which forces LaTeX to set the figure in the column, at the point you typed it in the source text. This is not always what you want, but a reasonable layout can usually be obtained by playing with the exact positioning of the table or graphic commands in the source.

The following, inserted into the preamble, defines two new environments, tablehere and figurehere, which insert tables and figures inline with the column text.

% Figures within a column...
\makeatletter
\newenvironment{tablehere}
{\def\@captype{table}}
{}
\newenvironment{figurehere}
{\def\@captype{figure}}
{}
\makeatother

In addition to this, you'll probably need to scale graphics to the column width. This can be achieved with the \resizebox command inside a figure or our new figurehere environment. For example,

\begin{figurehere}
\centering
\resizebox{\columnwidth}{!}{\includegraphics{gf-graphs.eps}}
\caption{\label{gf-graphs}Graph showing applied RF frequency against magnetic field from the sweep coils, the gradient of the lines reveals information about $g_f$ for each isotope. The left gradient, corresponding to \chem{^{85}{Rb}} is $7.73 \times 10^9 (\pm 1.5 \times 10^8)$.  The right gradient, corresponding to \chem{^{87}{Rb}}, is $5.24 \times 10^9 (\pm 6.00 \times 10^7)$ }
\end{figurehere}

\columnwidth is defined by the multicol package to be the width of a text-column!

Blog Comments

2008-09-14T00:09:00.003+01:00

In the last couple of months I have been receiving a lot more comments on posts to this blog. This is, I think, a good thing. Thanks to everyone who took the time to leave a comment.

Some of the posts on this blog are meant to inform others, to precipitate discussion. In such cases, I welcome all comments, but could I ask that you please take the time to read the post before commenting, and that you understand that a number of posts are informative, rather than attempting to be original work. Those which are original work are usually (but not always) prefixed by a tale of how I got interested in that line of thought. On those other posts, the informative ones, comments telling me that this has been done before are not constructive.

Here's how you can post a constructive comment on those informative posts:

Add something I forgot, missed out or didn't cover. Provide more details, discuss caveats and cases where the general principle doesn't work.
Provoke discussion; ask questions, pose new problems or extend existing ones.
Provide your own take on the issue, even if the Boost people, TR1, physicists on Mars, or whoever, have been there and done it before. Rehashing stuff never hurts, and it might improve your own understanding as well as that of any other readers!
Offer links to further details, other blog posts about the same topic, etc. (Note: Don't take this as an excuse for spamming, I read every comment posted to this blog and I will remove those containing links to irrelevant material)

Once again, I'd like to thank everyone who has posted a comment so far. It's nice to know that people read this stuff, and even nicer to know that they take the time to comment. One of my goals here is to communicate science and technology to others, and a part of that communication is high-quality discussion.

Writing Scientific Documents Using LaTeX

2008-09-13T23:33:00.002+01:00

Back in 2007, I wrote a short document describing some of the basic LaTeX commands and explaining how to go about creating a document using LaTeX. I did this because I was writing third year lab reports collaboratively with my lab partners (as per the regulations for the third year lab at Warwick) and I wanted us all to be using LaTeX, so the reports would look good and so that editing the sections together would be easier. In order to do this, though, I had to educate my lab partners in the use of LaTeX.

Since then, the document has been read (and requested!) by other physicists, scientists and mathematicians at Warwick and in the wider world. I've updated it several times to add content and to correct minor mistakes. It was not until this week, however, that someone suggested I put it online and make it freely available to anyone who wants it (and can find it).

So, I spent a while tonight polishing it up, adding a few more subsections, and creating what is now the 13-page Fourth Edition. It is available online at http://www.physical-thought.com/scidoc.pdf

LHC First Beam: Anti-clockwise Beam

2008-09-10T15:21:00.004+01:00

This morning, a clockwise beam of protons was passed all the way around the LHC's 27km beamline, making a total of 3 complete revolutions. This afternoon, an anti-clockwise beam completed 2 revolutions. (I believe those numbers are correct; I obtained them from the CERN live webcast).

This is great news for the LHC; everything seems to have worked out incredibly well, and to have made full traversals, both within an hour of each first attempt, is a huge success!

LHC First Beam: Full Circle

2008-09-10T09:26:00.005+01:00

A proton beam has completed a circuit of the 27km of the Large Hadron Collider (LHC) beamline!

Perhaps I'll make a post later about the physics going on at the LHC. Meanwhile, I'm getting some coffee.

LHC First Beam: 7km Left

2008-09-10T09:16:00.002+01:00

The proton beam has made it around 20km of the LHC, according to the latest BBC News coverage. They keep showing a CERN schematic showing the LHC beamline layout and the beam progress. Does anyone know if this is online anywhere?

LHC First Beam

2008-09-10T08:57:00.003+01:00

I'm sitting here watching the BBC News coverage of the LHC first beam. Currently, it seems that protons have made it half way round the 27km beamline.

Over the course of the day, the guys at CERN plan to get a proton beam around the entire 27km!

Groups and Matrix Fun

2008-08-22T16:29:00.012+01:00

In the process of working through Introduction to Elementary Particles (Griffiths), this summer, I've been trying most of the end-of-chapter exercises. Some are reminders of things I did a while ago, others are a different take on something I did fairly recently, and have proved interesting. Still other questions are topics I had previously taken for granted due to lack of time to work out for myself. This post is about one such question.

The question involved showing that sets of n*n matrices with various properties form groups. The conditions for a set to form a group under some operation (multiplication, here) are:

Closure: If A and B are in the set, then the product AB is also in the set, i.e. there exists some M such that AB = M
Identity: There is an element I such that IA = AI = A for all elements A in the set
Inverse: For every element A, there is an inverse, A^-1, such that AA^-1 = A^-1 A = I
Associativity: A(BC) = (AB)C

The fourth is pretty obvious for matrices anyway, so I'll ignore that. Equally, the identity matrix,

I = \left(\begin{array}{cccc}
1&0&\ldots & 0\\
0 & 1 & \ldots & 0\\
\vdots & \vdots & \ddots &\vdots \\
0 & 0  & \ldots & 1\end{array}\right)

provides the identity element in all of the cases considered here.

All we need to look at, then, are closure and inverse.

Set of all n*n unitary matrices
A unitary matrix is one where

U^{-1} = U^\dagger = (U^\star)^T

, i.e. the inverse of a unitary matrix is its transpose complex-conjugate.

To check closure, we require

AB = M ~,~ M^{-1}=M^\dagger

, where A and B are two unitary matrices (i.e. if A and B are members of the set, we need to check that their product M is also):


\begin{align*}
AB = M\\
M^\dagger = (AB)^\dagger = B^\dagger A^\dagger\\
\mathrm{If~M~is~unitary,~~~~}\\
MM^\dagger = ABB^\dagger A^\dagger\\
= AI_nA^\dagger = I_n \\
\Rightarrow MM^\dagger = I_n = MM^{-1}
\end{align*}

which provides closure of the set.

The other thing to check is the inverse. Since A is unitary,

AA^\dagger = I = A^\dagger A

implies that

A^\dagger

is unitary, and therefore

A^{-1} = A^\dagger

is unitary!

Thus, the set of all n*n unitary matrices forms a group, known as U(n).

Set of all n*n unitary matrices with determinant +1
To determine closure, we already know that the product of two unitary matrices is unitary, from above. All we need to check is that the determinant remains +1:


\begin{align*}
\mathrm{det}(AB) = \mathrm{det}(A)~\mathrm{det}(B) \\
\Rightarrow \mathrm{det}(AB) = 1\times1 = 1
\end{align*}

And for the inverse, since

A^{-1} = A^\dagger

, and

AA^\dagger = I

, and

\mathrm{det}(I) = +1

, then

\mathrm{det}(AA^\dagger) = \mathrm{det}(A)~\mathrm{det}(A^\dagger) = +1

. Since we know that

\mathrm{det}(A) = +1

as it is a member of the set of unitary matrices with determinant +1, then

\mathrm{det}(A^\dagger}) = +1 \Rightarrow \mathrm{det}(A^{-1}) = +1

as required.

This forms the group SU(n).

Set of all n*n orthogonal matrices
Orthogonal matrices are just real unitary matrices. Since I showed above that the set of all unitary matrices forms a group, the set of orthogonal matrices clearly forms a (sub) group too! If you want to be convinced, follow exactly the steps from above, replacing

A^\dagger

with

A^T

, the matrix transpose.

This forms the group O(n).

Set of all n*n orthogonal matrices with determinant +1
Exactly as above for the unitary (or orthogonal) matrices, and exactly as above for the matrices of determinant +1. This forms the group SO(n).

Noncopyable: A C++ class which cannot be copied

2008-08-21T14:58:00.007+01:00

Yes, I know this is reinventing the wheel [boost.org], but sometimes reinventing the wheel is a good way of explaining to people why it's round!

Consider the following problem: you have a handle on some resource (file, memory, socket, database connection, whatever...) and you want to have a single instance of a class responsible for managing that handle, freeing the resource when you're done with it. You want to be able to pass references to this one object around, but you don't want to encounter a situation where something made a copy of the object, so that the handle is still hanging around after the resource itself has been freed. This is, for example, useful in implementing smart pointers, socket libraries, etc. where you do not want to share ownership of the resource, but you want to share access to it, through the wrapper class which owns it.

In C++, a copy of an object can be made in a number of situations. The most obvious is direct assignment:

struct X {
  int foo;
};

X a, b;
a.foo = 2;
b = a;
// We now expect b.foo == 2

Less obvious, but equally valid, is a function call using pass-by-value:

int Square(X object) { return object.foo*object.foo; }

X a;
std::cout << Square(a) << "\n";

In the examples above, the copy is harmless (though passing by reference could be more efficient!) If, however, the object contained a pointer you don't wish to relinquish ownership of, or if copying doesn't make sense for a certain class, you can prevent it altogether!

Copying in C++ is handled in two ways. The first example above uses the copy assignment operator, and the second example uses a copy constructor. Both of these are provided implicitly by the compiler if a class does not declare its own. Here are some relevant quotes from the C++ standard:

12.1.10 A copy constructor for a class X is a constructor with a first parameter of type X& or of type const X&.
12.8.2 A non-template constructor for class X is a copy constructor if its first parameter is of type X&, const X&, volatile X& or const volatile X&, and either there are no other parameters or else all other parameters have default arguments.
12.8.4 If the class definition does not explicitly declare a copy constructor, one is declared implicitly...
12.8.9 A user-declared copy assignment operator X::operator= is a non-static non-template member function of class X with exactly one parameter of type X, X&, const X&, volatile X& or const volatile X&...
12.8.10 If the class definition does not explicitly declare a copy assignment operator, one is declared implicitly...

The standard tells us exactly what we need to know here: the copy constructor should take a form similar to

X::X(const X&);

and the copy assignment operator should take a form similar to

X& X::operator=(const X&);

If either of the two are not defined, the compiler provides an implicit version which is a public member function.

How does this help us to make an object which is noncopyable? Well, if we declare our own copy constructor and our own copy assignment operator, we prevent the compiler from generating its own. What is more, we can choose the visibility of our manually-declared versions. In particular, we can make them private (or protected, if it makes sense for derived classes to be able to implement copy semantics somehow...)

By making the copy constructor and the copy assignment operator private, we prevent their use outside of the class itself, and therefore prevent the rest of the program from being able to make copies of the class!

Now, if we want to make a class noncopyable, we can simply define its copy constructor and copy assignment operator as private members, documenting why we did that so that future maintainers don't make them public and screw up our best efforts to prevent copying, or we can do as they do in boost, and make a noncopyable base class!

If a class derives from the noncopyable base class, it too will be noncopyable, since the copy of an object of a derived class must necessarily invoke the copy constructor or copy assignment operator of the corresponding base class, and in this case they are private, so that can't be done!

Here's our noncopyable class:


class noncopyable {
 private:
  noncopyable(const noncopyable& x);
  noncopyable& operator=(const noncopyable& x);
 public:
  noncopyable() {}; // Default constructor
};

And here's a class which is derived from it, and a program which tries to make copies, first using the copy constructor as part of a function call, then using the copy assignment operator.


#include <iostream>

class noncopyable {
 private:
  noncopyable(const noncopyable& x);
  noncopyable& operator=(const noncopyable& x)};
 public:
  noncopyable() {}; // Default constructor
};

class X : noncopyable {
 public:
  int c;
};

void print(const X x)
{
 std::cout << x.c << "\n";
}

int main()
{
 X x1;
 x1.c = 3;
 print(x1);
 X x2;
 x2 = x1;
 return 0;
}

And finally, here's what GCC spits out when we try to compile this, not only informing us that what we're trying to do won't work, but also clearly stating why in the error message -- the class name noncopyable is clearly visible!


loki:~$ g++ -o nocopy nocopy.cpp 
nocopy.cpp: In copy constructor 'X::X(const X&)':
nocopy.cpp:5: error: 'noncopyable::noncopyable(const noncopyable&)' is private
nocopy.cpp:11: error: within this context
nocopy.cpp: In function 'int main()':
nocopy.cpp:24: note: synthesized method 'X::X(const X&)' first required here 
nocopy.cpp:24: error:   initializing argument 1 of 'void print(X)'
nocopy.cpp: In member function 'X& X::operator=(const X&)':
nocopy.cpp:6: error: 'noncopyable& noncopyable::operator=(const noncopyable&)' is private
nocopy.cpp:11: error: within this context
nocopy.cpp: In function 'int main()':
nocopy.cpp:26: note: synthesized method 'X& X::operator=(const X&)' first required here

By contrast, this version passes the object by reference, and compiles cleanly:


#include <iostream>

class noncopyable {
 private:
  noncopyable(const noncopyable& x);
  noncopyable& operator=(const noncopyable& x);
 public:
  noncopyable() {}; // Default constructor
};

class X : noncopyable {
 public:
  int c;
};

void print(const X& x)
{
 std::cout << x.c << "\n";
}

int main()
{
 X x1;
 x1.c = 3;
 print(x1);
 return 0;
}

Output:

loki:~$ ./nocopy 
3

The Plan: The Update

2008-08-21T14:45:00.003+01:00

In my second ever blog post, I outlined the reasons I started this blog. To recap, I planned to write about my ideas and projects in the hope that it would in some way document them, and that I might keep them alive long enough to make some actual progress.

Over the course of the last eight months, the blog has evolved in a very different direction. Some of it is about posting my ideas and projects. In fact, that constitutes a large proportion of the posts. What I didn't foresee, however, was that my interest in writing and explaining things would take over somewhat.

In the time this blog has been running, I've posted about electrodynamics, quantum mechanics, relativity and several other physics topics, as well as discussing fairly well-known but interesting aspects of a number of programming languages, often including small examples to demonstrate. In some cases, I've been playing with this stuff for the first time. In others, I've known about it for a while, but someone has asked me to explain it, or refreshed my interest, and I've felt like posting about it here.

In short, this blog has expanded beyond my original reason for creating it, and has become a vast collection of whatever I happened to be thinking about on any given day. I've also started to get comments on some of the posts, which can only imply that people actually bother to read this stuff. I'm amazed, and I'd like to thank everyone who has made a comment, especially if it has been to point out something I forgot!

As for the success of the original goal, reminding me about my projects and ideas, it has been of mixed worth. Some of the things I've returned to numerous times, others I simply haven't had time to look over again. At least I now have a written record of what I was thinking the first time around, now, so that should I ever decide to take up one of these projects, I can refer back to my original thoughts (and probably laugh at how wrong they were!)

Slight modifications to Python/C sizeof program

2008-08-21T14:27:00.002+01:00

I've made a few small changes to the sizeof.py program I posted here. The first is in response to a comment on the original post, noting that I should change


os.system(binfilename) # old
os.system('./' + binfilename) # new

so that the program works on systems where the current working directory, ".", is not in the $PATH environment variable.

The second change was to make the behaviour a bit better when you pass an undefined type to the sizeof program. The behaviour of the old version was to try to compile, showing the gcc errors resulting. This was not a problem, but the unlink(binfilename) line was before the unlink(cfilename) line, and if gcc failed, the binfile did not exist, so this unlink() failed, and the deletion of the source file never happened. I've now changed the order so that the .c file is deleted first, and sizeof.py should now leave no temporary files lying around, even if there's a compile error!

The third change was to merge two lines into one. The old version had


types = sys.argv[1:]
for type in types:

and the new version merges these into


for type in sys.argv[1:]:

In all, not big changes, but the second change at least is fairly important behaviour for handling cleanup of temporary files!

Here's the full source again:


#!/usr/bin/env python

import sys, time, string, os

code_pre = '#include <stdio.h>\nint main() {\n'
code_post = '\treturn 0;\n}'
code_mid = '\tprintf("TYPE: %d\\n", sizeof(TYPE));\n'

if len(sys.argv) > 1:
 binfilename = 'pyctypes_' + str(time.time())
 cfilename = binfilename + '.c'
 out = open(cfilename, 'w')
 out.write(code_pre)
 for type in sys.argv[1:]:
  out.write(code_mid.replace('TYPE', type))
 out.write(code_post)
 out.close()
 os.system('gcc -o ' + binfilename + ' ' + cfilename)
 os.system('./' + binfilename)
 os.unlink(cfilename)
 os.unlink(binfilename)

C++ Function Objects

2008-08-21T02:23:00.003+01:00

I had a conversation with a computer scientist this evening about object-oriented stuff you can do in C++. Unfortunately, it seems that most comp. sci. teaching is through Java, which I find severely lacking in many respects when it comes to good object-oriented programming (multiple-inheritance and operator overloading, to name two things C++ has which Java doesn't).

Anyway, one of the things that came up was the use of function objects, something which can't be done (at least in the same way) in Java because of the lack of operator overloading. Here is an example:

You have a vector of numbers and you want to get the total. You could do something involving iterators and a for() or while() loop, or you could create a nice generic class for future use. The basic idea is simple... the constructor initialises some private member variable to zero (or whatever makes sense for the type in question), then you provide an operator() which allows objects of that class to be called as though they were functions. Define operator() such that it adds its argument to the private member variable representing the total, and call the whole thing in a for_each().

The only point to consider is that std::for_each() takes a function-like-object as its fourth-argument, but does so by value. This means that if you create an accumulator object and pass it to for_each(), then a copy is made, for_each() acts using that copy, and the copy is destroyed when for_each() is finished. The original never gets its status updated.

There are two ways around this. The first is to pass the function object by reference, which means explicitly stating the template arguments to for_each(), something like:


std::for_each<std::vector<int>::iterator, Accumulator<int>&>(vec.begin(), vec.end(), acc);

The second method relies on the fact that for_each() returns its function object! None of the other standard algorithms do this, so you can only apply this neat trick with for_each(), but it is what I used in the code below, which represents a simple yet fairly comprehensive solution to the problem of adding each item in a container, regardless of the actual type of the item.


#include <algorithm>
#include <functional>
#include <iostream>
#include <vector>

template<typename T> class Accumulator {
 T total;
public:
 Accumulator() : total() {};
 void operator()(T& x) { total += x; };
 T& getTotal() { return total; };
};

int main()
{
std::vector<int> foo;
foo.push_back(1);
foo.push_back(2);
foo.push_back(3);

Accumulator<int> acc;
acc = std::for_each(foo.begin(), foo.end(), Accumulator<int>() );
std::cout << acc.getTotal() << std::endl;
}

The output of this program is, as expected, the number 6.

The template class Accumulator constructs its total variable using the default constructor, which should provide a sensible default "zero" value. The getTotal() and operator() use references to avoid copying potentially large objects.

The call to for_each() assigns its result to an Accumulator object and uses a temporary Accumulator as the function passed in.

This is a fairly simple example of C++ function objects, and there are much more powerful (and much more complicated) possibilities. I wanted to present this example since it should be fairly easy to follow.

Please note that at no point do I claim to have invented this technique. I'm discussing it here to inform others, and for no other reason. You see, I've always thought it was important to pass on knowledge, to help others to understand things. That is all I'm doing here. If you have a problem with that, I don't care.

FeynMF: Feynman Diagrams in LaTeX

2008-08-20T16:38:00.004+01:00

When I posted the seventeen fourth-order Compton Scattering Feynman Diagrams, I was asked in a comment to show the LaTeX code responsible for generating them. I can't find the file for those diagrams, but I've drawn a QCD diagram here; quark-antiquark annihilation and production, with a gluon propagator.

Here is the diagram:

And here is the LaTeX code:


\documentclass{article}
\usepackage{feynmf}
\begin{document}
 \setlength{\unitlength}{1mm}
 
 \begin{fmffile}{feyn}
  \begin{fmfgraph*}(40,25)
   \fmfleft{q1,q2}
   \fmfright{q3,q4}
   \fmf{fermion}{q1,v1,q2}
   \fmf{fermion}{q3,v2,q4}
   \fmf{gluon}{v1,v2}
   \fmflabel{$q$}{q1}
   \fmflabel{$\bar{q}$}{q2}
   \fmflabel{$q^\prime$}{q4}
   \fmflabel{$\bar{q}^\prime$}{q3}
   \fmfdot{v1,v2}
  \end{fmfgraph*}
 \end{fmffile}
 
\end{document}

In order for this to work, MetaFont and FeynMF need to be installed into the LaTeX distribution on your system. The commands to generate the (PDF, in this case) output are:


pdflatex quark
mf '\mode=localfont; input feyn;'
pdflatex quark
pdflatex quark

Python / C Program to display size of C types

2008-08-20T16:02:00.004+01:00

The code below is a Python program which takes a number of command-line arguments, each representing a C type-name. It constructs a C program which prints the size (in bytes) of each of the given data types, compiles (using gcc) that program, and runs it, then deletes the temporary source and binary files. It is useful if you want to know the size of a particular type on a given architecture but don't want to manually write out a C file that looks something like:


#include <stdio.h>
int main() {
    printf("int: %d\n", sizeof(int));
    return 0;
}

every time you want to know! In fact, the code above is exactly what this Python script produces and executes, but it inserts a printf line for each type you tell it you want to know the size of!

Here's the code:


#!/usr/bin/env python

import sys, time, string, os

code_pre = '#include <stdio.h>\nint main() {\n'
code_post = '\treturn 0;\n}'
code_mid = '\tprintf("TYPE: %d\\n", sizeof(TYPE));\n'


if len(sys.argv) > 1:
 binfilename = 'pyctypes_' + str(time.time())
 cfilename = binfilename + '.c'
 out = open(cfilename, 'w')
 out.write(code_pre)
 types = sys.argv[1:]
 for type in types:
  out.write(code_mid.replace('TYPE', type))
 out.write(code_post)
 out.close()
 os.system('gcc -o ' + binfilename + ' ' + cfilename)
 os.system('./' + binfilename)
 os.unlink(binfilename)
 os.unlink(cfilename)

And here's a sample run:


loki:~$ ./sizeof.py "signed char" "unsigned long long" double float short
signed char: 1
unsigned long long: 8
double: 8
float: 4
short: 2

Estimating the volume of water in a shower

2008-08-19T14:45:00.006+01:00

I noticed this morning (actually, afternoon... I was lazy today and didn't get up until almost 13:00) that in the summer I have showers with the water at about 30C, and in winter I turn it up to about 38C. Ignoring the obvious explanation that I just like to be warmer when it's colder outside, I decided to try to work out the approximate volume of water in the shower at any moment in time from those two temperatures and the air temperatures: around 25C on a warm summer day, and around 20C with the central heating on in winter.

Assuming that I like a constant heat transfer (and ignoring the actual direction of it... for the values given the heat is actually flowing out of my body, not in) and using

Q = mc\Delta T

for the energy transfer, where m is mass, c is specific heat and dT is temperature difference, we can write:


\begin{align*}
Q_\mathrm{Winter} = c_A\rho_{A20}V_A(20-37) + c_W\rho_{W38}V_W(38-37) \\
Q_\mathrm{Summer} = c_A\rho_{A25}V_A(25-37) + c_W\rho_{W30}V_W(30-37)
\end{align*}

where

\rho

is the density of air (A) and water (W) at the given temperature, V is the volume of air or water and c is the specific heat.

Assuming that Q(Winter) = Q(Summer),


\begin{align*}
-17c_A\rho_{A20}V_A + c_W\rho_{W38}V_W = -12c_A\rho_{A25}V_A - 7c_W\rho_{W30}V_W \\
\left(7\rho_{W30} + \rho_{W38} \right)c_WV_W = \left(17\rho_{A20} - 12\rho_{A25}\right)c_AV_A
\end{align*}


\frac{c_W V_W}{c_A V_A} = \frac{17\rho_{A20} - 12\rho_{A25}}{7\rho_{W30}-\rho{W_38}}

The specific heats are:

c_W = 4.18~\mathrm{kJ~kg^{-1}~K^{-1}} ~~,~~ c_A = 1.005~\mathrm{kJ~kg^{-1}~K^{-1}}

, giving:


V_W = 0.24 V_A \left( \frac{17\rho_{A20} - 12\rho_{A25}}{7\rho_{W30}-\rho{W_38}} \right)

And finally, the densities:
Air at 20C: 1.200870572 kg/m^3
Air at 25C: 1.180721737 kg/m^3

Water at 30C: 995.7 kg/m^3
Water at 35C: 992.4 kg/m^3

The total volume of the shower I stand in is about 1.27 m^3, and ignoring the volume of my body, this gives a volume of water as:


V_W = 2.39\times 10^{-4} ~\mathrm{m^3}

1 m^3 is 1000 litres, so this is a volume of 0.23 litres of water in the shower at any given time.

Considering that a shower is a fairly fine spray of continuously falling water, this is not a bad estimate. I don't know the actual water delivery rate of the average shower, but I am willing to bet it isn't far off this value.

However, there are a number of problems with this model:

The densities quoted for air are valid only for dry air. Water vapour reduces the air density, and in a shower there's going to be a fair amount of waver vapour around!

I have neglected the volume, mass and specific heat of the human body

Not all of the air in the shower is in contact with the body, similarly not all of the water is

The water tends to form a film around the body, insulating it from the air, so the effect of water temperature is likely to be much more exaggerated than that of air temperature

I'm sure there are more flaws with this model, those were just the obvious ones. In any case, it was an interesting exercise and the number I got out of it sounds vaguely sensible, if perhaps a little on the low side! I'd say it was more likely that there was perhaps getting on for a litre of water in the shower, at any given time.

Multiplying a number by -1: the NEG instruction

2008-08-19T00:35:00.002+01:00

I completely forgot about the NEG instruction available on IA32 CPUs... it, well, negates things. gdb was more than happy to remind me about it when I wrote this test code:


int f(int x) {
 return 2*x;
}

int main() {
 int x = 2;
 int y = -f(x);
 return 0;
}

I compiled it using gcc -ggdb and then fired up gdb. I instructed gdb to disassemble main, and the relevant lines are quite clear:


0x0804837c <main+30>:   call   0x8048354 
0x08048381 <main+35>:   neg    %eax

The compiler translates the negation of the result of a function call into a call followed by a negation of its return value using the NEG instruction. This is so obvious I don't know why I didn't realise it immediately, but hey, now I know!

Multiplying a number by -1: A bitwise approach

2008-08-19T00:23:00.005+01:00

In my previous post, I mentioned that when the compiler encounters a line such as


fx = -f(x);

it has to multiply the result of f(x) by -1. This can be interpreted in a number of ways, for example:


// multiply by -1
fx = -1 * f(x);
// or subtract from 0
fx = 0 - f(x);

Since compilers are tricky beasts, and processors are trickier still, it's entirely possible a bitwise approach is taken. Indeed, the processor may have a single instruction for "multiply by -1". I can't be bothered to check the Intel manual right now, but in any case, here's a program for bitwise multiplication by -1, and the output it produces:


#include <stdio.h>

int main()
{
 int x;
 for(x = -5; x < 5; ++x)
  printf("x: %d\t-x: %d\n", x, ~x + 1);
 return 0;
}


loki:~$ ./neg-trial
x: -5   -x: 5
x: -4   -x: 4
x: -3   -x: 3
x: -2   -x: 2
x: -1   -x: 1
x: 0    -x: 0
x: 1    -x: -1
x: 2    -x: -2
x: 3    -x: -3
x: 4    -x: -4

It would be interesting to find out exactly how a line like y = -f(x); is compiled... perhaps I'll dig around a little more and post again about this!

Profiling a Minimum / Maximum finding algorithm

2008-08-18T23:47:00.004+01:00

Today I extracted the Minima/Maxima finding algorithm from Pebble v1 and re-implemented the necessary parts to make it work in a more-or-less standalone manner. The min/max class itself requires only some functor class representing the function it operates on. For this I created a base class Function with a virtual double operator()(double x);

In addition, the min/max algorithm in Pebble throws a Pebble::Exception on certain error conditions, so I had to provide an Exception class as well.

The reason for doing this is twofold. First, it's the best algorithm I've found for finding the minima of a function, so it'll be useful in Pebble v2. Second, Pebble v1 uses a minimum finder to locate maxima. It does this because if a function f(x) has a maximum at some x, then the function -f(x) has a minimum at the same x. In order to use the minima-finder, then, we had to pass "negated" versions of the functions to it. Pebble v1 does this using a Negater class which derives publically from the Pebble::Math::Calculus::Integrand class, a functor representing the mathematical functions handled within Pebble.

In the code I extracted today, I replaced Pebble::Math::Calculus::Integrand by Pebble::Math::Function. As I said above, this is an abstract base class with virtual public member functions (I don't want to say methods, it sounds too Java-esque). Anyway, in the new Pebble::Math::Minimiser, I provided two equivalent functions, getMinimum() and getMaximum(). The code is exactly the same between the two, apart from the calls to the function itself. In getMinimum(), the function is called as-is, so a function wrapped in a Negater object will have its maxima found through the getMinimum() function.

However, getMaximum() changes the calls to -f(x), eliminating a whole level of indirection and virtual member function lookup. I expected this version to be faster, but to check this I profiled a test-case, finding the maximum of a probability distribution function used in the original Pebble code. The test-case performs one million constructions and maximum-finding operations on this distribution, and valgrind --tool=callgrind was used to profile the execution for a version using a Negater and for a version directly using getMaximum().

Here are the call graphs, first for the getMinimum(Negater(fn)):

and then for getMaximum(fn):

The numbers in the graphs are the absolute cost of execution, measured in terms of instruction fetch (for those who know what that means... I'm only partly sure myself!) As I expected, the version which uses -f(x) is faster than the version that uses a Negater to wrap f(x), and the execution cost shows this. The absolute cost of getMinimum() is 5,539,000,000 with the Negater, and 4,947,000,000 without. Looking at the source-code profiling information, the lines in which f(x) are evaluated actually have a higher direct cost for the -f(x) version. This is to be expected, since there's a multiplication by -1 going on (or a subtraction from zero, or some cunning bitwise trick to make a number negative). Both versions involve a virtual function call, so the direct overhead of that is the same.

The difference comes from the overhead introduced by the Negater object. The cost of a line containing the source:


fx = (*f)(x);

has a cost of 874,000,000 for 19,000,000 calls to Pebble::Math::Negater::operator().

By contrast, the equivalent line in the getMaximum() version:


fx = -(*f)(x);

only has a cost of 551,000,000 for the same 19 million calls directly to Mode1::operator().

If you didn't understand a word of the above, I'll summarise the key points here:

Obviously, eliminating a level of indirection makes things faster, with fewer overheads
Duplicating the code, changing only one or two lines, in this way, is generally considered to be a bad thing
Having the Negater class around could be useful for other situations where one cannot easily just put a negative sign into the code, especially if I decide to use dynamic libraries, where one will not be able to simply change the code at leisure
In this case, considering that in my previous profiling of Pebble v1, the Minimiser and Integrator are the biggest bottlenecks, the speed increase benefit outweighs the code-duplication drawback, so I'm going to leave it in there!

Googling MD5 Hashes

2008-08-18T12:52:00.008+01:00

It occurred to me last night that it's likely there are a lot of MD5 hashes floating around on web pages with the associated plaintext right next to them, and that Google, by indexing most of the web, would perhaps be one of the world's most comprehensive lookup tables for MD5 (or other cryptographic function) data.

This morning, I set out to generate the MD5 hashes of a number of "commonly used passwords". Personally, I can't believe anyone uses these in the real world, but there's plenty of evidence to suggest they do!

Once I'd generated the hashes, I Googled each hash and made a note of the number of search results returned. Additionally, I checked that within the first page of Foogle results, the plaintext appeared clearly in such a way that it is obviously the plaintext corresponding to that hash value. Every single hash I tried (except those for which no results were found) had its plaintext within the first page of Google results.

Here's the word list, hash for each word, and number of Google results as of about 11:30 this morning.

And here's a graph representing these results. The last 'word' is a character sequence I made up on the spot, and reassuringly it returned no results!

(Click image for larger version)

It would appear that I'm not the only person (and certainly not the first) to have had this idea. In fact, there exists a proof-of-concept web site at Google Hash ( http://www.nth-dimension.org.uk/utils/ghash.php ) where the aim is to get Google to index the hash and its associated plaintext!

For those that care about such things, the graph was produced from columns 1 and 3 (plaintext and search result count) of the data I posted above in a file with the hash column removed. The gnuplot command was:

plot 'counts.txt' using ($0-0.2):2:xticlabels(1) with boxes title '',
'counts.txt' using 0:2:2 with labels center offset 0,1 notitle

Of course, there were also titles, axis labels and PostScript output settings, but they're not important here.

Surf Photography: Rip Curl Boardmasters 2008

2008-08-17T14:43:00.003+01:00

I took this on Thursday August 7th, 2008 at Fistral Beach, Newquay Cornwall. The Rip Curl Boardmasters is a surf competition and festival which has been at Fistral for as long as I can remember, and it's a great atmosphere.

This was taken at the 300mm end of my 70-300mm lens with a shutter speed of 1/640 s, aperture f/8.0 and ISO 200. It was subsequently adjusted in Photoshop to improve the contrast a little, then cropped to remove a bit of empty sea.

Pebble v1 Call Graph

2008-08-17T13:12:00.003+01:00

Following up on my recent post about a C++ implementation of the Mersenne Twister pseudo-random number generator (PRNG), I decided to find out just how much the performance of the PRNG actually affects a Monte Carlo system such as Pebble. Clearly, the random number generation is a critical part of the procedure, but as I'm about to show, it is definitely not the slowest part of Pebble, by a long way.

The valgrind tool allows a programmer to check their programs for many things, including memory leaks (the original goal of the valgrind project). One of the neat things it does is function call profiling and graphing. The idea is simple, you run a program inside valgrind using the callgrind option, and it tracks every function call made within the program. This lets you figure out how many times a function is being called, and how much of the program execution time is spent within a given function, allowing you to easily profile the code and determine functions or classes which are good candidates for optimisation.

Here is the call graph for Pebble v1, which uses a C++ implementation of the Mersenne Twister which I wrote last summer, rather than the new one I have recently posted about.

(Click for larger version)

The random number generation, performed by Pebble::Math::Random::MersenneTwister::getDouble and getLong, account for a mere 3.33% of the time spent in the event generator function Pebble::EventGenerator::generate. In contrast, finding minima of the probability distribution functions (using Pebble::Math::Utility::Minimiser::getMinimum) accounts for 31.29%, and integrating the PDFs for normalisation purposes (using Pebble::Math::Calculus::RombergIntegrator::integrate) accounts for 50.82% of the time spent generating events.

It is clear from this that, if I want to improve on the speed of generation in "Pebble v2", I need to focus my efforts on efficient minimia-finding and numerical integration functions.

There is also 7.47% of the time taken to generate events spent adding particles to an Event object, using the Pebble::Event::addParticle function. A lot of this time is accounted for by the memory allocation, and as I mentioned in my previous post, eliminating such object-oriented extravagances, while maintaining a sensible object-oriented design for the important parts of the program, should help to boost the speed of any new version.

Mersenne Twister

2008-08-15T16:21:00.004+01:00

A key part of the Pebble program was the ability to generate lots of good-quality pseudo-random numbers for Monte Carlo (MC). The Mersenne Twister (MT) algorithm is ideally suited to this, and is indeed the algorithm we used. The original author, Makoto Matsumoto, has published C source for the algorithm, which I converted into a C++ class for the project, using the same interface as several other random number generators (RNGs) provided within Pebble.

So, following my post yesterday about considering building a "Pebble v2" or similar, I've spent the afternoon looking back into Mersenne Twister. The C code is fairly small, and reasonably easy to understand (surprisingly for something as complicated as a PRNG!) and I had no difficulty converting it to C++ the first time around. This time, I've ignored a lot of the complicated interface I built in the first time, since I no longer intend to provide other RNGs -- Mersenne Twister is more than good enough for any MC purposes.

Anyway, I modified the MT code to encapsulate it in a C++ class (code below), then ran some timing tests against the C version. In order to get useful timing data, I had to generate 1 million integers and 1 million real numbers, timing the program run for the entire thing. A real number is generated from a single integer, so this is roughly equivalent to two million integer number generations, plus a little overhead in converting half of them to a number between 0 and 1.

For the purposes of this test, I removed all the program output code from both the original C code and my new C++ version, timing the programs as they generated and threw away two million random numbers. The reason for removing output code was to prevent any skew in the results from the different output systems used by C and C++.

The test itself consisted of running a program to generate 1M integers and 1M reals, using the C or C++ versions of the MT algorithm, respectively. Each test was run ten times, using the bash time command to provide timing details. I discarded the real time taken, which refers to what is known as wall time -- the amount of actual elapsed time according to a clock on the wall -- and focused only on the user and system times, which are the number of seconds the CPU actually spent running the program in user-space and in system-calls, respectively.

The graph below shows the results. Both are incredibly fast, with the slowest taking a mere 0.156 seconds of CPU time. If I can generate two million random numbers in 0.2 seconds, this will certainly not be a bottleneck in any MC system!

(Click image for larger version)

The comparison, however, is slightly discouraging. The system times are incredibly low, 0.002 and 0.001 seconds for C and C++ respectively, averaged over ten runs. This is almost certainly because I removed all of the I/O, so the only time spent in the "system" would be the initial memory allocation and the final deallocation. In fact, I'd expect the C++ variant to take longer here, as the MersenneTwister object must be constructed. Still, these times are so close to zero it's hard to talk about them meaningfully.

The "user" times tell a different story. These are also averaged over ten runs, and show that the C++ version took almost twice as long as the C version. I expected there to be a little overhead from dynamically allocating the array used to store the internal state of the MT generator, but that should be a one-time thing, and the rest of the code is almost identical. I can only assume that the discrepancy arises from the function calls being mapped to member functions in C++ rather than global functions in C.

Many would argue that the factor of 2 slowdown is completely unacceptable, even if the C++ version does hide the details of MT (no defines, no global variables, just a few member variables). I agree that the reduction in speed is quite dramatic here, but compared to other components of a functioning Monte Carlo system, this overhead is negligible, and for any program written predominantly in C++, encapsulating the functionality and hiding the details is a good idea. This is where the hardcore C fans would argue that no system which requires speed or efficiency should be written in C++, but I disagree, backed up by plenty of evidence for C++ development in the physical sciences, and in particle physics in particular!

MersenneTwister.hpp


/*
 A C++ program for MT19937, with initialization improved 2002/1/26.
 Original C Coded by Takuji Nishimura and Makoto Matsumoto.
 Converted to a C++ class by Andrew J. Bennieston

 Before using, initialize the state by using init(seed)
 or initArray(init_key, key_length).

 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
 All rights reserved.
 C++ Class modifications copyright (C) 2008, Andrew J. Bennieston.

 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions
 are met:

   1. Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.

   2. Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.

   3. The names of its contributors may not be used to endorse or promote
      products derived from this software without specific prior written
      permission.

 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


 Any feedback is very welcome.
 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
 email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)

 Feedback on the C++ modifications to andrew @ physical - thought . com (remove spaces)
*/


class MersenneTwister {
  private:
      const int N;
      const int M;
      const unsigned long MATRIX_A;
      const unsigned long UPPER_MASK;
      const unsigned long LOWER_MASK;
      int mti;
      unsigned long* mt;
  public:
      // Constructors
      MersenneTwister()
          : N(624), M(397), MATRIX_A(0x9908b0dfUL),
          UPPER_MASK(0x80000000UL), LOWER_MASK(0x7fffffffUL),
          mti(N+1)
          { mt = new unsigned long[N]; };

      // Destructor
      ~MersenneTwister() { delete[] mt; };
   
      // Initialisation
      void init(unsigned long s);
      void initArray(unsigned long init_key[], int key_length);

      // Random number generation
      unsigned long genInt32(); // Random integer on [0,0xffffffff]
      long genInt31(); // Random integer on [0,0x7fffffff]
      double genRealClosed(); // Random real on closed range [0,1]
      double genReal(); // Random real on half-open range [0,1) (i.e. not including 1)
      double genRealOpen(); // Random real on open range (0,1) (i.e. not including 0 or 1)
      double genReal53(); // Random 53-bit real on [0,1)
};

MersenneTwister.cpp


#include "MersenneTwister.hpp"

/*
* Initialise mt[N] with a seed
*/
void MersenneTwister::init(unsigned long s)
{
mt[0] = s & 0xffffffffUL;
  for (mti=1; mti<N; mti++) {
      mt[mti] =
    (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
      /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
       * In the previous versions, MSBs of the seed affect
       * only MSBs of the array mt[].                   
       * 2002/01/09 modified by Makoto Matsumoto
  */
      mt[mti] &= 0xffffffffUL;
      /* for >32 bit machines */
  }
}

/*
* Initialise by an array with length given by key_length
*/
void MersenneTwister::initArray(unsigned long init_key[], int key_length)
{
int i, j, k;
  init(19650218UL);
  i=1; j=0;
  k = (N>key_length ? N : key_length);
  for (; k; k--) {
      mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
        + init_key[j] + j; /* non linear */
      mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
      i++; j++;
      if (i>=N) { mt[0] = mt[N-1]; i=1; }
      if (j>=key_length) j=0;
  }
  for (k=N-1; k; k--) {
      mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
        - i; /* non linear */
      mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
      i++;
      if (i>=N) { mt[0] = mt[N-1]; i=1; }
  }

  mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
}

/*
* Generate random number on [0,0xffffffff]
*/
unsigned long MersenneTwister::genInt32()
{
unsigned long y;
  static unsigned long mag01[2]={0x0UL, MATRIX_A};
  /* mag01[x] = x * MATRIX_A  for x=0,1 */

  if (mti >= N) { /* generate N words at one time */
      int kk;

      if (mti == N+1)   /* if init() has not been called, */
          init(5489UL); /* a default initial seed is used */

      for (kk=0;kk<N-M;kk++) {
          y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
          mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
      }
      for (;kk<N-1;kk++) {
          y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
          mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
      }
      y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
      mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];

      mti = 0;
  }

  y = mt[mti++];

  /* Tempering */
  y ^= (y >> 11);
  y ^= (y << 7) & 0x9d2c5680UL;
  y ^= (y << 15) & 0xefc60000UL;
  y ^= (y >> 18);

  return y;
}

/*
* Generate a random number on [0,0x7fffffff]
*/
long MersenneTwister::genInt31()
{
return (long)(genInt32() >> 1);
}

/*
* Generate a random number on [0,1]-real
*/
double MersenneTwister::genRealClosed()
{
return genInt32() * (1.0 / 4294967295.0);
  /* divided by 2^32 - 1 */
}

/*
* Generate a random number on [0,1)-real
*/
double MersenneTwister::genReal()
{
return genInt32() * (1.0 / 4294967296.0);
  /* divided by 2^32 */
}

/*
* Generate a random number on (0,1)-real
*/
double MersenneTwister::genRealOpen()
{
return (((double)genInt32()) + 0.5) * (1.0 / 4294967296.0);
  /* divided by 2^32 */
}

/*
* Generate a random number on [0,1)-real with 53-bit resolution
*/
double MersenneTwister::genReal53()
{
unsigned long a = genInt32() >> 5;
unsigned long b = genInt32() >> 6;
  return(a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
}

Thoughts on a Pebble... and the 100th post!

2008-08-14T17:47:00.004+01:00

This is my 100th blog post here at Physical Thought. To mark the occasion I've decided to write what promises to be a fairly lengthy textual post, in contrast to the large number of diagrams, images and photographs which have made their way onto the blog recently.

I've been thinking about software development recently, particularly in the field of high energy physics (hep). In the summer of 2007, I developed an event generator for double beta-decay events with colleague Sam Smith. It was written in C++ and the work was carried out as part of a summer research project with the particle physics group at Warwick.

Now, a year on, I've been thinking about that code and how it could be improved, retrospectively. With these thoughts in mind, I intend to write not only this blog post, but also to look into improving or rewriting the event generator, which we called Pebble.

I'm not saying that I'll definitely rewrite it, or that if I do I will come up with anything better than the original, but for the time being, here are some thoughts on the architecture of "Pebble v1", including what worked and what didn't.

The 'Mode' Functions
At the core of Pebble, probability distribution functions (PDFs) were sampled using Monte Carlo techniques. We implemented one- and two-dimensional PDFs as realisations of abstract interfaces for each dimensionality, and passed them through various mathematical algorithms to integrate (for normalisation) and find the minima or maxima.

On the whole, this approach worked well and provided the flexibility to support nine different decay modes with a single event generator function. The class hierarchy allowed us to separate mode selection from event generation, and the event generation code didn't have to care about the exact PDF being used.

There were also some problems (one or two of which were unavoidable). Notably, delta functions and straight lines caused problems which we had to build in workarounds for. Better core support for discontinuities or other weirdness in a mathematical function could have helped out here. Another issue was that, for the two-dimensional function, one could not normalise and find maxima until the first parameter had been fixed, which means that normalisation and maximum finding had to be performed once-per-event rather than once-per-run c.f. the one-dimensional PDF. This is unavoidable; you need to fix the first value before you can evaluate the PDF with a randomly chosen second parameter, but the overheads in initialisation, normalisation (i.e. integration) and finding maxima made this a bottleneck in Pebble v1.

The bottleneck mentioned above is perhaps the single biggest point of failure of Pebble v1, and any upated version ought to consider very carefully the instantiation overhead of the classes used in event generation, moving as much work as possible to an initialisation step which occurs just once, before event generation begins.

Mathematical Facilities
The basic mathematical facilities provided by Pebble v1 include integration of one- and two-dimensional functions, minimisation and normalisation. The original plan was to allow these to be templated and chained together such that one might be able to obtain a (pointer to a) function representing the normalised version of a raw PDF. This didn't work out so well in the original design, but the mathematical facilities themselves worked just fine, and carried relatively low overhead! Since normalisation of a 2D PDF has to happen for every event, the integration function really could have done with some more optimisation.

Input and Output
Program input was through command-line parameters while output was to tab-separated text files. Originally, we planned to create an Event object for each generated event, containing a Particle object for each particle involved in the event. In turn, the Particle objects would contain information about their energy and momentum (as a four-vector), particle type, etc. The creation of so many objects proved to be a huge time-sink and we eventually resorted to a flat data format, but the interface had already been geared up for this highly OO approach and it could have been handled better.

The PDFs form part of the C++ code and are compiled in directly. This means that, in Pebble v1, to add a PDF or modify its behaviour, the code must be completely recompiled. Unfortunately, it is hard to get around the need to compile the PDFs if one wishes to retain the event generation speed, but a library interface for adding additional PDFs as shared objects could have made things easier here without losing much (or any!) speed.

As for input, the eventual plan was to separate the interface and the event generation sufficiently to create an event generator library which could be tied to a large number of interfaces. This would have been possible, but we ran out of time. A better approach would have been to develop the library from the outset, adding interfaces as standalone programs which link (dynamically, ideally) to the library. The same point stands for output, for example providing a graphical interface or an interface to the CERN GEANT4 software for particle tracking, which was planned but time did not permit.

Thoughts on a possible Pebble v2
These are random bits and pieces that are currently swimming through my head, so I'll present them as a bulleted list rather than in full paragraph form.

Separate functionality from interface, completely. This should be done right from the beginning, and everything should be replacable without losing any speed or power. The key here is the ability to not only switch interface, but to switch out the PDFs, the integrator or the random number generator (for example) to obtain completely different behaviour from the same core codebase.
Don't overdo the object orientation; data is best left as data, not as Particle objects
Provide Particle objects "just in case"; some interfaces may require them. Ideally make these compatible with known software (e.g. GEANT4) so that there is minimal interface overhead in translating between disparate object formats.
Make extensive use of modern C++ features including namespaces, templates and the standard library containers / algorithms. These not only help the core programmer but also provide a standard set of tools for anyone working with the library.
Provide helpers for repetetive or complicated tasks, similar to the ModeMaker program in Pebble v1, which generated Pebble Mode Functions (PDFs) from a mathematical expression suitable for use in C++. Perhaps a Python (or similar) program to automate the generation of C++ functions from a mathematical expression, or even an interface with a Python interpreter allowing mode functions to be compiled into a shared library and dynamically linked to the currently running Pebble v2 backend. This would be slower (on first run) than the built-in mode functions, but would allow a fully interactive form of event generation.
Emphasise computational speed throughout. Spitting out a lot of events in a short time is a good thing!
Keep track of every last byte of memory. Valgrind is an excellent tool and the feeling you get when it tells you your program isn't leaking a single byte of memory is good too!

That is about all I can think of for now. I may play about with a bit of code over the next few days and see if I can come up with a basic framework on which to build something like Pebble v2... maybe.

J/Psi Decay Again

2008-08-13T21:17:00.004+01:00

So I haven't posted in a while... that is because I've been away sans computer (well, sans internet, anyway) for the past three weeks!

In the meantime, I've been reading Introduction to Elementary Particles (Griffiths) and I've found a Feynman diagram which really encapsulates the name of the J/Psi (well, the Psi part anyway) particle (charm-anticharm). These diagrams (although looking slightly like they should belong in an anatomy book...) are by far the nicest presentation of the decay of this particle that I've seen, made possible simply by using curved lines rather than straight ones!

NOTE: In order for these diagrams to look like the letter Psi, I had to adopt the convention used in the Griffiths book, which is for time to progress upwards. This is in contrast to my usual (and preferred) format for Feynman diagrams, in which time runs from left to right!

The top diagram is the OZI-suppressed decay of a ground-state J/Psi to three pions. The bottom one is the decay of an excited state particle to a D+ and a D-. Admittedly my drawing skills are somewhat lacking, but I didn't want to scan / photograph directly from the book, so this'll have to do!

(Click on the image for a larger version)

Some more complicated Feynman diagrams

2008-07-21T17:47:00.003+01:00

I decided to test FeynMF with a few more complicated diagrams (not necessarily because of the process involved, but because they're difficult to lay out nicely).

The first diagram is the strong decay of an excited state of the J/Psi meson. It has to be an excited state otherwise the c/anti-c doesn't have enough energy (mass) to produce the extra d/anti-d needed for the final state:

This diagram turned out to be easy enough for FeynMF to draw. The next one, however, is the decay of the ground-state of the J/Psi. This has to go via c/anti-c annihilation then production of more stuff from the resulting gluons. The ground-state of the J/Psi has a longer lifetime than typical strongly-decaying particles because this decay mode is OZI suppressed; the QCD coupling of gluons to fermions is a function of the momentum transfer. At low momentum transfer (equivalently, large separation) the coupling is large, but at high momentum transfer (or small separations) the coupling is very small. This leads to the asymptotic freedom that quarks enjoy within nucleons (where they are essentially free particles) but also leads to the relatively small amplitude for the following Feynman diagram, where all the mass energy of a c/anti-c pair is transferred through gluons to the final state. This is a high momentum-transfer process, so the QCD coupling is small, and we have three gluons, so it's smaller still!

Qualitatively, the OZI rule states that Feynman diagrams which can be cut in half between incoming and outgoing particles by cutting only internal gluon lines, and without cutting quark lines, are suppressed with respect to those which cannot.

FeynMF had trouble with this diagram. Well, that isn't true - it drew it just fine, but it didn't look very good. I had to change various tensions of the particle lines to get it to look like this, and it still isn't great. I am sure I could play some more and get it looking better, but since I wanted to test FeyMF's layout capabilities here, I didn't want to do too much over and above the "default" rendering it provides.

Finally, our old friend: Beta Decay! This is a simple process, down quark changing to an up quark, emitting a W- which turns into an electron and a neutrino. The difficulty arises in putting in the spectator quark lines for the other two nucleons present.

To do this I had to resort to the FeynMF "Immediate" mode, after freezing the d -> e + nu part of the diagram. With a bit more manipulation I could probably make this look much better!

Compton Scattering Feynman Diagrams

2008-07-21T15:04:00.010+01:00

I've started working through the problems in Griffiths' Introduction to Elementary Particles. One such problem is to draw the 17 fourth-order (that is, 4-vertices) Feynman diagrams for Compton Scattering.

The basic (2-vertex) diagram is shown below. An electron absorbs a photon and re-emits it at a later time.

NOTE: In all of my Feynman diagrams, time runs from left to right. The vertical axis does not, however, correspond to "space". Time runs left to right, nothing else matters.

I set about drawing all of the fourth-order diagrams, and eventually decided that I'd use them as an opportunity to learn FeynMF properly. FeynMF is a package for drawing Feynman diagrams in LaTeX (as the above diagram was, and several that I've posted previously). Until now, I'd been copying snippets from various examples to produce the simple diagrams I needed.

So I set about drawing all seventeen diagrams using FeynMF. And now I'm going to post them here, explaining what is special about each one (or rather, each related set). If you can see any reason that one of these diagrams isn't allowed, let me know. As far as I'm aware they're all valid, but I could be wrong. Seventeen is a lot of diagrams to keep track of, and although I don't think I've done the same thing twice, it is entirely possible!

The first set are related by the self-interaction between charged particles; QED loops. Here, each charged fermion can emit a photon and re-absorb it at a later time:

The next set involve interactions between the incoming or outgoing fermions and the propagator fermion:

And, since we can equivalently have Z0 exchange instead of photon exchange, the previous five diagrams with the photon replaced by a Z0:

Next,we must consider pair production from the incoming or outgoing photons. We can construct production/annihilation loops for fermions (either electron/positron, or quark/antiquark):

We also have to consider the fact that the W bosons are charged, so they can couple to the photon also. This allows us to have a boson loop in the pair production/annihilation:

And finally, we need to consider coupling between the initial and final state particles. The initial and final electrons can interact by exchange of a photon or a Z0, and the photons can interact by exchange of either a W+ or a W-:

Autostitch

2008-07-17T20:33:00.003+01:00

After reading about Autostitch on Pierre's blog, I decided to try it. I took around 27 photos of one of my bookshelves and let Autostitch loose on them. The result (after cropping the rough edges where it stitched the end photos together): pretty damn impressive.

I can't wait to try it on some landscapes!

Mathematica Notebook: Fourier Transform of a Cat

2008-07-16T21:07:00.004+01:00

UPDATE (April 19, 2009):
The Mathematica notebook is no longer available online. If you want it, you'll have to e-mail me and ask. It is left as an exercise for the reader to find my e-mail address.

I was asked in a comment on the Fourier Transform of a Cat post whether I was planning on uploading the Mathematica notebook. Since I can't upload it directly to Blogger, I have had to find an alternative hosting location. It's currently hosted in the physical-thought.com domain, which doesn't hold much else right now.

The notebook is now available in a ZIP file at http://www.physical-thought.com/files/fourier-cat.zip and the original post has been edited to reflect this. I stripped the output from the notebook to keep it small, so anyone who wants to use it will have to re-evaluate the cells. The original tiger.jpg is provided in the ZIP file too.

I also found the following image in the Mathematica notebook; it's a combination of Fourier transforms with some image data subtracted from it, and it's a rather striking visual effect.

Book: Introduction to Elementary Particles

2008-07-15T22:54:00.002+01:00

The other book I ordered, Introduction to Elementary Particles (Griffiths) arrived today. In what appears to be a consistent trend on the Amazon Marketplace, the seller claimed that the book was shipped from the UK in 3 to 5 working days, yet the packaging it arrived in bore postage stamps and markings from Korea. Nevermind, at least it did arrive, and arrived before I leave for Cornwall later this month.

From my preliminary flick through the pages, it looks like a very detailed, yet very readable book. The thing that impressed me most was the particle data tables printed inside the front and rear covers. Many particle physics books gloss over the exact values of quark masses, particle constituents, spin and various other important quantum numbers, forcing you to look online for a more comprehensive reference for anything beyond up and down quarks, electrons and the W or Z bosons!

The tables presented here seem to be very useful, and will prove valuable when I'm sitting on a beach working through the book, without access to a computer to google the mass of the something meson, or the quark content of a [Greek letter here] baryon.

Cover image below.

Close-Up Photography

2008-07-12T22:59:00.004+01:00

The company I ordered my step-up ring adapter from, SRB-Griturn, sell adapters which I'd never even considered before. Amongst them are adapters to reverse-mount a lens to a camera body, supposedly for close-up photography. The idea is that you take your lens off the camera body, and reverse mount it to the body using the lens filter thread into one side of this adapter, and the other side fits into the camera where the lens would normally go.

And they recommend you "try it first", by holding the lens in reverse up against the body and taking a photo of something.

Now, here is a photo taken using my Canon 18-55mm lens mounted normally. This is about the closest focus I can get to an object:

After taking that photo, I removed the lens, turned it round, and put it up against the hole in the camera. There are a number of obstacles to overcome here.

You lose autofocus functionality, meaning you have to manually focus the lens
The lens focussing isn't really sufficient in this configuration, so you end up moving the camera backwards and forwards anyway
While holding the lens against the body with one hand, and leaving the other hand poised over the shutter button, you need a third hand to change the lens focus
With the camera so close to the subject, light is an issue. There's not enough of it. I switched on a spot-light to get enough light in

After all that, I managed to take the photo below. This was with the lens rather close to the phone, and the results speak for themselves:

There's a fairly small area which is in focus. This is to be expected, given this configuration, and could probably be improved if I had enough hands. I suspect that is what these reverse-mounting adapters are for; freeing up the hand that is holding the lens, so you can focus with it instead!

Even so, I'm amazed by how close this gets you, and of the image quality resulting from it (fuzzy focus aside!)

Audiosurf

2008-07-12T18:42:00.002+01:00

As I'm sure most people already know, Audiosurf is a game which can be installed through Steam. It allows you to "ride your music" by generating a track from your MP3 files. The tempo and beat of the music alter the gradient of the track, and thus how fast you go, while other characteristics affect other aspects of the track generation. There are a number of different play modes, though my personal favourite is Mono Pro, in which you have to avoid grey blocks, and hit coloured blocks to get points. It's all fairly self-explanatory.

Anyway, a couple of months ago (i.e. in the middle of revision) my flatmates all started playing it, and I downloaded the demo and tried it on my somewhat old computer. It didn't work so well, so I gave up. Since then, however, I reinstalled Windows, as my previous installation was quite broken (to say the least). Since then, Audiosurf seems to work passably, albeit not in full-screen mode. As a result, I shelled out the $11 or so to buy Audiosurf, so I could continue to play beyond four tracks, and it's definitely great fun!

Here's a screenshot of the 'Freeride' mode, which lets you play music with the visualisation of the Audiosurf track (and even move the vehicle around) but you don't get points, and there is nothing to avoid.

Starting the Day

2008-07-11T11:42:00.004+01:00

I appear to have fallen into a well-defined routine every morning. It starts with getting up, whenever that happens to be, then showering of course.

After that, I get a bowl of bran flakes and cut up some fresh strawberries from the garden onto them, apply milk, and pour a glass of orange juice. I take the cereal and orange juice up to my computer, which I switched on after showering, and log in.

The next ten or fifteen minutes are a careful balance of eating and the following, roughly in order:

1. Start MSN
2. Start Xchat and log into IRC
3. Start Audacious (or Winamp, if I'm on Windows) and get some music playing
4. Open Firefox
5. Open the following websites:
- If (Monday || Wednesday || Friday) xkcd
- Gmail
- Facebook
- Questionable Content
6. Read xkcd and QC
7. Check email
8. Return pokes on Facebook, Decline being turned into a zombie, seeing how my friends rate me or other such crap, and finally close Firefox (albeit usually only briefly).

About this time, I manage to finish my cereal and go back downstairs to boil the kettle for tea or (if it's after 11am) coffee.

Upon returning to the computer with my tea or coffee, I'm ready to actually do something useful. At the moment I'm doing some paid programming work, which means working in Windows on Visual Studio 2008. Otherwise, I might play with some C++ or Python code, write a new class or plugin for Cerenkov or do something else I thought up.

I somehow manage to waste a lot of time during the day, so even though I have lots of ideas for projects, many of them are never finished and some don't even get started. This is a bit of a problem, since I'd quite like to see some of these things through to completion, but there are just far too many distractions around these days.

Google Analytics

2008-07-10T10:51:00.002+01:00

I own and operate two domains, one of which I don't particularly use (but plan to in the future) -- physical-thought.com. The other is nmap-tutorial.com, home to a comprehensive guide to the Nmap port scanner, and the number 1 result on google for the words nmap tutorial. It's also generally in the first page of results for the term nmap itself!

A while ago (a year or two) I added Google Analytics to the site, allowing me to track visitor trends, page views, unique visitors, their web browser and other similar data. nmap-tutorial.com receives an average of 3,000 unique visitors a month, and around 4,000 page views. 70% of that traffic comes from search engines, around 20% from referring sites (it is linked from the nmap documentation page, as well as several other locations that older versions of the tutorial have been posted), and around 10% direct hits.

Most visitors come from the US (1,206 in the past month), Canada (120) or the UK (223). India also features prominently with 160 visitors in the past month, and the rest of the world also played their parts.

A massive 58% of my visitors use Firefox, compared to only 30% on Internet Explorer. Opera falls into third place with 5% of the visitors, and Safari takes a mere 3%. Other browsers account for the rest but no single browser seems to dominate in the obscure end of the market!

A more interesting report is the Browser/Operating System combination. Over 34% of my visitors use Firefox on Windows. 30% use Internet Explorer on Windows, and 20% use Firefox on Linux. The rest aren't really worth mentioning. This was interesting mainly because I would have expected far more Linux users on a site like nmap-tutorial.com. Perhaps the "tutorial" nature attracts those just getting into security, and hence still clinging to their Windows operating systems!

Anyway, today I added Google Analytics code to this blog. There are two reasons for doing this; firstly I want to keep track of how many readers I have. My current mechanism for doing this is to use the Adsense page hits count, which is not a particularly accurate indicator since I may visit a page several times in the course of writing a new post.

Secondly, I'm intrigued as to the browser/OS combination of readers of this blog, given its fairly technical nature incorporating programming, physics and other such stuff, I'm expecting Firefox/Linux to feature prominently, but we'll see.

Filter Ring Step Adapters

2008-07-09T20:12:00.005+01:00

This morning I was running round Nottingham trying to find somewhere that sold filter ring adapters; these are devices which attach to the filter thread of a camera lens and provide another filter thread of a different size, allowing you to use smaller or larger filters than the lens ordinarily permits.

In my case, I wanted to be able to use my 62mm circular polarising filter on my 58mm thread lens, so I needed a 58-62mm step-up ring. This proved impossible to find. Jacobs had some adapters, but most were step-down (which surely would obscure part of the lens and thus be of very little use at all). The few step-up adapters they did have were the wrong sizes, so I moved on to London Camera Exchange, who informed me they did not keep adapters in stock, but could order them in.

Onwards to Jessops, who also do not stock them. It was here, however, that I received the most useful response; that a company called SRB, based in Luton, produced every kind of adapter imaginable, and took online orders!

So, after returning home, I Googled SRB Luton, and found the most useful photographic adapter producer in the world: http://www.srb-griturn.com/. They stock step-up and step-down ring adapters as well as a whole range of other more exotic things; adapters to make a macro lens out of two lenses, adapters for telescopes and microscopes... the list is almost endless.

I placed my order this afternoon, selecting 1st class delivery. Hopefully the adapter should arrive soon, and I can certainly recommend the website if you're looking for an adapter. In fact, if they don't have what you need, you can ask them to quote you a price for making it! It really doesn't get any better.