Christmas posts

Christmas time - blog time.
Today I will switch finally to Fityk and it's many uses.
One of them is spectroscopy, Mossbauer spectroscopy. A very fascinating and useful method.
I will skip large part of introduction. For people using Mossbauer for the first time I highly recommend to get a copy of [1] or [2]. Well, yes, one can use Fityk for Mossbauer, but mostly for the simple cases and the following restrictions apply.

• The effective thickness t_eff of the sample must be small, so that the fitting of the spectra can be done with Lorentzian peakshapes.
  
• The values you estimate are not too presize.
  

The Effective thickness is defined in the textbooks like:

t_eff = d * n_a * s₀ * f * a_a,

where:
d is the thickness of the studied material in cm;
n_a - the number of the Mossbauer atoms, for example, iron per unit cell volume (1/cm³);
a_a - abundance of the Mossbauer isotope (natural abundence of ⁵⁷Fe is 0.02119);
s₀ - the maximum resonance crossection, it is constant for every Mossbauer isotope (depending on excitations of nuclear spin) and has value of 256.6*10^-20 cm² for iron;
f - recoiless fraction or Lamb-Mossbauer factor, in general it has value ~0.7, and this value changes depending on the isotope, crystallographic environment, etc. With heating f->0, and f->1 under pressure or with cooling.

Effective thickness is unitless. If t_eff<<1, than we can use Lorentz functions to extract quantitative information from Mossbauer spectra. The basic desctiption of the lineshapes is:

• singlet(one Lorentz function)
  
• duplet (two Lorentz functions)
  
• sixtet (six Lorentz functions)
  

Singlet is a single Lorentz function, it has been already builtin. Doublet and sixtet shapes one adds to the Fityk commands like this:

define DubletLorentzUA(areaL=0.1, areaR=0.1, fwhm=0.19, center=0, qua=0)=
LorentzianA(areaL, center-qua/2, fwhm)+
LorentzianA(areaR, center+qua/2, fwhm)+
define SixtetLorentzA(area16=3, area25=2, area34=0.1, width16=1, width25=1, width34=0.19, center, qua=0, bhf=10.59)=
LorentzianA(area16*area34, center+qua/2-0.5*bhf, width16*width34)+
LorentzianA(area25*area34, center-qua/2-0.289*bhf, width25*width34)+
LorentzianA(area34, center-qua/2-0.079*bhf, width34)+
LorentzianA(area34, center-qua/2+0.079*bhf, width34)+
LorentzianA(area25*area34, center-qua/2+0.289*bhf, width25*width34)+
LorentzianA(area16*area34, center+qua/2+0.5*bhf, width16*width34)

This is all folks. Couple of comments in the end. Hyperfine magnetic field is used in Tesla units in literature, however, one measures this parameter in mm/s, velocity (this is why bhf=10.59). The conversion factor from mm/s to T is 3.097.
What about folding spectra? The simple recipy is below, if we assume channel 256 as folding point:

Y[0..256]=y[256-n-1]+y[256+n] in @0
a[256..512]=false
delete(not a) in @0
plot [:] [:]

At first line we fold, at second and at the third line we delete channels above the folding point. We do not need them after folding anymore. And the last thing to mention is the recallibration of the velocity limits of the Mossbauer experiment:

#calibration using known values of speeds a=true
$xconv=-10.0955 #here we provide calibration velocity
X=-2*n/(M-1)*$xconv+$xconv
plot [] []

I hope this post was not too technical and some can find something usefull.

gnuplot, cygwin and gnuplot custom points

There are not many free or commersial graphical packages out there which allow to use custom, user-defined points. Amazing Gnuplot allows you to use images for points. This post is based on gnuplot demo.

There are many ways one can create custom shapes/points in gnuplot. One of them is to provide gnuplot special postscript file. Another brilliant example is described here. However, I have not seen any example using images as points. For example .png images. That is a trick I would like to share.

One can find all the files used in example in my Skydrive folder. Today we play with the birthrate data available from WolframAlpha engine. For instance, I would like to plot the the probable number of people in Germany, France, Italy and Spain starting from 2008 to 2010. And I would like to use images for points to separate different countires. Here is the image I would to make in Gnuplot:

