Menu
▾
▴
#469 voigt function
this patch add the voigt function (convolution of a gaussian and a lorentian) widely used in spectroscopy:
R.J. Wells Rapid Approximation to the Voigt/Faddeeva Function and its Derivatives JQSRT 62 (1999), pp 29-48.
http://www.atm.ox.ac.uk/user/wells/voigt.html
p.s. It's my first patch... I hope it's in the right format!
I am just having a first look at this patch now.
It looks basically OK, and I don't see any reason we can't add this after cleanup.
I have a couple of comments and questions:
- Is there a reason that all the variables in the routine are declared as static?
- Not all C compilers accept C++ style comments (// rather than /* ... */), so gnuplot code should not use them.
- Most people, certainly including me, have no idea what this function is used for. Please provide an explanation for the documentation. I went to the URL you quote, but even the JQSRT paper it links to does not bother to explain what the function is used for! (Yeah, I know. Mathematics journals and articles are notorious for failing to include background/significance in the introduction. We practicing scientists are appalled. ) At the very least it needs documentation telling what the parameters are and what are the domain/range limitations.
- How do I test that it is working? It might be best to provide a demo, voigt.dem, so that the output can be verified on different platforms.
- I get two compilation warnings:
specfun.c:722: warning: ISO C90 forbids mixed declarations and code
specfun.c:804: warning: ISO C90 forbids mixed declarations and code
The Voight function is a convolution of a gaussian and a lorentian and it takes into account two phyical mechanism of absorption (or emission) line spectra.
It's widely used in astrophysics.
An instroduction can be found here: http://en.wikipedia.org/wiki/Voigt_profile
The implemented fit widely used is explained here
http://www.eso.org/sci/publications/messenger/archive/no.80-jun95/messenger-no80-37-41.pdf
but this implementation is rather old and very rough, the humlicheck version that I copied int gnuplot (coming from the last post link) is much better and has a lower error (1e-4).
Nowadays there are even better approximations of the Voight function but it's more expensive and complicated. IMHO Humicheck it's a good bet.
I will add a figure showing different voigt profiles dove with the gnuplot implementation and the ESO (midas FitLyman package) showing that this implementation is much better.
voigt curves differences between hulicheck and ESO-MiDAS-FITLYMAN implementations
Sorry I forgot to answer to the other questions more about coding:
the static declaration it's just because I reused the ibeta and confrac just below were I put the voigt function...
a simple plot would be
plot [-10:10] for [n=1:5] voigt(x,n) title sprintf("Voigt (x,%d)",n)
and yes it's important because when I did a bibliografic seach of the voigt function I found a lot of times the voigt(y,x) rather than voigt(x,y)
the latex formula is
H(x,a)=\int{\frac{e^{y^2}}{(x-y)^2+a^2}}dy
... and there is no need for a static declaration of all the variables....
static double humlik(double x, double y)
{
const double c[6] = { 1.0117281, -0.75197147, 0.012557727,
0.010022008, -2.4206814e-4, 5.0084806e-7 };
const double s[6] = { 1.393237, 0.23115241, -0.15535147,
0.0062183662, 9.1908299e-5, -6.2752596e-7 };
const double t[6] = { 0.31424038, 0.94778839, 1.5976826,
2.2795071, 3.020637, 3.8897249 };
const double rrtpi = 0.56418958; // 1/SQRT(pi)
double a0, d0, d2, e0, e2, e4, h0, h2, h4, h6,
p0, p2, p4, p6, p8, z0, z2, z4, z6, z8;
double mf[6], pf[6], mq[6], pq[6], xm[6], ym[6], xp[6], yp[6];
double old_y = -1.;
bool rg1, rg2, rg3;
double xlim0, xlim1, xlim2, xlim3, xlim4;
double yq, yrrtpi;
if (y != old_y) {
old_y = y;
yq = y * y;
yrrtpi = y * rrtpi;
rg1 = true, rg2 = true, rg3 = true;
if (y < 70.55) {
xlim0 = sqrt(y * (40. - y * 3.6) + 15100.);
xlim1 = (y >= 8.425 ? 0. : sqrt(164. - y * (y * 1.8 + 4.3)));
xlim2 = 6.8 - y;
xlim3 = y * 2.4;
xlim4 = y * 18.1 + 1.65;
if (y <= 1e-6)
xlim2 = xlim1 = xlim0;
}
}
double abx = fabs(x);
double xq = abx * abx;
if (abx >= xlim0 || y >= 70.55) // Region 0 algorithm
return yrrtpi / (xq + yq);
else if (abx >= xlim1) { // Humlicek W4 Region 1
if (rg1) { // First point in Region 1
rg1 = false;
a0 = yq + 0.5; //Region 1 y-dependents
d0 = a0 * a0;
d2 = yq + yq - 1.;
}
return rrtpi / (d0 + xq * (d2 + xq)) * y * (a0 + xq);
}
else if (abx > xlim2) { // Humlicek W4 Region 2
if (rg2) { //First point in Region 2
rg2 = false;
// Region 2 y-dependents
h0 = yq * (yq * (yq * (yq + 6.) + 10.5) + 4.5) + 0.5625;
h2 = yq * (yq * (yq * 4. + 6.) + 9.) - 4.5;
h4 = 10.5 - yq * (6. - yq * 6.);
h6 = yq * 4. - 6.;
e0 = yq * (yq * (yq + 5.5) + 8.25) + 1.875;
e2 = yq * (yq * 3. + 1.) + 5.25;
e4 = h6 * 0.75;
}
return rrtpi / (h0 + xq * (h2 + xq * (h4 + xq * (h6 + xq))))
* y * (e0 + xq * (e2 + xq * (e4 + xq)));
}
else if (abx < xlim3) { // Humlicek W4 Region 3
if (rg3) { // First point in Region 3
rg3 = false;
//Region 3 y-dependents
z0 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y
+ 13.3988) + 88.26741) + 369.1989) + 1074.409)
+ 2256.981) + 3447.629) + 3764.966) + 2802.87)
+ 1280.829) + 272.1014;
z2 = y * (y * (y * (y * (y * (y * (y * (y * 5. + 53.59518)
+ 266.2987) + 793.4273) + 1549.675) + 2037.31)
+ 1758.336) + 902.3066) + 211.678;
z4 = y * (y * (y * (y * (y * (y * 10. + 80.39278) + 269.2916)
+ 479.2576) + 497.3014) + 308.1852) + 78.86585;
z6 = y * (y * (y * (y * 10. + 53.59518) + 92.75679)
+ 55.02933) + 22.03523;
z8 = y * (y * 5. + 13.3988) + 1.49646;
p0 = y * (y * (y * (y * (y * (y * (y * (y * (y * 0.3183291
+ 4.264678) + 27.93941) + 115.3772) + 328.2151) +
662.8097) + 946.897) + 919.4955) + 549.3954)
+ 153.5168;
p2 = y * (y * (y * (y * (y * (y * (y * 1.2733163 + 12.79458)
+ 56.81652) + 139.4665) + 189.773) + 124.5975)
- 1.322256) - 34.16955;
p4 = y * (y * (y * (y * (y * 1.9099744 + 12.79568)
+ 29.81482) + 24.01655) + 10.46332) + 2.584042;
p6 = y * (y * (y * 1.273316 + 4.266322) + 0.9377051)
- 0.07272979;
p8 = y * .3183291 + 5.480304e-4;
}
return 1.7724538 / (z0 + xq * (z2 + xq * (z4 + xq * (z6 +
xq * (z8 + xq)))))
* (p0 + xq * (p2 + xq * (p4 + xq * (p6 + xq * p8))));
}
else { // Humlicek CPF12 algorithm
double ypy0 = y + 1.5;
double ypy0q = ypy0 * ypy0;
int j;
for (j = 0; j <= 5; ++j) {
double d = x - t[j];
mq[j] = d * d;
mf[j] = 1. / (mq[j] + ypy0q);
xm[j] = mf[j] * d;
ym[j] = mf[j] * ypy0;
d = x + t[j];
pq[j] = d * d;
pf[j] = 1. / (pq[j] + ypy0q);
xp[j] = pf[j] * d;
yp[j] = pf[j] * ypy0;
}
double k = 0.;
if (abx <= xlim4) // Humlicek CPF12 Region I
for (j = 0; j <= 5; ++j)
k += c[j] * (ym[j]+yp[j]) - s[j] * (xm[j]-xp[j]);
else { // Humlicek CPF12 Region II
double yf = y + 3.;
for (j = 0; j <= 5; ++j)
k += (c[j] * (mq[j] * mf[j] - ym[j] * 1.5)
+ s[j] * yf * xm[j]) / (mq[j] + 2.25)
+ (c[j] * (pq[j] * pf[j] - yp[j] * 1.5)
- s[j] * yf * xp[j]) / (pq[j] + 2.25);
k = y * k + exp(-xq);
}
return k;
}
}
It's a great idea to put the Voigt function into gnuplot. It's really very often used for peak fitting in spectroscopy.
Can you please update the code including the patches you have put into the comment? Please the reference where the code comes from and its precision -- I've found this very similar:
http://www.sfr-fresh.com/unix/misc/grace-5.1.22.tar.gz:a/grace-5.1.22/src/humlik.c
Please add also a patch to gnuplot.doc about this function. Mention that under the same name voigt() are both Voigt and Fadeeva functions (real and complex argument).
Finally, I'm pleased to mention that prof. Josef Humlíček is the chief of the Department of Condensed Matter Physics I work in. His paper "Analysis of spectrum line profiles" from 1984 is available from our homepage http://www.physics.muni.cz/ufkl/Publications/ufkl_publ.shtml#odkaz1984 .
patch version 2 without static variables
I've uploaded the patch.
The code probably comes from grace (or at least from the same source), I really do not remember...
I first tried the fortran routine given in the first link, than I stripped out the routine from the ESO package FITLYMAN and I've also looked to the fityk package (pretty good fitting gui....) but there they use a more refined version.
thanks to prof. Josef Humlíček!
> The code probably comes from grace (or at least from the same source), I
> really do not remember...
Comments in the code itself indicate that it was auto-translated into C from a Fortran program by R.J. Wells. But Wells' ftp site is no longer accessible, so I can't see what it says there. As best as I can make out, the proper attribution of the code is the following. Is this correct?
/*
* Calculate the Voigt/Faddeeva function with relative error less than 10^(-4).
* (see http://www.atm.ox.ac.uk/user/wells/voigt.html\)
* arguments:
* x, y - Faddeeva/Voigt function arguments
* return value -- voigt
*
* Algorithm: Josef Humlíček JQSRT 27 (1982) pp 437
* Fortran program by J.R. Wells JQSRT 62 (1999) pp 29-48.
* Translated to C++ with f2c program and modified by Marcin Wojdyr
* Minor adaptations from C++ to C by E. Stambulchik
* Adapted for gnuplot Tommaso Vinci
*/
I could still use some help understanding the domain and arguments to this function.
Is it true that x and y are the real and imaginary components of a single complex number?
If so, since gnuplot supports complex numbers, should this be made into a single-argument function taking a complex number as input?