RFI

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

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

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

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!

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

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

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

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

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

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

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

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

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

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.

Edit as of 2013 August 29: OpenSSH now has the ControlPersist option, which, in combination with ControlMaster and ControlPath, can be used to make the controlling SSH process background itself. You still can't (as far as I can tell) add new forwardings later on, but at least the problem of having the first connection be magical has gone away. I use this combination for one particular machine that accepts only password authentication (don't get me started) but that I use only as a gateway machine. Now the first time I try to connect through it I get prompted for a password but later connections just reuse the link. And because all it's doing is forward connections using the -W option, I don't care that I can't add port forwardings. 

Toxic waste

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

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

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

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

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:

1. 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.


2. 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.


3. 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.


4. 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.


5. 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

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.