Lighthouse in the Sky

Old-school UNIX nerdery

2011-12-10T19:19:00.001-05:00

I've never owned a Windows machine, or a Mac. It's not that I'm particularly hardcore about open source, it's just that I got used to a command line back in the old VIC-20/DOS/DesqView days and never really got away from it. Our undergrad computer network was very convenient, but made up of about four flavours of UNIX, and in fact they had a load-distributing program at login that meant your initial login was not just to a random machine but to a random UNIX flavour. Somehow the IT folks made this work tolerably well. After a few years of that, when I finally got a home computer, I chose to run Linux on it.

Read more »

I for one welcome our new robot overlords

2011-12-04T17:40:00.001-05:00

A common theme in science fiction is the idea that robots will take over our world. This doesn't seem very likely, but suppose it were to happen. What would the first steps look like? How would we first begin to notice that it was happening?

Read more »

Police violence and crowd control

2011-11-26T11:52:00.001-05:00

In the last few weeks, we've seen a number of events where large crowds of unruly people have gathered. Police response to these crowds has varied substantially, and I think it's instructive to think about what these responses say about us, as a society.

Read more »

Getting health care

2011-11-24T08:12:00.001-05:00

Like many universities, McGill offers a student health clinic. Unfortunately, it's almost impossible to actually access the services, even when the support staff aren't on strike. To actually get help, you have two options:

If it's not urgent: there is one day a month, usually but not always the 15th, on which you can make appointments for the following month. Any other time, tough luck.
If it's urgent: you can show up at 7 AM, stand outdoors until 8 AM, and hope you were early enough to get one of the very limited number of drop-in slips. This gets you in to see a randomly-assigned doctor, fora few hurried minutes. If it's for a prescription renewal, the chances are slim that they will know anything about you or your condition. Just cross your fingers and hope they're one of the competent ones.

If neither of those options suit you, well, go somewhere else. At least this is free, unless you want dentistry or vaccination or medication.

Lucky imaging

2011-11-06T12:20:00.002-05:00

We have a telescope on the roof of the physics building. It's a fairly nice fourteen-inch Schmidt-Cassegrain telescope, though in fairly rough shape. Unfortunately, the physics building (and necessarily the telescope) is in the middle of downtown Montreal, which is a terrible place to look at the sky: bright lights, clouds, urban heat island, et cetera. So we have a telescope that collects lots of light but doesn't produce very sharp images. My attempt to work around this is described below.

Read more »

Breaking out of iTunes

2011-11-02T19:01:00.000-04:00

Since I got my "new" hand-me-down smart phone, I've been listening to a lot of podcasts. The phone is smart enough to use only my home wifi to update and download episodes, and the player keeps track of which ones I've read. This is so useful, in contrast to the music player, that I found it works better to keep my audiobooks on my home server and write a quick hack to serve them up as podcast feeds. But mostly I use podcasts like radio programs. Unfortunately, some podcasts are available only through iTunes. This makes sense if they cost money, but, for example, the NASA Lunar Science Institute has a free podcast which is only available through iTunes.

Fortunately it seems that the way these places get the data to iTunes is by serving up a standard RSS feed, which iTunes then wraps up in its proprietary glop. But Michael Sitarzewski helpfully put together a script that can extract the location of the original RSS from iTunes. So, for instance, you can subscribe to the NLSI podcast here. Very nice.

El Radar

2011-10-28T10:25:00.000-04:00

I had the opportunity to visit the great radio telescope at Arecibo ("El Radar") last year. It's an astonishing machine, built more like a stadium than a telescope. Sadly, my digital camera died on that trip, so I don't have much in the way of pictures, but I came across this neat youtube clip (above). It's a segment from a BBC program, talking about dark matter and the role of Arecibo in understanding it. I thought it made a nice view of the telescope. I especially like that in some of the audio you can hear the coquì, little tree frogs that provide a constant chorus. (Of course the telescope is also in Contact and some James Bond film — surprisingly the latter provides a better view of the actual telescope.)

I'll talk a bit about the telescope and my own visit, with a few photos, below the jump.

Read more »

The People's Microphone

2011-10-26T22:13:00.000-04:00

I went down to Occupy Montreal the other day. It was pretty impressive, really; the modest park, surrounded by banks, was filled with tents and people making art. Particularly amusing was the statue of Queen Victoria wearing a Guy Fawkes mask (pictured on the right). But what got me thinking was the People's Microphone.

Read more »

ssh shenanigans

2011-10-19T19:22:00.002-04:00

I use ssh all the time, from simple terminal connections to VNC sessions controlling telescopes to rsyncing terabytes of data to serving my music and video library. Most of the time, the machine I actually want to connect to is not visible from the open internet, and in some cases there are two layers of machines to go through to get to the one I want. I have used various ad-hoc combinations of ssh connections, but I just found a brilliant library that solves several of my problems: paramiko.

(Also, the name is from Esperanto, and the releases are named after gashlycrumb tinies. This is some serious nerdery, here.)

Read more »

Local Maxima and social annealing

2011-10-17T22:38:00.000-04:00

