With the help of Michael Croucher of the walkingrandomly weblog, who graciously spent some time optimising one of the slowest components of the Matlab code, we were able to come up with an approximation to the Bessel function which was about two times faster than Matlab’s internal function. However, then I had to smoothen out the crossover point using error functions, and now it has lost some of the speed it had gained in favour of smoothness. It is now about 40% faster than Matlab’s Bessel function, when approximated using only 5 polynomial terms.

It looks like this, with the crossover gradient centre at 4.5:

So there remains a small discrepancy for such low polynomial numbers, but it should be possible to obtain a good speed increase with this in cylinder and disk scattering functions. Please grab yourself a copy of the code, and I’ll be happy (albeit perhaps slow in the coming two weeks) to include more improvements!

bessel_li_gross

p.s.: I forgot to mention how to compile the code for your computer. At this point, there’s a compiled version available for intel Macs. If you do not have that one, but you do have the matlab (”mex”) compiler, I suggest you follow the following steps:

1. make sure your mex compiler can find the following include files: mex.h, math.h, omp.h (in my case, I copied a suitable version of omp.h into the bessel_li_gross directory

2. go to the bessel_li_gross directory and type “mex blg_component.c”, or if your compiler complains about not being able to find the omp.h file that you just put in the directory: “mex -I./ blg_component.c” (adds the current path to the search path for files under unix)

3. Play. You can also send me your .mex file so I can add it to the package.

Binning of data is often done to SAXS data in order to reduce the size of the data, “improve statistics” or otherwise retransform the data to better suit the subsequent plotting or fitting of the data. Common 2D and 1D binning methods are “linear” binning, where the q-range is divided into bins of equal size, and “logarithmic” binning, where the q range is divided just so, that the bins are equidistant on a logarithmic q-scale.

Depending on your binning method, you may be weighting your data differently in subsequent fits of no weighting is done by uncertainty or error of the intensity. If you would like to weigh your data points in the fitting by their errors, you need to compute the error for each pixel, which is often given as the square root of the number of counts on your detector. This does imply that you either need a single photon counting detector or that you need to scale your detector output to photon counts.

An alternative to this is to rebin your data just so, that all points have the same intensity, meaning that you have wider bins at higher q than at low q. Since each bin then contains the same relative amount of intensity, the error in each bin should be the same (similar solutions can be achieved with non-position sensitive detectors such as the cyberstar, by counting in each position until a certain threshold value has been reached).

I have written a small Matlab program which achieves just this, and is fast enough to do this for large datasets (i.e. 2D images consisting of 2048×2048 pixels) within seconds. I invite you all to take a look at it, and, as always, suggestions or alternative methods are more than welcome. The input is self-explanatory, but the output isn’t necessarily. The output consists of four vectors:

• qbin: the mean centres of the bins
• dqbin: the width of the bins
• Ibin: the intensity in each bin
• Ibdq: the binned intensity divided by the width of the bins. This is given for plotting purposes (since plotting qbin vs. Ibin would give you more or less a straight line).

Good luck!

equal intensity binning routine:

 Read the rest of this entry »

Herewith a small program I wrote a while ago and some documentation on how it’s built up.

The Matlab program can calculate the scattering pattern of a superellips of revolution. I hope it is of some use to some of you, if only for satisfying curiosity. The program runs in Matlab, opens a GUI, which requires the input of five variables.

The first two dictate the limits in q that one is interested in. The minimum q value determines the maximum size of the real-space box, and the distance between the minimum q and the maximum q determines the pixel resolution in real space. Good values for these are, for example, 0.01 for the minimum q and 2 for the maximum q. 

The three other variables determine the size and shape of the superellips-of-revolution. The superellips is generated  with width a (in Angstrom), height b, and curvature r. The ellips is then revolved around the vertical axis to generate the ellips-of-revolution. A nice value for these are, for example, a=20, b=100, r=4. 

Please see the attached documentation and have fun running the program. The program is distributed under the GPL V3 license. geomsuperell.pdfsuperellips-of-revolution scattering pattern simulator