fityk: /home/wojdyr/fityk/src/voigt.cpp Source File

00001 
00002 // fastest_humlik.for and humdev.for - from Bob Wells Voigt Function Page
00003 // http://www.atm.ox.ac.uk/user/wells/voigt.html
00004 // Translated to C++ with f2c program and modified by M.W.
00005 // It can be slower than original, I haven't compared the speed.
00006 
00007 #include <math.h>
00008 #include "voigt.h"
00009 
00010 ///     Calculates the Faddeeva function
00011 ///     and partial derivatives of the Voigt function for y>=0
00012 ///     (from http://www.atm.ox.ac.uk/user/wells/voigt.html)
00013 /// arguments:
00014 ///  x, y - Faddeeva/Voigt function arguments
00015 ///  k - voigt              -- output
00016 ///  l - Imaginary part     -- output
00017 ///  dkdx - dVoigt/dx       -- output
00018 ///  dkdy - dVoigt/dy       -- output
00019 void humdev(const float x, const float y,
00020             float &k, float &l, float &dkdx, float &dkdy)
00021 {
00022     static const float c[6] = { 1.0117281f,     -0.75197147f,      0.012557727f,
00023                                 0.010022008f,   -2.4206814e-4f,    5.0084806e-7f };
00024     static const float s[6] = { 1.393237f,       0.23115241f,     -0.15535147f,
00025                                 0.0062183662f,   9.1908299e-5f,   -6.2752596e-7f };
00026     static const float t[6] = { 0.31424038f,     0.94778839f,      1.5976826f,
00027                                 2.2795071f,      3.020637f,        3.8897249f };
00028 
00029     static const float rrtpi = 0.56418958f; // 1/SQRT(pi)
00030     static const double drtpi = 0.5641895835477563f; // 1/SQRT(pi)
00031 
00032     static float a0, b1, c0, c2, d0, d1, d2, e0, e2, e4, f1, f3, f5,
00033                  g0, g2, g4, g6, h0, h2, h4, h6, p0, p2, p4, p6, p8,
00034                  q1, q3, q5, q7, r0, r2, w0, w2, w4, z0, z2, z4, z6, z8,
00035                  mf[6], pf[6], mq[6], mt[6], pq[6], pt[6], xm[6], ym[6],
00036                  xp[6], yp[6];
00037 
00038     static float old_y = -1.f;
00039 
00040     static bool rgb, rgc, rgd;
00041     static float yq, xlima, xlimb, xlimc, xlim4;
00042 
00043     if (y != old_y) {
00044         old_y = y;
00045         rgb = true, rgc = true, rgd = true;
00046         yq = y * y;
00047         xlima = 146.7f - y;
00048         xlimb = 24.f - y;
00049         xlimc = 7.4f - y;
00050         xlim4 = y * 18.1f + 1.65f;
00051     }
00052 
00053     float abx = fabs(x);
00054     float xq = abx * abx;
00055 
00056     if (abx > xlima) {  //  Region A
00057         float d = (float)1. / (xq + yq);
00058         d1 = d * rrtpi;
00059         k = d1 * y;
00060         l = d1 * x;
00061         d1 *= d;
00062         dkdx = -d1 * (y + y) * x;
00063         dkdy = d1 * (xq - yq);
00064     }
00065 
00066     else if (abx > xlimb) { //  Region B
00067         if (rgb) {
00068             rgb = false;
00069             a0 = yq + .5f;
00070             b1 = yq - .5f;
00071             d0 = a0 * a0;
00072             d2 = b1 + b1;
00073             c0 = yq * (1.f - d2) + 1.5f;
00074             c2 = a0 + a0;
00075             r0 = yq * (.25f - yq * (yq + .5f)) + .125f;
00076             r2 = yq * (yq + 5.f) + .25f;
00077         }
00078         float d = 1.f / (d0 + xq * (d2 + xq));
00079         d1 = d * rrtpi;
00080         k = d1 * (a0 + xq) * y;
00081         l = d1 * (b1 + xq) * x;
00082         d1 *= d;
00083         dkdx = d1 * x * y * (c0 - (c2 + xq) * (xq + xq));
00084         dkdy = d1 * (r0 - xq * (r2 - xq * (b1 + xq)));
00085     }
00086     else {
00087 
00088         if (abx > xlimc) {  // Region C
00089             if (rgc) {
00090                 rgc = false;
00091                 h0 = yq * (yq * (yq * (yq + (float)6.) + (float)10.5)
00092                         + (float)4.5) + (float).5625;
00093                 h2 = yq * (yq * (yq * (float)4. + (float)6.) + (float)9.)
00094                         - (float)4.5;
00095                 h4 = (float)10.5 - yq * ((float)6. - yq * (float)6.);
00096                 h6 = yq * (float)4. - (float)6.;
00097                 w0 = yq * (yq * (yq * (float)7. + (float)27.5)
00098                         + (float) 24.25) + (float)1.875;
00099                 w2 = yq * (yq * (float)15. + (float)3.) + (float)5.25;
00100                 w4 = yq * (float)9. - (float)4.5;
00101                 f1 = yq * (yq * (yq + (float)4.5) + (float)5.25)
00102                     - (float) 1.875;
00103                 f3 = (float)8.25 - yq * ((float)1. - yq * (float)3.);
00104                 f5 = yq * (float)3. - (float)5.5;
00105                 e0 = y * (yq * (yq * (yq + (float)5.5) + (float)8.25)
00106                         + (float)1.875);
00107                 e2 = y * (yq * (yq * (float)3. + (float)1.) + (float) 5.25);
00108                 e4 = y * (float).75 * h6;
00109                 g0 = y * (yq * (yq * (yq * (float)8. + (float)36.)
00110                             + (float)42.) + (float)9.);
00111                 g2 = y * (yq * (yq * (float)24. + (float)24.) + (float)18.);
00112                 g4 = y * (yq * (float)24. - (float)12.);
00113                 g6 = y * (float)8.;
00114             }
00115             float u = e0 + xq * (e2 + xq * (e4 + xq * y));
00116             float d = (float)1. / (h0 + xq * (h2 + xq * (h4 + xq
00117                                                             * (h6 + xq))));
00118             k = d * rrtpi * u;
00119             l = d * rrtpi * x * (f1 + xq * (f3 + xq * (f5 + xq)));
00120             float dudy = w0 + xq * (w2 + xq * (w4 + xq));
00121             float dvdy = g0 + xq * (g2 + xq * (g4 + xq * g6));
00122             dkdy = d * rrtpi * (dudy - d * u * dvdy);
00123         }
00124 
00125         else if (abx < (float).85) {  // Region D
00126             if (rgd) {
00127                 rgd = false;
00128                 z0 = y * (y * (y * (y * (y * (y * (y * (y * (y *
00129                         (y + (float)13.3988) + (float)88.26741) + (float)
00130                         369.1989) + (float)1074.409) + (float)2256.981) +
00131                         (float)3447.629) + (float)3764.966) + (float)
00132                         2802.87) + (float)1280.829) + (float)272.1014;
00133                 z2 = y * (y * (y * (y * (y * (y * (y * (y * (float)5.
00134                         + (float)53.59518) + (float)266.2987)
00135                         + (float)793.4273) + (float)1549.675) + (float)2037.31)
00136                         + (float)1758.336) + (float)902.3066) + (float)211.678;
00137                 z4 = y * (y * (y * (y * (y * (y * (float)10. + (float)80.39278)
00138                         + (float)269.2916) + (float) 479.2576)
00139                         + (float)497.3014) + (float)308.1852) + (float)78.86585;
00140                 z6 = y * (y * (y * (y * (float)10. + (float)53.59518)
00141                         + (float)92.75679) + (float)55.02933) + (float)22.03523;
00142                 z8 = y * (y * (float)5. + (float)13.3988) + (float)1.49646;
00143                 p0 = y * (y * (y * (y * (y * (y * (y * (y * (y *
00144                         (float).3183291 + (float)4.264678) + (float)
00145                         27.93941) + (float)115.3772) + (float)328.2151) +
00146                         (float)662.8097) + (float)946.897) + (float)
00147                         919.4955) + (float)549.3954) + (float)153.5168;
00148                 p2 = y * (y * (y * (y * (y * (y * (y * (float)
00149                         1.2733163 + (float)12.79458) + (float)56.81652) +
00150                         (float)139.4665) + (float)189.773) + (float)
00151                         124.5975) - (float)1.322256) - (float)34.16955;
00152                 p4 = y * (y * (y * (y * (y * (float)1.9099744 + (
00153                         float)12.79568) + (float)29.81482) + (float)
00154                         24.01655) + (float)10.46332) + (float)2.584042;
00155                 p6 = y * (y * (y * (float)1.273316 + (float)4.266322)
00156                         + (float).9377051) - (float).07272979;
00157                 p8 = y * (float).3183291 + (float)5.480304e-4;
00158                 q1 = y * (y * (y * (y * (y * (y * (y * (y * (
00159                         float).3183291 + (float)4.26413) + (float)27.6294)
00160                          + (float)111.0528) + (float)301.3208) + (float)
00161                         557.5178) + (float)685.8378) + (float)508.2585) +
00162                         (float)173.2355;
00163                 q3 = y * (y * (y * (y * (y * (y * (float)1.273316 +
00164                         (float)12.79239) + (float)55.8865) + (float)
00165                         130.8905) + (float)160.4013) + (float)100.7375) +
00166                         (float)18.97431;
00167                 q5 = y * (y * (y * (y * (float)1.909974 + (float)
00168                         12.79239) + (float)28.8848) + (float)19.83766)
00169                         + (float)7.985877;
00170                 q7 = y * (y * (float)1.273316 + (float)4.26413)
00171                         + (float).6276985;
00172             }
00173             float u = (p0 + xq * (p2 + xq * (p4 + xq * (p6 + xq * p8))))
00174                                                         * 1.7724538f;
00175             float d = 1.f / (z0 + xq * (z2 + xq * (z4 + xq
00176                                                 * (z6 + xq * (z8 + xq)))));
00177             k = d * u;
00178             l = d * 1.7724538f * x * (q1 + xq * (q3 +
00179                     xq * (q5 + xq * (q7 + xq * .3183291f))));
00180             dkdy = 2 * ((double) x * (double) l
00181                                     + (double) y * (double) k - drtpi);
00182         }
00183 
00184         else {     // Use CPF12
00185             float ypy0 = y + 1.5f;
00186             float ypy0q = ypy0 * ypy0;
00187             k = (float)0.;
00188             l = (float)0.;
00189             for (int j = 0; j <= 5; ++j) {
00190                 mt[j] = x - t[j];
00191                 mq[j] = mt[j] * mt[j];
00192                 mf[j] = (float)1. / (mq[j] + ypy0q);
00193                 xm[j] = mf[j] * mt[j];
00194                 ym[j] = mf[j] * ypy0;
00195                 pt[j] = x + t[j];
00196                 pq[j] = pt[j] * pt[j];
00197                 pf[j] = (float)1. / (pq[j] + ypy0q);
00198                 xp[j] = pf[j] * pt[j];
00199                 yp[j] = pf[j] * ypy0;
00200                 l += c[j] * (xm[j] + xp[j]) + s[j] * (ym[j] - yp[j]);
00201             }
00202             if (abx <= xlim4) {  // Humlicek CPF12 Region I
00203                 float yf1 = ypy0 + ypy0;
00204                 float yf2 = ypy0q + ypy0q;
00205                 dkdy = (float)0.;
00206                 for (int j = 0; j <= 5; ++j) {
00207                     float mfq = mf[j] * mf[j];
00208                     float pfq = pf[j] * pf[j];
00209                     k += c[j] * (ym[j] + yp[j]) - s[j] * (xm[j] - xp[j]);
00210                     dkdy += c[j] * (mf[j] + pf[j] - yf2 * (mfq + pfq))
00211                                + s[j] * yf1 * (mt[j] * mfq - pt[j] * pfq);
00212                 }
00213             }
00214             else {               //  Humlicek CPF12 Region II
00215                 float yp2y0 = y + (float)3.;
00216                 for (int j = 0; j <= 5; ++j) {
00217                     k += (c[j] * (mq[j] * mf[j] - ym[j] * 1.5f)
00218                              + s[j] * yp2y0 * xm[j]) / (mq[j] + 2.25f)
00219                           + (c[j] * (pq[j] * pf[j] - yp[j] * 1.5f)
00220                             - s[j] * yp2y0 * xp[j]) / (pq[j] + 2.25f);
00221                 }
00222                 k = y * k + exp(-xq);
00223                 dkdy = 2 * ((double) x * (double) l
00224                                   + (double) y * (double) k - drtpi);
00225             }
00226         }
00227         dkdx = 2 * ((double) y * (double) l
00228                                         - (double) x * (double) k);
00229     }
00230 }
00231 
00232 
00233 // ========================================================================
00234 
00235 
00236 
00237 ///   Calculates the Faddeeva function with relative error less than 10^(-4).
00238 ///     (from http://www.atm.ox.ac.uk/user/wells/voigt.html)
00239 /// arguments:
00240 ///  x, y - Faddeeva/Voigt function arguments
00241 /// return value -- voigt
00242 float humlik(const float x, const float y)
00243 {
00244 
00245     static const float c[6] = { 1.0117281f,     -0.75197147f,      0.012557727f,
00246                                 0.010022008f,   -2.4206814e-4f,    5.0084806e-7f };
00247     static const float s[6] = { 1.393237f,       0.23115241f,     -0.15535147f,
00248                                 0.0062183662f,   9.1908299e-5f,   -6.2752596e-7f };
00249     static const float t[6] = { 0.31424038f,     0.94778839f,      1.5976826f,
00250                                 2.2795071f,      3.020637f,        3.8897249f };
00251 
00252     const float rrtpi = 0.56418958f; // 1/SQRT(pi)
00253 
00254     static float a0, d0, d2, e0, e2, e4, h0, h2, h4, h6,
00255                  p0, p2, p4, p6, p8, z0, z2, z4, z6, z8;
00256     static float mf[6], pf[6], mq[6], pq[6], xm[6], ym[6], xp[6], yp[6];
00257     static float old_y = -1.f;
00258     static bool rg1, rg2, rg3;
00259     static float xlim0, xlim1, xlim2, xlim3, xlim4;
00260     static float yq, yrrtpi;
00261     if (y != old_y) {
00262         old_y = y;
00263         yq = y * y;
00264         yrrtpi = y * rrtpi;
00265         rg1 = true, rg2 = true, rg3 = true;
00266         if (y < 70.55) {
00267             xlim0 = sqrt(y * (40.f - y * 3.6f) + 15100.f);
00268             xlim1 = (y >= 8.425 ?  0.f : sqrt(164.f - y * (y * 1.8f + 4.3f)));
00269             xlim2 = 6.8f - y;
00270             xlim3 = y * 2.4f;
00271             xlim4 = y * 18.1f + 1.65f;
00272             if (y <= 1e-6)
00273                 xlim2 = xlim1 = xlim0;
00274         }
00275     }
00276 
00277     float abx = fabs(x);
00278     float xq = abx * abx;
00279 
00280     if (abx >= xlim0 || y >= 70.55)         // Region 0 algorithm
00281         return yrrtpi / (xq + yq);
00282 
00283     else if (abx >= xlim1) {            //  Humlicek W4 Region 1
00284         if (rg1) {  // First point in Region 1
00285             rg1 = false;
00286             a0 = yq + 0.5f;  //Region 1 y-dependents
00287             d0 = a0 * a0;
00288             d2 = yq + yq - 1.f;
00289         }
00290         return rrtpi / (d0 + xq * (d2 + xq)) * y * (a0 + xq);
00291     }
00292 
00293     else if (abx > xlim2) {  // Humlicek W4 Region 2
00294         if (rg2) {  //First point in Region 2
00295             rg2 = false;
00296             // Region 2 y-dependents
00297             h0 = yq * (yq * (yq * (yq + 6.f) + 10.5f) + 4.5f) + 0.5625f;
00298             h2 = yq * (yq * (yq * 4.f + 6.f) + 9.f) - 4.5f;
00299             h4 = 10.5f - yq * (6.f - yq * 6.f);
00300             h6 = yq * 4.f - 6.f;
00301             e0 = yq * (yq * (yq + 5.5f) + 8.25f) + 1.875f;
00302             e2 = yq * (yq * 3.f + 1.f) + 5.25f;
00303             e4 = h6 * 0.75f;
00304         }
00305         return rrtpi / (h0 + xq * (h2 + xq * (h4 + xq * (h6 + xq))))
00306                  * y * (e0 + xq * (e2 + xq * (e4 + xq)));
00307     }
00308 
00309     else if (abx < xlim3) { // Humlicek W4 Region 3
00310         if (rg3) {  // First point in Region 3
00311             rg3 = false;
00312             //Region 3 y-dependents
00313             z0 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y
00314                     + 13.3988f) + 88.26741f) + 369.1989f) + 1074.409f)
00315                     + 2256.981f) + 3447.629f) + 3764.966f) + 2802.87f)
00316                     + 1280.829f) + 272.1014f;
00317             z2 = y * (y * (y * (y * (y * (y * (y * (y * 5.f  + 53.59518f)
00318                     + 266.2987f) + 793.4273f) + 1549.675f) + 2037.31f)
00319                     + 1758.336f) + 902.3066f) + 211.678f;
00320             z4 = y * (y * (y * (y * (y * (y * 10.f + 80.39278f) + 269.2916f)
00321                     + 479.2576f) + 497.3014f) + 308.1852f) + 78.86585f;
00322             z6 = y * (y * (y * (y * 10.f + 53.59518f) + 92.75679f)
00323                     + 55.02933f) + 22.03523f;
00324             z8 = y * (y * 5.f + 13.3988f) + 1.49646f;
00325             p0 = y * (y * (y * (y * (y * (y * (y * (y * (y * 0.3183291f
00326                     + 4.264678f) + 27.93941f) + 115.3772f) + 328.2151f) +
00327                     662.8097f) + 946.897f) + 919.4955f) + 549.3954f)
00328                     + 153.5168f;
00329             p2 = y * (y * (y * (y * (y * (y * (y * 1.2733163f + 12.79458f)
00330                     + 56.81652f) + 139.4665f) + 189.773f) + 124.5975f)
00331                     - 1.322256f) - 34.16955f;
00332             p4 = y * (y * (y * (y * (y * 1.9099744f + 12.79568f)
00333                     + 29.81482f) + 24.01655f) + 10.46332f) + 2.584042f;
00334             p6 = y * (y * (y * 1.273316f + 4.266322f) + 0.9377051f)
00335                     - 0.07272979f;
00336             p8 = y * .3183291f + 5.480304e-4f;
00337         }
00338         return 1.7724538f / (z0 + xq * (z2 + xq * (z4 + xq * (z6 +
00339                 xq * (z8 + xq)))))
00340                   * (p0 + xq * (p2 + xq * (p4 + xq * (p6 + xq * p8))));
00341     }
00342 
00343     else {  //  Humlicek CPF12 algorithm
00344         float ypy0 = y + 1.5f;
00345         float ypy0q = ypy0 * ypy0;
00346         for (int j = 0; j <= 5; ++j) {
00347             float d = x - t[j];
00348             mq[j] = d * d;
00349             mf[j] = 1.f / (mq[j] + ypy0q);
00350             xm[j] = mf[j] * d;
00351             ym[j] = mf[j] * ypy0;
00352             d = x + t[j];
00353             pq[j] = d * d;
00354             pf[j] = 1.f / (pq[j] + ypy0q);
00355             xp[j] = pf[j] * d;
00356             yp[j] = pf[j] * ypy0;
00357         }
00358         float k = 0.f;
00359         if (abx <= xlim4) // Humlicek CPF12 Region I
00360             for (int j = 0; j <= 5; ++j)
00361                 k += c[j] * (ym[j]+yp[j]) - s[j] * (xm[j]-xp[j]);
00362         else {           // Humlicek CPF12 Region II
00363             float yf = y + 3.f;
00364             for (int j = 0; j <= 5; ++j)
00365                 k += (c[j] * (mq[j] * mf[j] - ym[j] * 1.5f)
00366                          + s[j] * yf * xm[j]) / (mq[j] + 2.25f)
00367                         + (c[j] * (pq[j] * pf[j] - yp[j] * 1.5f)
00368                            - s[j] * yf * xp[j]) / (pq[j] + 2.25f);
00369             k = y * k + exp(-xq);
00370         }
00371         return k;
00372     }
00373 }
00374 
00375