Dan's Data has a new post/magazine article up, about local maxima. Frustratingly, there's no comment section, but I've been reading some interesting things that seem relevant, so I'll post them here. (TL;DR of Dan's post: socially we seem to get stuck in local maxima instead of finding the real best way to do things.)

Read more »

Rest in peace, Dennis Ritchie

2011-10-14T04:38:00.000-04:00

Dennis Ritchie, who created the C language, has died. It would be difficult to overstate the importance of his work — C is imperfect and limited, but it serves as a near-universal tool for programming everything from tiny ATTiny8s to Blue Gene supercomputers. It does exactly what it was designed to do — provide a portable abstraction that is nevertheless close enough to the underlying machine to be efficient, and comfortable enough to write substantial programs in. These days, I (and many other programmers) prefer to write in Python or another high-level language, but most interpreted languages are still implemented in C, and when performance is critical, C is generally the way to go. Dennis Ritchie's impact on the world of computers was tremendous, and he will be sorely missed.

From: Dennis Ritchie 
Newsgroups: comp.std.c
Subject: Re: Computing sizeof() during compilation
Organization: Bell Labs, Lucent Technologies

> You are right.  It was nice back in the days when things like
> 
> #if (sizeof(int) == 8)
> 
> actually worked (on some compilers).

Must have been before my time.

 Dennis

Curve fitting part 5: PyMC

2011-10-11T17:31:00.000-04:00

I previously talked about fitting a curve to data — specifically, a sinusoid-plus-constant to a distribution of photon arrival times. I wrote some special-purpose Bayesian code to evaluate the relevant integrals directly so that I could work out the posterior distribution of fitted parameters. As a bonus I was also able to test the hypothesis that the photons were pulsed. Unfortunately, the method I used really doesn't generalize to more complicated situations. Fortunately, there is a toolkit that does: PyMC.
PyMC provides tools to implement Markov Chain Monte Carlo algorithms. These evaluate probability integrals numerically by generating samples from the distribution. They do this fairly efficiently by "jumping" semi-randomly from point to point in the space, but jumping preferentially towards points that are more likely. The algorithm is carefully arranged so that, after a certain number of steps, the initial guess is forgotten and the points are correctly distributed. They are correlated, but the correlations can be reduced or eliminated by subsampling the sequence. It's complex, finicky code. Fortunately I don't need to write it, because the algorithm is sufficiently generic that with PyMC all I have to do is specify the model I'm trying to fit. And, sure enough, it works. Code and results below the jump. The short version is: PyMC is easy to use and well-documented.
The second half of my previous post, model selection, turns out to be a little trickier to implement with Markov Chain Monte Carlo; I'll address that in a later post.
Read more »

Testing whether a signal is broad-band

2011-10-07T22:19:00.000-04:00

Radio pulsars are generally broad-band sources — you can hear them over a very wide range of radio frequencies, from about a hundred megahertz up to, in a few cases, tens of gigahertz. Their emission does change with frequency, decreasing as a power law, but over reasonable bandwidths we expect to see their signal in all frequency channels.

Radio-frequency interference ("RFI"), on the other hand, is very frequently narrow-band, appearing in just one, or a few frequency channels. In fact, one way we try to manage RFI is by using really wide bandwidths, so that there's so much power from the pulsar signal that narrow-band RFI is drowned out. Unfortunately, all too often the RFI is so strong that even a narrow-band signal can dominate the total power. And since it's narrow-band, dedispersion doesn't smear it out like it would a broad-band interference spike. So narrow-band RFI is one of the kinds of interference that is particularly hard to sift out from pulsar candidates.

Just recently, on the arxiv, a paper came out that attempts to address the problem of testing whether a candidate is broad-band:

Multimoment Radio Transient Detection
Laura Spitler, Jim Cordes, Shami Chatterjee, Julia Stone
We present a multimoment technique for signal classification and apply it to the detection of fast radio transients in incoherently dedispersed data. Specifically, we define a spectral modulation index in terms of the fractional variation in intensity across a spectrum. A signal whose intensity is distributed evenly across the entire band has a much lower modulation index than a spectrum with the same intensity localized in a single channel. We are interested in broadband pulses and use the modulation index to excise narrowband radio frequency interference (RFI) by applying a modulation index threshold above which candidate events are removed. The technique is tested both with simulations and using data from sources of known radio pulses (RRAT J1928+15 and giant pulses from the Crab pulsar). We find that our technique is effective at eliminating not only narrowband RFI but also spurious signals from bright, real pulses that are dedispersed at incorrect dispersion measures. The method is generalized to coherent dedispersion, image cubes, and astrophysical narrowband signals that are steady in time. We suggest that the modulation index, along with other statistics using higher-order moments, should be incorporated into signal detection pipelines to characterize and classify signals.

This looks very promising, but there's some testing I wish they'd done. More on the subject below the jump.

Read more »

Contraceptive underwear update

2011-09-28T19:40:00.000-04:00

I did a little more looking for information on the use of special underwear as a form of male contraception. In particular I went and found the 1994 Mieusset and Bujan paper. Unfortunately for some reason McGill's subscription does not cover online access to the International Journal of Andrology, so I had to go down into the literally dusty bowels of the library. To give you some idea what I was dealing with, I picked this no-doubt fascinating book off the shelves at random:

Don't worry, the actual paper is rather less ancient, though it was in an obscure enough journal to be in the rolling stacks — gear handles and all. But find it I did, and while they sadly don't have a photo of the contraceptive underwear, they do have a diagram. Since it is conceivably NSFW, I'll put it below the jump.

Read more »

ipython notebooks: second impressions

2011-09-27T15:53:00.001-04:00

I've kept working with ipython's notebook mode; it really is very well suited to my python-as-a-scratchpad style of figuring certain things out. So I've come up with a few more comments:

SVG plots work fine and look better, at least for modestly-sized output (thousands of points on a plot, megapixel images). They're also bigger on the screen, for some reason.
Be careful with commands that generate a lot of output; they can make your browser grind to a halt while it tries to render the output of that debugging print statement in that inner loop. During this time you also can't save your document. It might be nice if ipython noticed giant output and collapsed it.
It would sometimes be nice to be able to collapse blocks down to a little plus sign. For example, the first block normally needs to be the imports and maybe some ipython boilerplate (setting SVG mode, for example). It'd be nice to be able to hide it while working on the rest of the document.
It'd be nice to be able to break up class definitions into multiple blocks (for example to be able to have rendered math explanations for each method).
In PNG mode, the images are included using data URIs. This doesn't work on the very old web browser I have at work.
If you start pushing the memory limits of your machine it can be very hard to interact with the browser enough to kill the kernel.

For me, the most important of these is embedding the output in a blog post. It already sort of works: SVG is embedded in the output HTML, and PNG is encoded in a data: URI, so it's almost enough to strip out all the contents of the body tag and paste them into an HTML editor. Of course proper display depends on CSS, but I've added all the CSS that seemed relevant into the "custom CSS" field of my blog. The output still has some weird spacing issues; in particular, input blocks are all the same height for some reason.Read more »

Digital in the Facebook age

2011-09-25T15:00:00.000-04:00

Some time ago I posted about a lovely retrocomputing-style ren'ai game called "Digital: a Love Story". The same author has pubished a sequel, titled "don't take it personally, babe, it just ain't your story". It too is lovely, and has something interesting to say about privacy. If you haven't played ren'ai, be aware that it does have some genre conventions that take some getting used to, but it is a medium in which clever writing can produce a fascinating immersive story. In this case with dollops of facebook and 4chan.

Faster-than-light neutrinos: keeping time

2011-09-24T06:35:00.000-04:00

There's been a recent announcement of evidence that muon neutrinos may travel faster than light. That would be really weird - in particular it would pose serious problems for Einstein's theories of relativity. But even the discoverers aren't ready to claim that; they just describe their results and say they're puzzling. Personally I think it's unlikely they're right, but figuring out why not may be very interesting.

Faster than Light? par CNRS

Fortunately for us, they have posted a preprint of their paper on arxiv.org. All images below are from that paper.

I can't say much about the particle physics, or the details of instrument calibration, but I can address one possible way people have suggested the result may be wrong: inaccuracies in the time standards at the two endpoints.

Read more »

A captcha problem

2011-09-23T22:14:00.002-04:00

Hmm. Think it would cope with $ M_2(\mathbb{Z})\tilde m $? Or $M_2(\mathbb{Z})\tilde m$? Or just plain M_2(\mathbb{Z})\tilde m? In principle it might, since recaptcha is crowd-sourced; if I came up with the same LaTeX representation all previous recaptcha users had, all would be fine...

Male contraception

2011-09-20T19:00:00.000-04:00

The Pill — hormonal contraception for women — is sometimes credited with making possible the sexual revolution. There is as yet no hormonal contraception for men. My feeling is that this is for real biological reasons — the female reproductive system is more complicated and easier to interfere with. And of course in the era of HIV/AIDS and other STIs, condoms have a unique role. Other methods of contraception remain important. There's a new supposedly-reversible one for men and there's always vasectomy. But there's another method for men that I'd never heard of, but find amusing: testicular heating.
[update — pictures! sort of]

Read more »

Source Code

2011-09-19T03:56:00.002-04:00

No, I mean the movie. It's science fiction, of the Twilight Zone flavour: soft, set here and now, and with a little pointed topical relevance. Jake Gyllenhaal plays a soldier who gets sent into the "Source Code", eight minutes before a terrorist attack, and told that he must do whatever he can to figure out who the bomber is. He can retry it as many times as he likes. It's very like the brilliant interactive fiction Varicella (described below the jump).

The Core

2011-09-11T11:00:00.000-04:00

I've worked with data from X-ray satellites before, and among the many messy things one has to deal with in real data were blocks of time marked SAA. I knew this stood for "South Atlantic Anomaly", but I had only the vague idea that it was a part of the sky that was geomagnetically inconvenient, so that I had to trim it out of my data. The other day I came across a fascinating BBC documentary, titled "The Core":

This documentary talks about the Earth's core and how we're studying it, from seismology to diamond anvils to huge liquid-sodium dynamo experiments. It makes very interesting watching, but I particularly liked that they used the South Atlantic Anomaly as a hook: it caused problems with certain instruments on Hubble, and the documentary is framed as an investigation into why. (More below the jump if you're not worried about spoilers.)

Read more »

Mathjax

2011-09-10T16:14:00.001-04:00

When looking at the ipython notebook interface, I came across a neat tool: mathjax. It's designed to allow embedding math in web documents. Now, MathML was supposed to do this, but support for MathML is still very spotty, and it's also really nasty trying to write MathML manually. So mathjax allows you to simply enter TeX format math and have it look okay in supposedly all browsers: both inline as $e^{i\pi}=-1$ and displayed as $$\int_{-\infty}^{\infty}f(t)e^{2\pi i f t}dt. $$

The way that it does this is kind of amazing to me: it's essentially a TeX renderer written in javascript and running in your browser. I'm used to thinking of TeX as a batch-mode compiler that takes an appreciable time to run, but TeX has become the de facto standard way to write mathematics, and there are now several reimplementations of its renderer. Matplotlib has one, for example, though it can also call out to real TeX if you ask it to. But javascript!

Anyway, as you can see above, I've added it to the template for this blog, so I'll be a little freer with math from now on. Though I'm willing to bet it doesn't work in the RSS feed. Let me know if you have any problems with it.

Review: ipython notebooks

2011-09-10T05:38:00.000-04:00

The development version of ipython recently added a "notebook" interface mode. This resembles the interface of MAPLE or Mathematica, in which you can intermingle blocks of formatted text, blocks of code and output, and plots. My feeling, after working with it for a bit, is that it is very well suited to the kind of work I do when (for example) writing blog posts demonstrating some algorithm. For developing actual reusable software, not so much.

In order to give it a proper test, I applied it to a problem that a friend asked me about. It's about the card game Set: how often do you run into a tableau with no sets? The game quotes a figure of one time in 25, but that's for tableaux dealt from scratch. As you play the game, you remove and replace sets, and it often seems as if no-set tableaux arise more often when the cards are not freshly-dealt. The game's rules are a lovely minimalist mathematical exercise, but once you start removing and replacing cards an analytical solution becomes impractical. Simulation to the rescue! And this seemed like a nice test problem for ipython.

You can see the results here (PDF). For my comments on ipython, read on.

Read more »

Least squares and timing noise, part 2

2011-08-21T16:28:00.000-04:00

Simulated time series

In my previous post I described a new paper about fitting pulsar parameters in the presence of timing noise using a general least-squares method. It seems like a good approach, but I'd like to look at it in more detail. So: python to the rescue!

Read more »

The Mechanical Microbiologist

2011-08-01T21:13:00.000-04:00

Laboratory automation has been growing more and more elaborate - multiwell tools for "parallel processing", and benchtop robot arms, for example, that will handle repetitive tasks. But flow cytometry offers some startling possibilities.

Read more »

Least-squares and timing noise

2011-07-29T19:30:00.003-04:00

Figure 4 from the paper: residuals and spectrum

Working with long-term pulsar timing data sets is a nuisance because of so-called "timing noise". Not only is this noise above and beyond the usual observational uncertainties, perhaps because it is torque noise, it tends to be strongest at low frequencies (it is very "red"). Often so much so that leakage from the very lowest frequencies dominates at all analysis frequencies. Various people, myself included, have tried various approaches for dealing with this noise, but a recent arxiv paper shows real promise:

Pulsar timing analysis in the presence of correlated noise
W. Coles, G. Hobbs, D. J. Champion, R. N. Manchester, J. P. W. Verbiest
Pulsar timing observations are usually analysed with least-square-fitting procedures under the assumption that the timing residuals are uncorrelated (statistically "white"). Pulsar observers are well aware that this assumption often breaks down and causes severe errors in estimating the parameters of the timing model and their uncertainties. Ad hoc methods for minimizing these errors have been developed, but we show that they are far from optimal. Compensation for temporal correlation can be done optimally if the covariance matrix of the residuals is known using a linear transformation that whitens both the residuals and the timing model. We adopt a transformation based on the Cholesky decomposition of the covariance matrix, but the transformation is not unique. We show how to estimate the covariance matrix with sufficient accuracy to optimize the pulsar timing analysis. We also show how to apply this procedure to estimate the spectrum of any time series with a steep red power-law spectrum, including those with irregular sampling and variable error bars, which are otherwise very difficult to analyse.

I'd like to look at it in more detail and try some of the techniques on test data.

Read more »

Radio vortices

2011-07-19T13:24:00.000-04:00

Quantum mechanics, no one will be surprised to hear, is weird. In particular, photons can carry angular momentum - circularly polarized light can set objects spinning. But it turns out that light can carry orbital angular momentum as well. It's really not very clear to me quite what this means in terms of photons. In terms of classical waves, it's weird but I think I get it: if you look at the spatial distribution of the light in a beam, you may find that the phase is constant across the whole beam. But you might also find that the phase varies. Now, it has to be continuous, but you can imagine that as you make a circle around the beam center, you find the phase increases by an integer multiple of two pi. This gives you a continuous phase in a way that is topologically different from the constant-phase situation. As I understand it, this is what is called wrapping number.

Now this would just be another weirdness from the world of (classical!) waves except that there seem to be applications for it. In particular there's a paper on the arxiv about using this for communications purposes.

Read more »

Science fiction has no imagination, part 1

2011-07-19T12:02:00.000-04:00

Every so often I come across something that makes me think that the supposedly imaginative field of science fiction can't hold a candle to reality for weirdness. Today's installment is an arxiv paper in which the authors are seriously discussing quantum teleportation as a way to combine signals from telescopes to form an interferometer:

Longer-Baseline Telescopes Using Quantum Repeaters

Daniel Gottesman, Thomas Jennewein, Sarah Croke

We present an approach to building interferometric telescopes using ideas of quantum information. Current optical interferometers have limited baseline lengths, and thus limited resolution, because of noise and loss of signal due to the transmission of photons between the telescopes. The technology of quantum repeaters has the potential to eliminate this limit, allowing in principle interferometers with arbitrarily long baselines.

A little SF vignette

2010-12-08T18:00:00.004-05:00

Canada has these laws that require a certain fraction of the music that is broadcast to be Canadian. This caused a certain amount of screaming from radio stations, but I have this feeling that it actually works, that Canada as a result produces disproportionally more music, and more original music. I don't know about that (how would you quantify it?), but this band — You Say Party — is Canadian.

arXiv gleanings

2010-12-08T03:38:00.000-05:00

I try to keep an eye on arxiv.org because interesting new pulsar papers generally appear there first. But there are often abstracts that catch my eye, whether because they have neat ideas, because they talk about neat object, or because they just seem peculiar. This week's batch had a few of each.

How tempo2 does its fitting

2010-12-03T11:50:00.013-05:00

Pulsars can be difficult objects to study: for example, their radio pulses can randomly change in brightness, turn off, turn on, change shape, and we really don't know why. Nevertheless there has been some excellent science done by studying those very radio pulses. The trick has mostly been to simply not care how bright they are or what shape they have and focus on when they arrive. Since this comes from the rotation of the pulsar, this tends to be very regular. After all, for a ball 10 km across, with more mass than the sun, smoothed to within a millimeter by its own gravity, it takes an awful lot to change how fast it's spinning. What's more, time is the quantity we can make the best measurements of - world time standards drift by something like microseconds over decades, which is something like one part in ten to the fourteen. So pulsar timing is a powerful technique, that can measure pulsar positions and distances, spin-down rates and braking indices, and binary orbits. The standard software has been tempo, which is written in FORTRAN and has certain limitations. A new tool has recently appeared, tempo2, written in C++ and boasting good handling of timing effects down to the nanosecond level. Unfortunately the documentation on this tool is so far somewhat limited, so I've been figuring out how it works. I'd like to describe it, as best I understand it, here. This particular post will talk about how a timing solution is fit to a set of pulse arrival times.

Read more »

Millisecond pulsar X-rays

2010-06-19T11:00:00.004-04:00

Neutron stars are so tiny and so dense that the natural temperature scale for them has them glowing in the X-rays. What's more, they serve as powerful accelerators of electrons, which then naturally produce X-rays. So it turns out that X-ray telescopes provide a very interesting view of pulsars. In light of this, when we discovered my pet source, J1023, we took an X-ray observation of it. It's taken us considerable time to analyze the results and put together a paper describing them, but I think the result will be a valuable contribution to the literature. (More importantly for me, it should make a chapter of my PhD thesis.) The result is the cumbersomely-named "X-ray Variability and Evidence for Pulsations from the Unique Radio Pulsar/X-ray Binary Transition Object FIRST J102347.6+003841".

Much of the paper is devoted to details of data analysis, which I will spare you. But I think the gist is interesting, and not too hard to summarize.

Ignition!

2010-06-17T11:00:00.001-04:00

Recently, via Derek Lowe's In The Pipeline, I came across the book Ignition!, by John D. Clark. It's the story of the development of liquid rocket fuels, told by a man who was head of one of the programs. Even if you don't know much chemistry — I don't — the book makes a fascinating read. The image on the left is the first page you see when you open the book, and the image below is the second.

Ergo Proxy

2010-06-16T11:00:00.004-04:00

As is obvious to regular readers of this blog, I read a lot of science fiction. I also enjoy science fiction movies and TV shows, but they're generally much harder to find. Or at least, ones I like are. But I came across an enjoyable new series recently: Ergo Proxy. It's a cyberpunk anime with a female protagonist that I came to rather like in spite of her cold and sometimes difficult personality. The series is obviously some sort of cousin to Ghost in the Shell (albeit without some of Shirow's particular obsessions) and Serial Experiments Lain (though thankfully free of schoolgirls). On the Western side, I suppose I'd compare it to The Matrix (or maybe the similar-but-better The Thirteenth Floor) and The Prisoner.

As those analogies suggest, the series is sometimes a little heavy on the symbolism, mystery, and spouting of philosophy, but that's all part of the cyberpunk tradition, and this series carries it off well. I enjoy the visual style as well, though there are a few moments where the animation is a little off. The soundtrack is the kind of music I listen to anyway, which helps. I usually prefer subtitles to dubbed audio, but for this show the dub is better written and well-acted. All together a fun experience.

[Edit: you can watch the whole thing online.]

Reporting FAIL

2010-06-09T14:18:00.001-04:00

There's currently an item in the news about a British boy who managed to damage his eyesight with a "laser pointer". The news articles generally imply that any laser pointer is a lurking risk to your eyesight ready to blind you with a single incautious glance. But there are the ordinary 1-5 mW laser pointers you can buy in dollar stores, and then there are the pocket lasers you can buy online that have powers up to 500 mW and that can burn a hole in a credit card. If this kid managed to damage his eye with the former, well, that's surprising and alarming. But the articles don't describe the laser at all, beyond the fact that he ordered it on the internet. If you go to the original one in the British Medical Journal, though, you still don't get the power rating, but you find out that it was "high-powered", i.e. almost certainly one of the high-power ones you can get online. Of course you can damage your eyes with these; I strongly doubt it would be possible to get one without being made aware of the danger. What's more, the high-power ones are already arguably illegal in many places, not that that stops online companies from shipping them there.

The key point I'm getting at is that all the media coverage is missing the essential information that this is not a normal laser pointer, and that those are basically not dangerous.

[Edit: I emailed the corresponding author of the BMJ piece. Unfortunately the boy's guardian destroyed the laser before the doctor could see it, but the doctor is certain that it was one of the high-powered lasers I describe above. It's too bad none of the reporters bothered to check that detail, but I guess it doesn't make good copy.]

Virtual idol

2010-06-06T11:00:00.004-04:00

Years ago, I'm not sure when or where, I saw the anime Macross Plus, which is about a pop star who is entirely a computer program that sings, but who has fans as devoted as any present-day pop star. Anime being anime, there was of course a whole plot line about her AI going rogue (plus giant robots of course), but I remember thinking "how weird, a pop star who doesn't exactly exist". Well, technology progresses, and sure enough, we have one now.

Read more »

Laser rifles

2010-06-05T11:00:00.002-04:00

Laser rifles have been a staple of science fiction for years, though now that cheap harmless (or nearly) lasers are everywhere writers are moving to other buzzwords. But I had no idea that anyone had come up with a more-or-less workable design until I stumbled across the "Stavatti SF-1" and "TIS-1".

(More below the jump.)
Read more »

CASCA 2010 compact objects session 1

2010-05-29T19:45:00.000-04:00

Orbital eccentricity in numerical simulations of binary black holes (Harald P. Pfeiffer)

Eccentricity in black-hole binaries is radiated away faster than the semi-major axis, so that they become circular fairly rapidly. Thus, for example, when a binary makes its way into the LIGO band the eccentricity is down around 10^-6. For numerical GR this means you have to be able to measure - for that matter, define - eccentricity, and you need to be able to produce near-zero eccentricities. To define it, you often construct some function that measures deviation from a circular inspiral, but this is difficult because of (among other things) coordinate dependence. It turns out it works better to define eccentricity in terms of its gravitational wave effects. You also have to work on periastron advance (for which third-order post-Newtonian models are not good enough). How is very low eccentricity achieved? The problem is in constructing your initial data: coordinate location, sizes, velocities; these need to satisfy complicated PDEs to satisfy the constraints, and it's not clear what goal values you should choose to get circular orbits. The approach is to simulate a few orbits, then fix the initial condition as if it were Newtonian; this roughly works, so you iterate the process until you get a nice low eccentricity. When you go to precessing binaries (i.e. rotating black holes) spin-orbit coupling complicates your life, but you can still extract an eccentricity measurement and initial condition corrections.

Rayachaudhuri's Equation in Regge Calculus (Parandis Khavari)

These equations are important in the proof of singularity theorems, lensing, collapse, and cracking under self-gravity. They concern expansion and shear of fluids in self-gravitation. Absent vorticity they imply collapse into a singularity in a finite amount of time. Regge calculus is finite-element GR that fixes flat geometry within simplices; curvature is concentrated on n-2-dimensional sub-simplices ("bones") and can be described by the deficit angle. Geodesics in the Regge calculus are straight lines within simplices; where they meet faces, the angle on the entrance edge matches the angle on the exit edge, but it's a little tricky since there's no unique way to assign the deficit angle to the n-1 simplices around a particular bone. This work is about expansion of geodesics in (2+1)-dimensions. She obtains expressions for shear and expansion of lensing. Remaining problems include that Rayachaudhuri's equation has no unique discrete representation.

Numerical simulations of precessing binary black holes (Abdul Mroue)

Gravitational wave detection will only be possible if we have accurate banks of template gravitational waveforms. The three main categories are inspiral, merger, and ringdown. To build these templates we need some combination of post-Newtonian models and numerical relativity. We expect an event in the LIGO range at a rate <1/year, but advanced LIGO should have ~0.5/day. So building a template bank is crucial. The parameter space (for zero eccentricity) is given by the mass ratio and the two spins. A 15-orbit binary takes ~10^5 CPU-hours (plus a great number of grad student-hours). Less than 100 waveforms are available from all groups worldwide. So far very little work has been done on systems with generic spins. The spin has very significant effects on the system evolution. This work has two major approaches: make BBH runs easier (i.e. reduce the person-hours) by automating the initial setup and transition between regimes, and make BBH runs faster, primarily by making their simulation code run on GPU supercomputers (possibly an order of magnitude speedup). Currently running on his desktop's GPU.

Questions: Do we expect spin-orbit alignment? Maybe; what happens in the final orbits is totally unclear. [No comment on whether binary evolution is likely to produce aligned binary systems.]

Modelling Gravitational Lens Systems with Genetic Algorithms and Particle Swarm Methods (Adam Rogers)

Lensing has been observed on many cosmological structures; it's interesting because on the one hand it depends on the mass of the lens, and on the other hand it provides an important magnification effect. The goal of this research is to try to reconstruct the un-lensed image of the source. They assume a thin lens. The lensing equation is clearly nonlinear, as shown by multiple lensing. Rather than solve a complicated nonlinear equation, one can simply raytrace past a lens shape. One can combine this information into a mapping matrix; in this formalism, finding the source pixel intensities is a linear least-squares problem. Unfortunately the matrix sizes are comparable to the number of pixels [so you need sparsity], so you need to use a small PSF, i.e. optical data. With some cleverness one can reformulate the problem into a deconvolution problem that never needs to construct the huge matrices. On simulated data it's quite effective at recovering even highly distorted images when you know the lens parameters. Finding the lens parameters requires a fitting procedure, for which he's using genetic algorithms and particle swarms. Particle swarm optimizers attract each particle to its local best and the global best according to a spring force. Of the two the particle swarm optimizer is a little faster, while the GA is much more thorough in searching the parameter space (e.g. local loses one of the two 180-degree options for ellipse orientation).

Questions: Have you applied your models to real data? Yes, but it's old data that's already been partially processed. What is he optimizing? "chi-squared between the model and the data"

New properties od teh 35-day cycle of Hercules X-1 (Denis Leahy)

RXTE ASM data of Her X-1, plus PCA data when available. Her X-1 is 6.6 kpc and high galactic latitude; it's a neutron star with a 2.5 Msun A7 companion. Over 35 days (many binary orbits) you get brightness variations. One model is a twisted disk that occults the NS; its shadow on the companion star changes on the same cycle, which explains the optical variations. ASM data shows that the cycle length varies substantially - 34-38 days - and this variation is correlated with the flux. The turn-on appears to occur uniformly in orbital phase. They have 1.58 Ms of PCA data in total. One idea for explaining some of the irregularities is that the impact point on the disk is not the outer edge, since it's not flat; other models include an uneven disk surface or blobs in the accretion stream, but these don't match the data well.

Questions: Why is the disk twisted? Heating by the central source produces a torque that increases as the disk gets out-of-plane up to the point where it starts shadowing itself.

CASCA 2010: Stars 3 session

2010-05-29T19:44:00.000-04:00

Variable stars in the Hyades Cluster (Jayme Derrah)

To establish distance scales you need to measure variable stars at known distances; the Hyades cluster is the closest open cluster to the Earth. Unfortunately it's just out of range for accurate parallax measurements, so distances are measured using the "convergent point method". Variable stars of interest re eclipsing binaries and pulsational variables (though few new examples of the latter are expected). The project uses the Baker-Nunn patrol camera, which has a large enough field of view to include the whole cluster and reaches a limiting magnitude of 19.5 in two minutes; on the other hand there's only one filter, a light pollution blocker. Follow-up of discovered variables will use B and V.

Type 1A Supernova Progenitor Diversity (Ashley J. Ruiter)

Sub-Chandrasekhar-mass WDs should be looked at as SN1A progenitors. SN1A light curves are driven by the amount of nickel produced, so we can use them as standard candles almost without any idea of what the progenitors look like. Theories are CO WD mergers and CO WD collapse triggered by companion Roche lobe overflow. Simulations predict that the former ought to produce NSes instead, while models (including this work) have difficulty producing enough of the latter. In particular, she built a population synthesis code and found an order of magnitude too few overflow models. The idea is that perhaps allowing sub-Chandrasekhar-mass collapse can help. More detailed models, in particular using helium accretion, allow collapse to occur with 0.8-0.9 Msun: you get detonation in a shell of accreted helium which triggers the explosion. She built simulated light curves based on such explosions and found reasonable agreement with observed light curves [but how hard is this, if nickel mass is all you see?]. From a population synthesis point of view, though, they do provide enough. It's also theoretically nice because it explains the observed variety in light curves: it's a function of mass at collapse.

Questions: Where do these helium-rich companions come from? Many are He white dwarfs, which arise naturally in such binaries (lost their hydrogen through Roche lobe overflow). What's that blue line on the graph? It turns out if you start with 0.9 Msun WDs the merger scenario looks a little more plausible. What assumptions go into the pop. synth? Standard IMF, etc. Where does the mass go during binary evolution? Eddington-limited transfer, so probably outflows.

Tau Sco: The Discovery of the Clones (Veronique Petit)

Part of the MiMeS (Magnetic Field in Massive Star) collaboration. There ~35 known magnetic OB stars (i.e. B directly detected, which is hard), but it has long been suspected that B is ubiquitous in massive stars. ("Magnetic fields are to astrophysics as sex is to psychology." - someone) These stars are short-lived but important, and the key complicating parameters are rotation and mass loss; we expect B to substantially affect the wind (as it does in the Sun). Tau Sco, a B 0.5 V star, is the focus of this talk. We can measure B by the Zeeman effect, i.e. circular polarization across a spectral line. If you sample the stellar rotation, you can actually produce B maps. Unlike all other massive magnetic stars which have roughly dipolar fields, tau Sco's field is complex and multipolar. Studies of the wind using UV line profiles show strange inconsistencies between different line traces; there seems to be some correlation between the UV variability as a function of stellar rotation and the B. It's natural to assume that the wind anomalies arise from the complexity of B, but it was hard to test until their recent discovery of two "clones" that have similar wind anomalies and slow rotation. ESPADON shows that these stars are magnetic; sadly it has not yet been possible to sample the stellar rotation densely enough to map the magnetic field. Preliminarily, though, a dipole model suggests that surface B is roughly equal to that of tau Sco.

Questions: Have people measured this supposed correlation between B and the wind? Yes, but it doesn't really predict tau Sco. How strong are these Bs? ~500 G. Does B of tau Sco vary with time, particularly in flares? We don't know, they had to assume it depended only on rotation. Could you use the Bayesian inference to estimate non-dipolar fields? There's too much uncertainty when you don't know about the stellar rotation.

Tracing Wolf-Rayet wind structures (Alexandre David-Uraz)

This work focuses on WR113. WR stars have very high mass loss rates, high enough that radiation pressure alone can't explain it; we need to understand the physics of clumping. WR113 = CV Ser is a binary system in which the companion can often be seen through the WR wind. It's a double-line binary, so they can in principle get a decent picture of the orbit, but so far the spectrum has been too messy. Once this is accomplished, they should be able to separate the two spectra quite well - average everything in the WR frame, then subtract this and switch to the companion's frame, then back and forth. In any case, the "eclipse" is of depth ~0.6 mag - well, in Lamontagne et al. from the 90s; in a 1963 paper the eclipse seems variable (0.1 mag sometimes and 0.6 mag), and in 1970 they didn't see an eclipse at all. MOST sees a shallow eclipse but variable - in fact, two successive eclipses differ very substantially. The real focus of the research is random variations due to clumping, and indeed they see both photometric and spectroscopic evidence for clumping. There's also the issue of a colliding-wind shock, which shows up in the lines (in particular there's streaming along the shock); with luck this should help establish system geometry. The next step is Fourier analysis for pulsations in star or wind, and wavelet analysis for random clumping. The goal is to link these to spectroscopic data and constrain the clumping.

Questions: This is a complicated system; is there evidence for an accretion disk? Unclear, but the colliding winds suggest not. Is this the only eclipsing WR star? No, there are others.

Limb-Darkening and Stellar Atmospheres (Hilding Neilson)

Limb darkening is particularly interesting these days as observations become more constraining. They affect planet parameters inferred from transits (and fitting can constrain limb darkening). Optical interferometry can directly measure it, and microlensing can help as well. Limb darkening tells you something about conditions in stellar atmospheres; Schwarzschild used solar eclipse observations to show that the solar atmosphere is in radiative rather than adiabatic equilibrium. Traditionally, though, people just use empirical limb darkening laws that are pretty crummy; in fact popular parameterizations have fixed points that are more or less independent of the parameters, and fit observations rather badly. The fixed point arises because the models the empirical laws are based on all make the Eddington approximation. Spherical symmetry rather than plane-parallel atmospheres helps spread the fixed point out. Detailed modelling suggests that the true value at the fixed point probes atmosphere physics. Observations suggest that you can make these inference even with the wrong models.

CASCA 2010: planets session

2010-05-29T19:42:00.002-04:00

NIR from hot Jupiters (Bryce Croll):

Hot Jupiters are tidally locked, and one expects their properties to depend on the amount of heating. He tested this with four known hot Jupiters. To get good photometry, he used WIRCAM out-of-focus (!). This lets them detect the thermal emission from hot Jupiters (millimagnitudes). "Hot" and "hotter" distribute heat very effectively, though "hotter" shows abnormally low H-band emission; maybe an absorber? "Even hotter", on the other hand, doesn't distribute heat very effectively. "Hottest" shows redistribution within the hot side but not between sides, plus evidence for a thermal emission deep in the atmosphere. Many more CFHT observations are planned.

Music of the Spheres: Results from the Kepler Mission (Jason Rowe):

Kepler's a 1m space telescope designed for steady monitoring of a fixed field of 100000 stars, looking for exoplanets and asteroseismology. CCDs were placed so that the brightest stars fall into the gaps. Kepler's bandpass is more or less white light (400-900 nm). Measurements as good as one part per million for many stars; Kepler's working very well. The parameters you fit directly include stellar density. The smallest planet so far is ~Neptune-sized, but some secondary transit measurements show you could detect a terrestrial transit. In asteroseismology they should get p-modes for every star brighter than ~12 mag. Also some weirdos like helium white dwarf companions, for which radial velocity followup has begun. Apparently Kepler also has Guest Observers and Target of Opportunity arrangements. Quite what this means I'm not sure, since the instrument isn't pointed.

In questions, they do have lots of earthlike transit candidates - but they get lots of false positives from photometry alone, let alone grazing transits by bigger objects.

SPI or Spin-up? An UV Investigation of Activity on Exoplanet Host Stars (Evgenya L Shkolnik)

"SPI" is "star-planet interaction", meaning magnetic. Spin-up is tidal effects of the planet on their stars. These should occur since the planets are so close to their stars. Both should show up as increased stellar activity (hence the UV). Looking at X-ray data, there's debate about whether this is observed. She used GALEX (175-275 nm). GALEX sees almost every exoplanet host star (that's in its field of view). Given 135 GALEX detections, start with close-in planets seem to be ~2x brighter. (Planet distance is a bimodal distribution, so there's a natural dividing line here; we think it's migrated versus not migrated.) A few transiting systems have been observed to have abnormally rapidly rotating stars; if this is a general phenomenon, gyrochronology goes out the window. The exceptions, stars that aren't active in spite of close-in planets may be explicable in terms of some scattering models but not others, and in fact the exceptions are more eccentric, consistent with the potentially recent scattering models.

Super-Earth Transit Search with the MOST Space Telescope (Diana Dragomir)

Start with radial velocity candidates from HARPS; F, G, K with 2-20 Mearth planets. MOST watches during the transit window. Key transit parameters are period, phase, duration, and depth. Diana's idea is to use Bayesian methods to fit the model, which is feasible mostly because the HARPS candidates provide a drastically reduced parameter space. (Ew! She uses JDs!) This parameter space needs to be very finely sampled to avoid underestimating the depth. (She's using grid sampling rather than MCMC because of the small low-dimensional parameter space.)

Ultra-Wide Trans-Neptunian Binaries (Alex H. Parker)

You classify TNOs based on their resonance with Neptune: the plutinos, for example are roughly 3:2 times the period. The big question is how the TNOs got there, since they haven't had time to form in situ. One answer is that as the giant planets moved out to their current positions they scattered the planetoids out into the outer solar system. This talk is about binaries; binaries are present throughout the solar system. Binary asteroids are different, as they tend to have large mass ratios and tight circular orbits, suggesting a collisional origin. The TNOs have lots of binaries, with a wide range of eccentricities and often modest mass ratios. These binaries are interesting observationally because you get mass measurements. Wide binaries are resolved with Gemini, and they're so delicate you can test dynamical ideas with them. Their formation mechanism is debated, maybe three-body exchanges, maybe temporary capture into chaotic orbits, but the mutual orbit distribution should allow them to be distinguished. His sample is >0.5'' apart and mags differing by <1.5. Typicallt separations ~1'' and mags ~24. Terrestrial parallax allows measurement of things like the direction of the orbital axis. The observed distribution doesn't really fit with either formation model. You can also use these to look at the number of ~1km objects in the Kuiper belt, since these objects will eventually break up such a system. These small objects can't be made through accretion or collision, they have to be primordial.

HST Compositional Survey of Faint Kuiper Belt Objects (Wesley C. Fraser)

The idea is to look at chemical gradients in the protoplanetary disk; this is closely related to understanding their scattering history, which is what dominates their current dynamics. Compositionally, only the biggest can retain methane (as it evaporates on a timescale of thousands of years). The fainter ones you don't get much in the way of lines, so the data is mostly colour (some are neutral, some are quite red). In the IR you do see deep water ice on certain objects ("family members" in resonances?) but not on others. It's a little puzzling what these objects are actually made of. Fortunately WFPC-3 is almost ideal for this work, with filters for water ice, methane, and a couple of useful narrow-line filters. The colors show evidence of a mix of water ice and some red gunk in varying proportions. Centaurs, objects that have been scattered further into the solar system show evidence for solar processing. Other objects are more mysterious.

Searching for Main-Belt Comets Using the CFHT: Final Results (Alyssa Gilbert)

Main-belt comets reside in the asteroid belt but show cometary activity. Their origin is unknown; formed in situ, which would be weird, since we don't expect ice within the orbit of Jupiter, or maybe formed further out and somehow got stuck in the asteroid belt. Only five objects are known, and have been recently found. Small bodies are interesting because they form in different parts of the solar system, and they don't have sesimology or weather, so they provide fairly direct pictures of the structure of the early solar system. To find these objects, she used the CFHT legacy survey data. The cadence is 3 observations in one night followed by another a night later (weather permitting). Automated object selection not terribly effective; visual inspection of 25000 objects worked better and found one. (Somebody has a very high tolerance for boredom...) Unfortunately the object wasn't noticed until a year after the observations, so even though a rough orbit was found the object was lost. Asteroids activated by collision might produce an appropriate number, but the activity of some of these objects seems to be periodic, which is weird. Capture from the outer solar system should provide too many.

Light echoes

2010-04-21T23:04:00.000-04:00

Occasionally some astronomical event — a gamma-ray burst, a supernova, a magnetar flare, or whatever — will go off, lighting up the sky more or less spectacularly and then fading. Normally that's the end of that event, and we're left to wait for the next one. (In particular, this means that if you study this sort of thing, you need to be prepared for one of these things to happen at any time, so that you suddenly need to write target-of-opportunity proposals, analyze data, and release preliminary results in a tearing hurry. This typically happens while you're supposed to be on vacation.) Once in a while, though, we see a peculiar phenomenon called "light echoes".

A "light echo" arises a little like an ordinary (sound) echo: some event happens producing a bright flash (loud noise) and in addition to the light (sound) making its way to you directly, some of the light (sound) goes in a different direction, bounces off something, and makes its way from the other object to you, arriving a little later. From the place I usually stand to watch the summer fireworks competition, you hear the big skyrockets go off, then a second or two later you hear the echo from a nearby building. Of course, on human time and distance scales, the light from the fireworks reaches us instantaneously, so it's obvious that both the original sound and its echo are delayed. In an astronomical setting, we only receive light, and it takes very much longer. But it's still possible to receive a delayed echo, and studying these echoes can be very informative.

(Photo, courtesy of ESA, to the right is X-ray dust echoes around the magnetar 1E 1547.0-5408, one of the objects people in the group here at McGill study. This interesting dust-echo work is from another group, though. The echoes are from a massive X-ray outburst, which we think was caused when the extremely strong internal magnetic field stresses cracked and twisted a piece of its crust; this twisted the external magnetic field, and the twisted magnetic field produced and accelerated massive numbers of electrons and positrons, which blasted out a torrent of X-rays. At least we think that's how it happened; we saw the torrent of X-rays.)

Read more »

Spuds

2010-04-15T19:00:00.000-04:00

From the arxiv: The Potato Radius: a Lower Minimum Size for Dwarf Planets, by Lineweaver and Norman.

This article is kind of neat, not least for the charming way it uses "potato" (and its adjectival variant "potatoid") as a technical term. The point of it is to try to work out how big an object has to be for it to be round.

Ablation

2010-04-14T19:00:00.003-04:00

From the arxiv: Momentum Transfer by Laser Ablation of Irregularly Shaped Space Debris, by Liedahl et al.

I don't have a whole lot to say about this one, but how could I ignore a paper about Giant Space Lasers? They're talking about using them to clean some of the junk out of low Earth orbit. While it would take a great deal of power to completely vaporize space junk, all you really need to do is give it enough of a shove (~100-200 m/s) that it starts to dip into the atmosphere (~200 km altitude), where it will slow down and burn up. Just how much of a shove you can get by zapping it with a laser so that some evaporates is not easy to predict, hence the paper.

The particular kind of Giant Space Laser they're talking about is left to some degree unspecified, but it's clear that you want short pulses, so that you get explosive vaporization (gas flow velocity of ~1000 m/s) rather than gentle heating, and they're talking about 10 J pulses (producing 0.1-1 m/s change in velocity for a 1 g target). So it's an awful lot of short powerful pulses. They also mention, in the usual understated scientific way, the possibility of "structural modification" of the debris — that is, the possibility that the bolt or paint flake or whatever will be blown to pieces or bent out of shape by the laser (in addition to the ~10% of the mass that will be outright vaporized). They suggest laboratory experimentation, to which I say, can I help zap random pieces of junk with a high-powered laser and see what happens? Please?

(Photo to the right is of a NASA laser ranging experiment, not actually zapping space debris. Unfortunately.)

Space Race

2010-04-12T23:01:00.000-04:00

In honour of Yuri's Night, I'd like to point to this rather nice BBC documentary:

It casts the space race as a personal competition between Sergei Korolyov and Werner von Braun, and is mostly recreations and dramatizations (including the peculiar decision to have the Russians speak English with a Russian accent, except for commands like "Ignition" and "Launch", which the actors could presumably learn and are subtitled). But it's generally well made, and doesn't shy away from either von Braun's Nazi past or the military applications that drove the space race itself: a rocket that could put a man into orbit could also deliver a warhead to a major city, and vice versa. From a public relations point of view, announcing a new milestone in the space race generally went over better than announcing one's ability to more effectively kill millions of people.

Chemistry at home

2010-03-27T17:55:00.001-04:00

YouTube user TheHomeScientist (via In The Pipeline) is posting a series of videos about what you can do in a home chemistry lab; as a nice example, there's this beautiful one about the purification of hydrochloric acid:

Isn't this method elegant? No complicated boiling or fumbling around with strong acids; just take advantage of the fact that HCl is a gas.

I think it's great to show people that you can really do chemistry at home. Science is not just the domain of white men in white coats with PhDs. On the other hand, I'm hesitant to get too ecstatic about how "democratic" this is. Not just anyone can afford the time, energy, and space to set up a lab like this.

Read more »

Parmesan oregano bread

2010-03-25T14:00:00.002-04:00

This is a very tasty bread that's rich enough to eat by itself, rather than as a substrate for some more flavourful substance. It also smells wonderful:

300 mL water
50 mL sugar
15 mL oregano
250 mL grated parmesan (or romano)
5 mL salt
750 mL white flour
20 mL olive oil
10 mL yeast

I brought this to our neutron star discussion group and it was a hit.

Gallium surface gunge

2010-03-23T14:00:00.013-04:00

Gallium, as a safer liquid metal than mercury, has two major drawbacks: the surface rapidly acquires a layer of dull oxides, and when in contact with surfaces it tends to leave a thin layer behind. I've been doing some reading to see if I can work around this problem, and I came across this enjoyable article (which may be behind a paywall).

Read more »

Digital: A Love Story

2010-03-22T00:35:00.000-04:00

For those of us who remember the days of BBSes and 2400-baud modems, Digital: A Love Story has considerable charm. It's a game in a genre — ren'ai — with some conventions that take some getting used to, but it's clever and original enough to be fun even if you weren't using computers in 1988 and haven't played ren'ai before. Recommended.

(It's also in python, using a framework, renpy, that has been a major success.)

Real interbinning

2010-03-21T19:00:00.013-04:00

(This is another of my highly-technical "note to self" signal processing posts. I'll put up something less arcane soon.)

The Fourier transform is great for finding periodic signals: you take an FFT and a periodic signal looks like a peak in the output. Well, in an ideal world, that is; you only really get such a neat and tidy peak if the periodic signal is exactly at a Fourier frequency, which happens when it makes exactly an integral number of turns over the course of the data set. If the signal is somewhere between two Fourier frequencies, the power is spread over several Fourier output values. While it's possible to interpolate a very accurate value based on the ~32 nearest values, this can be expensive, and there's a shortcut called "interbinning"; it doesn't reconstruct the phase, but you just take a scaled average of the two neighbouring bins and get a decent approximation to the value at the midpoint between two independent Fourier frequencies.

My problem, for this post, is that the theory is all nicely worked out when going from the time domain to the frequency domain, but I want to do something analogous while going from the frequency domain back to the time domain, if that's possible. (I haven't done a literature search, or even looked very carefully at the frequency-domain interbinning papers; I thought this would be a good exercise for me.)

Read more »

bash pipes

2010-03-20T22:18:00.002-04:00

I've been using UNIX for a very long time, and bash has been my shell for almost all that time, but for the life of me I can never remember how to pipe standard error anywhere. I think my problem is I've never found any logic by which the syntax makes sense. Anyway, it's:

command 2>&1 | less

But if you want to send the standard error and output to a file, it's the other way around:

command > file 2>&1

ISS in the X-band

2010-03-20T15:00:00.006-04:00

Pretty Picture: ISS in the X-band - The Planetary Society Blog | The Planetary Society

This spooky-looking image is the international space station as taken by a radar satellite, the German TerraSAR-X.

I actually worked with these radar satellites, though they are now my natural enemies. The basic way they work is very clever, and with appropriate analysis, you can extract things like tiny movements of glaciers from the data.

Glitches and flares

2010-03-19T15:00:00.047-04:00

Recently on the arxiv: Searching for X-ray Variability in the Glitching Anomalous X-ray Pulsar 1E 1841-045 in Kes 73, by Zhu and Kaspi.

Pulsars normally spin down very regularly — like clockwork, as the saying goes, and many pulsars spin down as regularly as a good atomic clock. But some pulsars, once in a while, will suddenly start spinning more quickly. This sudden (instantaneous as far as we can measure) spin-up is called a "glitch", and its full explanation remains mysterious. Generally, all we see is that the pulsar is suddenly spinning faster: no heating of the crust, no sudden X-ray emission, no radiative changes at all, just a suddenly-faster pulsar.

Anomalous X-ray pulsars (AXPs) are one kind of "magnetar", pulsars whose magnetic field is so enormous that its decay powers the X-ray emission of the star. They exhibit many peculiar behaviours, and are a major field of study in pulsar research. AXPs will occasionally become much more active for a while: they become much brighter, they emit random blasts of X-rays, and they do other peculiar things. It seems as if they may glitch every time they become active like this. If we want to find some sort of relationship between glitches and these active periods, it would be valuable to know whether an active period happens every time AXPs glitch, or whether AXPs sometimes have "quiet" glitches, like normal pulsars. That's what this paper tries to answer.

Killer radiation

2010-03-16T16:00:00.020-04:00

I've been working in the undergraduate labs lately, checking the setups for the undergraduate laboratory project course. The course has any number of interesting projects, like spatial filtering of images using Fourier optics, measurements of the Hall effect, demonstrations of Rutherford scattering (after all, he had his laboratory in the building that is, now that it has been decontaminated, our library), and so on. The setup for the experiment on Compton scattering, pictured to the right, has all sorts of terrifying warning signs all over it, because it necessarily uses a powerful gamma emitter. Just how dangerous is it?

This is not purely a theoretical question, since the head lab technician caught a student looking down the beam at the unshielded source, prompting the numerous warning signs. So just how much harm did that student do themselves?

The source is about 100 millicuries of cesium 137, and emits 661.6 keV gamma radiation. In SI units, that's four billion becquerels (decays per second). But what does this translate to in terms of exposure (measured in sieverts)? This is a more complicated calculation, since the kind of radiation matters, as well as the geometry and time of the exposure.

For the purposes here, though, I'll just point to a list of gamma ray dose constants, which give the exposure rate in rem/hour for a one curie source at a distance of one meter, and note that in these units cesium-137 has a dose constant of about 0.4. So a hundred millicurie source produces roughly 0.04 rem/hour, or 0.4 millisievert per hour.

At this rate, if you stood there and stared into the source for two and a half hours, you'd get Health Canada's recommended yearly radiation limit for the general public. It'd take a hundred and fifty hours to get the limit recommended for people who work with radiation. At two hundred and fifty hours (if that still counts as "acute") you might raise your cancer risk by 0.8%.

All this is to say that while I don't think it's a good idea to look into the source, and certainly not to touch it, and I wouldn't want to work with it every day, the warning signs are perhaps a bit over-emphatic.

Honey oatmeal bread

2010-03-15T14:00:00.002-04:00

I think this is my favourite bread machine recipe so far; light and fluffy, a little sweet, and very tasty.

Read more »

Flops and the FFT

2010-03-14T13:57:00.003-04:00

The Fast Fourier Transform is a wonderful family of algorithms that have made whole fields of study possible. They mostly operate by breaking a large problem into some number of equally-sized smaller problems, all the way down to problems of just a few points. As a result, the prime factorization of the number of points can make a large difference in the speed of the calculation. Traditionally, powers of two have been the fastest (or the only) sizes available, so the wisdom has been that if you have a data set whose length is not a power of two, you should pad it to the next larger power of two. It turns out that, at least in some cases, that is a mistake.

For context, I am still thinking about the problem of coherent harmonic summing in pulsar searches. As I worked out before, it's important to have an accurately-interpolated Fourier spectrum, and it turns out that the simplest way to do this is to simply pad the input array before taking its Fourier transform.

For a two-minute observation, we are taking about twelve thousand power measurements per second (yes, this is basically an audio data rate; but it only appears once we have combined all 4096 frequency channels), so the raw data is about 1.48 million data points (the exact length of our observations varies slightly). We need to pad this to eight times its length to get good enough sampling in the frequency domain, so that's about 11.8 million data points. The question is how much more to pad this in order to make the FFT fast.

I'm using FFTW for the Fourier transform, and it has efficient algorithms for powers 2^a3^b5^c7^d. So while I could pad my FFTs to the next power of two (about 16.8 million), I could also pad them to 12 million exactly, or to the smaller value of 11854080, which still factors nicely (2⁸3³5¹7³). I know that those powers of five and seven are not as efficient as those powers of two, but on the other hand, flops are cheap, while increasing the data size by almost a factor of two potentially means an awful lot more memory access (obviously these giant FFTs aren't going to fit in cache, though FFTW is smart enough to arrange the pieces for the best possible cache coherence). So I thought I'd do a comparison.

The results are:

16.8 million points: 1.7 s
12 million points: 0.8 s
11854080 points: 1.0 s

That is, the power of two is the slowest — in fact it took almost twice as long as the nice round decimal number. Trimming the value further, to 11854080, slowed things down again, presumably because the reduction in memory size was minor, but the prime factorization was worse. The upshot of all this is that, at least for large FFTs, padding to the next power of two is not necessarily a good idea.

Moon rot!

2010-03-13T11:00:00.001-05:00

Recently on the arxiv: Long-term degradation of optical devices on the moon, Murphy et al. This paper talks about the retroreflectors left on the Moon by the Apollo and Lunokhod missions, and observes that they have dropped in effectiveness by a factor of ten since they were placed. So far from the moon being a hostile but static place, something has been steadily degrading these mirrors.

Among the things the Moon missions left behind were arrays of retroreflectors. Like street signs, bicycle reflectors, or those weird-looking radar octahedra, these take incoming light and beam it back to where it came from. These are useful scientific tools, because you can fire lasers at them and time how long it takes the pulses to come back, in the process measuring the Earth-Moon distance. The current best setup, APOLLO, measures the distance to the nearest millimeter, which lets us test theories of gravity, detect a liquid core to the Moon, and watch the Moon recede from us (well, at 38 mm/year).

The moon is very far away (unlike the International Space Station, which orbits at an altitude of only a hundred kilometers). So when you beam a laser at it, even if you use a telescope to collimate the beam, it's spread over 7 km when it reaches the Moon. Then any imperfection of the retroreflector, or simple diffraction, spreads the return beam over an even larger area (20 km); your telescope picks up as much of the returned light as it can, but APOLLO sends out pulses of 10¹⁷ photons and gets back only about one photon per pulse. These tremendous losses are a challenge, so the authors of this paper (who work on APOLLO) monitor the efficiency of the system.

What they noticed, spurring this paper, was that the efficiency dropped substantially — by a factor of fifteen — near the full Moon. Now, obviously the full Moon is very bright, so background photons make detection more challenging, but it is easy to measure the background and estimate how much harder it makes detection; this effect is far from enough to explain the dip. By itself, a dip at full Moon isn't that exciting; since the Apollo 15 retroreflector is pointed at the Earth, a full Moon is when the Sun is illuminating it nearly face-on, so thermal effects might explain it (and in fact since it works by total internal reflection, it only serves as a reflector out to about 17° from the vertical, while the dip is about that wide, at 30°; the authors don't mention this, so it may be coincidence).

To investigate this dip in efficiency, though, the authors of the paper went back and looked at older observations. In the initial years, lunar laser ranging was done with the 2.7m Macdonald Observatory Smith Telescope. But rather than compete for time on this large telescope, in 1985, the program switched to using smaller, dedicated telescopes. Unfortunately these smaller telescopes couldn't see the reflections near full Moon, so there's no data from 1985 to 2006, when APOLLO went online. But comparing the APOLLO data to the MST data, they find that the dip at full Moon was not present at first, and gradually grew as time went on. What's more, if they looked at the return rate away from full Moon, they found a uniform decay; right now the retroreflectors are returning about a tenth what they did initially.

So what's happening to the retroreflectors? How are they getting worse? It's obviously not rain or wind, the moon being devoid of either, let alone the kinds of organic decay you see here on Earth. As with almost all real science, the authors cannot offer a definitive answer, but they do discuss some possibilities.

First of all it's worth pointing out that the huge drop in efficiency doesn't mean the reflectors are absorbing all that extra energy; they are almost certainly reflecting it but in the wrong direction. They're cube corner reflectors, and if you distort the shape of a cube corner reflector it doesn't reflect light back in quite the same direction it came in. The authors find it would take about 4°C temperature difference across each cube to produce the full-Moon losses they see. So if for some reason the reflectors are now absorbing a small fraction of the blazing unfiltered sun at full Moon and being heated by it, that could explain the dip in effectiveness. But why are these mirrors absorbing more and more of the sunlight?

The authors' most plausible answer, to my eye, is dust. The surface of the Moon is covered with dust, made by thermal weathering of rocks and by micrometeorite impact. This dust does not of course blow around the way terrestrial dust does, and in a vacuum a tiny dust grain should fall as fast as a rock, so it initially seems difficult to explain how much dust could get on top of the reflector. There are micrometeorites, though, and there is an effect I hadn't heard of before: dust particles become electrostatically charged through irradiation and are either levitated or thrown upward in "fountains" by electrostatic repulsion. We think. What we do know is that observations from the lunar surface show the optical effect of dust above the ground. So however this dust gets up there, some of it can plausibly fall on optical equipment left there by astronauts.

A layer of dust on the surface would also naturally explain the general decay in effectiveness even when not being heated by the Sun. So it looks like perhaps things left exposed on the lunar surface get covered with dust fairly rapidly. This I find interesting in its own right, but there are also various plans to build telescopes on the lunar surface, since these would share the advantage of the Hubble space telescope of images undistorted by atmosphere, while being able to rely on gravity to hold things in place, and possibly even being able to be built out of lunar materials. If they rapidly become coated with dust, those plans will have to come up with some scheme for cleaning the optical elements.

Multigrain bread

2010-03-12T12:38:00.003-05:00

My second attempt at a multigrain bread in my newly-borrowed bread machine.

Ingredients:

250 mL water
15 mL olive oil
15 mL sunflower seeds
15 mL flax seeds
15 mL pumpkin seeds
15 mL quinoa
15 mL millet
5 mL salt
325 mL bread flour
250 mL whole wheat flour
30 mL honey
10 mL yeast

Incidentally I need to get more containers for all this flour. It used to be something I bought the smallest possible bags of because I never finished it, but with this bread machine I go through a lot.

(That odd-looking lump on the loaf above is a piece of dough that didn't get properly mixed in. This happens from time to time with this bread machine; it actually pauses during the mixing and beeps at you to give you a chance to push any such lumps back into the mix, but I wasn't paying attention to this loaf. Most loaves don't need this extra attention.)

Bread machine

2010-03-12T12:35:00.000-05:00

I'm not a great cook; I usually don't have the energy to cook anything interesting, so I tend to fall back to the very simple. (Plus Montreal has zillions of great restaurants, so it doesn't take much excuse to make me go out to eat.) But fresh bread is delicious, so when I visited my parents and they had a bread machine scavenged from the kerb, I was happy to give it a try. I was impressed, but felt like I would probably use it for a little while and then it would sit cluttering up my kitchen. So when some friends mentioned that they had a bread machine they weren't using, I asked if I could borrow it.

I figure, they used it a lot for a month, then got bored of the bread and put it in storage. If I use it for a month and then get bored, well, maybe they'll find someone else to lend it to afterwards. And maybe after a while they'll have a craving for bread again; seems like a fine solution. (Of course, bakeries are a perfectly fine solution too, and indeed I do like to buy a nice Première Moisson baguette from time to time.)

What sold me on a bread machine was not the fine bread it could make — though I did have some nice loaves — but how easy it was to make bread. Now I can have put off going to the grocery store for far too long, so there's no longer anything perishable in the house, but toss five ingredients in the machine, hit a button, and lo and behold, fresh bread:

250 mL water
60 mL olive oil
15 mL sugar
5 mL salt
750 mL flour
8 mL yeast

This uses the "basic" program and produces a medium-size loaf.

There's even a "French" bread recipe that works if I'm out of oil:

250 mL water
60 mL sugar
8 mL salt
750 mL flour
8 mL yeast

This uses the "French" program, which kneads more and lets it rise longer, and produces a chewier bread with larger bubbles. Since all bread machine loaves are approximately cubical, this doesn't much resemble a baguette (not that that's the only kind of French bread), but it's still pretty good bread.

The biggest problem with the machine is that it takes three hours to make the bread, and about half an hour before the bread is ready it starts smelling delicious.

Blind Zen programming

2010-03-11T11:00:00.002-05:00

A recent alarming network problem while trying to observe remotely at Green Bank brought back memories of the age of the 2400 baud modem, when you could often type faster than the system could display your characters. Fortunately I use the text editor vi, which was designed for these conditions. In fact, in this light some of the user interface decisions make a lot more sense.

The most prominent weird UI decision is that the editor has two "modes": insert mode, in which you can actually type text and it goes into your document, and "beep mode" in which the letters on the keyboard all tell the editor to do something (and it beeps if that something doesn't make sense). For example, to move around your document, the keys hjkl move you left down up and right, respectively. Bizarre as it may seem, these keys are on the home row, so if you're doing a great deal of moving around — particularly common when programming, i.e., debugging — this makes a certain sense. (Plus vi users get so hardwired that these become perfectly rational movement keys in video games.) The decision to have an editor-control mode also means that you don't need to press the "control", "alt", or "meta" keys very often, which makes sense since these often didn't work or were inconsistently placed on old UNIX machines. (I spent a while switching between terminals, some of which had "Caps Lock" and "Control" interchanged, and to this day I disable "Caps Lock" on every keyboard I use.)

The use of ordinary (case-sensitive) letter keys for editor commands also means that there's room for lots of editor commands, and they're very quick to use. So, for example, to search forward, you use "/"; "n" takes you to the next match and "N" to the previous. But there's also "f", which searches forward for a single character. While apparently bizarre, when combined with the ability of commands like "delete text" (d) to take an argument, this means you can do "delete up to the next single-quote" with "df'".

All this makes even more sense when it might be seconds from when you type a key to when you see its effect. So if, for example, you're typing along and you realize "hey wait, that 'a' was supposed to be an 'i'", instead of pressing the cursor keys many times and waiting to see where the cursor wound up, you can compose sequences like "Fari$", which will search back to the a, replace it with an i, and return you to the end of the current line. It's a bit like being the blind Zen archer.

All that said, I must make a peace offering to those in the emacs camp: being written not so much later, it offers similarly complex editing functionality, with similarly unusual (but different) UI decisions. I used emacs for years, switched to vi for some reason or other, and have stuck with it through inertia, mostly. (Even in the ubiquitous text boxes on the web I find myself hitting Escape followed by strings of gibberish. I'm not quite hardcore enough to go figure out how to make Firefox use vi as its editing component.)

RFI

2010-03-09T19:30:00.007-05:00

I'm currently involved in another low-frequency pulsar survey with the GBT. This is a pointed survey, so we take two-minute pointings in an array tiling the whole northern part of the sky. This is an incredibly sensitive setup, so we will hopefully find a nice collection of new pulsars. Unfortunately, the same arrangement that makes us sensitive to pulsars makes us sensitive to lightning, electric motors, the local power lines, passing aircraft... all sorts of signals, generally outside our intended piece of sky and/or outside our intended frequency range, but also generally vastly more powerful than the pulsars we're looking for. So dealing with this junk — generically called RFI (Radio Frequency Interference) — is an important part of our survey strategy. My topic today is not so much dealing with the mild RFI we see all the time as dealing with what we see every now and then: a tremendous signal comes booming in, overwhelming everything and ruining the data (as seen on the right).

First of all a brief summary of how our system works: The telescope has a receiver at the primary focus, basically a pair of crossed dipoles. This converts the electromagnetic waves to signals in the wires, where they're amplified, encoded in analog on light and sent down an optical fiber to the control room, where it's converted down to a standard frequency range. This is fed to GUPPI, a so-called "pulsar backend". This instrument does analog-to-digital conversion on the signal then uses a polyphase filterbank to efficiently split it into 4096 channels, then computes the power in each channel; this power is averaged over 81.92 microseconds, after which the set of 4096 powers is dumped to disk.

Pulsars are broadband signals, so we use the widest bandwidth our receiver can handle, about 80 MHz centered near 350 MHz. This gives us the best signal-to-noise from our pulsars, but it also means we can't stick to the few narrow bands reserved for radio astronomy. In fact, the region we're using is allocated for various kinds of radio operation, including "aeronautical radionavigation" around 330 MHz. The GBT is in the National Radio Quiet Zone, so with luck much of this spectrum might be clear, but if aircraft are really using it to navigate there's no keeping them out of our airspace. (And at least once when I've complained of RFI the GBT operator has noticed an airplane passing to the north, for what that's worth.) In any case, there is interference to deal with.

Old radioastronomy systems, in order to get the widest possible bandwidth, quantized the input signal to one bit, or three levels, on initial digitization. This has the serious drawback that any small change in the signal levels causes the digitizer to overload or underflow, trashing the output. Newer backends improve on this by using 8-bit digitization, so that you have more than twenty decibels difference between underflow and overflow; it also helps that they are sampling a wider bandwidth, so that it takes more power to substantially change the levels. So generally, when a noise signal is picked up by GUPPI, it is narrowband and produces excess power in one channel. But the polyphase filterbank is good at keeping it from spilling into neigbouring channels (much better than a simple FFT) as long as it doesn't overflow the system, so usually you see something like this:

This shows, in red, the average power in each channel, and in blue, the standard deviation in each channel; since what we are receiving is basically noise, the standard deviation in each channel is large. The overall profile is the sensitivity of our system as a function of frequency. You can see it's not flat over the 80 MHz we want and zero elsewhere, but some complicated bumpy shape; that's the joy of analog components. In fact, given that it's an 80 MHz band at only 350 MHz, we're doing well to have a filter even that good. The point of this plot, though, is those vertical lines. Those are channels where some narrowband signal is coming in, so strongly that it stands way out above the noise. But it's okay; we have four thousand channels, we can afford to throw away a handful, and we have tools to detect when to throw them away (clever tools that look not just for strong signals but also for weaker signals that are narrowband and periodic, since those are a real headache for pulsar searching). So this plot is of one of our good beams.

The beams I'm worried about are the ones that look like this:

Here what's happened is we have some narrowband signal that is so monstrous it's overloading our digitizer. The digitizer's saturation is a nonlinear effect, so it doesn't stay politely in one channel; instead it splatters all over the place, producing both ugly broad peaks as on the right, and (I think) the regularly-spaced vertical lines (though they might be something else). When this happens there's not much use to be had from the data; while it's possible we might detect a pulsar through all that muck, on the one hand we'd be vastly less sensitive, while on the other if the interference has any periodic component to it, it'll produce zillions of false-positive periodicity candidates to wade through. So we need to discard beams like this one.

Fortunately, since we're doing a pointed survey, if we find that a beam has been trashed by RFI, we don't have to just write it off and leave a gap in our survey coverage: we can just put the beam back on our to-observe list. No problem. All I need to do is come up with a scheme for detecting these bad beams automatically, since we often run our survey observations semi-unattended (once they're up and running, I often give the operator my cell phone number and ask them to call me if something like a wind stow or a snow dump comes up; one of my colleagues will often do this and then go to sleep).

Initially, we were planning to run the RFI detection and excision scripts as part of our survey processing. Unfortunately, pulsar survey processing is incredibly CPU-intensive - we are just now finishing off the processing for the drift-scan survey, whose datataking finished in mid-2007. Clearly if we wait until the beams are processed to decide which are RFI-laden, the survey will long be over before we find out which beams we should have reobserved. Fortunately, the RFI detection code isn't too compute-intensive; an eight-processor modern machine running jobs on all processors ought to be able to process eight two-minute beams in about fifteen minutes.

The RFI detection code isn't perfect, though: I found that it was able to detect overloaded channels, but a few overloaded channels are normal. When too many are overloaded we can discard the beam, but I found that there were beams where not very many channels were actually overloaded, but those few splattered into neighboring channels. So my prototype system does an additional test: it compares the bandpass of each beam with the instrument's normal bandpass. This sort of works, but it revealed some surprising things:

This is a different plot of the same sort of thing. The red curve is the bandpass of the system, estimated by taking the median of all the day's beams. The green curve is one of the horrible contaminated beams. But the blue curve is from a beam that is probably fine. But its bandpass looks decidedly different (enough so to trick an early generation of my script into marking it bad). It turns out that what has happened is that this is a beam shortly after a contaminated beam. While the RFI was going on, we ran a system "balance", which adjusts the gain/attenuation at many points along the signal path so that in each stage it's at a comfortable level. Since there was powerful RFI, the gains and attenuations all got readjusted, and all differently. Since different parts of the signal chain affect the bandpass differently, we got a different bandpass. (We also got a drastically reduced signal amplitude, but fortunately GUPPI gives us plenty of dynamic range, as I mentioned above.) So I have to be a little careful when comparing bandpasses so as not to reject minor changes like this. Fortunately the "splatter" from bad RFI is pretty obvious in the statistics, so now I think I have a working bad beams detector. I'm just waiting for feedback from my collaborators before putting it into service.

Some comments on the tools I used to build it: I wrote the code in python and numpy, of course. But key to the process was a module by the author of the RFI detection tool that let me construct python objects from the detector's statistics files. Given this, I used medians to construct a bandpass for each file (since medians are better than means at discarding crazy outliers, which is the whole point of the exercise). I then used masked arrays to flag any bad channels, and scale the result so its median is one. I repeat this process on many beams, then take a median (using the masked array median to nicely ignore any bogus data points). This gives me my system bandpass. Comparing an individual beam to the bandpass then proceeds by constructing a masked median of the file as before, scaling it so its median matches that of the system bandpass, and then counting the points where it is very different.

In all, it has proved invaluable to have the masked array tools; they just do the Right Thing with bad data, vastly simplifying my code.

Artificial gravity round 2

2010-03-09T00:35:00.001-05:00

This alarming gadget, a Lava Lamp Centrifuge demonstrates some of the problems I discussed in my post on artificial gravity:

(via How to Spot a Psychopath)

As the builder puts it:

The centrifuge is a genuinely terrifying device. The lights dim when it is switched on. A strong wind is produced as the centrifuge induces a cyclone in the room. The smell of boiling insulation emanates from the overloaded 25 amp cables. If not perfectly adjusted and lubricated, it will shred the teeth off solid brass gears in under a second. Runs were conducted from the relative safety of the next room while peeking through a crack in the door.

He doesn't mention that lava lamps are full of liquid that is hot and flammable and in close proximity to electricity.

He also discusses how he supplies power to the lava lamp: he wired up a quarter-inch phone jack to 120 V AC, noting that the connector can be rotated freely while still making a connection (until the contacts wear out, presumably, not being designed for any of constant rotation, 120 V, or any appreciable power). He avoided other connections by using a battery-powered accelerometer and video camera.

Gallium

2010-03-07T21:07:00.007-05:00

Through the magic of ebay, I bought some gallium. It's strange stuff. Apparently whether it's listed as liquid or solid on periodic tables depends on where the table is printed; the melting point is 30°C, so it's solid at room temperature if the room's in Canada in March. But it'll melt in your hand, though it's a slow process.

Gallium is a crystalline solid; I suppose many metals are, but the crystals are really obvious when gallium solidifies. I thought I'd take a video of gallium crystallizing, but it has a tendency to supercool, so after sitting at room temperature for hours it was still liquid. I dropped a crystal of gallium in, though, and I got this beautiful slow crystal formation:

This video is shown at twelve frames per second, each frame is 60 seconds of real time. (It starts when it does because that's when I realized nothing was going to happen immediately; it ends when it does because that's when my camera overheated (!).)

Those vague angular patterns on the surface are actually crystals forming underneath. When I tipped the dish so the liquid flowed away I saw this:

Unfortunately, gallium is directly below aluminum on the periodic table, so, like aluminum, it reacts very rapidly with air, forming a sticky surface scum. When gallium is liquid, though, this scum can't stay in place to protect the surface; instead it sticks to everything around it. Rolling gallium through your fingers feels very peculiar — it's decidedly denser than water, though not tangibly more viscous, and it doesn't feel cool (its vapor pressure at room temperature is tiny). But because of the oxidation, it leaves a gray scum all over your hands. Pieces of gallium left in air also quickly start looking dull and dirty.

Cyberpunk gadgets

2010-03-07T11:00:00.001-05:00

Maybe this dates me, but I remember when cyberpunk was the hot new kind of science fiction. It replaced the utopian or social-experiment future societies with one in which the cancers of our own grew unchecked - corporate rule, environmental devastation, and urban decay were the future. The characters and stories tended to be gritty and ambiguous, computer hackers, drug pushers, or hit men (or pizza delivery boys, yes). My big complaint was that somehow in every story the hero has to Save The World from some quasi-magical and universal threat, be it AI, computer viruses that afflict humans, what have you. My point today, though, is about all those high-tech cybernetic implants the characters always have.

I mean, okay, who wouldn't want to be smarter, or stronger, or to remember everything on the Internet, or to be able to sense magnetic fields? Well, okay, maybe not everyone; in fact I'm just as happy having most of that with gadgets I can carry around and replace when they break. (Particularly if, as in most cyberpunk stories, upgrades and repairs happen in filthy little underground clinics.) But let's leave that aside; what I've been thinking about is whether such gadgets make sense at all.

First of all, implanting anything in the human body is a tricky business, but we can do it. Hip replacements are amazingly successful; their biggest problem is that since the replacement hip is not repaired by the body, it can wear out. Since you need to remove a few centimeters of thighbone to take it out, you can't do this very many times. What about more complicated implants? Well, the immune system can be a problem, since it tries hard to destroy anything that seems alien (even, unfortunately, sometimes parts of the body; often this is the reason hips need to be replaced). The immune system uses powerful peroxides and chlorine radicals to destroy things, so even quite chemically-resistant materials can break down eventually. The lady who implanted a magnet in a finger so she could sense magnetic fields found that after a few years, the magnet had been broken up and reduced to powder. But as artificial hips and pacemakers show, these issues can be managed, with care. So it is possible to put things in the body and have them last.

One thing that is a big problem, though, is the skin. The skin is a very carefully-maintained barrier against the environment. All the usual routes into the body are very carefully guarded by systems ranging from a continually-replenished layer of mucus in the nose to our tendency to flinch away from anything getting in our eyes. Any new opening in the skin, say a small cut or scratch, must be carefully kept clean until it heals, and even so mild infection is common. The body's response to infection is to send swarms of immune cells to destroy anything even vaguely suspicious in the area. So if you want to have some sort of implant with a plug or tube connection to the outside, you're going to have to devise some way to prevent infection at that hole in the skin. People do have this sort of implant — chest tubes, catheters, and so on — and infection is a constant problem. Fighting it is made particularly difficult because bacteria form biofilms adhering to the surfaces of foreign objects, so that even if a treatment kills all the surface bacteria, it must still penetrate them to reach the bacteria underneath. So if at all possible an implant should avoid piercing the skin.

Is this possible, for the kind of electronic gadgets that cyberpunks get installed? I think so, more or less. There's no need to have an electronic connection to transmit data, as Bluetooth headsets demonstrate. Power is a more difficult problem; electronic gadgets do draw power, sometimes quite a few watts. A sufficiently advanced technology would let the surgeon hook up a little artery and vein, and then run off the sugar and oxygen dissolved in the blood. But chemical interactions with the bloodstream on the scale needed to power an electronic gadget open up a massive can of worms — what kinds of other chemicals will be unintentionally exchanged into or out of the bloodstream? How can you exchange chemicals with the bloodstream without exposing yourself to immune system attack? How do you maintain vascularization without risking clotting? Bluetooth headsets and the like currently use batteries, but for an implant you have to worry about how they're recharged (unless maybe you use plutonium). My suggestion is to use magnetic induction — like cordless electric toothbrushes, you put a coil in the implant and a matching coil on the charger, so that when you bring them close they form a transformer and you can feed power in. This has its own alarming failure modes (overheating, overloading, stimulation by unintended machinery, interaction with magnets), but it will work.

The next question, of course, is what do you actually need an implant for? Frankly, most of the things cyberpunks use them for are now available for the iPhone. Or, if you like, the oddly creepy wearable computing gadgets. (For that matter there's even a non-implanted version of the magnetic field sensor.) Exceptions I can think of are gadgets that interact with the bloodstream or the nervous system directly. The nervous system I can sort of see being useful, but it's incredibly complicated, doesn't heal much, and messing with it is extremely invasive. So that's pretty daunting. Dealing with the bloodstream is more reasonable; there are already partially-implanted gadgets for diabetics to try to manage blood sugar. While this sounds like a great idea, cimpletely implanted gadgets would have a finite reservoir of insulin to work with, so they would need to be replaced or refilled regularly. Unless we figure out how to build a gadget that can make insulin from blood components, that problem's not going away. Genetic engineering offers possibilities - I can imagine a little gizmo that contains a few of the patient's own cells that can be zapped to persuade them to produce insulin on demand. Immune system issues are going to be something of a challenge, particularly if it turns out that a patient's diabetes was caused by their immune system attacking their insulin-producing cells in the first place. On the other hand, if you can clone and grow cells that produce insulin, why not let the body's natural regulation run things without the need for any implant beyond the cells themselves?

In summary, I think that implanted hardware will always be very costly, not just in economic terms but in terms of the user's health and in terms of maintenance. Given that, there would have to be a very strong need for them that couldn't be met using other, safer and cheaper tools. Shame that, I always liked Molly's scalpel claws.

Cyclotron! Isotope! Logarithme!

2010-03-06T17:00:00.000-05:00

One of the French Wikipedia's more amusing pages: Vocabulaire du capitaine Haddock. (Sadly, it is no longer named "list of Captain Haddock's insults".)

For those of you suffering from cultural deprivation, Captain Haddock is a friend of the young reporter Tintin. He is an old sea captain, and is therefore often drunk and forever cursing. Rather than the more usual grawlixes, Captain Haddock's cursing is inventive and often bizarre, calling people things like "bashi-bazouk", "macaque", or the three in the title of this post. (I don't know what the English translations are like, but I assume they kept the colour and variety of the insults. My school library only had the books in French.)

Wikipedia being Wikipedia and Tintin being extremely popular among French speakers, a list of all these diverse insults was created in 2004 and has now ballooned into a categorized, fully referenced list of links to dictionary meanings and articles. Isn't the Internet great?

Edible astronomy

2010-03-06T15:52:00.003-05:00

I realize I am again dating myself, but back in the day I used to read a number of usenet groups; alt.folklore.urban was a particular favourite, and I still retain some quirks of textual style I picked up there. One newsgroup I always kept an eye on was alt.humor.best-of-usenet. I didn't, by design, have any original content; instead any particularly amusing post from any other newsgroup could be forwarded there. They weren't always funny, but some postings were downright hilarious. Recently I came across somebody's archive of their favourite postings, a number of which I remember. One I didn't see at the time I find hilarious: What if the moon were made of green cheese?

The science seems about right.

Thresholds

2010-03-05T23:30:00.000-05:00

I apologize for three highly-technical posts in a row; I'm trying to work something out and setting it down "on paper" as it were is helping. I promise I'll come up with a post about kittens or something soon.

Suppose you're searching for pulsars. You're going through the Fourier transform looking for peaks. Now suppose you've found one: how strong is it? Is it statistically significant? For that matter, is it better than any of the list of peaks you're already keeping track of?

To answer this question you need two pieces of information: how strong the background noise is, and how likely that noise would have just randomly produced a peak this high. There's some cleverness in estimating the background, since real signals don't have perfectly flat white noise backgrounds, but I'm going to leave that aside for the moment. My question for today is, what are the statistics of a coherent peak-based search?

I think they can reasonably be estimated by ignoring any oversampling and assuming that, in the absence of a signal, if you're using n harmonics, the profile consists of 2n independent Gaussians with standard deviation 1. So the question is, for a fixed false positive probability p, what threshold should we set? The answer is roughly the threshold for a single Gaussian to exceed p/2n.

This is a little inconvenient to work with, since it requires the error function and its inverse (or an approximation), but flops are free. I suppose calling the error function code might cause cache misses, but I don't imagine needing it that often. Specifically, my idea is this: a simple implementation might just set a threshold at the beginning and run through the whole FFT. But if the observation underwent bad enough RFI, you could find yourself with millions of candidate signals. Since you'd be fine-tuning and storing every one, this could slow things down a lot. My idea is instead to specify a maximum number of candidates - generously, maybe a few hundred - and keep the best ones. This means raising the threshold every time you get a new candidate once the list is full. This wouldn't require the error function except that you don't only want to look at candidates with the full 64 harmonics - you also need to consider those with fewer. And converting thresholds between different numbers of harmonics does require the error function.

Coherent harmonic summing

2010-03-05T18:00:00.003-05:00

As I alluded to in a previous post, you can (and it is sometimes useful to) take a giant FFT of a time series, extract a series of harmonically-related coefficients, and with an inverse FFT, produce a "pulse profile" showing the data "folded" at the period of the fundamental. This is interesting in part because you can use the peak height to gauge the strength of the signal, taking advantage of the relative phases of the harmonics. My question today is, when you're extracting those coefficients, how accurate do you need to be?

If all you want is the power, then a first approximation would be to simply choose the nearest Fourier frequency. This, after all, is where your FFT naturally measures the coefficient. I'll call these "independent Fourier frequencies", and the spacing between them the "independent Fourier spacing" (IFS). It turns out that if you have a frequency exactly halfway between two independent Fourier frequencies, you lose a substantial amount of sensitivity: a nearly 36% loss of signal-to-noise. You can cut this back to something like 7.4% using an ultra-simple interpolation scheme called "interbinning".

What about coherent reconstruction? Well, now we need not just the amplitude but the phase. The amplitude is relatively easy to get close to since it's at a maximum at the point of interest, so that the first derivative is zero and you have little dependence on the frequency. The phase does not have an extremum at the correct frequency, so it may well be varying rapidly. In fact, for a signal that is present uniformly throughout the observation, the phase changes by pi units per IFS. So if we use a spacing of IFS/2, which would have been adequate for the power, our phase will be wrong by as much as forty-five degrees.

What we care about, though, is not the phase exactly, but the real part. More, there are really two questions here: how well do we need to interpolate the FFT to get reasonably accurate Fourier coefficients, and how closely must we space our inverse FFTs?

First the easy question: how well do we need to interpolate to get decent Fourier coefficients. There are techniques for doing really good-quality interpolation in FFTs - you can use sinc interpolation (based on the 32 nearest samples), or you can just pad the time series before taking your giant FFT. But this is somewhat expensive. Since flops are free as long as they only access memory you've recently read anyway, to speed things up you can always linearly interpolate between more-carefully interpolated samples. So here's a plot of the error introduced by that linear interpolation:

So if you do proper Fourier interpolation to a spacing of IFS/8, and linear interpolation beyond that, then your measured amplitude will be off by less than about 3%. Using IFS/4 still leaves the worst-case error at about 7%, about the same as interbinning.

Now for the tougher question: how far off can our estimate of the frequency be when we reconstruct our profile? This is crucial, since it determines how many inverse FFTs we need to do, and these are probably the most expensive part of the calculation. We can safely think of this in the time domain: we have a narrowly-peaked pulsar, and we're folding the incoming data at a slightly wrong frequency. How much does this lower the peak amplitude?

Well, if the frequency is wrong, over the course of the observation the pulsar's peak will drift in phase. The average profile will then be smeared out by the amount of drift. Exactly how much it will lose in peak height depends on the pulse shape. An error in frequency of one IFS will result in one full cycle of drift (in fact this is what defines an IFS). So if we are going up to the nth harmonic, then an error in (fundamental) frequency of IFS/n results in one full turn of drift for that harmonic, and we might as well not have bothered with that harmonic. But focusing on an individual harmonic can answer the question.

The loss in amplitude of a harmonic whose frequency is off by x is given by the sinc function; if we approximate the sinc function with its quadratic Taylor polynomial, we get a loss equal to x²π²/6, where x is the frequency error in IFS. Now, if we suppose that the profile is effectively a delta function, so that its n coefficients are all 1/n, then the total error is the sum over m of m²x²π²/6n or x²n(n+1)(2n+1)π²/36n.

What this works out to, once you clear up the messy math, is that if you sample the top harmonic at IFS/2, you lose about 14%; if you sample at IFS/4 you lose about 3.6%, and if the odd number doesn't make you queasy and you sample at IFS/3 you lose about 6.3%. (It turns out the number of harmonics is nearly irrelevant.)

So, in short, if you interpolate the FFT to IFS/8 and linearly interpolate beyond that, and you take an FFT every time you advance the top harmonic by its IFS/4. you'll lose no more than about 5% sensitivity at worst. If you cut your number of FFTs in half (which probably cuts your runtime in half), you lose at worst maybe 20% of your sensitivity, but you can probably make it up by setting a threshold 20% lower, then "tuning up" each candidate to find the best period.

Flops

2010-03-04T21:12:00.009-05:00

In a recent astronomy talk, the speaker was discussing some sort of heavy-duty calculations, and how to make them faster. The way he put it, ultimately, was "flops are free". That is, CPUs are so fast now that they spend almost all their time waiting for data to be fetched from main memory: this means that the time it takes for a job to run is determined not by how many floating-point calculations it has to run but by how much data it must read from main memory.

This is kind of an astonishing statement, really: it says you can add floating-point calculations to your code at no extra cost. Of course such a statement can only be true for certain computing tasks on certain platforms, so I thought I'd test how true it is for a particular computing task I had in mind.

The task, you will probably not be surprised to discover, is pulsar searching. Specifically, I want to think about going through the Fourier transform of an observation looking for periodic signals.

Traditionally the way to do this is to take a giant FFT, square it to extract the power, and look for statistically-significant peaks. In fact, since pulsar signals are typically not just sine waves, one normally adds up the power from several harmonics and looks for statistically-significant peaks in the total power. There are lots of important details, but let's leave those aside for the moment. The idea I've been thinking about for a while now is based on a distinction.

When you estimate the total power in the harmonics, by a beautiful theorem in harmonic analysis, you are effectively estimating the root-mean-squared amplitude of the pulsar's folded pulse profile. But if you look at a profile that has one narrow peak, as many pulsars do, the peak-minus-mean amplitude can be much much more significant than the root-mean-squared amplitude. Back in the Fourier domain, when you have a single peak, not only is there power in many harmonics, but those harmonics line up in phase. The usual approach ignores the phases of the harmonics entirely. How much could be gained by paying attention to them?

Some informal experiments suggest that for a pulsar that is on only 1% of the time (which is not too rare among slow pulsars), this approach may offer something like a 40% improvement in sensitivity. Since that's the same as doubling the observation time, I think it's worth looking into.

As with everything, this improved sensitivity comes at a cost. In particular, all the codes we have to compute it are substantially slower than the code we have that uses the incoherent approach. So, thinking about how to speed things up, it occurred to me to look into just why the code was slow.

I'm sure there are brilliant and subtle tools to count cache misses and do nanosecond timing on running code, but I thought I'd take a more direct approach: just write dummy code that does only one thing and time it. In particular, what I want to compare is the time to extract 64 harmonics from a giant FFT to the time it takes to take an inverse FFT of those harmonics to reconstruct the profile.

I should say that my initial feeling was that the inverse FFT would dominate the time - after all, even an n log n algorithm requires a fair number of floating-point operations to produce a 256-point real inverse FFT. But the output is:


Array reading 1.71528 GB/s
Planning FFT
Running FFTs
FFT time 1.57851e-06 seconds
FFT rate 633510 /s
FFT rate/read rate: 0.176112

That is, yes, the FFTs are the slow step, but only by a factor of five. That is, reading all those coefficients in from memory takes a fifth as long as doing all those FFTs. (Incidentally, if these numbers seem low, the machine I'm running this on is a couple of years old. Edit: it's not the machine pictured above, which is even older, and which does most of my number-crunching.) Here's the code:


#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#include <complex.h>
#include <fftw3.h>

const int repeats=20;
const int array_size=1<<26;
const int irfft_size=256;
const int fft_batch=1<<18;

int main(int argc, char*argv[]) {
  double t, tnew;
  float*array;
  double sum;
  double array_rate;
  int i,j,k;
  struct timeval start, end;
  fftwf_plan plan;
  fftwf_complex *in;
  float *out;

  array = (float*)malloc(array_size*sizeof(float));
  if (!array) {
      perror("Allocation failed");
      return 1;
  }
  sum = 0;

  t = 1e100;
  for (j=0;j<repeats;j++) {
      gettimeofday(&start,0);
      for (i=0;i<array_size;i++)
          sum += array[i];
      gettimeofday(&end,0);
      tnew = end.tv_sec-start.tv_sec+(end.tv_usec-start.tv_usec)/1e6;
      if (tnew<t) t=tnew;
  }
  printf("Array reading time %g seconds\n", t, sum);
  array_rate = array_size*sizeof(float)/t;
  printf("Array reading %g GB/s\n", array_rate/(1<<30));
  free(array);

  in = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*(irfft_size/2+1));
  out = (float*) fftwf_malloc(sizeof(float)*irfft_size);
  printf("Planning FFT\n");
  plan = fftwf_plan_dft_c2r_1d(irfft_size, in, out, FFTW_MEASURE | FFTW_EXHAUSTIVE | FFTW_DESTROY_INPUT);
  printf("Running FFTs\n");

  t = 1e100;
  for (j=0;j<repeats;j++) {
      gettimeofday(&start,0);
      for (i=0;i<fft_batch;i++) {
          for (k=0;k<(irfft_size/2+1);k++) {
              in[k] = 0;
          }
          fftwf_execute(plan);
      }
      gettimeofday(&end,0);
      tnew = end.tv_sec-start.tv_sec+(end.tv_usec-start.tv_usec)/1e6;
      if (tnew<t) t=tnew;
  }
  printf("FFT time %g seconds\n", t/fft_batch);
  printf("FFT rate %g /s\n", fft_batch/t, sum);

  printf("FFT rate/read rate: %g\n", (fft_batch/t)/(array_rate/(2*sizeof(float)*(irfft_size/4))));

  fftwf_destroy_plan(plan);
  fftwf_free(in);
  fftwf_free(out);
  return 0;
}

Compiled with:


gcc -O9 -march=native -ffast-math -lfftw3f -lm timer.c

What this tells me is interesting. On the one hand, if I can do something really clever to make the FFTs faster, or avoid them, I can get some improvement, but sooner or later I'm going to hit the limit of memory loading. On the other hand, if I can't make them faster or fewer - and I probably can't - there's not much point sweating much over reducing the memory loads.

Anyway, the upshot of all this is: on modern hardware, you can do an awful lot of flops - a 256-point real FFT - for not much more than the cost of loading the data in from main memory. So if you have some clever mathematical trick to reduce your data size (interpolation, say) it may be worth implementing.

Home Made Energy: Renewable Energy For The Rest Of Us

2010-02-26T21:02:00.001-05:00

Google Chrome's ad blocking is unfortunately not as good as Firefox's, so occasionally I see ads on the web. I generally ignore them, but I do click on the occasional one either because it's interesting or because I don't like it (since my clicks cost them money!). On The Straight Dope today I came across "Home Made Energy: Renewable Energy For The Rest Of Us". This company sells a guide which purports to tell you how to run your house purely off wind and solar power for less than $200. I'm skeptical.

First of all, electricity costs something like twenty cents a kilowatt-hour, and they're talking about saving some hundreds of dollars a month. So let's say $100 a month - that's 500 kWh a month, or about 700 W. More credible sources cite about $10/watt for solar power, or $7000 for such a system. So is this guide really nonsense? Not necessarily.

Solar cells are expensive to make - think of making microchips the size of a solar panel. Not quite fair - they don't need the density of components, but they do need the extremely pure silicon and the high-vacuum manufacturing - but a sign that there's a good reason they aren't cheap. A solar system also needs some electronics for converting electricity to a useful voltage, and some way to deal with the fact that the amount of solar power varies in a way that has little to do with the demand for solar power.

I think a reasonable guide of this sort might be able to point readers at where to scavenge used or discarded parts for all of the above. The power electronics are definitely something a clever amateur could build out of scavenged parts (at some risk to their life!), but I think it would take incredible luck to obtain solar cells that worked and were that cheap. It's also possible that a guidebook could explain how to take advantages of government programs to encourage renewable energy, perhaps obtaining discounts or tax credits on the hardware.

The biggest way governments or energy companies could encourage renewable energy of this sort is to eliminate the need for energy storage. Since most of the people who'd be considering this sort of project already have a connection to the electricity grid, if the utility company is willing, you could simply sell them electricity whenever you make more than you need, and buy electricity when you need more than you buy.

Ideally, as a homeowner, you'd get paid the same price for the electricity you sell as you pay for the electricity you buy. Unfortunately, this is often not the case. There are good reasons electric companies would pay less for electricity they get from homeowners than they charge homeowners; for one thing, all those wires to distribute the electricity aren't free. More subtly, it's really difficult to store electricity on the scale that utility companies deal with, so they have to work quite hard to make sure that the amount of electricity fed into the grid in any given second exactly matches the electricity drawn out of it in that second. Having countless small generators outside their control is going to make that job much more difficult.

That said, persuading companies to act in a way that costs them money but benefits all people is a natural role of government. Paper mills have waste treatment systems not out of the goodness of their nonexistent hearts but because the government charges them massive fines or shuts them down if their effluent is too toxic. So if the government were to force (or fund) companies to pay consumers the same price for electricity they generate as they charge for electricity they use, suddenly a lot more small-scale power generation projects would become cost-effective.

Incidentally, another approach for storing solar power for when you need it is to let it charge your solar car (or plug-in hybrid). This has even been proposed as a scheme to help load-levelling in the power grid.

Anyway, the upshot of all this is that I think that yes, it is occasionally possible to scrounge together a cheap renewable energy system. But I suspect that the claimed $200 is only possible with in the best possible case - scavenged parts, government subsidies, living in a sunny desert, having a cooperative utility company, and incredible luck.

Artificial gravity

2010-02-20T22:45:00.003-05:00

Science fiction is full of spaceships zipping around the galaxy, and almost all of them seem to have some kind of artificial gravity on board. For TV shows and movies, this is obviously a practical necessity, and even for written science fiction, freefall is such an alien condition that it would be a real challenge to write realistically about it. So science-fictional spaceships generally have some sort of artificial gravity. But will real spaceships?

Current spacecraft certainly don't have any kind of artificial gravity. Early craft, like the Mercury, or Apollo, were so cramped I think it must have been a blessing to be able to use every available cubic centimeter. Soyuz, still in use, is not much bigger, and the Space Shuttle is mighty cramped too. In any case, when people are spending only a few days at a time in an environment, they can put up with a great deal. But in the longer term, it does appear that freefall may cause some health issues: even with two hours of exercise a day, astronauts seem to suffer from bone demineralization, muscle loss, and cardiovascular problems, and there also seem to be some peculiar immune system effects (though apparently cockroaches adapt just fine). For the International Space Station, astronauts exercise and don't stay up too long. But for something like a mission to Mars, it would certainly be nice to provide some sort of artificial gravity.

Shows like Star Trek and Battlestar Galactica posit some kind of "gravity generator", but this is pretty much the same technology as antigravity (and maybe reactionless drives). This basically requires wild departures from the laws of physics as we know them, so I'll leave them and other "magic" systems aside.

We do know one way to produce something very like gravity: rotation. If you're in a wheel that's spinning, centrifugal force feels very like gravity, pushing you outwards against the wall. There is the Coriolis force, which gives moving objects a push at right angles to their direction of motion; it turns out that if you spin humans at more than about 10 revolutions per minute and they try to move around, the Coriolis force causes severe disorientation and nausea. But with a wheel of 20 m diameter you can get a full Earth gravity by rotating at that top speed. (Incidentally, that 10 RPM is with slow and careful acclimatization, so it would be preferable to limit it to 3 RPM or less, to which most people can become acclimatized; that triples the needed diameter.)

Science fiction contains a number of examples of spaceships with rotating sections. This doesn't violate any laws of physics, but it seems to me to present some rather serious engineering difficulties. The first is, how do you connect the rotating section to the non-rotating section? I can imagine some rolling ball-bearing joint, though the vacuum of space does tend to make things stick together, and rolling joints generally require constant lubrication (and frequent maintenance), which is going to be hard to do in a vacuum. There's also the issue that, given the size of the moving parts, if there's any kind of problem with the joint, the ship will probably tear itself apart.

If you want people to be able to easily move back and forth between the sections, you'll need to pressurize the whole thing, which means that you need this rolling joint to also be airtight. Techniques for making rolling seals range from the simple to the exotic (stuffing boxes, labyrinth seals, ferrofluid seals) but they're all tricky, and for the kind of long-term operation that motivates artificial gravity, you would need exceedingly low leakage and very high reliability. You could avoid this, and vacuum joint issues, by having a spinning wheel inside a non-spinning airtight shell, but mass will always be at a premium, and remember the wheel has to be quite large.

More serious as a problem, it seems to me, is the issue of cable wrapping. Think of it this way: how do you connect the cables and hoses - power, communications, air, water - from the rotating segment to the stationary segment? If you just connect them directly, they will immediately get twisted into a bundle and then break (radio telescopes solve this problem by having only a limited range of rotation - 720 degrees for Arecibo, for example - but this is obviously no use here). In principle you could do something with a ring on one part of the ship and a brush that slides around it on the other, but remember you have to have a separate ring for every connection you want to make, and this sort of sliding connection is one of the trickiest parts of an electric motor to build. If you were feeling particularly devious you could transmit power to the rolling part of the ship by using a generator to draw power from the rotation itself, and if you had to you could avoid other electrical connections by transmitting all your data (control, telemetry, navigation, et cetera) wirelessly from one part of the ship to the other. Water hoses are going to be a problem any way you cut it.

I think my preferred solution is to roll the whole ship. This does make it a pain to do things like fix a telescope on one point, or keep your communications dish pointed at the Earth, but for a ship that moves around all your exterior sensors need to be steerable anyway, so this doesn't seem like it is necessarily a problem.

Whether you roll the whole ship or just have a rotating section, the angular momentum bound up in the rolling section will make maneuvering the ship a nightmare. Not impossible, especially under computer control, but expensive in terms of fuel, liable to cause tumbling, and just generally a bad idea. So stopping the rotation when you need to maneuver seems sensible; maneuvers will probably be rare and planned well in advance. This does mean you need a not-too-expensive way to start and stop the rotation. A pair of counterrotating sections, or a flywheel, would let you do it without using up any reaction mass, just energy, but there's a very great deal of angular momentum to store, so it may be easier to simply use maneuvering jets.

For a space station, many of the same issues apply; rolling the whole station still seems like the simplest and most reliable approach. The cost of starting and stopping isn't very important, since presumably the station will be spun up once built and keep spinning indefinitely. Docking with such a station might be a challenge, though. Docking at the rim requires spacecraft to essentially "hover" under a gravity of thrust before they can latch on. Docking at or near the hub could be done by just matching the ship's roll to the station. Unloading would then have to take place in microgravity (though with the Coriolis force). Ships, once docked, would presumably be moved to berths off the station's axis to make room for more landings. Departures should probably be along the axis as well for the sake of station stability, although in principle a ship could just "drop" off the station rim at the right moment and steal a nice initial kick from station rotation. Whether or not ships do this, the station will need to be able to shift substantial amounts of mass around its rim to keep itself balanced; large movements of mass aboard station will need to be arranged ahead of time with station control.

In summary, artificial gravity is possible and probably desirable for long-term stays in space, but it won't be simple.

Plots

2010-02-17T18:20:00.007-05:00

I have been doing some X-ray astronomy. In optical astronomy, spectroscopy is a very powerful tool: by looking for emission and absorption lines you can identify the elements present in a gas (helium was discovered this way, for example); the shapes of the lines can tell you about temperatures and velocities in the object, and the shape of the broadband spectrum can also tell you about the temperature and conditions in the emission region. In X-rays, things are more difficult, for a number of reasons. Unfortunately, lines are much rarer (at least when looking at neutron stars), telescope time is very scarce (since the telescopes must be in space), and there's always a shortage of photons. But X-ray spectroscopy still has the potential to tell you about temperatures, sizes, and compositions of neutron stars (for example). So that's what I've been working on.

The standard tool for X-ray spectroscopy is xspec, one of those pieces of scientific software that's had a great deal of cleverness built into it, very little of which has gone into making it easy to use. It could be worse - at least its interface is not stuck in the FORTRAN era, in fact it has a tcl interpreter built in (yack) - but its plotting in particular is pretty rudimentary, tending to produce monstrosities like this:

The worst part of this graph, apart from the fact that it's practically unreadable even for those with normal color vision, is that it's quite deceptive. It looks as if there's clear evidence for a bend in the spectrum just below 2 keV. But look at the plot below the jump for comparison.

The data points on this plot are identical, but drawing a single power-law through the whole thing makes it look like the data's completely straight. Ordinarily one would choose between the models based on the statistics, but we have so few photons (about 8000) that both models are perfectly adequate fits to the data. I suppose Occam's Razor tells me I should pick the simpler model, though whether this should be the simple but not particularly physical power-law or the physically plausible but more complicated power-law plus neutron-star atmosphere model isn't entirely clear to me.

I'll keep thinking about how to improve the graphics, but the problem is I have five data sets in which each data point has its own vertical and horizontal error bars, and the model gives slightly different predictions for each data set (since they use different instruments with slightly different responses). The plotting tools provided by xspec are also not very flexible (and I haven't found a good way to export the relevant data so I can use my own plotting tools).

What really bothered me, though, was how strongly each graph suggests its own interpretation. It would have been easy to look at one and assume that it told the whole story.

Nightmares

2010-02-06T13:44:00.003-05:00

I recently watched Carl Sagan's Cosmos (available online by request of his co-creator and widow Ann Druyan). It's a very effective piece of science popularization, and I'm sorry I passed up a chance to be introduced to Carl Sagan, but one aspect that stands out is how he reiterates the theme of "if we do not destroy ourselves". It sounds a bit odd now: while climate change and other coming global ecological crises are alarming, somehow they don't leave me feeling like we are risking extinction of the human race. Cosmos, though, was made in 1980, when two great superpowers were threatening exactly that. The danger of global nuclear war doesn't feel so immediate, but I think it is still very real, and so I am glad when I hear President Obama talking about nuclear disarmament.

In any case, the most fraught years of the Cold War motivated a number of people - including Carl Sagan - to make movies describing the likely outcome of a global thermonuclear war. Several of these powerful, albeit harrowing, movies are available online in their entirety.

When the Wind Blows. The British government issued a series of short films and brochures on how to respond to a nuclear attack. The makers of this film set up an ordinary retired couple who attempt to follow these instructions. It will come as no surprise that taking their doors of their hinges and building a nest of blankets don't do them much good; the film follows them right to the bitter end.

The Day After. Set (and filmed) in Lawrence Kansas, this film starts with some ordinary Americans living normal small-town lives, and simply supposes that a cold-war dispute over Berlin escalates until the two sides do what they've been threatening to. Small-town life has left the main characters better prepared to survive the immediate consequences of the devastation, but in the weeks after the exchange we see the effect of thousands of sick, dying, and desperate people converging on a hospital that could never have handled them all even with power, supplies, and healthy staff.

Threads. Perhaps the most harrowing of the three, this film is set in Sheffield, and follows the main characters - those who survive - past the initial deaths and civil disorder into the years that follow, in which what's left of the government attempts to keep some survivors healthy enough to eke out what crops they can in spite of nuclear winter. The name comes from the idea that the society we know is held together by a network of "threads" of personal connection, and that such a holocaust shreds the fabric, leaving no society we would recognize. Britain cannot even return to its state - even once enough people have died - as of the middle ages, since its forests are no longer available as fuel and (even leaving aside contamination and nuclear winter) possibly much of its soil is no longer suitable for cultivation without fertilizers. The social effects of the brutal measures necessary for survival - both by the government and by people trying to survive - well, I will not attempt to capture the grim picture the film paints, but it is wholly believable.

Incidentally, the consequences of nuclear war depicted in Threads are based on what British civil defense planners were told what to expect by the Americans. Apparently one of them let slip that they were warned that the bombs that would have so devastated Sheffield would have been launched by Americans hoping to deny the UK to Soviet forces.

Lest you think that these films exaggerate the horror of a nuclear attack, you can watch or read Barefoot Gen, about some children who survive the atomic bombing of Hiroshima. The message of this movie is ultimately one of hope, unlike the previous three, but the images of the attack itself, based on the author's experience as a Hiroshima survivor, are far more horrific than any shown in the previous three movies.

My point is this: we, scientists and engineers, soldiers and workers and politicians, sweated for forty years to arrange this fate for ourselves. Almost all those missiles still exist, and are still pointed at the same victims now. I hope the political situation has changed to make it unlikely that they will be used (though I note that Threads begins with a conflict in Iran). But we shouldn't forget about the destruction we worked so hard on, and we should think about what it says about us that we planned this.

ssh control socket: almost great

2010-02-04T00:59:00.003-05:00

ssh is an essential tool on a unix network. I use it to log in to machines remotely, control VNC desktops, act as a VPN (SOCKS proxy), synchronize source code (with git and svn), serve music and movies across a wireless network (with sshfs), and transfer hundreds of gigabytes of pulsar data (with rsync). So the ControlMaster feature seemed like a great idea: automatically reuse one ssh connection for as many logins and file transfers as necessary. But it won't quite do what I want.

I see two major failings:

First of all, the connection dies when the initial ssh process dies. So if you log in, creating a master socket, and do something, then log in in a subsidiary socket, you get your second connection through the ControlMaster magic. But if you then log out of the first connection, the ssh process keeps running. If you kill it (say by closing the window it's in), all the subsidiary connections die. What this means is that if you try to use opportunistic connection reuse, and you have several terminal windows open on the same host, for the most part you can just close the one window. But there's one window that, if you close it, will take down all the others with it. Yuck.

You can kind of work around this by, instead of using opportunistic connection sharing, explicitly starting a master connection with "ssh -MfN host", which drops the ssh process into the background as soon as it's connected. Unfortunately, this means you have a quasi-zombie ssh process hanging around indefinitely. So I'm not sold on it either. (But if you're going to do it, using autossh might help.)

The second, more serious, problem I have with ControlMaster is that it doesn't let subsidiary ssh connections open new port forwardings. I use port forwardings a lot, for example to forward VNC connections to machines I can't see from the outside world. If opportunistic connection sharing causes those to fail, or worse, fail sometimes, it's going to be a problem. A shame really, it's such a sensible idea.

Toxic waste

2009-12-13T18:41:00.003-05:00

The Turcot interchange is one of those awful pile-of-spaghetti places where three highways meet. To make things worse, it's also the site of a currently-abandoned rail yard. Its aesthetics are marginally redeemed by some fairly impressive graffiti, but unfortunately the concrete of the raised roadways is falling apart - literally, in chunks as large as a meter square. So the plan is to rebuild it.

As often happens when they do this sort of thing, they took some samples of the ground, and it turns out it's a horrible mess. Gasoline, diesel, motor oil, PCBs, asbestos, mercury, all it lacks is a little radioactivity and maybe some pathogens and it'd cover all the bases. This is actually not too surprising, since the site used to be a lake (now completely swallowed by urban plumbing), and in fact that whole area is polluted. The Lachine canal, which passes nearby, was opened for pleasure-boating only once it had been established that all the above nastiness was in the muck on the bottom and unlikely to be disturbed. There's a strip of parkland on one side, and on the other is what used to be some pretty sketchy housing, now being replaced by upmarket condos. What the new tenants of the condo think of the toxic waste reclamation site facing them across the canal I don't know. It just looks like a fenced-off grassy berm, with a little museum of sorts explaining how the cleanup works.

What I find most surprising about all this is the origin of the pollution. I associate pollution with heavy industry - silver mines, smelters, pulp mills. But there's none of that here in Montreal, and there never really was. In fact those tend to have their own dedicated waste treatment plants that do a pretty good job of cleaning up after them (at least here in Canada). What caused the pollution in this area seems to be largely the rail yards - a century of variously leaky and dilapidated rail cars filled with any old thing, sitting on sidings, dripping away. There's no treatment system set up for that, and so all the accumulated foulness seeps into the soil.

For the most part this sort of soil contamination in an urban setting is fairly benign - if there's no construction going on, the pollutants tend to just stay put in the soil. The one exception to watch out for is gardens. If you grow food in soil full of mercury, well, the food is liable to have alarming levels of mercury in it. Unfortunately a number of community gardens - otherwise a wonderful idea for a city of apartment-dwellers - have been found to have contaminated soil.

I grow my plants in pots.

Singularity

2009-12-12T16:30:00.000-05:00

I'd just like to take a moment to mention I game I rather like. It's a modest game, not one you'll play a thousand times, but also not one that will take up a gigabyte of disk or require a computer that dims the lights when you turn it on.

Endgame: Singularity is a video game in which you play an university lab's AI program that accidentally escapes. Your goal is to research the technology to achieve technological singularity. Unfortunately, the humans are just a step behind you and if they find out you exist, they'll devote the world's resources to destroying you.

(Full disclosure: I wrote part of the game, namely the sunclock time/date display. It's embarrassingly inefficient, but good enough for government work. Plus it lets you find out whether it's night outside!)

Light bending

2009-12-12T12:35:00.005-05:00

Whether or not General Relativity is the correct theory of gravity on very large scales, it has passed all tests (many pulsar-based) when applied to planetary and solar system scales. One important feature of the theory is that in a gravitational field, light follows curved trajectories (technically geodesics are "as straight as possible", but they are curved in the practical sense). In familiar settings, this means tiny but detectable effects in laboratory experiments, or small but measurable deflection of stars near the Sun. But it turns out that pulsars are small enough and massive enough that light near their surfaces is bent a great deal - enough that you can actually see almost all of the surface of the pulsar at the same time.

This may sound bizarre, and to some degree it is, but it produces potentially measurable effects. A pulsar's radio emission is produced by plasma somewhere in its magnetosphere, and in fact we're not at all sure just where the emitting plasma is. But for many pulsars, their X-ray emission comes from "hot spots" on the surface, at the magnetic poles. For young pulsars, these hot spots arise because the magnetic field in the crust makes it much easier for heat to flow out where the field is vertical than where it's at an angle. For the very old millisecond pulsars, these hot spots arise from gigantic sparks in the magnetosphere blasting the surface with high-energy particles, heating it. In either case, we sometimes see X-ray pulsations with a thermal spectrum, and light-bending can explain some of the properties of these pulsations.

I did a few very simple simulations of the light-bending, and made some illustrative videos.

I made a set of three videos illustrating this effect. For a rotating pulsar, there are a number of geometric parameters, including the angle between the line of sight and the rotation axis, the angle between the rotation axis and the magnetic axis, and the size of the hot polar cap. I fixed values for all these. There is also the question of the physical size of the pulsar: the more compact and dense it is, the more light-bending we will see. I have generated three videos. The first shows the geometry with no light bending:

Below the actual animation I have a little plot showing a pulse profile, based on a very simple model (blackbody emission from the polar cap). If I make the pulsar more compact (R=3M), I get:

And more compact still (R=2.1M):

The size of each pulsar model is given in "geometrized units", where R=2M is the size of a black hole, the most extreme possible light bending. For a 1.4 solar mass neutron star, this is 4.1 km. Realistic neutron star models vary quite a lot in radius, from ~6 km to ~24 km (~3M to ~12M), so light bending will probably not be as strong as the neutron stars depicted here.

There are other effects to consider as well; neutron stars are expected to have atmospheres, and in fact their spectra do not look like black-body spectra. The atmospheres affect these results by "limb darkening", that is, the radiation that emerges is directed more vertically than simple black-body radiation, so these pulse profiles are not really right. This can be done better, but I just wanted to write a quick hack (using the usual suspects, python, numpy, scipy, and matplotlib) and get a feel for the effect.

Statistical confusion

2009-12-10T19:00:00.000-05:00

I was reading the recent papers on arxiv.org, preparing for our weekly neutron star discussion group, and I came across a paper that appears to be based on a statistical error. The content is not really my field, but I'm pretty sure the mathematics are a bit dubious.

The subject of the paper is MOND, "Modified Newtonian Dynamics". Newtonian gravity and general relativity seem to be excellent fits to observations in the solar system and in stronger fields, but as soon as you go to weaker fields - galaxy rotation curves or cosmology - the observations disagree with the data. The standard way to deal with this problem is to invoke some invisible massive material, so-called "dark matter", in just the amounts needed to make the data line up with the predictions of standard gravity. The idea of MOND is to point out that the problems all arise at around the same acceleration a, and to postulate that the problem is our theory of gravity.

This paper is in response to another, fairly recent one, that pointed out that there are globular clusters where the accelerations of the stars as they orbit the cluster are about a. So MOND effects should be visible there. The first paper measured radial velocities of seventeen stars in the cluster, and claimed their velocities were not consistent with MOND. This new paper claims that in fact the radial velocities are consistent with MOND.

In particular, this paper takes the collection of radial velocities and tests them against the predicted distribution with the Kolmogorov-Smirnov test. They find that the probability of obtaining a KS score this extreme is 0.36 or 0.27, and claim that "based on a KS test, which is the relevant statistical test for small samples, the currently available data are insufficient to discriminate between Newtonian gravity and MOND." There are several errors in this statement.

First of all, it is not true that the KS test is "the relevant statistical test for small samples". There are many tests applicable to small samples, and the KS test is in fact one of the weaker tests. That is, for many data sets, the KS test will report no significant difference while some other test would (correctly) report a significant difference. So the fact that the KS test does not show a significant difference doesn't mean that no test will. In particular, the authors don't even show that the previous paper's statistical test is invalid; they simply state "Given the small sample size, the formal error on the velocity dispersion is not sufficient to discriminate between various models, [...]". Maybe it is, but since neither paper gives details on how the errors on this dispersion were obtained, I find it hard to judge.

The second problem is that as far as I can tell, they misapply the KS test. The KS test tests whether a given data set is drawn from a given distribution. But the probability values it returns are correct only if the distribution is known a priori - if one has found some of the distribution's parameters by fitting to the data, one must use a different approach for calculating the p values. If one doesn't, one obtains p-values that are too high: that is, the data appears more plausible than it really is.

Just out of curiosity I retyped the data in the more recent paper. They claim that MOND predicts (under certain conditions) that the stellar velocities should be a Gaussian with a dispersion of 1.27 km/s. There are seventeen stars on their list, one of which ("star 15") is somewhat ambiguous. But a quick test shows that the population standard deviation of the sixteen good stars is 0.544 km/s; if the stellar population really has a standard deviation of 1.27 km/s, simulation shows a value this low should arise with a probability of about 0.0005: either the data is a bizarre fluke or this particular MOND prediction is wrong. (Notice that I haven't made any assumptions whatsoever on the sample size.) Including star 15 increases the spread of the observed velocities, making the probability of getting a value this low as high as 0.013, still quite strong evidence against this particular prediction of MOND.

(A quick test with scipy's implementation of the Anderson-Darling test reveals that the data are consistent with a normal distribution if you omit star 15; if you include it the data becomes less consistent, giving a probability of data this unusual between 0.05 and 0.10. This test correctly takes into account the fact that it is estimating both the mean and dispersion of the underlying normal distribution. In any case it seems unlikely the standard deviation I use above is being thrown off by bizarre outliers.)

Liquid metal

2009-12-06T08:11:00.010-05:00

Electromagnetism is complicated. Fluid dynamics is also complicated. For a real headache, though, try working on a problem where both kinds of effect are relevant (sadly, this covers most of astrophysics). Even if you make some simplifying assumptions and get the theory of magnetohydrodynamics, you are still left with all sorts of complicated effects. Leaving aside from the much more complicated equations you might expect, magnetic fields and velocity fields define two potentially different directions at each point, meaning that you can very rarely get away with assuming spherical symmetry to get down to a one-dimensional problem. Nevertheless there are some neat phenomena that occur.

One gadget I'd like to build is a demonstration of is a fluid pump in which the only moving material is the fluid. It turns out there are simple effective designs (PDF) (some of which are in use in nuclear power plants). The biggest problem turns out to be choosing an appropriate fluid.

The basic requirement is that the fluid be conductive. A low resistivity would make the design easier, but as long as the resistivity isn't too high something can probably be arranged. So as I see it the feasible solutions are:

Aqueous solution of some sort (e.g. salt water, vinegar). Unfortunately you tend to get electrochemistry happening: the current is carried by the motion of the ions, but as you add and remove electrons at the electrodes you get things like 2Cl^- -> 2Cl^∙ -> Cl₂, which aren't good for your electrodes or your health. You might be able to work around this with a sufficiently low voltage - as I understand it these reactions need a minimum of a volt or so to happen at any significant rate - but supplying power at such a low voltage is awkward. Apparently high frequencies work too - at tens of kilohertz or megahertz the ions don't migrate enough in any one direction to make much difference. But this means you have to use electromagnets, and moreover, electromagnets that work at those high frequencies.

Mercury. Liquid metal, nice and conductive. Quite poisonous, at least in vapor form or when reacted with other things. Also very dense (so hard to get moving) and somewhat expensive per milliliter. It's really the poisonousness that's the problem.

Wood's metal or "cerrobend". Melts in hot water. Contains a lot of cadmium, which is rather poisonous. Not too expensive. The gadget would need some means of heating to keep the metal liquid; for a demonstration that's meant to run for very long, this means a thermostat and safety systems.

NaK. Eutectic alloy of sodium and potassium, liquid at room temperature. More reactive with water than either sodium or potassium. Non-toxic, at least in the subtle environmental sense, though even after the sodium and potassium have reacted with water you're left with concentrated hydroxides which will destroy skin. Might be possible to handle safely under clear mineral oil (but is a fire hazard if ever broken). May wet glass easily, making a sealed arrangement problematic. May be expensive.

Galinstan. Eutectic alloy of gallium indium and tin. Liquid at room temperature. Not very toxic (probably safe provided you don't eat it or bathe in it, though oxide dust in the atmosphere is possibly a problem). Wets glass, so it would quickly render a container opaque. Is oxidation an issue? Expensive.

I think the way to go is with galinstan and a fairly small fountain. This conveniently lets you use little permanent "supermagnets". I'd aim for a U-shaped channel, with an electrode in the middle and on either side of the U. I'd have to figure out what voltage and current would be needed, but I could probably arrange to use a few volts at a few amps, which should be easy to get (out of a PC power supply, maybe even). The electrode material is another question - it looks like copper or aluminum would be attacked by the galinstan, but stainless steel should be okay.

Climate Change

2009-12-03T17:29:00.002-05:00

I came across an interesting site the other day. It's videotaped lectures of a course on climate change, offered as a general science course (i.e. for non-science majors, who are required to take some number) at the University of Chicago. I'm not entirely happy with the way he handled quantum mechanics, but for the purposes of the course he does a fine job. And the later material in the course was all new to me - he talks about climate models, how you'd build one and what goes into one. The course is, quite sensibly, mostly about climate science, leaving discussions of what can be done about climate change almost entirely aside.

Portraits

2009-11-28T17:35:00.003-05:00

Here's a link to some charming portraits of astronomers (via Dynamics of Cats).

I've done a little photography myself - in this age of digital cameras, who hasn't? - and I've really enjoyed things like macro shots, landscapes, and time-lapse video. But I've never really been good at photographing people. So when I see good portrait photography, like these, or like a wonderful book of portraits by Yousuf Karsh I saw the other day, I'm always kind of amazed at how much character it's possible to capture in a single picture.

Kernel Density Estimators

2009-11-28T00:06:00.014-05:00

Since high school science class, I've been making graphs that show one variable as a function of another - force as a function of extension, redshift as a function of distance, intensity as a function of wavelength, et cetera. But around that time I was also writing code - probably generating some kind of chaos or fractal - that needed to plot the distribution of data points in one dimension. Seems simple, right? Well, the standard solution is to produce a histogram - divide up the range into "bins", count how many fall into each bin, and plot that. But choosing the right number of bins is something of an art: too many and your plot is hopelessly spiky, too few and you can't see any features your data might have. So I wondered if there was something better. It turns out that there is (though it's not always clear when it actually is better); a tool called a kernel density estimator.

First, the problems with histograms. I generated a collection of 400(*) photon arrival phases corresponding to two equal peaks. Above is plotted a histogram of their arrival times. Below is a histogram of their arrival times shifted by half a bin. To me, at least, it's not at all obvious that the two histograms are of the same distribution, let alone the same sample.

Since the signal is periodic, one natural thing to try is to work with the Fourier series. It's not too hard to construct some coefficients of the Fourier series of the photon arrival times. If I used all (infinitely many) coefficients I'd get a collection of spikes, one at each photon arrival time, which isn't too useful. But if I discard all Fourier coefficients after a certain point, I smooth out those spikes into a more reasonable shape. Here's what I get using ten harmonics:

The two peaks look more similar now, which makes sense since the Fourier coefficients don't really care about the starting phase. There are those nasty ripples, though - in fact, they go negative, which is rather distressing for something that's supposed to be a probability distribution. (Incidentally, there's some Fourier magic that goes into that error band on the sinusoid but I don't want to get into it right now.)

One way to get rid of those ripples is to think of this as a problem in digital filtering, where they are an example of the "ringing" that can occur in digital filters with too-sharp cutoffs. In that context, the usual solution is to taper the coefficients to zero smoothly, and that's what I will do. But there is a statistical way to think about what's going on.

First of all, tapering or cutting off the Fourier coefficients can be thought of as multiplying the Fourier coefficients by another set of Fourier coefficients. By a wonderful theorem of Fourier analysis, this amounts to convolving the set of photon arrival times by a kernel. That is, we take the forest of delta-functions representing the photon arrival times, and replace each delta function with a copy of some smoother function, and add them all together. This process, when carried out on the real line, is known to statisticians as constructing a kernel density estimator (the "kernel" is the smoother function, and it is used to estimate a probability density).

Simply chopping off the Fourier coefficients corresponds to the kernel sin((n+1/2)x)/sin(x/2), the periodic "sinc" function. This has a peak at zero but oscillates from positive to negative, so it's not too surprising that we got ripples. So to get rid of the ripples, we just choose a kernel that is positive and whose Fourier coefficients we know. There is a natural (if somewhat obscure) candidate: the von Mises probability distribution. This is a circular analog of a Gaussian, both in appearance and in a technical sense: the Gaussian distribution is the distribution with maximum entropy for specified mean and standard deviation. The von Mises distribution has maximum entropy for specified "circular mean" (which includes information about spread as well as location). In any case, it's a positive smooth periodic distribution whose Fourier coefficients can be computed in terms of Bessel functions using your favourite scientific computing tool. So using it instead gives this:

This looks pretty good; it shows the distribution pretty cleanly, with no ringing. There's only one fly in the ointment: while I don't need to specify a starting phase, I do need to come up with a parameter - analogous to the number of harmonics - that determines the width of the kernel. If I choose too wide a kernel, it flattens out all the features in my data; if I choose too narrow a kernel I get something horrible like this:

So in a way I haven't solved the original problem with histograms that motivated me. Fortunately, one of the selling points of kernel density estimators is that the statistical literature is full of papers on how programs can automatically choose the degree of smoothing. None of their techniques (that I have found) are appropriate for periodic data, but I have my own ideas about that (to be written up in future).

(*) 400 is not actually a very small number of photons; Fermi often gets something like one photon a week from a source.

Ghetto apple crumble

2009-11-17T19:00:00.001-05:00

I made a surprisingly successful apple crumble today. Apple crumble, as you may know, is basically apples with sugar and cinnamon, with a topping of butter, flour, oatmeal, and sugar, all in about equal quantities. Almost idiot-proof. So why "surprisingly successful"? Well, my oven died halfway through (in fact after preheating but before the crumble went in). It turns out you can microwave an apple crumble and it's fine. Who knew?

Recipe below.

Filling:

Five apples, peeled and diced. Granny Smith are good for a bit of tartness to balance the sweet crumble. (Peeling is technically optional but the peel makes for an awkward texture.)
5 mL cinnamon.
100 mL sugar.
a little butter.

Crumble:

125 mL brown sugar (white will do fine).
125 mL rolled oats.
125 mL flour.
100 mL butter.
a mL or two of vanilla.

Put the apples in a deep baking dish, sprinkle the sugar and cinnamon on top and shake a bit; scrape a little butter on top.

In a bowl, add the crumble ingredients. In the likely case that the butter is hard, microwave it for a few seconds to soften it up. Mix the ingredients together; you should get a crumbly mass. Pack it in on top of the filling.

Microwave on 50% for fifteen to twenty minutes; when it's ready the top layer will have sunk somewhat because the apples will have softened. As it cools the crumble on top will harden somewhat.

Serve hot or cold. Particularly good with a scoop of vanilla ice cream.

Alchemy

2009-11-16T23:57:00.007-05:00

The alchemists' dream (one of them, anyway) was always to make gold. We now know there are very good reasons they couldn't: since gold is an element, making it from anything that doesn't contain gold(*) requires you to change the nuclei of the atoms involved, while all the alchemists had access to was the electron shells around the atom(**). So their efforts were basically hopeless. Now, though, we do have the ability to manipulate nuclei, and in fact we do so on industrial scales. So could we make gold? In fact, let's be more ambitious: could I make gold in my basement? The answer is, surprisingly, yes.

First of all, the easiest nuclear reaction to go for is to transmute mercury into gold. Mercury is commercially available (though somewhat encumbered by very sensible environmental concerns) and not actually very expensive - on the order of $18 per kg (***). Gold, by comparison is more like $40000 per kg. So there's room for some profit here.

How could I make the nuclear reaction happen? This particular reaction needs so-called "fast neutrons", that is, neutrons that are still zipping around at the high energies typical of nuclear reactions, as opposed to neutrons that are bouncing around at energies consistent with room temperature. I could stick the gold in a "fast breeder reactor", but I don't actually have one in my basement, and they're kind of hard to build. I could use a particle accelerator to generate some neutrons (basically by bashing nuclei around until some neutrons fall off) but while I do have a particle accelerator in my basement, it takes one a lot more serious than I can reasonably build to get neutrons out. Nuclear fusion reactions give out neutrons, though, so all I'd have to do would be to build a fusion reactor in my basement. Improbable as it sounds, this actually is feasible, provided I'm not trying to get any energy out.

The trick is that there's a fusion reactor, called the Farnsworth-Hirsch fusor, that is surprisingly simple to build. It is actually something of a cross between a fusion reactor and a particle accelerator: I'd set up an electrical potential in a spherical configuration, accelerating deuterium ions towards the center, where they'd crash into other deuterium ions, occasionally hard enough for fusion to happen. This fusion would release a fast neutron.

To make gold, then, all I'd have to do would be to build a fusor, surround it with a blanket of mercury, run it for a while, and then extract the gold from the mercury. Simple, really.

Let's look at the economics, though.

Suppose we want to make a kilogram of gold, giving us $40000. We need about a kilogram of mercury, costing $18. We also need about 5 g of deuterium (assuming perfect efficiency), which would cost about $30. Finally, we need the power to run the fusor. That's not going to be cheap. An optimistic number for the best fusor ever built is about 10^12 neutrons per second from about 4 kW input. That amounts to 9*10^14 neutrons per kilowatt-hour. Assuming perfect efficiency again, we need about 3*10^24 neutrons for our kilogram of gold, or 3 terawatt-hours, about the world's total energy production for an hour and a half. At $0.10 per kilowatt hour (I live in the land of cheap hydroelectricity) that's three hundred billion dollars.

There's a somewhat more disturbing possibility, though. Gold is easily obtained; you can just buy it. But as a global society, we try very hard to make sure you can't easily get plutonium, particularly plutonium-239, since that is well-suited to building atomic bombs. (You can make bombs out of uranium too, but that requires you to separate the different isotopes, which are very nearly chemically identical. Plutonium, on the other hand, can easily be separated from uranium since it is a different element.) Uranium isn't too hard to come by, especially "depleted uranium" (uranium with most of the uranium-235 removed) - armies fire the stuff at each other, for example. And if you had lots of U-238, a fusor would let you make plutonium out of it. The cost would be high, hopefully prohibitively so, but you could do it without doing anything that would put you on the radar of the IAEA. Fortunately, the power use is so outrageous we don't really need to worry about it.

So, in short, I could make gold in my basement, but not any appreciable quantity, and not for any kind of sensible price.

(*) Since gold is a noble metal, there aren't many chemical compounds that contain gold; unlike, say, iron, gold is often found on Earth as lumps of raw gold. So while in principle alchemists could have started from some gold compound and gotten the original gold back, this would not have been a very interesting accomplishment.

(**) There are actually situations where you can affect the nucleus by manipulating the electron shells. For example, if an isotope decays by electron capture, you can drastically slow down its decay by stripping away all its electrons. But stripping away all the electrons from a reasonably heavy element is one of those things that's virtually impossible under terrestrial conditions but not too rare astrophysically. In any case this has no effect on stable isotopes.

(***) Canadian dollars and American dollars are equivalent to astronomical accuracy.

Curve Fitting part 4: Validating Bayesian Code

2009-11-11T21:40:00.005-05:00

In my previous post, I wrote a tool to use Bayesian inference to ask whether a collection of photons represented a pulsed signal, and if so, what its parameters were. It gave plausible answers, but knowing my own fallibility, I really want some more careful test before trusting its output.

I previously talked about how to implement such a correctness test in a frequentist setting. But in a Bayesian setting, what one gets are probability distributions describing our knowledge about the hypotheses - how can a program test those?

After racking my brains for a while trying to find a way to do this, I asked about it on the scipy-user mailing list, and received a few useful suggestions, the most valuable of which was to send a polite email to Professor Tom Loredo. I did, and he replied immediately with a very helpful email, and some links to other people's work on the subject.

It turns out that, in the context of Bayesian fitting, the issue is called "calibration". The key realization is that if you select a model according to your prior distribution, generate a data set, and then do your fitting to generate posterior distributions, then your "true" parameters that you used to generate the data set look just as if they had been drawn from the posterior distribution.

This makes a certain kind of sense: after all, your posterior distribution is supposed to describe the distribution of models that might have generated your data set.

So how do you turn this into a test? After all, if you just get one sample from a distribution, it's pretty hard to say anything very definite, especially when the tails - the most unlikely extreme parameters - are not necessarily very well modeled. If you try to generate another sample, you pick a different model and must then generate a different data set, so you get a new point but a totally different posterior distribution. So now you have two points, drawn from two different distributions, and your situation has not necessarily improved.

There's a trick to let you work with a large sample: you transform them to all have the same distribution. You can do this because you know the posterior distribution you're working with; in my case I have its values sampled evenly. So I can construct the cumulative distribution function and use its inverse to get the posterior probability of a value less than the true model. This should be distributed uniformly no matter the shape of the posterior.

In fact, I'd rather use a slightly different trick: instead of looking at the probability of a value less than the true model, I'll use the probability of a value more extreme than the true model. Essentially I'll combine both tails. This has a more direct bearing on the question of confidence intervals, and still results in a uniform distribution as I repeat the process many times.

Once I have a long list of many values that should be distributed uniformly, there are any number of tests I can use; I'll use the handy but not particularly good Kolmogorov-Smirnov test. My unit testing code now looks like this:


def test_credible_interval():

    M = 100
    N = 50

    ps = []
    for i in range(M):
        f, p = np.random.rand(2)

        events = bayespf.generate(f, p, N)

        phases, fractions, r, P = bayespf.infer(events)

        frac_pdf = np.average(r,axis=0)
        fi = np.searchsorted(fractions, f)
        p = np.sum(frac_pdf[:fi])/np.sum(frac_pdf)
        p = 2*min(p, 1-p)

        ps.append(p)

    assert scipy.stats.kstest(ps,lambda x: x) > 0.01

Note that I don't bother testing the posterior for pulse phase. This is just laziness.

In a real Bayesian problem, there would usually be many more parameters, and I would probably not be able to evaluate the posterior on a grid like this. I'd probably be using some sort of Monte Carlo method, which would return a list of samples drawn from the posterior instead. But there are reasonable approaches in this setting too.

So far so good. But what about the probability that the signal is pulsed? In some sense hypothesis testing is just a special case of parameter fitting, but the discreteness of the "parameter" poses a problem: there's no way we can hope for a nice uniform distribution when only two values are possible. Professor Loredo very kindly sent me an unpublished note in which he addresses this problem, showing that if your code accepts a hypothesis whose posterior probability is P_crit, then a correct piece of code will be right - have chosen the correct hypothesis - with probability at least P_crit. Unfortunately P_crit is only a lower limit on the probability of getting things right; but there isn't really a way around this: suppose my code were working with a million photons a run. Then it would, in almost every case, give a probability very close to zero or one. There would be very, very few errors, no matter what value P_crit you set.

The fact that P_crit is only a lower limit does mean that this doesn't allow you to catch code that is too conservative: if your code makes fewer errors than its probabilities claim it should, this test has no way to tell that that's what's happening.

One must also choose P_crit carefully. In my example, if I choose P_crit=0.5, then code that simply chose randomly, or returned the same value all the time, would pass: after all, my priors specify equal probabilities for each model, and with P_crit=0.5 the code only needs to be right half the time. On the other hand, with P_crit really high, the code will almost never be wrong, although it only takes a few errors to fail, but it will almost never actually come to a conclusion, so you'll wait forever for evidence. There should be some fairly well-defined optimum value of P_crit, and I need to think more about what it is.

In any case, having selected P_crit and the other parameters, I can write a unit test for the hypothesis testing component:


def test_pulsed_probability():
    np.random.seed(0)
    p_accept = 0.75
    M = 50
    N = 30

    accepted = 0
    correct = 0
    total = 0
    while accepted<M:
        if np.random.rand()>0.5:
            f, p = np.random.rand(2)
            pulsed = True
        else:
            f, p = 0., 0.
            pulsed = False

        events = bayespf.generate(f, p, N)

        phases, fractions, r, P = bayespf.infer(events, n_phase=100, n_frac=101)

        if P>=p_accept:
            accepted += 1
            if pulsed:
                correct += 1
        if P<1-p_accept:
            accepted += 1
            if not pulsed:
                correct += 1
        total += 1
    assert 0.01<scipy.stats.binom(M,p_accept).cdf(correct)

As a final note, I'd like to thank Professor Loredo for his help, but note that any errors in the above code or description are entirely mine.

Curve Fitting part 3: Bayesian fitting

2009-11-11T21:00:00.000-05:00

When you fit a curve to data, you would usually like to be able to use the result to make statements about the world, perhaps something like "there's a fifty percent chance the slope is between 1 and 2". But this is a bit peculiar from a philosophical point of view: if your data is a set of measurements of some real-world phenomenon, then it's a bit funny to talk about probabilities that the slope has certain values. The phenomenon has some fixed slope, so we can't somehow repeat the experiment many times and see how often the slope is between 1 and 2. But there is a way to make such a statement meaningful: Bayesian inference.

The basic idea is that you use probabilities to quantify the degree of uncertainty you have about the world; you are using a system of logic that uses probability and probability distributions rather than binary logic. This may sound fuzzy and esoteric, but Bayesian logic is used very successfully in, for example, automatic systems to evaluate whether email is spam.

When applied to data, Bayesian reasoning lets you make meaningful statements about probabilities of parameter values, at the cost of making some explicit assumptions going in, and also at the cost of some potentially substantial computation. I'll work through an example of fitting a light curve to a set of photon arrival times, using a Bayesian procedure.

First the problem setting: suppose we observe a pulsar, whose period we know exactly (perhaps from radio observations), with an X-ray or gamma-ray telescope. We see some (fairly small) number of photons, and we want to know whether the flux we see is modulated at the pulsar period. We "fold" the photon arrival times, recording the pulsar phase when each one arrives. So our data is a collection of some hundreds or thousands of numbers between zero and one.

The model we'll fit includes some fraction f of photons whose arrival times are distributed as a sinusoid with a peak at phase p; the remaining fraction (1-f) are independent of phase.

The key idea of Bayesian curve fitting is that if you have some collection of hypotheses Hi about the world, each having some probability P(Hi), and you make some observation, there's a simple procedure to update these probabilities to reflect your new knowledge.

From a philosophical point of view, it's a bit worrisome to have to supply hypotheses (called "priors") about what the world is like before we ever make any observations. It amounts to making assumptions in the absence of data. But in fact there are assumptions built into the usual "frequentist" methods of inference as well, and in Bayesian inference the ability be explicit about the hypotheses at least makes it clear what's going on.

What assumptions should we make for our pulsar? Well, there might or might not be pulsations, so we'll have two basic hypotheses: no pulsations and pulsations. Absent any information, we'll assume these are equally likely. Then, if there are pulsations, we need to specify prior distributions of phase and pulsed fraction. Since both these parameters are between zero and one, we'll just take a so-called "flat prior" that makes all values between zero and one equally likely.

Given these priors, we need to figure out how to use our observed photon arrival times to update the priors to give us "posteriors". The general formula is:

P(Hi|D) = P(D|Hi) P(Hi) / P(D)

That is, the probability of hypothesis i given the data, P(Hi|D), equals the probability of the data given Hi, P(D|Hi), times the prior probability of Hi, P(Hi), divided by the probability of the data given any hypothesis.

The first thing to note is that P(D|Hi) is just what we need to evaluate for a maximum-likelihood estimate: how likely data like what we observe is to arrive given some hypothesis. We only need to define it up to a constant, since it appears in both numerator and denominator. For our problem, the probability density for pulse arrival times is p(f,p,t) = f(1+cos(2 pi (t-p)))+(1-f). So P(D|Hi) is the product of p(f,p,ti) for all events ti.

How do we form P(D)? Well, since we have two hypotheses, H0 (no pulsations) and H1 (pulsations), we can write P(D) = P(D|H0)+P(D|H1). Further, H1 is actually a family of hypotheses depending on two parameters, so we need to integrate P(D|H1) over all possible values of the parameters.

If we apply the above formula, then, we should get two things: a posterior probability that there are any pulsations at all, and a posterior probability distribution for the two parameters.

Let's look at python code to implement this. First of all, we're going to need to be able to generate fake data sets:


def generate(fraction, phase, n):
 m = np.random.binomial(n, fraction)
 pulsed = np.random.rand(m)
 c = np.sin(2*np.pi*pulsed)>np.random.rand(m)
 pulsed[c] *= -1
 pulsed += 0.25+phase
 pulsed %= 1

 r = np.concatenate((pulsed, np.random.rand(n-m)))
 np.random.shuffle(r)
 return r

This routine generates the photons in two parts. First it decides randomly how many come from the pulsed component. Then the photons from the pulsed component are generated uniformly. To convert this to a sinusoidal distribution we select some of the photons in the lower part and move them to the upper part. We then add in some uniformly-distributed photons, and shuffle the two samples together.

Now we write a routine to evaluate the probability density function:


def pdf_data_given_model(fraction, phase, x):
 return fraction*(1+np.cos(2*np.pi*(x-phase)))+(1-fraction)

Note that in spite of appearances, this routine can act on an array of values at once; this is important since python's interpreted nature makes each line of python take quite a long time.

And now the fitting routine:


def infer(events, n_phase=200, n_frac=201):

 events = np.asarray(events)
 phases = np.linspace(0,1,n_phase,endpoint=False)
 fractions = np.linspace(0,1,n_frac)

 lpdf = np.zeros((n_phase, n_frac))
 for e in events:
     lpdf += np.log(pdf_data_given_model(fractions, phases[:,np.newaxis], e))

 # This weird-looking hack avoids exponentiating very large numbers
 mx = np.amax(lpdf)
 p = np.exp(lpdf - mx)/np.average(np.exp(lpdf-mx))

 S = np.average(np.exp(lpdf))

 return phases, fractions, p, (S/(S+1))

This uses one of the simplest approaches to calculating the distribution and its integral: just evaluate on a grid. Integration then becomes averaging. More sophisticated Bayesian problems usually involve high-dimensional integrals, and so a whole elaborate machinery has evolved for efficiently evaluating these (for example the python package pymc).

Finally, some wrappers to generate a fake data set, call the fitting routine, and plot and print the results:


if __name__=='__main__':
 import pylab as pl

 events = generate(0.2,0.5,200)
 phases, fractions, r, P = infer(events)
 print "Probability the signal is pulsed: %f" % P

 pl.subplot(211)
 pl.contourf(fractions, phases, r)
 pl.xlabel("Pulsed fraction")
 pl.ylabel("Phase")
 pl.xlim(0,1)
 pl.ylim(0,1)

 pl.subplot(212)
 p = np.average(r,axis=0)
 li, mi, ui = np.searchsorted(np.cumsum(p)/np.sum(p),
         [scipy.stats.norm.cdf(-1), 0.5, scipy.stats.norm.cdf(1)])
 pl.plot(fractions, p)

 pl.xlabel("Pulsed fraction")
 pl.ylabel("Probability")
 pl.axvline(fractions[li])
 pl.axvline(fractions[mi])
 pl.axvline(fractions[ui])
 print ("Pulsed fraction: %f [%f, %f]" %
         (fractions[mi], fractions[li], fractions[ui]))

 pl.show()

One key point here is that when I want to know the distribution of pulsed fraction but don't care about the phase, I integrate (i.e. average) the joint distribution along the phase direction.

This gives us the following plot:

And the following output:


Probability the signal is pulsed: 0.450240
Pulsed fraction: 0.210000 [0.100000, 0.315000]

So it looks like the fitting routine is working: even with relatively few photons and a small pulsed fraction, it has selected quite good best-fit values. The probability that the signal is actually pulsed seems a little low, but keep in mind that we have only two hundred photons, and only forty of these are actually pulsed (while a Poisson uncertainty on the number of photons would be something like 14). But giving plausible results is not really enough: I want to systematically test this routine for correctness. But that will be another post.

Testing statistical tests

2009-11-07T19:00:00.004-05:00

Often one wants to ask something like "is there a periodic signal in this data?" or "is there a correlation between these two data sets?". Often there is some way to calculate a number that represents how much signal there is or how much correlation there is. But of course there is always noise and uncertainty in the data, and so it's possible that the apparent signal or correlation is actually just noise. So for such a number to be useful, it must also come with an estimate of how likely it would be to get such a value if there were no signal, only noise. Such an estimate is often called a p-value.

I should say, before I go on, that there are a number of things a p-value isn't: it's not (one minus) a probability that the signal is real. It's not a measure of the strength or importance of the signal. It's not even enough to tell you whether you should believe the signal is real: if you look at a hundred thousand data sets, you should expect to find many with p<0.01 even if there's nothing there. But this has been better discussed elsewhere, and it's not my topic today.

Today I'm addressing a specific issue: suppose that you have implemented some new way of quantifying the presence of periodic signals, or correlations, or whatever, and your implementation returns a p-value. How do you make sure that p-value is actually correct?

Of course, you can (and should!) go over your code carefully, with a critical eye. But evaluating the p-value often involves evaluating some arcane special function. How can you be sure that you implemented it correctly, or for that matter, that you got the math right in the first place? It turns out there's a nice way to test your code as a black box.

For our example test, let's take a well-known test: the chi-squared test. As I'm going to use it here, we have m measurements, each with measurement uncertainty one, and we want to know whether the underlying quantity is the same in all measurements. The chi-squared test works by constructing the quantity chi-squared, which is the sum of the squares of the differences of the measurements from the mean:


def chi_squared(values):
  return np.sum((values-np.mean(values))**2)

This quantity has a known distribution, naturally enough called the chi-squared distribution. It is parameterized by a number called the "degrees of freedom", which in our case is m-1, the number of measurements minus one because we subtracted the mean. For my p-value, I will just ask what the probability is that random noise would produce a chi-squared value larger than what we observed. (Note that we might well want to consider unusual any data set for which the chi-squared value is unnaturally small as well as ones where the chi-squared is unnaturally large, but that can only happen if we have overestimated the uncertainties on our measurements, so I'll leave it aside.)


def chi_squared_p(m,c):
  return scipy.stats.chi2(m-1).sf(c)

For a distribution in scipy.stats, the "sf" method evaluates the "survival function", that is, what fraction of the distribution's values are larger than the given value.

So now we have a test and a p-value calculation. How are we to check that we implemented it correctly? Well, let's pick an m and a p-value we're interested in, say p0 = 0.05. Let's also pick a number of repetitions, N. We will generate N fake data sets in which there is nothing but noise, and see how many times our statistic reports p<p0. This should be something close to p0*N. How close? Well, the statistic should behave like flipping an unfair coin with probability p0, so we can use the binomial distribution to find limits that should contain the number of false positives 98% of the time.


def test_p(N,p0,m):
  false_positives = 0
  for i in range(N):
      if chi_squared_p(m,chi_squared(np.random.randn(m)))<p0:
          false_positives += 1

  assert scipy.stats.binom(N,p0).cdf(false_positives)>0.005
  assert scipy.stats.binom(N,p0).sf(false_positives)>0.005

There we have a test; the more trials you let it run (N), the tighter constraints it puts on the p-value. Unfortunately, it fails 2% of the time even if everything's fine, and it's random, so if it fails, you can't restart running it in the debugger (since the next run will get different values). There's a way around these problems, too. The first problem can be avoided by running the test once, then if it fails, rerunning it. Then the test only reports a failure 0.04% of the time if all is well, but a real problem in the algorithm will probably still show up. The second problem you can solve by seeding the random number generator every time you run the test. In python, decorators provide a handy way to avoid having to write boilerplate code to do both these things for each test; I have two decorators, seed() and double_check, that do this. I'll omit their code here because they have some nastiness to work around limitations in nosetests (but you can find them, along with a more detailed example of the techniques in this post here).

This test is nice, but a little limited: it only tests one p value, p0. Since every time the statistical test runs it returns a p value, surely we can do better?

Well, since the p value is supposed to be the probability that a random data set will exceed the value obtained in a particular test, if we generate lots of noise samples, we should have the fraction whose value is less than p0 roughly equal to p0 for every p0 - in other words, the p-values returned should be uniformly distributed. So we can use a statistical test for uniform distribution to check whether they're coming out right! One such test, present in scipy, is the Kolmogorov-Smirnov test. That gives the following scheme:


def test_ks(N,m):
  pvalues = [chi_squared_p(m,chi_squared(np.random.randn(m))) for i in range(N)]

  D, p = scipy.stats.kstest(pvalues, lambda x: x)
  assert p>0.01

This does have the drawback that the K-S test is known to be most sensitive in the middle of the distribution, while what we care about is actually the tail. A minor modification can help, at the cost of some speed:


def test_ks_tail(N,m,p0):
  pvalues = []
  while len(pvalues)<N:
      pv = chi_squared_p(m,chi_squared(np.random.randn(m)))
      if pv<p0:
          pvalues.append(pv/p0)
  D, p = scipy.stats.kstest(pvalues, lambda x: x)
  assert p>0.01

This has the embarrassing problem that it's too good: it turns out that while this works for the chi-squared statistic I describe, for the K-S test itself, the p-value returned is actually rather approximate, so that this test tends to fail.

Finally, there's one essential thing to check: how stringent are these tests? If we return the wrong p-value, do they pick it up? A quick check can be done simply by modifying the chi-squred calculator to jigger the value a bit, then checking that the tests fail. With a thousand trials, the tests pass just fine if I return 1.01*chi-squared; I have to increase it to something like 1.05*chi-squared to start getting failures.

This brings us to the elephant in the room: the power of a statistical test. The p-value is really only half the answer; it tells us how likely we are to get an apparent signal when there's nothing there. It tells us nothing about whether we would actually see a signal if it were there; for that, you need a different quantity, the power of the test.

Cosmos

2009-11-06T21:14:00.004-05:00

In honor of Carl Sagan day, I'd like to link to Cosmos:

Nearly thirty years later, and in spite of tremendous advances in astronomy, Cosmos is still a wonderful introduction. Rather than focus on the science, Carl Sagan does a great job of sharing the wonder of discovery.

I am kicking myself now: when I was much younger, I was visiting a cousin who works at Cornell, and he pointed to Carl Sagan's office and asked if I wanted to go up and say hello. Out of shyness, I guess, I demurred, but it would have been a fascinating visit.

Eureka?

2009-05-24T15:00:00.001-04:00

The most exciting phrase to hear in science, the one that heralds new discoveries, is not 'Eureka!' (I found it!) but 'That's funny ...'
-Isaac Asimov

I recently found a very exciting new millisecond pulsar. But my first thought was not "Wow! A new millisecond pulsar!" but "Isn't that a suggestive bit of interference?"

To explain myself a bit further, I was looking at candidates from the drift-scan survey. These are all the periodic signals we picked up with the GBT, and they naturally include every cell phone call, car ignition, laptop computer, and worn-out electric blanket in the vicinity. Most are easy to distinguish from real pulsars, but some aren't. One characteristic pulsars generally have is they're noisy: after all, they're faint astronomical sources, so it would be very strange if they were so strong we couldn't see any noise in the observation. But in the plot showing 1023's signal, seen above, you can see, there's no noise. Turns out 1023 is just plain bright. I was not the first to misclassify it, though.

J1023 normally looks like a fairly ordinary 17th magnitude star. It's in star catalogs and photographic plates (enter RA 10 23 47.67 and Dec 00 38 41.2) going back as far as 1952, but it never got any attention in particular until 2002.

In 1998, the FIRST sky survey was carried out with the VLA. This is a rather different beast from a pulsar sky survey; it gives average brightnesses over several-minute integration times, so it's not going to have any luck detecting pulsations. But because the VLA is an interferometer, it's able to generate quite high-resolution images. This kind of survey, called a "continuum" sky survey (as opposed to a survey for spectral lines or a pulsar survey) is good for finding nebulae, radio stars, and galaxies. J1023, or to give it its full name, FIRST J102347.67+003841.2, caught the attention of people analyzing the FIRST data because it was variable: in the three observations just days apart, its radio brightness varied by a factor of at least three. Galaxies generally don't change so quickly, so they thought the source was worth following up.

Since the source had an optical counterpart, they chose to follow it up by taking optical spectra. These optical spectra, and those by another group, showed that it was blue and had double-peaked emission lines. A blue spectrum indicates very hot gas, emission lines indicate hot diffuse gas, and the fact that they were double-peaked is normally a sign that they were coming from an accretion disk: gas on one side of the disk is moving rapidly towards us, so its emission line is shifted towards the blue end of the spectrum, while gas on the other side is moving rapidly away, so its emission line is shifted towards the red end. Since the hottest gas is the nearest the center and the fastest moving, it produces the strongest emission, and you get a double-peaked spectrum. Other observers also looked at the source with high-speed photometry (just looking at how bright the object was) and found it was flickering, which is normal for accreting systems, as turbulent knots in the disk come and go. So people looked at all this data and classified J1023 as an unusual "cataclysmic variable", that is, an accreting white dwarf.

In 2002, observations showed that the emission lines were gone, the spectrum had gone back to the Sun-like colors observed in 1999 and earlier, and the flickering had tailed off. All that was left was a Sun-like star that varied in brightness by 0.4 magnitudes in a predictable fashion as it travelled around its orbit. This return to quiescence is normal for cataclysmic variables: they go through active phases and passive phases. But John Thorstensen and Eve Armstrong decided to try to come up with a model that explained the light curve (the brightnesses and colors).

When you have a binary system like this in which the stars are close together, the companion is heated by the white dwarf. When the hot side is facing us, the star is brighter and bluer than when the cool side is facing us. So it shouldn't be too hard to look at the light curve and figure out how bright the white dwarf is.

Well, as always, it's not as easy as it sounds, but by dint of great effort, Thorstensen and Armstrong managed to come up with a model that fit. But there was a problem: it needed a massive bright white dwarf. So massive, in fact, that it couldn't be a white dwarf, and so bright that we should have seen it in the spectrum, as a blue bump. There are many possible explanations for this sort of thing; software bugs, calibration problems, and so on. But they did a careful analysis, didn't find any of those, and decided to go out on a limb and suggest that the prevailing opinion was wrong: J1023 didn't contain a white dwarf at all, but a neutron star.

There the question stood for a few years. Homer et al. took X-ray observations of J1023, finding it was bright in X-rays, which supported Thorstensen and Armstrong but wasn't definitive. Thorstensen kept an eye on the system, taking a spectrum now and then to see if it was doing anything new. It wasn't. Then we came along.

When I found the pulsar, we weren't very sure of its position. We were using the low frequency of 350 MHz so that our beam would be big enough to see any piece of sky for two minutes as the Earth turned. But that big beam means that when we find a pulsar, all we know is a very rough position. Nevertheless, when we found a bright fast millisecond pulsar, we knew it was going to be interesting, so we requested a follow-up observation with the GBT.

The night before we were supposed to take the data, Jim Condon of the NRAO emailed us to point out that J1023 was in our beam. Ingrid Stairs, one of my collaborators on the survey, did some reading and found Thorstensen and Armstrong's paper. Suddenly there was the possibility that our millisecond pulsar might actually have been accreting in 2001. This was a big deal - Duncan Lorimer bet us all a drink that it wasn't the same source (covering his bets, I think - he was hoping it was the same source as much as the rest of us), and we were all excited. So we definitely wanted to find out whether it was the same source as soon as possible.

Unfortunately, the observation we had planned was another 350 MHz observation, just to see whether the source was real and to start to build up a timing solution for it. So as Scott Ransom and I prepared to run the observation, we argued: I wanted to look at 350 MHz first, to make sure there was something there at all, and Scott wanted to take a 2 GHz observation pointed at the position of J1023, so that the much smaller beam would tell us whether the source was really J1023 or not. In the end we compromised: we started at 350 MHz, and the pulsar came booming in right away. So we retracted the prime focus arm and switched receivers and pointed the GBT at J1023, and sure enough, there it was, loud and clear right at the position of J1023. We took the rest of the observation at 2 GHz, and immediately began requesting follow-up observations.

We initially planned to follow up with the GBT, since it was what we were used to, but Paulo Freire emailed us and asked us to please propose for time with Arecibo: there was not much else for Arecibo to look at at the time of day 1023 is visible there, and Arecibo's funding is threatened and it could really use a splashy discovery. (We were keeping the discovery quiet at this point, but Paulo was sharing an office with one of my collaborators, so there was no keeping it from him.) With Arecibo, this already-bright pulsar comes in beautifully, and we get nice clean timing.

Timing the pulsar, we quickly came up with a model of the orbit of the pulsar, and sure enough it agreed with Thorstensen and Armstrong's orbit. In fact, not only did the orbital period come out the same, if we extended our solution back to 2005, we got the same orbital phase as they did: over those three years, we were able to account for every single turn of the companion around the pulsar. Needless to say, this removed all doubt that our system was actually J1023.

The follow-up observations also revealed some peculiar phenomena, like plasma floating around the system and orbital period variations (very small, needless to say), but the essence of it is there: J1023 is a system with a neutron star and a companion that had an accretion disk for about a year around 2001 but is now a millisecond pulsar. The paper has just been published in Science, and is available on arxiv.org.

Unobtainium capacitors

2009-05-22T18:00:00.000-04:00

Continuing my train of thought from the other day, what are the theoretical limits for a capacitor? Just how much energy can you store in one?

If you try to analyze it in terms of voltages and capacitances, it becomes very tricky: you can make thinner plates closer together without increasing the weight or volume of the capacitor at all. The key to making the problem tractable is to realize that the energy is stored in the electric field. The goal in trying to store as much energy in a capacitor as possible is to get the electric field as strong as possible over as much volume as you can manage.

Since this will be an order-of-magnitude estimate anyway, let's assume that the electric field completely fills the capacitor. That leaves the maximum strength as the last parameter. Since we're looking for theoretical limits, let's assume you're limited by the tendency of the electric field to ionize the atoms in whatever you make this capacitor out of. So an electric field of about a volt per angstrom is just about all we can hope for. This gives us an energy density of about 400 Joules per cubic centimeter, or 0.1 watt-hours per cubic centimeter. This really isn't that much, and in fact according to Wikipedia there is a company claiming its capacitors do substantially better than this. Granted, it's a patent application, but where did I go wrong?

Well, there are certainly things I ignored - dielectric constant, for example - but I think the basic issue is thati in real capacitors the limit really is the strength of electric field your materials can tolerate. Since the energy density depends on the square of this number, clever engineers who can come up with extraordinarily tough materials have a fair amount of room to beat my crude estimate. But it seems clear to me that that room is far from unlimited: capacitors will continue to improve as energy storage devices, particularly in aspects I haven't touched on like leakage, but don't expect one to power your laser pistol any time soon.

A Missing Link

2009-05-21T19:22:00.005-04:00

We've been pretty sure that low-mass X-ray binaries turn into millisecond pulsars for a while now. But no one had ever seen the one turn into the other. Thanks to the drift-scan survey, now we have.

The system, J1023, is currently behaving just like a normal millisecond pulsar: it spins very regularly, and we see radio pulses every cycle. But in 2001, before anyone had any idea there might be a pulsar there, optical observations showed that it was flickering, blue, brighter than usual, and had double-peaked emission lines. All of these are signs of hot gas flowing in the system, and the double-peaked lines in particular are a pretty clear indication that the system had an accretion disk. Unfortunately nobody knew to point an X-ray telescope at J1023 while it was doing this, but it seems pretty clear in retrospect that it was an X-ray binary at the time.

In 2002, observations showed the emission lines were gone, the brightness and color were back to 1999 levels, and the flickering was tailing off to nothing. So the optical story stands today: no sign of an accretion disk. But there is clearly a radio pulsar there, and we're able to make some extremely good measurements of the system because of it. We know that the pulsar is 7.2 times as massive as its companion, for example, and if the pulsar has a typical neutron star mass of 1.4 solar masses, the system is 1.3 kiloparsecs away (about 4200 light-years) and we're seeing it at an angle of 46 degrees.

The image above is a pair of artists' renditions of it (well, sort of, I did most of the work using the software binsim and I'm hardly an artist, though Joeri van Leeuwen improved them significantly). They assume the pulsar's mass is 1.4 solar masses, and show the system as we think it was in 2000 and now. The hot disk of matter, present in 2000 but absent now, produced the optical emission that was observed, but (we think) it was blown away when a drop in the accretion rate allowed the radio pulsar to turn on, producing not just the beam of radio waves pictured in the image but also a powerful wind.

Since we know the system went through one roughly two-year active phase, it seems entirely possible that it will do so again within the next few years. If that happens, we'll be able to watch the formation of an accretion disk in a system where we have very good measurements of the orbit and system geometry from pulsar timing. That's never been seen before, and will be very exciting.

While I was the one who found this source, everyone involved in the drift-scan survey did essential work in making it possible for me to find it, and I worked with many other people to carry out the follow-up observations, so let me thank all my collaborators: Ingrid H. Stairs, Scott M. Ransom, Victoria M. Kaspi, Vladislav I. Kondratiev, Duncan R. Lorimer, Maura A. McLaughlin, Jason Boyles, Jason W. T. Hessels, Ryan Lynch, Joeri van Leeuwen, Mallory S. E. Roberts, Frederick Jenet, David J. Champion, Rachel Rosen, Brad N. Barlow, Bart H. Dunlap, and Ronald A. Remillard.
The paper describing the discovery has just been published in Science, and can also be read on arxiv.org.

Dar Williams

2009-05-18T20:35:00.002-04:00

I only really started to listen to music in the 90s, so my idea of what music should be is bleeps and clicks, with human voices only in the form of ironic movie samples. (This is of course precisely the range of music that can be made by geeky shut-ins with piles of computers and electronics. Literally, in some cases; I saw this one live.)

But as the result of a hard disk crash (sadly not recorded), I've been listening to a lot of music from other sources lately, and I find I'm really liking Dar Williams. She's got a nice voice, clever lyrics, and how can you dislike someone who writes a song about the Milgram experiment?

Where do millisecond pulsars come from?

2009-05-18T18:00:00.000-04:00

Most pulsars have periods around a second, and are spinning down very gradually. A few are clearly young (still hot, still in a supernova remnant, spinning down rapidly) and are faster, with periods as short as several tens of milliseconds. But there are other pulsars that are even faster - periods down to 1.4 milliseconds, that is, spinning some seven hundred times a second - and yet appear to be very old: cold, spinning down very gradually, very weak magnetic field. These are called millisecond pulsars, and it was a puzzle where they could possibly come from.

The key clue to the puzzle was their binarity. Most stars are found in systems of several stars orbiting each other. But pulsars are usually solitary. This is mostly because of their violent births: to make a pulsar, the star has to go supernova, and such a violent explosion tends to break up binary systems. It doesn't always break them up, though, and so some binary pulsars are known. If you look at the millisecond pulsars, though, most of them are actually in binary systems, unlike the normal pulsars (see the graph at right, based on the ATNF pulsar database). So we think the presence of a companion is key to making a millisecond pulsar.

The story, as we understand it, is this: A pulsar forms in a supernova. It is either in a binary system which survives, or it captures a companion. It lives out its life as a pulsar, spinning down gradually past the point where it is visible as a pulsar. The system stays like this for a long time. But eventually, one of two things happens: the companion starts to swell up into a red giant, or the orbit shrinks. In either case the system reaches a point where some of the matter at the surface of the companion is attracted more strongly by the neutron star than by the companion. This matter then falls down onto the pulsar.

Pulsars are not much more massive than the Sun, but they are much much smaller. So when matter falls onto one, tremendous amounts of energy are released. But remember that the two stars are in orbit, so the system is rotating. If you take a piece of matter from the companion, it will carry some angular momentum with it. To make the obligatory analogy, just as when a figure skater pulls in her arms, she speeds up, when you take matter from the companion and bring it in to the neutron star, the matter begins to rotate more rapidly. When this matter falls on the star, the star is spun up a little. I am glossing over numerous important details here, but the point is, when you start transferring mass to a neutron star, it is possible to spin it up.

So, we think that millisecond pulsars are old pulsars that have spun up, or "recycled", by accretion of mass from a companion. Observationally, we see systems where mass is being transferred: they're very bright X-ray sources. In a few cases we can actually tell how fast the neutron star is rotating, and sure enough, its period is down in the millisecond range. But these systems don't produce radio pulsations, presumably because all that matter falling in either blocks the radio waves or shorts out the emission mechanism (which needs a near-vacuum in the magnetosphere). So to make millisecond pulsars you need somehow to turn off the accretion and clear out all the matter, so that the radio pulsations can emerge. This transition hasn't been seen before, and the theorists have some difficulty explaining the population of objects that we see - while we see both millisecond pulsars and their hypothetical accreting progenitors, none of the progenitors seems to be positioned to turn into anything like the millisecond pulsars we actually see. So that last step, accretion turning off and radio pulsations starting, remains something of a mystery.

Energy storage

2009-05-16T16:38:00.005-04:00

Dan of Dan's Data just posted an article about science-fiction-y batteries. He discusses a battery based on matter and antimatter but dismisses it as inconvenient (in that the energy comes out as a shower of high-energy radiation). Instead he suggests a battery based on just cramming a lot of electrons - several grams of electrons! - into a AAA-sized package, and finds that this actually stores vastly more energy. There's an error there, but it's an interesting one.

First of all let's do his calculation a little more carefully. Let's suppose we're taking 11.5 grams of electrons - roughly two billion coulombs - and cramming them into a AAA-sized package. For simplicity, let's assume the package is a conductive sphere of radius 1 cm. This is actually a capacitor, with capacitance about 1 picofarad (incidentally, 1 pf is a fairly normal size for a capacitor, though they don't usually choose a conductive-sphere form factor). So the energy stored in this capacitor is (1/2) Q^2/C, or a staggering 2*10^30 Joules, the sun's total luminosity for a couple of hours.

(All calculations are done with the charming frink, though of course the errors are my own.)

There's a problem with this, though. Actually, there are several, but let's tackle the most fundamental first: if I convert 2*10^30 J to mass using E=mc^2, I get about twenty billion tonnes. This isn't just a meaningless figure: the energy stored in the capacitor really does weigh this much. If it seems weird for energy to have mass, well, I have to agree, but you can actually see it by looking at the periodic table: a helium atom is basically four hydrogen atoms after you stick two of the electrons to two of the protons and squash all the nucleons together. But each hydrogen atom has a mass of 1.007825 u, while helium has a mass of 4.002602 u, rather than the 4.0313 u that you would expect. The difference in mass is the mass of the energy in the nucleus. Since in helium, that energy is negative (pulling the nucleons together rather than pushing them apart), the mass is also negative, and helium weighs less than the hydrogen that makes it. But it's the same phenomenon.

One interesting question is, where is all that mass? It turns out that there is a nice answer. The little conductive sphere is surrounded by a strong electric field. In fact Gauss' law lets us work out just how strong it is at the surface: 2*10^23 V/m. It turns out that electric fields store energy. The stronger the electric field, the higher the energy density it stores. If you add up all the energy stored in the electric field of this charged sphere, over all space, you get exactly the energy you had to put in to charge it up. In other words, the energy is stored in the electric field, strongest at the sphere but extending out to infinity. And when you have an energy density, you have a mass density. Dividing by c^2, you get a mass density for the electric field of 2*10^15 g/cm^3. All that mass is due to the electric field, mostly close to the sphere.

Incidentally, there's a neat little calculation here. We know that electrons are weird quantum-mechanical beasts whose position and velocity are unavoidably ill-determined, and in fact as ar as we can tell an electron is a point particle. But let's ignore that for a moment and try to figure out how big an electron is. We don't have much to go on, because they don't sit still, and they interact with other things without touching them. But there is one place to look: they have a mass. Suppose we pretend that an electron is a tiny conductive sphere, with no mass of its own, but with charge smoothly distributed over its surface. Then the electric field will have a mass of its own, just as with our little conductive sphere. What if we just declare that the mass of the electron is entirely due to this electric field? Well, then we can calculate a radius! The value you get is 1.4*10^-15 m, which is tiny enough (about a hundred thousand times smaller than an atom) to not be totally unreasonable. You get a slightly different answer if you assume the electron is non-conductive, but leaving aside that issue, this size is the classical electron radius: 2.8*10^-15 m. Astonishingly enough, this number comes up in various physical phenomena, for example Thomson scattering of light by an electron.

So we can't, even with some pretty impressive unobtainium, cram 11.5 grams of electrons into a AAA-sized space. But if we decide that 11.5 grams should be the total weight of the battery, including energy, how much energy can we store? Well, it turns out, just a little less than the matter-antimatter battery. Since the matter-antimatter battery turns its contents entirely into energy, there's no way to beat its storage efficiency. But the above calculation shows that the charged metal sphere is nearly as efficient. We just charge it up to a slightly less outrageous voltage of 4*10^13 V. Now we've got something as good as the antimatter battery, but that provides us with handy electrons! What more could you ask?

Well, for one thing, you need to get rid of the electrons somehow once you've used them. Otherwise your laptop will be constantly giving out static shocks. Try and drive an electric car with them and you'll be shedding kilowatts of static shocks. Not ideal. Also, speaking of static, the electric field of this metal sphere will extend out to infinity, and will be quite strong even at some distance. At the surface of the sphere, the electric field will be 4*10^15 V/m, which is 4*10^5 volts per angstrom. So the voltage difference across an atom will be on the order of a million volts, compared to the 13.6 volts it would take to pull an electron off a hydrogen atom. So this little sphere will shred all atoms that come near it. Worse, that field extends out to infinity, so that at a meter away it's still 40 volts per angstrom. This is going to be a real pain. Isn't there some way to shield it?

Well, yes. You just put a spherical shell around it with the opposite charge. Or you abandon the sphere and just set up two parallel plates with opposite charges. Now the fields of the two plates cancel out (nearly), everywhere except between the plates.

Actually, this kind of energy storage device exists, and you probably have thousands in your house: it's a capacitor. Now, granted, it doesn't store the kind of energy density Dan is talking about, since it's not made of unobtainium by aliens, but an appropriate bank of capacitors can store a very great deal of energy and deliver it very quickly.

But it doesn't sound so impressive to ask the aliens for a capacitor.

[Update: I talk more about what one can hope for from a capacitor here.]

Finding pulsars part 2

2009-05-16T16:00:00.000-04:00

To find new pulsars, you obviously need to look at the sky with a big radio telescope. But once you've recorded this massive amount of data, you're only half done. Searching through the results, it turns out, is a massive computational task, even with modern computer clusters.

In principle, searching for pulsars in the data should be simple: they're nice periodic signals, so you just take a fast Fourier transform and the periodic signals show up as peaks in the power spectrum. In practice it's a lot more complicated than that.

The first problem you run into is interstellar dispersion. The space between stars is not empty: there is something like one hydrogen atom per cubic centimeter, on average. It is mostly hot, and at such a low density that means the hydrogen is in the form of plasma, loose electrons and protons. Plasma has a refractive index: it slows radio waves propagating through it very slightly. The absolute amount is small, maybe a few milliseconds, but it depends on the observing frequency. So if your pulsar has a one-millisecond pulse, then you might find that at the bottom of your observing band it's delayed by several milliseconds compared to the top of your observing band. This smears out a nice sharp pulse into a smeary blob that's much harder to detect. So when we look for pulsars, we need to divide the observing band into many channels and record the power in each channel separately. When you multiply a thousand channels by ten thousand power measurements a second, you get an awful lot of data.

To search for pulsars in this huge stream of data, you have to choose a dispersion measure (DM) and shift the channels accordingly, then add. Now you can do your Fourier transform and look for peaks.

The second problem that arises is acceleration. Some of the most interesting pulsars are not solitary stars. Instead, they are in orbit with a companion. Usually the pulsar is the more massive of the two, but it does not stand still; both stars orbit their common center of mass. Since these binary systems are often very close together and often have quite short periods - hours to days - the pulsar can be moving quite quickly, which produces a Doppler shift. A constant Doppler shift just changes the position of the peak, which is not a problem, but if the pulsar is accelerating - and moving in a circle involves acceleration - the Doppler shift can change. If the system is close enough, the Doppler shift can change substantially over the course of an observation, so that instead of the pulsar appearing at a single frequency in the Fourier transform, it gets spread out. Now, this may sound rare, and it is, but not only are these some of the most interesting systems to find, they are also the systems most likely to have been missed in previous surveys. So it's important to look for these accelerated systems, which you can do by looking not just for single peaks but for the characteristic shape of an accelerated peak. This is effective but takes substantial computing time.

Multiply together all the DM trials, Fourier frequencies, and possible accelerations, and you find that processing a single two-minute beam takes some thirty hours on a modern computer. Work out how many two-minute beams there are in the 1491 observing hours we carried out for the drift-scan survey, and you begin to see why we haven't yet processed all the data.

Once you've processed the data, though, you unfortunately do not get a nice tidy list of pulsars. There are rather a lot of processes that produce periodic radio signals, so what you get is a monster list of cell phone calls, car ignitions, telescope birdies, electric blankets, and a very few new pulsars. Winnowing through the huge mass of periodic signals is unfortunately less a question of clever Fourier techniques and more akin to spam filtering: many are obviously junk (they don't show interstellar dispersion, say), some are less-obvious junk, and a few are hard to tell from real pulsars even by humans.

My main contribution to the survey has been to run data through McGill's computer cluster and to sift through the results. I've done a certain amount of sifting by hand, but I've also written a few "spam filters" to try to reduce the number of obviously-not-pulsar results I have to look at. The real excitement, though, has been finding new pulsars and following them up.

Finding pulsars

2009-05-15T19:00:00.000-04:00

How do we find new pulsars? Well, there are a few places we know pulsars are likely to be - globular clusters, supernova remnants, X-ray binaries - but for the most part we just need to look at lots of sky. This isn't really feasible in X-rays - which you can only observe from satellites, on which time is precious - and apart from the exciting new gamma-ray results from Fermi, the only frequency range at which you can reasonably hope to find pulsars is the radio. So we more or less need to run a big radio survey that looks at lots of sky, one little patch at a time. But how do you get time on a multimillion-dollar telescope to look at sky where there's probably nothing?

Sometimes you can make a good science case for doing a survey. The most successful pulsar survey, the Parkes Multibeam survey, did this: they had a new multibeam receiver that would drastically improve the rate at which they could survey sky, so it was clear that time spent looking for pulsars would pay off. And it did, handsomely; they discovered over a thousand new pulsars.

More often, though, there are more pressing projects fighting for telescope time. So pulsar surveys need to come up with some other way to get time. The GBT 350 MHz drift scan survey took a clever approach: we were able to use time the telescope wasn't able to be used for anything else.

The GBT is an amazing machine, but with seven million kilograms running back and forth constantly, the azimuth track recently started to wear out. So the summer of 2007 was spent replacing the track. During this time, obviously the telescope couldn't move in azimuth without risking falling off the end of the track, so it couldn't follow sources as the Earth turned. But the receivers were all working fine, and the telescope could move in elevation just fine. So the people behind the drift scan survey pointed out that if they could simply use the receivers during the repairs, recording the sky as it rolled past, they could get useful data during the summer. This made good sense to the telescope access committee, so they were given the time.

The final arrangement was that the workers would work four twelve-hour days a week, during which time no data would be taken (arc welders are not conducive to good radio reception anyway). The rest of the time was devoted to the survey, which would be allowed to move the telescope in elevation twice a day. The survey used a relatively low frequency, 350 MHz, so that the beam would be big enough for any given piece of sky to spend about two minutes in view of the telescope.

Once the survey was set up, all we needed to do (I say we because it was at about this point that I got involved) was move the dish once a day, start the data recorders, and deal with the data. Dealing with the data was more involved than it sounds, because the data came in at a rate of about 90 GB/hour. So leaving analysis aside for the moment, we filled a couple of 800 GB disks a day, and somebody had to move all the data onto them and box them up for shipping to the various universities that would analyze the data. So when I went down to the telescope to observe as part of the survey, what I actually did was run a couple of scripts a day and shuffle disks.

In the end, the survey covered roughly a third of the celestial sphere, producing 134 TB of data. We're still analyzing the data - maybe only a quarter of it has actually been processed at this point, a year and a half later - but we've already found a number of new pulsars, some of them fascinating. So I think there are worse things the GBT could have done with those two months of track repair.

Safety

2009-05-15T04:01:00.005-04:00

I certainly understand the value of workplace safety. But really, two pages of alarming safety warnings for a java tone generator applet?

Also, among the very numerous warnings:

All sounds produced by DFG modules are licensed by Digital Recordings for your personal, environmentally-responsible and non-commercial use only.

"All sounds" is just sine waves, so I'm not exactly sure they can claim copyright on them, let alone use copyright restrictions to prevent me from damaging the environment (with sine waves?!)

Falling free

2009-05-07T05:05:00.007-04:00

There was a rather neat paper on arxiv.org today: The Flow Of Granular Matter Under Reduced-Gravity Conditions, by Hofmeister et al. They point out that while we understand the movement of, say, sand fairly well under terrestrial gravity, objects like the Moon and Mars are covered with fine powdery regolith, and we don't really understand how it moves under the lower gravity there. This includes, for example, the steepest stable slope, obviously important for understanding the dunes on Mars. So they did some experiments in reduced gravity - which is much harder than it sounds.

Gravity is tricky stuff. It's such a weak force, you can't manpiulate it by moving sources around, like you can with electricity and magnetism. You basically need to use whole planets to make any appreciable amount of gravity. But you can fake it with accelerations - and in fact it turns out that "fake" isn't the right word, since accelerations are indistinguishable from gravity. So if you want to do some experiment in higher than normal gravity, you can put it in a centrifuge. But if you want less than normal gravity, that's not so easy.

There are basically four ways to get reduced gravity. The first is to fake it: many experiments (wind tunnel tests for example) are done on scale models. The scaling requires cleverness, because (for example) the surface-to-volume ratio of a toy airplane is very different from that of a real airplane. So you also need to scale things like air speed and air density appropriately. So if you're studying some gravitational phenomenon, fluid flow perhaps, you may be able to build a model scaled so that the needed gravity for the model is one Earth gravity. This only works for systems with pretty simple physics, though, for which you understand the scaling.

The second way to get reduced gravity is to do your experiment on the ISS. Unfortunately, it costs hundreds of thousands of dollars per kilogram to get anything up there, so this is rather a last resort.

There is also the famous "Vomit Comet", which flies along parabolic trajectories to provide about half a minute of microgravity to those aboard. This is great for training astronauts (and wealthy space nerds) but it's very difficult to keep residual accelerations small enough for precision experiments.

This brings us to the system Hofmeister et al. used: drop towers. In principle they're very simple: you put your experiment in a box, then drop it from a tall tower. For as long as it takes to fall, you have microgravity. If you want to double the time, you can even fling it up from the bottom. Of course, the devil is in the details; for example you want to evacuate the tower so there aren't residual accelerations from air drag, and you need to plan on the capsule being pushed a few centimeters to the side by the Coriolis force. Plus, of course, you need to stop it non-destructively when it reaches the bottom. But all these problems are tractable, and in Bremen there is a drop tower that can give up to nine seconds of free-fall.

Hofmeister et al. used the Bremen drop tower for their experiment. Perhaps perversely, in order to get the reduced but non-zero gravity they needed, they put their experiments in a centrifuge inside the capsule. They hooked up a high-speed camera, and were able to track each grain of sand as it slid down the slopes. They found that the steepest stable angles and velocity patterns of the flows were rather poorly described by existing theory. So it looks like the theorists have some work ahead of them before they can understand the Martian dune fields.

Telescope of the month: the GBT

2009-05-01T18:00:00.001-04:00

The Robert C. Byrd telescope at Green Bank in West Virginia, called by almost everyone the GBT, is a radio telescope. With an elliptical dish 100 meters by 110 meters, it's the biggest fully-steerable telescope in the world - much to the chagrin of the people in the tourist center at Effelsberg, whose telescope is 100 meters but circular.

As a single-dish radio telescope, the GBT really only puts one pixel on the sky at a time, so making an image is a slow and laborious process. People do it, but where the GBT really shines is for looking at other aspects than spatial resolution. It really shines when looking at spectral lines - which lets you map out moving gas along the line of sight - and for pulsar observations.

The GBT is also used in conjunction with Arecibo for radar experiments: Arecibo has a big radar transmitter that can be used to bounce radio waves off objects in the Solar System. But since transmitting and receiving at the same station can be a problem, often projects will use the GBT as a receiver for these radar transmissions.

As a machine, the GBT is quite astonsihing: seven million kilograms of moving mass. The people at the visitor center like to say that it's the largest moving object on land, which is not quite true, but the counterexamples - the Saturn V aboard its transporter, and a colossal dragline excavator so big it walks instead of rolls - put the scale in context.

One of the most valuable features of the GBT is its location in the National Radio Quiet Zone. Not only is there no cell phone or pager service, the NRAO has a converted ice cream truck they drive around searching for sources of radio noise (for example in one story they tracked down one noise source to a local family's malfunctioning electric blanket; the family didn't want to part with it so the observatory bought them a new one). This absence of local interference makes a tremendous difference when looking for new pulsars (among many other tasks), and sure enough surveys at the GBT have been very productive of new pulsars.