Archive | December 2011

AMI Small Array Movie

A time lapse movie of the Arcminute Microkelvin Image (AMI) small array, I made by grabbing the webcam feed:

Watch on youtube.

Casacore installation: Addendum

So I was installing Casacore (as described here) on another machine and found I needed to do:

sudo cp -r lib/libcasa_*.so /usr/lib
sudo cp -r /usr/local/lib/libwcs*.so* /usr/lib/

To allow me to use some of the libraries. Bit odd but all seems to work, thought I should take a note of this!

A review of my 2011

2011 will always be one of the most memorable years in my life – it better be as I did get married 🙂
The year started off with me being back in the UK for my birthday which was great but my time at home was all to short and I jetted back off to Calgary (Canada). It is always a sad time to end up back in an empty flat far away from family. Though saying that I do have some
epic friends out in Calgary too. Once back in Calgary it wasn’t long before I was trying out something new… yep I went snowboarding for the first time:
me Snowboarding at COP
Was painful and so far has turned out to be the last time. I just can’t believe how poor my balance is.
In March my best-man Steve came over to visit me in Calgary and that was a lot of fun. Amongst the highlights of his trip over where seeing ice sculptures (and standing on the lake) at Lake Louise:
Lake Louise Ice Sculpture
It wouldn’t have been a proper trip over for him without taking in some hockey, so we went and
watched the outdoor WHL heritage game (that was damn cold!) and a couple of Flames games at the Saddledome:
Flames vs San Jose Sharks
We also went out and saw some world class snowboarding:
Snowboarding World Cup @ COP
At the end of April I found it hard to get home, as stuff was falling off the building across the road from me!
Building falling apart on 10th st...
In May, I said goodbye to Calgary as I had been offered a position at the University of Cambridge, back in the UK. I have to say I’m very glad to have come home and can spend a lot more time with my beautiful wife but I’ll always have very fond memories of Calgary and miss lots about it – in particular the people. I wrote a detailed review of my time in Calgary. Though before I left Calgary I got to go to the wonderfully named small town of Vulcan:
DSCF6677
Before I left Calgary I had a bit of a shock with a health issue, all was fine but lets just say I’ve since changed my work all hours, cut out go out drinking lots and living off caffeine.
In June I started at the Cavendish. I have to say I’m finding Cambridge a little different to the big
cities I’m used to. Oh and I got probably my best scientific result of the year accepted in June.
July was a very busy time. We were heavily involved in Wedding preparation though I did manage to stop for a few bits of fun. This included going to the British Grand Prix for the first time. Yes, it rained – it was a great time though.
Wet Racing
July also involved my stag do, but we don’t talk about that. It was great fun… we went go karting, had burgers and drank. Epic times.I even won a laser quest:
Top at Laser Quest
DSCF6157
August… well this was a hectic month but all got more relaxed after the 19th. On the 19th we got married:
Me and my wife
most of the better pictures aren’t mine and well I feel bad about linking them so won’t
(I wrote lots about the preparations before)
My wedding speech even made it onto youtube:

We took a short break to Bath as a honeymoon part 1.
In September I got to see my first supernova! (that’s a star blowing up by the way)
SN 2011fe
and I took my first photo of a comet, Garrard:
Comet Garrard C/2009 P1
October saw me start to settle into my new surroundings and I took my first trip out to the MRAO. I also had a paper on Ultra Steep Spectrum sources published.
At the of October we took our proper honeymoon out to Croatia, in particular Rovinj. It was absolutely beautiful.
Sunset over Rovinj
The whole Istrian peninsula was full of interesting things, including lots of Roman sites. We also went inland to the Plitvice Lakes National Park – somewhere I’d love to spend a bit more time walking around.
In November I got to go to the University of Birmingham Vale Fireworks for the first time in 2 years so I was quite pleased – always
feels good to celebrate the saving of parliament by blowing up part of the country:
DSCF9224
I got to take a trip over to Paris for ADASS XXI and I presented a poster on magnetic fields in GALFACTS. Elizabeth joined me after
the conference and we had a great weekend exploring Paris. Lots of fun.
Eiffel Tower at night
At the end of the month we went out to Wast Hills and did some more of Messier catalogue but I went outside and tried out my camera at some night time constellation shots (Orion here above the dome):
Wast Hills and Orion
and now its December. December has been a hectic month, and we are moving down to Cambridge in the new year so busy busy but we did manage to go out observing the other week and capture the Christmas tree cluser. Merry Christmas!
Christmas Tree Cluster
I’m looking forward to 2012… should be lots of fun. So far I have lots of observing planned, trip to India in Feb and tickets for the Olympics (alas just football). As a friend of mine said to me recently einen guten Rutsch ins neue Jahr

