"Fossies" - the Fresh Open Source Software Archive

Member "fityk-1.3.1/fityk/voigt.cpp" (13 May 2016, 16388 Bytes) of package /linux/misc/fityk-1.3.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "voigt.cpp" see the Fossies "Dox" file reference documentation.

    1 
    2 // fastest_humlik.for and humdev.for - from Bob Wells Voigt Function Page
    3 // http://www.atm.ox.ac.uk/user/wells/voigt.html
    4 // Translated to C++ with f2c program and modified by M.W.
    5 // It can be slower than original, I haven't compared the speed.
    6 
    7 #define BUILDING_LIBFITYK
    8 #include "voigt.h"
    9 #include <math.h>
   10 
   11 ///     Calculates the Faddeeva function
   12 ///     and partial derivatives of the Voigt function for y>=0
   13 ///     (from http://www.atm.ox.ac.uk/user/wells/voigt.html)
   14 /// arguments:
   15 ///  x, y - Faddeeva/Voigt function arguments
   16 ///  k - voigt              -- output
   17 ///  l - Imaginary part     -- output
   18 ///  dkdx - dVoigt/dx       -- output
   19 ///  dkdy - dVoigt/dy       -- output
   20 void humdev(const float x, const float y,
   21             float &k, float &l, float &dkdx, float &dkdy)
   22 {
   23     static const float c[6] = { 1.0117281f,     -0.75197147f,      0.012557727f,
   24                                 0.010022008f,   -2.4206814e-4f,    5.0084806e-7f };
   25     static const float s[6] = { 1.393237f,       0.23115241f,     -0.15535147f,
   26                                 0.0062183662f,   9.1908299e-5f,   -6.2752596e-7f };
   27     static const float t[6] = { 0.31424038f,     0.94778839f,      1.5976826f,
   28                                 2.2795071f,      3.020637f,        3.8897249f };
   29 
   30     static const float rrtpi = 0.56418958f; // 1/SQRT(pi)
   31     static const double drtpi = 0.5641895835477563f; // 1/SQRT(pi)
   32 
   33     static float a0, b1, c0, c2, d0, d1, d2, e0, e2, e4, f1, f3, f5,
   34                  g0, g2, g4, g6, h0, h2, h4, h6, p0, p2, p4, p6, p8,
   35                  q1, q3, q5, q7, r0, r2, w0, w2, w4, z0, z2, z4, z6, z8,
   36                  mf[6], pf[6], mq[6], mt[6], pq[6], pt[6], xm[6], ym[6],
   37                  xp[6], yp[6];
   38 
   39     static float old_y = -1.f;
   40 
   41     static bool rgb, rgc, rgd;
   42     static float yq, xlima, xlimb, xlimc, xlim4;
   43 
   44     if (y != old_y) {
   45         old_y = y;
   46         rgb = true, rgc = true, rgd = true;
   47         yq = y * y;
   48         xlima = 146.7f - y;
   49         xlimb = 24.f - y;
   50         xlimc = 7.4f - y;
   51         xlim4 = y * 18.1f + 1.65f;
   52     }
   53 
   54     float abx = fabs(x);
   55     float xq = abx * abx;
   56 
   57     if (abx > xlima) {  //  Region A
   58         float d = (float)1. / (xq + yq);
   59         d1 = d * rrtpi;
   60         k = d1 * y;
   61         l = d1 * x;
   62         d1 *= d;
   63         dkdx = -d1 * (y + y) * x;
   64         dkdy = d1 * (xq - yq);
   65     }
   66 
   67     else if (abx > xlimb) { //  Region B
   68         if (rgb) {
   69             rgb = false;
   70             a0 = yq + .5f;
   71             b1 = yq - .5f;
   72             d0 = a0 * a0;
   73             d2 = b1 + b1;
   74             c0 = yq * (1.f - d2) + 1.5f;
   75             c2 = a0 + a0;
   76             r0 = yq * (.25f - yq * (yq + .5f)) + .125f;
   77             r2 = yq * (yq + 5.f) + .25f;
   78         }
   79         float d = 1.f / (d0 + xq * (d2 + xq));
   80         d1 = d * rrtpi;
   81         k = d1 * (a0 + xq) * y;
   82         l = d1 * (b1 + xq) * x;
   83         d1 *= d;
   84         dkdx = d1 * x * y * (c0 - (c2 + xq) * (xq + xq));
   85         dkdy = d1 * (r0 - xq * (r2 - xq * (b1 + xq)));
   86     } else {
   87 
   88         if (abx > xlimc) {  // Region C
   89             if (rgc) {
   90                 rgc = false;
   91                 h0 = yq * (yq * (yq * (yq + (float)6.) + (float)10.5)
   92                         + (float)4.5) + (float).5625;
   93                 h2 = yq * (yq * (yq * (float)4. + (float)6.) + (float)9.)
   94                         - (float)4.5;
   95                 h4 = (float)10.5 - yq * ((float)6. - yq * (float)6.);
   96                 h6 = yq * (float)4. - (float)6.;
   97                 w0 = yq * (yq * (yq * (float)7. + (float)27.5)
   98                         + (float) 24.25) + (float)1.875;
   99                 w2 = yq * (yq * (float)15. + (float)3.) + (float)5.25;
  100                 w4 = yq * (float)9. - (float)4.5;
  101                 f1 = yq * (yq * (yq + (float)4.5) + (float)5.25)
  102                     - (float) 1.875;
  103                 f3 = (float)8.25 - yq * ((float)1. - yq * (float)3.);
  104                 f5 = yq * (float)3. - (float)5.5;
  105                 e0 = y * (yq * (yq * (yq + (float)5.5) + (float)8.25)
  106                         + (float)1.875);
  107                 e2 = y * (yq * (yq * (float)3. + (float)1.) + (float) 5.25);
  108                 e4 = y * (float).75 * h6;
  109                 g0 = y * (yq * (yq * (yq * (float)8. + (float)36.)
  110                             + (float)42.) + (float)9.);
  111                 g2 = y * (yq * (yq * (float)24. + (float)24.) + (float)18.);
  112                 g4 = y * (yq * (float)24. - (float)12.);
  113                 g6 = y * (float)8.;
  114             }
  115             float u = e0 + xq * (e2 + xq * (e4 + xq * y));
  116             float d = (float)1. / (h0 + xq * (h2 + xq * (h4 + xq
  117                                                             * (h6 + xq))));
  118             k = d * rrtpi * u;
  119             l = d * rrtpi * x * (f1 + xq * (f3 + xq * (f5 + xq)));
  120             float dudy = w0 + xq * (w2 + xq * (w4 + xq));
  121             float dvdy = g0 + xq * (g2 + xq * (g4 + xq * g6));
  122             dkdy = d * rrtpi * (dudy - d * u * dvdy);
  123         }
  124 
  125         else if (abx < (float).85) {  // Region D
  126             if (rgd) {
  127                 rgd = false;
  128                 z0 = y * (y * (y * (y * (y * (y * (y * (y * (y *
  129                         (y + (float)13.3988) + (float)88.26741) + (float)
  130                         369.1989) + (float)1074.409) + (float)2256.981) +
  131                         (float)3447.629) + (float)3764.966) + (float)
  132                         2802.87) + (float)1280.829) + (float)272.1014;
  133                 z2 = y * (y * (y * (y * (y * (y * (y * (y * (float)5.
  134                         + (float)53.59518) + (float)266.2987)
  135                         + (float)793.4273) + (float)1549.675) + (float)2037.31)
  136                         + (float)1758.336) + (float)902.3066) + (float)211.678;
  137                 z4 = y * (y * (y * (y * (y * (y * (float)10. + (float)80.39278)
  138                         + (float)269.2916) + (float) 479.2576)
  139                         + (float)497.3014) + (float)308.1852) + (float)78.86585;
  140                 z6 = y * (y * (y * (y * (float)10. + (float)53.59518)
  141                         + (float)92.75679) + (float)55.02933) + (float)22.03523;
  142                 z8 = y * (y * (float)5. + (float)13.3988) + (float)1.49646;
  143                 p0 = y * (y * (y * (y * (y * (y * (y * (y * (y *
  144                         (float).3183291 + (float)4.264678) + (float)
  145                         27.93941) + (float)115.3772) + (float)328.2151) +
  146                         (float)662.8097) + (float)946.897) + (float)
  147                         919.4955) + (float)549.3954) + (float)153.5168;
  148                 p2 = y * (y * (y * (y * (y * (y * (y * (float)
  149                         1.2733163 + (float)12.79458) + (float)56.81652) +
  150                         (float)139.4665) + (float)189.773) + (float)
  151                         124.5975) - (float)1.322256) - (float)34.16955;
  152                 p4 = y * (y * (y * (y * (y * (float)1.9099744 + (
  153                         float)12.79568) + (float)29.81482) + (float)
  154                         24.01655) + (float)10.46332) + (float)2.584042;
  155                 p6 = y * (y * (y * (float)1.273316 + (float)4.266322)
  156                         + (float).9377051) - (float).07272979;
  157                 p8 = y * (float).3183291 + (float)5.480304e-4;
  158                 q1 = y * (y * (y * (y * (y * (y * (y * (y * (
  159                         float).3183291 + (float)4.26413) + (float)27.6294)
  160                          + (float)111.0528) + (float)301.3208) + (float)
  161                         557.5178) + (float)685.8378) + (float)508.2585) +
  162                         (float)173.2355;
  163                 q3 = y * (y * (y * (y * (y * (y * (float)1.273316 +
  164                         (float)12.79239) + (float)55.8865) + (float)
  165                         130.8905) + (float)160.4013) + (float)100.7375) +
  166                         (float)18.97431;
  167                 q5 = y * (y * (y * (y * (float)1.909974 + (float)
  168                         12.79239) + (float)28.8848) + (float)19.83766)
  169                         + (float)7.985877;
  170                 q7 = y * (y * (float)1.273316 + (float)4.26413)
  171                         + (float).6276985;
  172             }
  173             float u = (p0 + xq * (p2 + xq * (p4 + xq * (p6 + xq * p8))))
  174                                                         * 1.7724538f;
  175             float d = 1.f / (z0 + xq * (z2 + xq * (z4 + xq
  176                                                 * (z6 + xq * (z8 + xq)))));
  177             k = d * u;
  178             l = d * 1.7724538f * x * (q1 + xq * (q3 +
  179                     xq * (q5 + xq * (q7 + xq * .3183291f))));
  180             dkdy = 2 * ((double) x * (double) l
  181                                     + (double) y * (double) k - drtpi);
  182         }
  183 
  184         else {     // Use CPF12
  185             float ypy0 = y + 1.5f;
  186             float ypy0q = ypy0 * ypy0;
  187             k = (float)0.;
  188             l = (float)0.;
  189             for (int j = 0; j <= 5; ++j) {
  190                 mt[j] = x - t[j];
  191                 mq[j] = mt[j] * mt[j];
  192                 mf[j] = (float)1. / (mq[j] + ypy0q);
  193                 xm[j] = mf[j] * mt[j];
  194                 ym[j] = mf[j] * ypy0;
  195                 pt[j] = x + t[j];
  196                 pq[j] = pt[j] * pt[j];
  197                 pf[j] = (float)1. / (pq[j] + ypy0q);
  198                 xp[j] = pf[j] * pt[j];
  199                 yp[j] = pf[j] * ypy0;
  200                 l += c[j] * (xm[j] + xp[j]) + s[j] * (ym[j] - yp[j]);
  201             }
  202             if (abx <= xlim4) {  // Humlicek CPF12 Region I
  203                 float yf1 = ypy0 + ypy0;
  204                 float yf2 = ypy0q + ypy0q;
  205                 dkdy = (float)0.;
  206                 for (int j = 0; j <= 5; ++j) {
  207                     float mfq = mf[j] * mf[j];
  208                     float pfq = pf[j] * pf[j];
  209                     k += c[j] * (ym[j] + yp[j]) - s[j] * (xm[j] - xp[j]);
  210                     dkdy += c[j] * (mf[j] + pf[j] - yf2 * (mfq + pfq))
  211                                + s[j] * yf1 * (mt[j] * mfq - pt[j] * pfq);
  212                 }
  213             } else {               //  Humlicek CPF12 Region II
  214                 float yp2y0 = y + (float)3.;
  215                 for (int j = 0; j <= 5; ++j) {
  216                     k += (c[j] * (mq[j] * mf[j] - ym[j] * 1.5f)
  217                              + s[j] * yp2y0 * xm[j]) / (mq[j] + 2.25f)
  218                           + (c[j] * (pq[j] * pf[j] - yp[j] * 1.5f)
  219                             - s[j] * yp2y0 * xp[j]) / (pq[j] + 2.25f);
  220                 }
  221                 k = y * k + exp(-xq);
  222                 dkdy = 2 * ((double) x * (double) l
  223                                   + (double) y * (double) k - drtpi);
  224             }
  225         }
  226         dkdx = 2 * ((double) y * (double) l
  227                                         - (double) x * (double) k);
  228     }
  229 }
  230 
  231 
  232 // ========================================================================
  233 
  234 
  235 
  236 ///   Calculates the Faddeeva function with relative error less than 10^(-4).
  237 ///     (from http://www.atm.ox.ac.uk/user/wells/voigt.html)
  238 /// arguments:
  239 ///  x, y - Faddeeva/Voigt function arguments
  240 /// return value -- voigt
  241 float humlik(const float x, const float y)
  242 {
  243 
  244     static const float c[6] = { 1.0117281f,     -0.75197147f,      0.012557727f,
  245                                 0.010022008f,   -2.4206814e-4f,    5.0084806e-7f };
  246     static const float s[6] = { 1.393237f,       0.23115241f,     -0.15535147f,
  247                                 0.0062183662f,   9.1908299e-5f,   -6.2752596e-7f };
  248     static const float t[6] = { 0.31424038f,     0.94778839f,      1.5976826f,
  249                                 2.2795071f,      3.020637f,        3.8897249f };
  250 
  251     const float rrtpi = 0.56418958f; // 1/SQRT(pi)
  252 
  253     static float a0, d0, d2, e0, e2, e4, h0, h2, h4, h6,
  254                  p0, p2, p4, p6, p8, z0, z2, z4, z6, z8;
  255     static float mf[6], pf[6], mq[6], pq[6], xm[6], ym[6], xp[6], yp[6];
  256     static float old_y = -1.f;
  257     static bool rg1, rg2, rg3;
  258     static float xlim0, xlim1, xlim2, xlim3, xlim4;
  259     static float yq, yrrtpi;
  260     if (y != old_y) {
  261         old_y = y;
  262         yq = y * y;
  263         yrrtpi = y * rrtpi;
  264         rg1 = true, rg2 = true, rg3 = true;
  265         if (y < 70.55) {
  266             xlim0 = sqrt(y * (40.f - y * 3.6f) + 15100.f);
  267             xlim1 = (y >= 8.425 ?  0.f : sqrt(164.f - y * (y * 1.8f + 4.3f)));
  268             xlim2 = 6.8f - y;
  269             xlim3 = y * 2.4f;
  270             xlim4 = y * 18.1f + 1.65f;
  271             if (y <= 1e-6)
  272                 xlim2 = xlim1 = xlim0;
  273         }
  274     }
  275 
  276     float abx = fabs(x);
  277     float xq = abx * abx;
  278 
  279     if (abx >= xlim0 || y >= 70.55)         // Region 0 algorithm
  280         return yrrtpi / (xq + yq);
  281 
  282     else if (abx >= xlim1) {            //  Humlicek W4 Region 1
  283         if (rg1) {  // First point in Region 1
  284             rg1 = false;
  285             a0 = yq + 0.5f;  //Region 1 y-dependents
  286             d0 = a0 * a0;
  287             d2 = yq + yq - 1.f;
  288         }
  289         return rrtpi / (d0 + xq * (d2 + xq)) * y * (a0 + xq);
  290     }
  291 
  292     else if (abx > xlim2) {  // Humlicek W4 Region 2
  293         if (rg2) {  //First point in Region 2
  294             rg2 = false;
  295             // Region 2 y-dependents
  296             h0 = yq * (yq * (yq * (yq + 6.f) + 10.5f) + 4.5f) + 0.5625f;
  297             h2 = yq * (yq * (yq * 4.f + 6.f) + 9.f) - 4.5f;
  298             h4 = 10.5f - yq * (6.f - yq * 6.f);
  299             h6 = yq * 4.f - 6.f;
  300             e0 = yq * (yq * (yq + 5.5f) + 8.25f) + 1.875f;
  301             e2 = yq * (yq * 3.f + 1.f) + 5.25f;
  302             e4 = h6 * 0.75f;
  303         }
  304         return rrtpi / (h0 + xq * (h2 + xq * (h4 + xq * (h6 + xq))))
  305                  * y * (e0 + xq * (e2 + xq * (e4 + xq)));
  306     }
  307 
  308     else if (abx < xlim3) { // Humlicek W4 Region 3
  309         if (rg3) {  // First point in Region 3
  310             rg3 = false;
  311             //Region 3 y-dependents
  312             z0 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y
  313                     + 13.3988f) + 88.26741f) + 369.1989f) + 1074.409f)
  314                     + 2256.981f) + 3447.629f) + 3764.966f) + 2802.87f)
  315                     + 1280.829f) + 272.1014f;
  316             z2 = y * (y * (y * (y * (y * (y * (y * (y * 5.f  + 53.59518f)
  317                     + 266.2987f) + 793.4273f) + 1549.675f) + 2037.31f)
  318                     + 1758.336f) + 902.3066f) + 211.678f;
  319             z4 = y * (y * (y * (y * (y * (y * 10.f + 80.39278f) + 269.2916f)
  320                     + 479.2576f) + 497.3014f) + 308.1852f) + 78.86585f;
  321             z6 = y * (y * (y * (y * 10.f + 53.59518f) + 92.75679f)
  322                     + 55.02933f) + 22.03523f;
  323             z8 = y * (y * 5.f + 13.3988f) + 1.49646f;
  324             p0 = y * (y * (y * (y * (y * (y * (y * (y * (y * 0.3183291f
  325                     + 4.264678f) + 27.93941f) + 115.3772f) + 328.2151f) +
  326                     662.8097f) + 946.897f) + 919.4955f) + 549.3954f)
  327                     + 153.5168f;
  328             p2 = y * (y * (y * (y * (y * (y * (y * 1.2733163f + 12.79458f)
  329                     + 56.81652f) + 139.4665f) + 189.773f) + 124.5975f)
  330                     - 1.322256f) - 34.16955f;
  331             p4 = y * (y * (y * (y * (y * 1.9099744f + 12.79568f)
  332                     + 29.81482f) + 24.01655f) + 10.46332f) + 2.584042f;
  333             p6 = y * (y * (y * 1.273316f + 4.266322f) + 0.9377051f)
  334                     - 0.07272979f;
  335             p8 = y * .3183291f + 5.480304e-4f;
  336         }
  337         return 1.7724538f / (z0 + xq * (z2 + xq * (z4 + xq * (z6 +
  338                 xq * (z8 + xq)))))
  339                   * (p0 + xq * (p2 + xq * (p4 + xq * (p6 + xq * p8))));
  340     }
  341 
  342     else {  //  Humlicek CPF12 algorithm
  343         float ypy0 = y + 1.5f;
  344         float ypy0q = ypy0 * ypy0;
  345         for (int j = 0; j <= 5; ++j) {
  346             float d = x - t[j];
  347             mq[j] = d * d;
  348             mf[j] = 1.f / (mq[j] + ypy0q);
  349             xm[j] = mf[j] * d;
  350             ym[j] = mf[j] * ypy0;
  351             d = x + t[j];
  352             pq[j] = d * d;
  353             pf[j] = 1.f / (pq[j] + ypy0q);
  354             xp[j] = pf[j] * d;
  355             yp[j] = pf[j] * ypy0;
  356         }
  357         float k = 0.f;
  358         if (abx <= xlim4) // Humlicek CPF12 Region I
  359             for (int j = 0; j <= 5; ++j)
  360                 k += c[j] * (ym[j]+yp[j]) - s[j] * (xm[j]-xp[j]);
  361         else {           // Humlicek CPF12 Region II
  362             float yf = y + 3.f;
  363             for (int j = 0; j <= 5; ++j)
  364                 k += (c[j] * (mq[j] * mf[j] - ym[j] * 1.5f)
  365                          + s[j] * yf * xm[j]) / (mq[j] + 2.25f)
  366                         + (c[j] * (pq[j] * pf[j] - yp[j] * 1.5f)
  367                            - s[j] * yf * xp[j]) / (pq[j] + 2.25f);
  368             k = y * k + exp(-xq);
  369         }
  370         return k;
  371     }
  372 }
  373 
  374