How do we do that? First, we need images, I use icons from Mark James. Of course we need gnuplot. Linux users can get a recent version from their Linux distribution repository or compile it themselves, I assume that the same is true for Mac fans. The Windows users can obtain specially compled version from gnuplot site or use one provided in cygwin. I preffer the latter one because this way I have bash command completition possibilities. It makes my life definetely easier. Some a little bit outdated cygwin gnuplot tips can be found here. If you need advice, please do not hesitate to ask. For complete set, you can install linux virtual machine or famous andLinux distribution. It is your call.

Once we have everything and before we begin, here is how I use gnuplot. Simple data visualizations are done in command line like this:

plot [] [] "birthrates.txt" u 1:2 title "Germany" w p pt 6 ps 3

This command line will plot XY data, where X are years, Y is the population at year X with circles, pointsize 3. For quick and dirty data visualization it works. More commands and their description can be found in gnuplot manual. This manual corresponds to gnuplot version 4.4. If you use older versions of gnuplot, you might find differences in syntax/features.
For publication purposes, I highly recommend to write a small script/macro file with right settings of fonts, scales, etc. This way you will always be prepared and you don not have to type every time gnuplot commands. It makes life easier, gnuplot learning process faster. For instance, let us suppose, you made a beautiful plot and want to save it as picture. You have done it, it is simple:

gnuplot> set term png size 640x480 enhanced font "Tahoma,10"
gnuplot> set output "my_beautiful_figure.png"
gnuplot> replot
gnuplot> set term x11; # or wxt if you have it
gnuplot> replot

Well I can do it one time, two times, but if I am going to do it regular fashion, I will create a special file, containing the same commands. This implies - commands without 'gnuplot>' of course. Any time I need a picture, I will load this macro file ('gnu2png') with a command:

gnuplot> load "gnu2png"

This is all for the introduction to the gnuplot. Let us finally do the real stuff. :)

The file 'birthrates.txt' contains the following data:

#        germany          france              italy                  spain
2008   82300000.00    63900000.00    59600000.00      44500000.00
2009   82221650.40    64243143.00    59900384.00      44939660.00
2010   82143375.39    64588128.68    60202281.94      45383663.84

And here comes the first part of my gnuplot script. Full version can be found here.

1:	#general stuff
2:	unset key
3:	set title "Expected population for selected countries Data taken from www.wolframalpha.com" font "Tahoma,8"
4:
5:	set xtics font "Tahoma,8"
6:	set ytics font "Tahoma,8"
7:	unset ylabel
8:	set xlabel "Year" font "Tahoma,9"

We do the following. Line 2: deleting plot's legend. Line 3: setting the title for the plot. Line 5, 6, 7: setting the fonts for xtics, ytics and xlabel.

10:	#picture/screen size
11:	picx=400; picy=400;
12:
13:	#point size - depends on image size
14:	iw=16;ih=16; #real image size
15:	ixscale=1.3; iyscale=1.3; #scale
16:	iw=iw*ixscale;ih=ih*iyscale;
17:
18:	#x and y minimum values
19:	xi=2007.; xa=2021.;
20:	yi=4.0e7; ya=9.0E7;
21:
22:	#set terminal
23:	set term x11 size picx,picy
24:
25:	#ImageCalculated[Height,Width]:
26:	#not precize, but good first approximation
27:	icw=abs(xa-xi)/picx; ich=abs(ya-yi)/picy;
28:
29:	#helpful styles to avoid too many words
30:	set macros
31:	filestyle="binary filetype=png"
32:	filestyle1="with rgbimage notitle"

Here is the necessary tricky part. Line 11: setting the terminal or picture size. Yes, yes, it is the size we will use for the final image. Line 14: setting the sizes of real images, here, .png images we will use for each point. We assume all points will have initial .png images of same dimensions (16x16px). Line 15: setting the adjustable scale. These parameters control the x and y scaling of the points (~point size). Line 16: recalculating widths and heights of images. Lines 19,20: setting the limits for our figure. [xy]i are the minimum, [xy]a are the maximum values. Line 23: setting the visual terminal. It is important to see what we get. Line 27: providing necessary parameters, we need to convert plot coordinates to pixels, and recalculate pixels to plot coordinates. At this point we calculate the conversion terms for width (icw) and height (ich). Line 30, 31, 32: setting some helpful parameters, for more information, please read the manual.