Observing last night – Jupiter + Dumbbell

It wasn’t particularly clear during our observing session, which was mostly aimed at taking a look at the ISS and Jupiter with our own eyes, but we focused the Meade out at Wast Hills and took a couple of images…
Jupiter_clouds
If you change the contrast a little you can clearly see the Moons:
Jupiter_moons
It was rather cloudy but we still managed a short exposure of the Dumbbell Nebula:
dumbell_neb_rgb

Creating FITS files in c++

A simple example, very similar to that given in the CFITSIO guidebook, on how to create a FITS file using CFITSIO. In this case I’m also building against some casacore libraries, but these aren’t going to be used in this little code snippet but the idea is to use casacore todo further analysis. I’m hoping to post more here over time. Anyway the code (this can also be found as a filebuild_fits.cpp:

/* Create a FITS file, using cfitsio and some casacore libraries
by Samuel George
21-11-2011
Compile: g++ build_fits.cpp -o build_fits -lcasa_casa -lcfitsio
*/
#include // STL iostream
#include
#include
extern “C”{
#include
}
using std::cerr;
using std::cout;
using std::endl;
using std::string;
int main()
{
cout << "Create a FITS file" << endl;
int lenTime(10), status (0), lenFreq(20);
long naxis(2), naxes[2] = {lenTime,lenFreq};
long nelements (lenTime*lenFreq);
long fpixel (1), exposure (1500);
char comment[] ="Total Exposure Time";
cout << comment << endl;
float pixels[lenFreq][lenTime];
// create an array of pixels
try {
for (int ii(0); ii < lenTime; ii++){
for (int jj(0); jj <lenFreq; jj++){
pixels[jj][ii] = 10.0*(ii+jj);
}
}
} catch (string message) {
cerr << "Error creating pixel array" << message << endl;
}
try { // write the image to a fits file…
fitsfile *fptr;
fits_create_file (&fptr, "!output.fits", &status);
fits_create_img (fptr, FLOAT_IMG,naxis,naxes,&status);
/* Write a keyword – its the address you pass */
fits_update_key(fptr,TLONG,"EXPOSURE",&exposure,comment,&status);
//write an array to the image
fits_write_img(fptr, TFLOAT, fpixel, nelements, pixels[0],&status);
fits_close_file(fptr,&status);
status = 0 ;
} catch (std::string message) {
cerr << message << endl;
}
}

Save the code and compile like so:
g++ build_fits.cpp -o build_fits -lcasa_casa -lcfitsio
and then run with ./build_fits

Spectral Extraction in Python

I’ve wrote about how to read in FITS files in Python before, but I thought I’d readdress as I’ve been writing lots about fitting and wanted to build up to fitting properly calibrated data. So in this example I will add how to extract a spectra too. For this example I’m going to use a calibration lamp taken at the University of Birmingham Observatory. Both the script described here and the data will be given below. The output looks like:
Spectral Plotting
So let’s start with reading in a FITS file. Firstly ensure you have the most useful per-requisites ([Pyfits], [numpy], [scipy]) and import them:

import pyfits, numpy, scipy

Now read in the file:

input_file = “A.fits”
hdulist = pyfits.open(input_file)

Get the data:

img_data = hdulist[0].data

Get the header:

img_header = hdulist[0].header

(more on FITS)
Now lets use some calibration data (which can be downloaded):

import pyfits, numpy, scipy, pylab
from scipy import *
from pylab import *
input_file = “Neon.fits”
hdulist = pyfits.open(input_file)
img_data = hdulist[0].data

We know that the spectra, by looking at the FITS file with DS9 has data in a small aperture. For now we will ignore any averaging. Lets take line 144 as we know there is data.

data_use = img_data[144]

Plot the spectra:

fig = figure(figsize=(9, 9)) #make a plot
ax1 = fig.add_subplot(111)
ax1.plot(img_data[144]) #spectrum
setp(gca(), ylabel=”Un-normalised power”, xlabel=”Pixel Position”)
pylab.savefig(“plot_spectra.png”)

We have now extracted the data for a spectral column. I will follow this post up in the not to distant future with how to now identify the lines and run spectral calibration. Once you have this you can apply to a target spectra. Of course I’ve so far also ignored dark/bias counts but these can be easily dealt with.
The above can be downloaded as extract_spectra.py and the data Neon.fits (this is zipped, use gunzip).

Christmas Tree Cluster

Probably not the best picture I’ll ever take of this source, I blame the bright moon, but here is a picture of the Christmas Tree Cluster that we took last night….
Christmas Tree Cluster
… so “ho ho ho, merry Christmas”.

Multiple Gaussian Fitting in Python

Yesterday I showed you [how to fit a single Gaussian in some data]. Today lets deal with the case of two Gaussians. This came about due to some students trying to fit two Gaussian’s to a shell star as the spectral line was altered from a simple Gaussian, actually there is a nice P-Cygni dip in there data so you should be able to recover the absorption line by this kind of fitting. Anyway, fitting 2 Gaussian’s is basically the same thing as fitting one in python but with the added function. I’m thinking I’ll work up how to deal with a whole spectrum including source finding. In this case what you have to deal with is that there are two sources and so a rough estimation of the peak position of both is crucial to the fit (well in the way it is implemented).
You basically want to end up with something like this:
Fit 2 Gaussians in Python
Where the green points are the data, the blue dashed line the fit and the red line where the maximum is in the array. The fitting is done via a least-squares fitting routine.
In this example we will first start by generating some data, skipping the input from a file (but you can use the code from the 1 Gaussian example).

import scipy, numpy, pylab, asciidata
from numpy import *
from pylab import *
from scipy import *
from scipy.optimize import leastsq
#generate some data
gaussian = lambda x: 3*exp(-(10-x)**2/10.) + 1*exp(-(30-x)**2/10.)#change the parameters as you see fit
y_power = gaussian(arange(100))
x_pos = arange(100)
gauss_fit = lambda p, x: p[0]*(1/sqrt(2*pi*(p[2]**2)))*exp(-(x-p[1])**2/(2*p[2]**2))+p[3]*(1/sqrt(2*pi*(p[5]**2)))*exp(-(x-p[4])**2/(2*p[5]**2)) #1d Gaussian func
e_gauss_fit = lambda p, x, y: (gauss_fit(p,x) -y) #1d Gaussian fit
v0= [1,10,1,1,30,1] #inital guesses for Gaussian Fit. – just do it around the peaks
out = leastsq(e_gauss_fit, v0[:], args=(x_pos, y_power), maxfev=100000, full_output=1) #Gauss Fit
v = out[0] #fit parameters out
covar = out[1] #covariance matrix output
xxx = arange(min(x_pos),max(x_pos),x_pos[1]-x_pos[0])
ccc = gauss_fit(v,xxx) # this will only work if the units are pixel and not wavelength
fig = figure(figsize=(9, 9)) #make a plot
ax1 = fig.add_subplot(111)
ax1.plot(x_pos,y_power,’gs’) #spectrum
ax1.plot(xxx,ccc,’b–‘) #fitted spectrum
ax1.axvline(x=xxx[where(ccc == max(ccc))[0]][0],color=’r’) #max position in data
setp(gca(), ylabel=”power”, xlabel=”pixel position”)
pylab.savefig(“plotfitting.png”)
print “p[0], a1: “, v[0]
print “p[1], mu1: “, v[1]
print “p[2], sigma1: “, v[2]
print “p[3], a2: “, v[0]
print “p[4], mu2: “, v[1]
print “p[5], sigma2: “, v[2]

The full detailed listing is given in [multiplegaussfit.py]

Gaussian Fitting in python

I spend a lot of my time working on noise statistics and of course and an important part of this is how to fit signals. I was asked earlier for an example code on how to fit a Gaussian, in particular fitting well defined signals. In python this is quite straight forward as you already have a bunch of nice libraries to count on.
So the first thing todo is to ensure you have these libraries. These are scipy, numpy, matplotlib and asciidata. Of course, you don’t need matplotlib if you don’t want todo the plotting but or asciddata if you don’t want to import an ascii file. (note: to students in 3rd year obs lab at Brum – you have this on hydra)
If you want to create some test data you can do something like:

from pylab import *
gaussian = lambda x: 1*exp(-(3-x)**2/10.) #change the parameters as you see fit
data = gaussian(arange(100))
x_pos = arange(100)

Now on to the fitting, this is simplistic and works well with “clean” data.

import warnings
warnings.simplefilter(“ignore”,DeprecationWarning)
import matplotlib
matplotlib.use(‘Agg’) #avoid X11 issues.
import scipy, numpy, pylab, asciidata
from numpy import *
from pylab import *
from scipy import *
from scipy.optimize import leastsq
gauss_fit = lambda p, x: p[0]*(1/sqrt(2*pi*(p[2]**2)))*exp(-(x-p[1])**2/(2*p[2]**2)) #1d Gaussian func
e_gauss_fit = lambda p, x, y: (gauss_fit(p,x) -y) #1d Gaussian fit
data = “test_gauss.dat” #your input data – should be x,y
demo = asciidata.open(data)
x_pos = double(demo[0]) #this is expected to be pixel number
y_power = double(demo[1])
v0= [1.0,mean(x_pos),1.0] #inital guesses for Gaussian Fit. $just do it around the peak
out = leastsq(e_gauss_fit, v0[:], args=(x_pos, y_power), maxfev=100000, full_output=1) #Gauss Fit
v = out[0] #fit parammeters out
covar = out[1] #covariance matrix output
xxx = arange(min(x_pos),max(x_pos),x_pos[1]-x_pos[0])
ccc = gauss_fit(v,xxx) # this will only work if the units are pixel and not wavelength
fig = figure(figsize=(9, 9)) #make a plot
ax1 = fig.add_subplot(111)
ax1.plot(x_pos,y_power,’gs’) #spectrum
ax1.plot(xxx,ccc,’b–‘) #fitted spectrum
ax1.axvline(x=xxx[where(ccc == max(ccc))[0]][0],color=’r’) #max position in data
setp(gca(), ylabel=”power”, xlabel=”pixel position”)
pylab.savefig(“plotfitting.png”)
print “p[0], a: “, v[0]
print “p[1], mu: “, v[1]
print “p[2], sigma: “, v[2]

Here is the output plot (red line is the maximum, blue is the fit and green points are the data):
plotfitting
This code can all be found in a more useful format [simplegaussfit.py]
Download this, change the input data parameter to your asciifile and run as simplegaussfit.py – ensuring you have the correct file permissions to execute (you can change this with chmod).
Related to this page is another page I had on [fitting arbitrary functions to data].

Solar Spectrum

The Sun’s spectrum taken with a digital camera and a low resolution spectrograph, you can see some nice absorption lines:
Solar Spectrum