#469 voigt function

closed-accepted
nobody
None
5
2010-03-21
2010-01-19
Tommaso Vinci
No

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!

Discussion

  • Ethan Merritt
    Ethan Merritt
    2010-02-15

    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

     
  • Tommaso Vinci
    Tommaso Vinci
    2010-02-16

    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.

     
  • Tommaso Vinci
    Tommaso Vinci
    2010-02-16

    voigt curves differences between hulicheck and ESO-MiDAS-FITLYMAN implementations

     
    Attachments
  • Tommaso Vinci
    Tommaso Vinci
    2010-02-16

    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

     
  • Tommaso Vinci
    Tommaso Vinci
    2010-02-16

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

     
  • Petr Mikulik
    Petr Mikulik
    2010-02-16

    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 .

     
  • Tommaso Vinci
    Tommaso Vinci
    2010-02-17

    patch version 2 without static variables

     
    Attachments
  • Tommaso Vinci
    Tommaso Vinci
    2010-02-17

    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!

     
  • Ethan Merritt
    Ethan Merritt
    2010-02-17

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

     
  • Ethan Merritt
    Ethan Merritt
    2010-03-21

    • status: open --> closed-accepted