Grace, known also by its command xmgrace, is a program for creating scientific graphs and is also capable of common data analysis tasks like fitting, Fourier transformation, etc. Unfortunately, the built-in mathematical functions are in some cases not enough. The Grace authors added therefore a mechanism to facilitate the loading of external functions.

This post basically follows the Grace User’s Guide part about DL Modules, just using a different example.

In NMR it is common to describe the magnetization build-up or decay with a Kohlrausch-William-Watts (KWW) function. Essentially it is an exponential function with stretching parameter \(\beta\):

\[S(t)=e^{-\left(\frac{t}{T_1}\right)^{\beta}}\]

In \(test\) xmgrace syntax the fit function can then be described with:

exp(-(x/A0)^A1)

A two step magnetization decay can be described with the following function:

\[S(t)=M_0\left[ \left(1-f\right)\cdot e^{-\left(\frac{t}{T_{1,2}}\right)^{\beta_1}}+f\cdot e^{-\left(\frac{t}{T_{1,2}}\right)^{\beta_2}}\right]\]

In xmgrace syntax this becomes:

A0*( (1-A1)*exp(-(x/A2^A3) + A1*exp(-(x/A4)^A5) )

with A0 = \(M_0\) the magnetization amplitude, A1 = \(f\) the component part, A2, A4 = \(T_{1,1}\), \(T_{1,2}\) the time constants and A3, A5 = \(\beta_1,\,\beta_2\) the stretching parameters.

Now it would be nice to access such a long function directly instead of typing/copying it every time we need it. This is now the point were loading external functions comes in to play.

What we need (according to the Grace manual):

  • shared library compiled with position independent code (PIC)
  • the function you want to use must map to one of the following Grace types and return a double (see table below)
  • if the function does not correspond to any of the Grace types you need to create wrapper of the function.
Grace type Description
f_of_i a function of 1 int variable
f_of_d a function of 1 double variable
f_of_nn a function of 2 int parameters
f_of_nd a function of 1 int parameter and 1 double variable
f_of_dd a function of 2 double variables
f_of_nnd a function of 2 int parameters and 1 double variable
f_of_ppd a function of 2 double parameters and 1 double variable
f_of_pppd a function of 3 double parameters and 1 double variable
f_of_ppppd a function of 4 double parameters and 1 double variable
f_of_pppppd a function of 5 double parameters and 1 double variable

Simple case

Let’s assume for the beginning we want to use the digamma function from the GNU Scientific Library (GSL). It is generally defined as:

\[\psi(x) = \Gamma'(x)/\Gamma(x)\]

The actual GSL function is gsl_sf_psi (double x) and returns a double, and we want to call it directly in Grace.

Assuming the GSL library is installed in the standard way (apt-get or rpm) so that it is in the default search path for dynamically loaded libraries, there is nothing more to do than starting Grace with:

xmgrace -pexec \
  'USE "gsl_sf_psi" type f_of_d from "/usr/lib/libgsl.so.0" alias "psi"'

You can now access the \(\Psi(x)\) function with psi(x) anywhere in Grace. In this way you can get access to all of the special functions offered by the GSL, from the manual:

The library includes routines for calculating the values of Airy functions, Bessel functions, Clausen functions, Coulomb wave functions, Coupling coefficients, the Dawson function, Debye functions, Dilogarithms, Elliptic integrals, Jacobi elliptic functions, Error functions, Exponential integrals, Fermi-Dirac functions, Gamma functions, Gegenbauer functions, Hypergeometric functions, Laguerre functions, Legendre functions and Spherical Harmonics, the Psi (Digamma) Function, Synchrotron functions, Transport functions, Trigonometric functions and Zeta functions.

Of course you can use any other math library, or you can create your own, which we will do in the next section.

Our own code

The bi-exponential decay can be coded in C. Here is the content of an implementation of the KWW function and the bi-exponential function discussed above:

#include <math.h>
double kww(double signal, double T1, double beta, double t) {
	double y;
	y = signal*exp(-pow(t/T1,beta));
	return y;
}


double biexp(double f, double T11, double beta1, double T12, double beta2, double t) {
	double y;
	y = kww(1-f, T11, beta1, t) + kww(f, T12, beta2, t);
	return y;
}

As you can see we define first the KWW functions, and use them to define our bi-exponential function. Note also that I do not use the \(M_0\) parameter. This is due to the fact that there can be not more than 5 parameters and one variable (see the table about Grace types).

In order to use this library we need to compile it first as a shared library (-shared) with position independent code (-fPIC) and, because we use math functions, we also need to link it to libm.so (-lm):

gcc -O2 -shared -fPIC -lm xmgrace_functions.c \
	-o xmgrace_functions.so

Then we need to start Grace with the following parameters (in the same directory 1):

xmgrace -pexec \
 'USE "biexp" type f_of_pppppd from "xmgrace_functions.so" alias "biexp"'

This tells Grace to use the biexp function which takes 5 parameters and 1 variable (f_pppppd) from the xmgrace_fucntions.so library and provide it to the user as a function with name biexp. Finally, here is an example where I used the biexp function to calculate a decay:

Plotting a bi-exponential decay in Grace via an external function.

If you use this function on a regular basis it is better to put the Grace command in your ` ~/.grace/gracerc.user`. This file is read upon every start of Grace and will provide your custom function in every instance, ready to be used. Here is for example an excerpt of my grace user file:

DEVICE "PNG" DPI 600
DEVICE "PNG" OP "transparent:on,compression:6" 
USE "biexp" type f_of_pppppd from "xmgrace_functions.so" alias "biexp"

Do not forget to change/set your LD_LIBRARY_PATH.

Unfortunately this is not supported on Mac OS X, but it works nice on Linux. The example above is tested and worked on Debian Wheezy.

  1. You can start Grace anywhere, but you have to make sure it can find your library. One way is to set the LD_LIBRARY_PATH variable to the proper directory. How to do this depends on your shell (BASH, TCSH, DASH,…)