Krioma.net

Its been ages since I did any proper C/C++ so I had to remind myself of a few things, particularly the compilation of libraries. Libraries are really useful as it simplifies multiple use and sharing of components. I'm a Linux person so I'm going to look at this with a Linux hat on. In Linux there are two types; Static libraries (.a) and dynamically linked shared object libraries (.so). Let take a look at static libraries.

Lets create a library, firstly lets define some code, say in a file example1.c:

To create a library you can do:

cc -c example1.c

which compiles the code. Then you need to make the libary:

ar -cvq libexample.a example1.o

(you can of course include other object files here)

Say you call the function in you maincode.c to now compile this do:

cc -o maincode maincode.c libexample.a

Taken at the University of Birmingham's Observatory at Wast Hills, the Omega Nebula:

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

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):

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

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:

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]

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

... so "ho ho ho, merry Christmas".

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:

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

I had an odd CASA issue yesterday. It was complaining about missing leap second information for TAI_UTC! Seemed really quite odd to me as I hadn't done anything to my installation... anyway simple fix is to update the latest leap second information into CASA by running a simple command in the root casa install directory:

cd $CASA_ROOT/data rsync -avz rsync.aoc.nrao.edu::casadata .

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

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

If you change the contrast a little you can clearly see the Moons:

It was rather cloudy but we still managed a short exposure of the Dumbbell Nebula:

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:

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:

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:

We also went out and saw some world class snowboarding:

At the end of April I found it hard to get home, as stuff was falling off the building across the road from me!

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:

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.

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:

August... well this was a hectic month but all got more relaxed after the 19th. On the 19th we got married:

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)

and I took my first photo of a comet, Garrard:

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.

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:

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.

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):

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!

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

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!

At the end of each year, yes I know I'm a few days early but bah I have to move house at the start of the new year, I like to summarise the music I've been listening to. The last few years have been much easier as last.fm keeps hold of that info. I suspect that its missing about 20% of the music I listen to, which isn't scrobbled for one reason and another. It might be a little higher than that but I don't like to include the car journeys when we are using the CD player.

I was a little shocked to see the number 1 artist of the year, and indeed a few different names in the top 10. I'm putting this down to a Calgary scene influence... (number, artist, plays)

1 Daft Punk 1,163
2 Linkin Park 664
3 Avril Lavigne 552
4 My Chemical Romance 526
5 Tiësto 522
6 Rise Against 285
7 Angels & Airwaves 271
8 Dead by April 259
9 Lostprophets 239
10 Disturbed 238

I have to say I did love the Daft Punk Tron album so I'm guessing this is what dominates that. I also find Tiesto quite good for background music whilst I'm reading (in particular I really think Kaleidoscope fits in really well with the Garth Nix books I've been reading).

Not that much new for 2011 in this list though. Wonder if that will change in 2012.

Been handed a dead hard drive and want to try and recover some of the data? That's the problem I had the other day. So I turned to some good old linux tools to try and recover the data. One of my favourite recovery tools is ddrescue. It works really well and given enough time can do a pretty solid job of data recovery - of course a better solution is proper backups but that's not always the way people go. Anyway, this is as much a note for myself but here is how to do some basic recovery using it:

sudo ddrescue -r 3 /dev/sde2 imaging loging

This images /dev/sed2 and produces the image file "imaging" with the log "logging".

To try and extract files you could use something like foremost

sudo foremost -w -i imaging -o /recovery/foremost

To make an audit of the files that can be recovered, recovering them with:

sudo foremost -i imaging -o /recovery/foremost2

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

Watch on youtube.