At this point the preparations are over, now we need to plot the data. Bad news: we need to plot point by point. This command allow us to plot only one point of the graph

plot [xi:xa] [yi:ya]  "de.png" \
@filestyle \
origin=(2008-icw*iw/2,82300000-ich*ih/2) \
dx=icw*ixscale dy=ich*iyscale \
@filestyle1

We plot a figure with x axis starting at 'xi' and ending at 'xa'. The same - for the y axis. We use 'de.png' icon for one point. By using 'origin=' command we set the position of the point. In this case the center of the 'de.png' image should correspond to the value of the point we want to plot. The 'dx=', 'dy=' commands set the size of point. Remember that we need to convert pixels to plot coordinates? This is what we do here. The '@filestyle' and '@filestyle1' substitution parameters we have defined above. And once again - we have to provide this line/lines for every single point on our graph.

The Good news: well people, we can use bash (native or cygwin's) or even Excel's/OpenOffice Calc capabilities to quickly get and convert data to the format we need. Believe me, it is fast if one uses spreadsheets. Some of us use them any way. I will not describe it here due to simplicity of the task. The hint is: we need to construct a special line containing the x and y coordinates of the point. These coordinates are the only thing one needs to change for gnuplot description of the single point. I will provide the solution for Linux/cygwin users. Behold, this line outputs the lines necessary to plot Spain related data (5th column):

cat birthrates.txt |grep -v "#"|awk '{print " \"es.png\" @filestyle origin=("$1"-icw*iw/2,"$5"-ich*ih/2) \
dx=icw*ixscale dy=ich*iyscale @filestyle1, \\"}'

Let us take a closer look. What we see here? The first two commands:
cat birthrates.txt |grep -v "#"| output the whole data file 'birthrates.txt' line by line, excluding the lines which contain '#' symbol. Line by line the output is piped to awk command which selects the first ($1) and the fifth ($5) columns and prepares a line we need for gnuplot. Alhough, this solution looks ugly, it is fast and efficient.

This is all folks for today. Just in case: the data, scripts and images one can download from Skydrive.

WolframAlpha and why I love its widgets

I am trying out new Capabilities of the Wolfram Alpha engine. Well, it is not that complicated as Mathematica, but still very useful, especially the new Widgets feature.

Let's learn by example. I need to have some easily accessible tool to recalculate lattice parameters from hexagonal to rhombohedral setting. How do I do that?

Well, first of all, I go to developer.wolframalpha.com and register myself as a 'developer'. Next, I watch the tutorial demonstrating the process of creating widgets. Five minutes later I start to shape my own.
The principle is very simple. You create a formula, assign variables, shape your widget, make some style adjustments and .. it is ready!

If we take a closer look.. the formula behind the widget is simple. I take A and C lattice parameters for the hexagonal cell and recalculate them:

{ArcCos[(15*a^2-c^2)/(2*c^2+6*a^2)]/pi*180, Sqrt[(c^2+3*a^2)/9]} where a=[//number:${A}//], c=[//number:${C}//]

As a result I obtain two values: one of them is the rhombohedral angle (Degrees), the other is the rhombohedral lattice parameter. The only trick was to make WolframAlpha engine to reuse the variables c,a or to set the simultaneously. We set all of them by where a=[//number:${A}//], c=[//number:${C}//].
Is that simple? I believe it is!

Now, let us answer the question - why we need it. My reasons are simple. I need something online, so I do not have to crawl all over my folders trying to find the solution. It also saves me some amount of disk space and backup efforts, which are more annoying than lack of space.
WolframAlpha will never displace Mathematica, but this tool is indeed worth to use. The widget which performs a reverse transformation (rhombohedral to hexagonal) can be found here or below.
Good day to you!

simple fityk tricks

Finally, I got stuck with experimental data and have some free time.
Ok:) tips I would like to share.

Fityk is great because:

• tool is free and works on different platforms (Windows, Mac OS?, Linux)
  
• program allows use of user defined models
  
• fitted models can be saved
  

Of course I miss features of:

• saving color of lines, size of points for different functions used in the model
  
• custom formating of the output (i errors)
  
• using constrains
  
• changing of visual baseline (in Fityk it is fixed to 0, all shapes are plotted relative to it) to some user defined constant value.
  

But I still like it. Comparing other free software (Gnuplot) to Fityk, I must say that these programs do different things and are in principle supplementary. Today's simple tricks:

• estimating errors of the fitted model: s=max(y-F(x)). We set values for the standard deviation as a constant value for all points. This value is considered as maximum deviation of experimental points y[n] from the fitted model F(x).
  
• setting manual constrains, like limiting parameter $param to positive values. For this we need to create another variable, for instance, $dummy and set: $dummy=~{10}; $param=sqrt($dummy^2). This way we define that Fityk may vary the $dummy variable, and $param variable depends linearly on it and is always positive. This trick is very useful if you need to avoid undesirable subtraction/addition of different Gaussian or Lorentz peaks (Spectroscopy, X-ray diffraction). Need to mention, that you have to control both widths/heights of the Lorentz/Gauss.
  

I strongly recommend to read the Fityk's Users Guide for some additional tricks like summing up the datasets, etc. The summation of datasets can be used as datasets subtraction, for instance, to remove calculated background from the X-ray dataset. Let us assume that @0 and @1 contain the observed spectrum and the calculated background. A simple line y=-y in @1; @2= sum_same_x @0+@1; delete @1; delete @0; will give us a spectrum without the background.

GSAS vs Fullprof - what I like, what do you like?

I am far from being professional with modern crystallography tools. However, here comes my impression of what I like in these two software packages (GSAS, Fullprof). I do not include comments on JANA2006, another marvelous piece of software. My comments are subjective and there are more comments than tips and tricks to it.

Fullprof. I know, the Fullprof .pcr files are difficult to handle, but I find them more and more attractive with highlighting (my previous post).

• LeBail fitting is fast. First, you set DUM parameter to 1, Pcr=1, AUT=0 (skip auto assignment to the fitting variables), Number of refined parameters=0 (according to manual). Then, please set JBT=2 and IRF=0 for all phases under consideration. After that you run the .pcr file with the Fullprof, fit background (polynomial or by manual selection of the background points). Setting Pcr=1 is really crucial if you want to prevent messing up your fitting.
  
• I really like combination of highlighting of .pcr file syntax with manual numerical assignment of codewords for different fitted parameters. I usually set AUT=0 in the beginning, then fit the background(assigning CODEWORDS 11-61), then I fit the lattice + zero shifts (some other CODEWORDS), then the peak shapes and etc. If now I assign the Number of fitted parameters only to 6, I will fit only the background. Then, I can slowly increase the Number of fitted parameters and fit the lattice and so on. Only if I am sure I will not mess up the fitting, I can set the AUT=1. I find this way of selecting the parameters to fit is faster, specially with LeBail fitting, than that checkbox style of GSAS.
  
• Setting up the magnetic structure is way faster in Fullprof compared with the GSAS.
  
• I like the way how simulations work in Fullprof. However, there is an important trick to it. It looks like Fullprof does not like really large numbers (version 2000/2006 in peak intensities/heights). If you have a really strange simulation diffraction pattern with peaks cut off to 0 intensity in the middle, pay attention to: scale of the phase, ATZ of the phase (take a look in .out file), occupation numbers. These are the pictures, first shows bad simulation, second - good.
  
  
  
  
  
  
• RTFM:)
  

GSAS.

• Yes, GSAS has a way much better tutorials. The combination of CMPR+EXPGUI is nicely explained in the same tutorials.
  
• The simulation of the diffraction spectra (except the magnetic phase) works like a breeze in GSAS.
  
• RTFM is really important.
  
• The texture fitting of the GSAS is more sophisticated than in Fullprof. If I have strange textures in my material, I will prefer GSAS.
  
• Background fitting in GSAS is outstanding, different models, fitting of the manually selected background points.
  

Of course, there are many things I do not like, but they are minor. Usually, I have many complains about the user interface. But this is not the complain I would express to the authors of these packages. The guys have done a brilliant job. Thank you guys!!!
However, I think, that investing whole lot of money and research to the creation of even better tools with the cooperation with the IT guys should be at the same priority as investing in the International Tables of Crystallography, International Crystallography Database and their support. I would like to have better tools.