"Fossies" - the Fresh Open Source Software Archive

Member "fityk-1.3.1/fityk/bfunc.cpp" (13 May 2016, 33693 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 "bfunc.cpp" see the Fossies "Dox" file reference documentation.

    1 // This file is part of fityk program. Copyright 2001-2013 Marcin Wojdyr
    2 // Licence: GNU General Public License ver. 2+
    3 /// Built-in function definitions
    4 
    5 #define BUILDING_LIBFITYK
    6 #include "bfunc.h"
    7 
    8 #include <boost/math/special_functions/gamma.hpp>
    9 
   10 #include "voigt.h"
   11 #include "numfuncs.h"
   12 
   13 using namespace std;
   14 using boost::math::lgamma;
   15 
   16 namespace fityk {
   17 
   18 
   19 void FuncConstant::calculate_value_in_range(vector<realt> const&/*xx*/,
   20                                             vector<realt>& yy,
   21                                             int first, int last) const
   22 {
   23     for (int i = first; i < last; ++i)
   24         yy[i] += av_[0];
   25 }
   26 
   27 CALCULATE_DERIV_BEGIN(FuncConstant)
   28     (void) x;
   29     dy_dv[0] = 1.;
   30     dy_dx = 0;
   31 CALCULATE_DERIV_END(av_[0])
   32 
   33 ///////////////////////////////////////////////////////////////////////
   34 
   35 CALCULATE_VALUE_BEGIN(FuncLinear)
   36 CALCULATE_VALUE_END(av_[0] + x*av_[1])
   37 
   38 CALCULATE_DERIV_BEGIN(FuncLinear)
   39     dy_dv[0] = 1.;
   40     dy_dv[1] = x;
   41     dy_dx = av_[1];
   42 CALCULATE_DERIV_END(av_[0] + x*av_[1])
   43 
   44 ///////////////////////////////////////////////////////////////////////
   45 
   46 CALCULATE_VALUE_BEGIN(FuncQuadratic)
   47 CALCULATE_VALUE_END(av_[0] + x*av_[1] + x*x*av_[2])
   48 
   49 CALCULATE_DERIV_BEGIN(FuncQuadratic)
   50     dy_dv[0] = 1.;
   51     dy_dv[1] = x;
   52     dy_dv[2] = x*x;
   53     dy_dx = av_[1] + 2*x*av_[2];
   54 CALCULATE_DERIV_END(av_[0] + x*av_[1] + x*x*av_[2])
   55 
   56 ///////////////////////////////////////////////////////////////////////
   57 
   58 CALCULATE_VALUE_BEGIN(FuncCubic)
   59 CALCULATE_VALUE_END(av_[0] + x*av_[1] + x*x*av_[2] + x*x*x*av_[3])
   60 
   61 CALCULATE_DERIV_BEGIN(FuncCubic)
   62     dy_dv[0] = 1.;
   63     dy_dv[1] = x;
   64     dy_dv[2] = x*x;
   65     dy_dv[3] = x*x*x;
   66     dy_dx = av_[1] + 2*x*av_[2] + 3*x*x*av_[3];
   67 CALCULATE_DERIV_END(av_[0] + x*av_[1] + x*x*av_[2] + x*x*x*av_[3])
   68 
   69 ///////////////////////////////////////////////////////////////////////
   70 
   71 CALCULATE_VALUE_BEGIN(FuncPolynomial4)
   72 CALCULATE_VALUE_END(av_[0] + x*av_[1] + x*x*av_[2] + x*x*x*av_[3]
   73                                           + x*x*x*x*av_[4])
   74 
   75 CALCULATE_DERIV_BEGIN(FuncPolynomial4)
   76     dy_dv[0] = 1.;
   77     dy_dv[1] = x;
   78     dy_dv[2] = x*x;
   79     dy_dv[3] = x*x*x;
   80     dy_dv[4] = x*x*x*x;
   81     dy_dx = av_[1] + 2*x*av_[2] + 3*x*x*av_[3] + 4*x*x*x*av_[4];
   82 CALCULATE_DERIV_END(av_[0] + x*av_[1] + x*x*av_[2] + x*x*x*av_[3]
   83                                       + x*x*x*x*av_[4])
   84 
   85 ///////////////////////////////////////////////////////////////////////
   86 
   87 CALCULATE_VALUE_BEGIN(FuncPolynomial5)
   88 CALCULATE_VALUE_END(av_[0] + x*av_[1] + x*x*av_[2]
   89                        + x*x*x*av_[3] + x*x*x*x*av_[4] + x*x*x*x*x*av_[5])
   90 
   91 CALCULATE_DERIV_BEGIN(FuncPolynomial5)
   92     dy_dv[0] = 1.;
   93     dy_dv[1] = x;
   94     dy_dv[2] = x*x;
   95     dy_dv[3] = x*x*x;
   96     dy_dv[4] = x*x*x*x;
   97     dy_dv[5] = x*x*x*x*x;
   98     dy_dx = av_[1] + 2*x*av_[2] + 3*x*x*av_[3] + 4*x*x*x*av_[4]
   99                + 5*x*x*x*x*av_[5];
  100 CALCULATE_DERIV_END(av_[0] + x*av_[1] + x*x*av_[2]
  101                           + x*x*x*av_[3] + x*x*x*x*av_[4] + x*x*x*x*x*av_[5])
  102 
  103 ///////////////////////////////////////////////////////////////////////
  104 
  105 CALCULATE_VALUE_BEGIN(FuncPolynomial6)
  106 CALCULATE_VALUE_END(av_[0] + x*av_[1] + x*x*av_[2] + x*x*x*av_[3] +
  107                         x*x*x*x*av_[4] + x*x*x*x*x*av_[5] + x*x*x*x*x*x*av_[6])
  108 
  109 CALCULATE_DERIV_BEGIN(FuncPolynomial6)
  110     dy_dv[0] = 1.;
  111     dy_dv[1] = x;
  112     dy_dv[2] = x*x;
  113     dy_dv[3] = x*x*x;
  114     dy_dv[4] = x*x*x*x;
  115     dy_dv[5] = x*x*x*x*x;
  116     dy_dv[6] = x*x*x*x*x*x;
  117     dy_dx = av_[1] + 2*x*av_[2] + 3*x*x*av_[3] + 4*x*x*x*av_[4]
  118                 + 5*x*x*x*x*av_[5] + 6*x*x*x*x*x*av_[6];
  119 CALCULATE_DERIV_END(av_[0] + x*av_[1] + x*x*av_[2] + x*x*x*av_[3] +
  120                         x*x*x*x*av_[4] + x*x*x*x*x*av_[5] + x*x*x*x*x*x*av_[6])
  121 
  122 ///////////////////////////////////////////////////////////////////////
  123 
  124 void FuncGaussian::more_precomputations()
  125 {
  126     if (fabs(av_[2]) < epsilon)
  127         av_[2] = epsilon;
  128 }
  129 
  130 CALCULATE_VALUE_BEGIN(FuncGaussian)
  131     realt xa1a2 = (x - av_[1]) / av_[2];
  132     realt ex = exp(- M_LN2 * xa1a2 * xa1a2);
  133 CALCULATE_VALUE_END(av_[0] * ex)
  134 
  135 CALCULATE_DERIV_BEGIN(FuncGaussian)
  136     realt xa1a2 = (x - av_[1]) / av_[2];
  137     realt ex = exp(- M_LN2 * xa1a2 * xa1a2);
  138     dy_dv[0] = ex;
  139     realt dcenter = 2 * M_LN2 * av_[0] * ex * xa1a2 / av_[2];
  140     dy_dv[1] = dcenter;
  141     dy_dv[2] = dcenter * xa1a2;
  142     dy_dx = -dcenter;
  143 CALCULATE_DERIV_END(av_[0]*ex)
  144 
  145 bool FuncGaussian::get_nonzero_range(double level,
  146                                      realt &left, realt &right) const
  147 {
  148     if (level == 0)
  149         return false;
  150     else if (fabs(level) >= fabs(av_[0]))
  151         left = right = 0;
  152     else {
  153         realt w = sqrt(log(fabs(av_[0]/level)) / M_LN2) * av_[2];
  154         left = av_[1] - w;
  155         right = av_[1] + w;
  156     }
  157     return true;
  158 }
  159 
  160 bool FuncGaussian::get_area(realt* a) const
  161 {
  162     *a = av_[0] * fabs(av_[2]) * sqrt(M_PI / M_LN2);
  163     return true;
  164 }
  165 
  166 ///////////////////////////////////////////////////////////////////////
  167 
  168 void FuncSplitGaussian::more_precomputations()
  169 {
  170     if (fabs(av_[2]) < epsilon)
  171         av_[2] = epsilon;
  172     if (fabs(av_[3]) < epsilon)
  173         av_[3] = epsilon;
  174 }
  175 
  176 CALCULATE_VALUE_BEGIN(FuncSplitGaussian)
  177     realt hwhm = (x < av_[1] ? av_[2] : av_[3]);
  178     realt xa1a2 = (x - av_[1]) / hwhm;
  179     realt ex = exp(- M_LN2 * xa1a2 * xa1a2);
  180 CALCULATE_VALUE_END(av_[0] * ex)
  181 
  182 CALCULATE_DERIV_BEGIN(FuncSplitGaussian)
  183     realt hwhm = (x < av_[1] ? av_[2] : av_[3]);
  184     realt xa1a2 = (x - av_[1]) / hwhm;
  185     realt ex = exp(- M_LN2 * xa1a2 * xa1a2);
  186     dy_dv[0] = ex;
  187     realt dcenter = 2 * M_LN2 * av_[0] * ex * xa1a2 / hwhm;
  188     dy_dv[1] = dcenter;
  189     if (x < av_[1]) {
  190         dy_dv[2] = dcenter * xa1a2;
  191         dy_dv[3] = 0;
  192     } else {
  193         dy_dv[2] = 0;
  194         dy_dv[3] = dcenter * xa1a2;
  195     }
  196     dy_dx = -dcenter;
  197 CALCULATE_DERIV_END(av_[0]*ex)
  198 
  199 bool FuncSplitGaussian::get_nonzero_range(double level,
  200                                           realt &left, realt &right) const
  201 {
  202     if (level == 0)
  203         return false;
  204     else if (fabs(level) >= fabs(av_[0]))
  205         left = right = 0;
  206     else {
  207         realt w1 = sqrt (log (fabs(av_[0]/level)) / M_LN2) * av_[2];
  208         realt w2 = sqrt (log (fabs(av_[0]/level)) / M_LN2) * av_[3];
  209         left = av_[1] - w1;
  210         right = av_[1] + w2;
  211     }
  212     return true;
  213 }
  214 
  215 bool FuncSplitGaussian::get_area(realt* a) const
  216 {
  217     *a = av_[0] * (fabs(av_[2])+fabs(av_[3])) / 2. * sqrt(M_PI/M_LN2);
  218     return true;
  219 }
  220 
  221 ///////////////////////////////////////////////////////////////////////
  222 
  223 void FuncLorentzian::more_precomputations()
  224 {
  225     if (fabs(av_[2]) < epsilon)
  226         av_[2] = epsilon;
  227 }
  228 
  229 CALCULATE_VALUE_BEGIN(FuncLorentzian)
  230     realt xa1a2 = (x - av_[1]) / av_[2];
  231     realt inv_denomin = 1. / (1 + xa1a2 * xa1a2);
  232 CALCULATE_VALUE_END(av_[0] * inv_denomin)
  233 
  234 CALCULATE_DERIV_BEGIN(FuncLorentzian)
  235     realt xa1a2 = (x - av_[1]) / av_[2];
  236     realt inv_denomin = 1. / (1 + xa1a2 * xa1a2);
  237     dy_dv[0] = inv_denomin;
  238     realt dcenter = 2 * av_[0] * xa1a2 / av_[2] * inv_denomin * inv_denomin;
  239     dy_dv[1] = dcenter;
  240     dy_dv[2] = dcenter * xa1a2;
  241     dy_dx = -dcenter;
  242 CALCULATE_DERIV_END(av_[0] * inv_denomin)
  243 
  244 bool FuncLorentzian::get_nonzero_range(double level,
  245                                        realt &left, realt &right) const
  246 {
  247     if (level == 0)
  248         return false;
  249     else if (fabs(level) >= fabs(av_[0]))
  250         left = right = 0;
  251     else {
  252         realt w = sqrt(fabs(av_[0]/level) - 1) * av_[2];
  253         left = av_[1] - w;
  254         right = av_[1] + w;
  255     }
  256     return true;
  257 }
  258 
  259 ///////////////////////////////////////////////////////////////////////
  260 
  261 void FuncPearson7::more_precomputations()
  262 {
  263     if (fabs(av_[2]) < epsilon)
  264         av_[2] = epsilon;
  265     if (av_.size() != 5)
  266         av_.resize(5);
  267     // not checking for av_[3]>0.5 nor even >0
  268     av_[4] = pow(2, 1. / av_[3]) - 1;
  269 }
  270 
  271 CALCULATE_VALUE_BEGIN(FuncPearson7)
  272     realt xa1a2 = (x - av_[1]) / av_[2];
  273     realt xa1a2sq = xa1a2 * xa1a2;
  274     realt pow_2_1_a3_1 = av_[4]; //pow (2, 1. / a3) - 1;
  275     realt denom_base = 1 + xa1a2sq * pow_2_1_a3_1;
  276     realt inv_denomin = pow(denom_base, - av_[3]);
  277 CALCULATE_VALUE_END(av_[0] * inv_denomin)
  278 
  279 CALCULATE_DERIV_BEGIN(FuncPearson7)
  280     realt xa1a2 = (x - av_[1]) / av_[2];
  281     realt xa1a2sq = xa1a2 * xa1a2;
  282     realt pow_2_1_a3_1 = av_[4]; //pow (2, 1. / a3) - 1;
  283     realt denom_base = 1 + xa1a2sq * pow_2_1_a3_1;
  284     realt inv_denomin = pow(denom_base, - av_[3]);
  285     dy_dv[0] = inv_denomin;
  286     realt dcenter = 2 * av_[0] * av_[3] * pow_2_1_a3_1 * xa1a2 * inv_denomin /
  287                                                       (denom_base * av_[2]);
  288     dy_dv[1] = dcenter;
  289     dy_dv[2] = dcenter * xa1a2;
  290     dy_dv[3] = av_[0] * inv_denomin * (M_LN2 * (pow_2_1_a3_1 + 1)
  291                        * xa1a2sq / (denom_base * av_[3]) - log(denom_base));
  292     dy_dx = -dcenter;
  293 CALCULATE_DERIV_END(av_[0] * inv_denomin)
  294 
  295 
  296 bool FuncPearson7::get_nonzero_range(double level,
  297                                      realt &left, realt &right) const
  298 {
  299     if (level == 0)
  300         return false;
  301     else if (fabs(level) >= fabs(av_[0]))
  302         left = right = 0;
  303     else {
  304         realt t = (pow(fabs(av_[0]/level), 1./av_[3]) - 1)
  305                   / (pow(2, 1./av_[3]) - 1);
  306         realt w = sqrt(t) * av_[2];
  307         left = av_[1] - w;
  308         right = av_[1] + w;
  309     }
  310     return true;
  311 }
  312 
  313 bool FuncPearson7::get_area(realt* a) const
  314 {
  315     if (av_[3] <= 0.5)
  316         return false;
  317     realt g = exp(lgamma(av_[3] - 0.5) - lgamma(av_[3]));
  318     //in f_val_precomputations(): av_[4] = pow (2, 1. / a3) - 1;
  319     *a = av_[0] * 2 * fabs(av_[2]) * sqrt(M_PI) * g / (2 * sqrt(av_[4]));
  320     return true;
  321 }
  322 
  323 ///////////////////////////////////////////////////////////////////////
  324 
  325 void FuncSplitPearson7::more_precomputations()
  326 {
  327     if (fabs(av_[2]) < epsilon)
  328         av_[2] = epsilon;
  329     if (fabs(av_[3]) < epsilon)
  330         av_[3] = epsilon;
  331     if (av_.size() != 8)
  332         av_.resize(8);
  333     // not checking for av_[3]>0.5 nor even >0
  334     av_[6] = pow(2, 1. / av_[4]) - 1;
  335     av_[7] = pow(2, 1. / av_[5]) - 1;
  336 }
  337 
  338 CALCULATE_VALUE_BEGIN(FuncSplitPearson7)
  339     int lr = x < av_[1] ? 0 : 1;
  340     realt xa1a2 = (x - av_[1]) / av_[2+lr];
  341     realt xa1a2sq = xa1a2 * xa1a2;
  342     realt pow_2_1_a3_1 = av_[6+lr]; //pow(2, 1./shape) - 1;
  343     realt denom_base = 1 + xa1a2sq * pow_2_1_a3_1;
  344     realt inv_denomin = pow(denom_base, - av_[4+lr]);
  345 CALCULATE_VALUE_END(av_[0] * inv_denomin)
  346 
  347 CALCULATE_DERIV_BEGIN(FuncSplitPearson7)
  348     int lr = x < av_[1] ? 0 : 1;
  349     realt hwhm = av_[2+lr];
  350     realt shape = av_[4+lr];
  351     realt xa1a2 = (x - av_[1]) / hwhm;
  352     realt xa1a2sq = xa1a2 * xa1a2;
  353     realt pow_2_1_a3_1 = av_[6+lr]; //pow(2, 1./shape) - 1;
  354     realt denom_base = 1 + xa1a2sq * pow_2_1_a3_1;
  355     realt inv_denomin = pow (denom_base, -shape);
  356     dy_dv[0] = inv_denomin;
  357     realt dcenter = 2 * av_[0] * shape * pow_2_1_a3_1 * xa1a2 * inv_denomin /
  358                                                       (denom_base * hwhm);
  359     dy_dv[1] = dcenter;
  360     dy_dv[2] = dy_dv[3] = dy_dv[4] = dy_dv[5] = 0;
  361     dy_dv[2+lr] = dcenter * xa1a2;
  362     dy_dv[4+lr] = av_[0] * inv_denomin * (M_LN2 * (pow_2_1_a3_1 + 1)
  363                            * xa1a2sq / (denom_base * shape) - log(denom_base));
  364     dy_dx = -dcenter;
  365 CALCULATE_DERIV_END(av_[0] * inv_denomin)
  366 
  367 
  368 bool FuncSplitPearson7::get_nonzero_range(double level,
  369                                           realt &left, realt &right) const
  370 {
  371     if (level == 0)
  372         return false;
  373     else if (fabs(level) >= fabs(av_[0]))
  374         left = right = 0;
  375     else {
  376         realt t1 = (pow(fabs(av_[0]/level), 1./av_[4]) - 1)
  377                                                 / (pow(2, 1./av_[4]) - 1);
  378         realt w1 = sqrt(t1) * av_[2];
  379         realt t2 = (pow(fabs(av_[0]/level), 1./av_[5]) - 1)
  380                                                 / (pow(2, 1./av_[5]) - 1);
  381         realt w2 = sqrt(t2) * av_[3];
  382         left = av_[1] - w1;
  383         right = av_[1] + w2;
  384     }
  385     return true;
  386 }
  387 
  388 bool FuncSplitPearson7::get_area(realt* a) const
  389 {
  390     if (av_[4] <= 0.5 || av_[5] <= 0.5)
  391         return false;
  392     realt g1 = exp(lgamma(av_[4] - 0.5) - lgamma(av_[4]));
  393     realt g2 = exp(lgamma(av_[5] - 0.5) - lgamma(av_[5]));
  394     *a =   av_[0] * fabs(av_[2]) * sqrt(M_PI) * g1 / (2 * sqrt(av_[6]))
  395          + av_[0] * fabs(av_[3]) * sqrt(M_PI) * g2 / (2 * sqrt(av_[7]));
  396     return true;
  397 }
  398 
  399 ///////////////////////////////////////////////////////////////////////
  400 
  401 void FuncPseudoVoigt::more_precomputations()
  402 {
  403     if (fabs(av_[2]) < epsilon)
  404         av_[2] = epsilon;
  405 }
  406 
  407 CALCULATE_VALUE_BEGIN(FuncPseudoVoigt)
  408     realt xa1a2 = (x - av_[1]) / av_[2];
  409     realt ex = exp(- M_LN2 * xa1a2 * xa1a2);
  410     realt lor = 1. / (1 + xa1a2 * xa1a2);
  411     realt without_height =  (1-av_[3]) * ex + av_[3] * lor;
  412 CALCULATE_VALUE_END(av_[0] * without_height)
  413 
  414 CALCULATE_DERIV_BEGIN(FuncPseudoVoigt)
  415     realt xa1a2 = (x - av_[1]) / av_[2];
  416     realt ex = exp(- M_LN2 * xa1a2 * xa1a2);
  417     realt lor = 1. / (1 + xa1a2 * xa1a2);
  418     realt without_height =  (1-av_[3]) * ex + av_[3] * lor;
  419     dy_dv[0] = without_height;
  420     realt dcenter = 2 * av_[0] * xa1a2 / av_[2]
  421                         * (av_[3]*lor*lor + (1-av_[3])*M_LN2*ex);
  422     dy_dv[1] = dcenter;
  423     dy_dv[2] = dcenter * xa1a2;
  424     dy_dv[3] =  av_[0] * (lor - ex);
  425     dy_dx = -dcenter;
  426 CALCULATE_DERIV_END(av_[0] * without_height)
  427 
  428 bool FuncPseudoVoigt::get_nonzero_range(double level,
  429                                         realt &left, realt &right) const
  430 {
  431     if (level == 0)
  432         return false;
  433     else if (fabs(level) >= fabs(av_[0]))
  434         left = right = 0;
  435     else {
  436         // neglecting Gaussian part and adding 4.0 to compensate it
  437         realt w = (sqrt (av_[3] * fabs(av_[0]/level) - 1) + 4.) * av_[2];
  438         left = av_[1] - w;
  439         right = av_[1] + w;
  440     }
  441     return true;
  442 }
  443 
  444 bool FuncPseudoVoigt::get_area(realt* a) const
  445 {
  446     *a = av_[0] * fabs(av_[2])
  447               * ((av_[3] * M_PI) + (1 - av_[3]) * sqrt(M_PI / M_LN2));
  448     return true;
  449 }
  450 
  451 ///////////////////////////////////////////////////////////////////////
  452 
  453 void FuncVoigt::more_precomputations()
  454 {
  455     if (av_.size() != 6)
  456         av_.resize(6);
  457     float k, l, dkdx, dkdy;
  458     humdev(0, fabs(av_[3]), k, l, dkdx, dkdy);
  459     av_[4] = 1. / k;
  460     av_[5] = dkdy / k;
  461 
  462     if (fabs(av_[2]) < epsilon)
  463         av_[2] = epsilon;
  464 }
  465 
  466 CALCULATE_VALUE_BEGIN(FuncVoigt)
  467     // humdev/humlik routines require with y (a3 here) parameter >0.
  468     float k;
  469     realt xa1a2 = (x - av_[1]) / av_[2];
  470     k = humlik(xa1a2, fabs(av_[3]));
  471 CALCULATE_VALUE_END(av_[0] * av_[4] * k)
  472 
  473 CALCULATE_DERIV_BEGIN(FuncVoigt)
  474     // humdev/humlik routines require with y (a3 here) parameter >0.
  475     // here fabs(av_[3]) is used, and dy_dv[3] is negated if av_[3]<0.
  476     float k;
  477     realt xa1a2 = (x-av_[1]) / av_[2];
  478     realt a0a4 = av_[0] * av_[4];
  479     float l, dkdx, dkdy;
  480     humdev(xa1a2, fabs(av_[3]), k, l, dkdx, dkdy);
  481     dy_dv[0] = av_[4] * k;
  482     realt dcenter = -a0a4 * dkdx / av_[2];
  483     dy_dv[1] = dcenter;
  484     dy_dv[2] = dcenter * xa1a2;
  485     dy_dv[3] = a0a4 * (dkdy - k * av_[5]);
  486     if (av_[3] < 0)
  487         dy_dv[3] = -dy_dv[3];
  488     dy_dx = -dcenter;
  489 CALCULATE_DERIV_END(a0a4 * k)
  490 
  491 bool FuncVoigt::get_nonzero_range(double level,
  492                                   realt &left, realt &right) const
  493 {
  494     if (level == 0)
  495         return false;
  496     realt t = fabs(av_[0]/level);
  497     if (t <= 1) {
  498         left = right = 0;
  499     } else {
  500         // I couldn't find ready-to-use approximation of the Voigt inverse.
  501         // This estimation is used instead.
  502         // width of Lorentzian (exact width when shape -> inf)
  503         realt w_l = av_[3] * sqrt(t - 1);
  504         // width of Gaussian (exact width when shape=0)
  505         realt w_g = sqrt(log(t));
  506         // The sum should do as upper bound of the width at given level.
  507         realt w = (w_l + w_g) * av_[2];
  508         left = av_[1] - w;
  509         right = av_[1] + w;
  510     }
  511     return true;
  512 }
  513 
  514 /// estimation according to
  515 /// http://en.wikipedia.org/w/index.php?title=Voigt_profile&oldid=115518205
  516 ///  fV ~= 0.5346fL + sqrt(0.2166 fL^2 + fG^2)
  517 ///
  518 /// In original paper, Olivero & Longbothum, 1977,
  519 /// http://dx.doi.org/10.1016/0022-4073(77)90161-3
  520 /// this approx. is called "modified Whiting" and is given as:
  521 ///  alphaV = 1/2 { C1 alphaL + sqrt(C2 alphaL^2 + 4 C3 alphaD^2) }
  522 ///   where C1=1.0692, C2=0.86639, C3=1.0
  523 /// Which is the same as the first formula  (C1/2=0.5346, C2/4=0.2165975).
  524 ///
  525 /// Voigt parameters used in fityk (gwidth, shape) are equal:
  526 ///   gwidth = sqrt(2) * sigma
  527 ///   shape = gamma / (sqrt(2) * sigma)
  528 /// where sigma and gamma are defined as in the Wikipedia article.
  529 static realt voigt_fwhm(realt gwidth, realt shape)
  530 {
  531     realt sigma = fabs(gwidth) / M_SQRT2;
  532     realt gamma = fabs(gwidth) * shape;
  533 
  534     realt fG = 2 * sigma * sqrt(2 * M_LN2);
  535     realt fL = 2 * gamma;
  536 
  537     realt fV = 0.5346 * fL + sqrt(0.2165975 * fL * fL + fG * fG);
  538     return fV;
  539 }
  540 
  541 bool FuncVoigt::get_fwhm(realt* a) const
  542 {
  543     *a = voigt_fwhm(av_[2], av_[3]);
  544     return true;
  545 }
  546 
  547 bool FuncVoigt::get_area(realt* a) const
  548 {
  549     *a = av_[0] * fabs(av_[2] * sqrt(M_PI) * av_[4]);
  550     return true;
  551 }
  552 
  553 const vector<string>& FuncVoigt::get_other_prop_names() const
  554 {
  555     static const vector<string> p = vector2(string("GaussianFWHM"),
  556                                             string("LorentzianFWHM"));
  557     return p;
  558 }
  559 
  560 bool FuncVoigt::get_other_prop(string const& name, realt* a) const
  561 {
  562     if (name == "GaussianFWHM") {
  563         realt sigma = fabs(av_[2]) / M_SQRT2;
  564         *a = 2 * sigma * sqrt(2 * M_LN2);
  565         return true;
  566     } else if (name == "LorentzianFWHM") {
  567         realt gamma = fabs(av_[2]) * av_[3];
  568         *a = 2 * gamma;
  569         return true;
  570     } else {
  571         return false;
  572     }
  573 }
  574 
  575 ///////////////////////////////////////////////////////////////////////
  576 
  577 void FuncVoigtA::more_precomputations()
  578 {
  579     if (av_.size() != 6)
  580         av_.resize(6);
  581     av_[4] = 1. / humlik(0, fabs(av_[3]));
  582 
  583     if (fabs(av_[2]) < epsilon)
  584         av_[2] = epsilon;
  585 }
  586 
  587 CALCULATE_VALUE_BEGIN(FuncVoigtA)
  588     // humdev/humlik routines require with y (a3 here) parameter >0.
  589     float k;
  590     realt xa1a2 = (x - av_[1]) / av_[2];
  591     k = humlik(xa1a2, fabs(av_[3]));
  592 CALCULATE_VALUE_END(av_[0] / (sqrt(M_PI) * av_[2]) * k)
  593 
  594 CALCULATE_DERIV_BEGIN(FuncVoigtA)
  595     // humdev/humlik routines require with y (a3 here) parameter >0.
  596     // here fabs(av_[3]) is used, and dy_dv[3] is negated if av_[3]<0.
  597     float k;
  598     realt xa1a2 = (x-av_[1]) / av_[2];
  599     realt f = av_[0] / (sqrt(M_PI) * av_[2]);
  600     float l, dkdx, dkdy;
  601     humdev(xa1a2, fabs(av_[3]), k, l, dkdx, dkdy);
  602     dy_dv[0] = k / (sqrt(M_PI) * av_[2]);
  603     realt dcenter = -f * dkdx / av_[2];
  604     dy_dv[1] = dcenter;
  605     dy_dv[2] = dcenter * xa1a2 - f * k / av_[2];
  606     dy_dv[3] = f * dkdy;
  607     if (av_[3] < 0)
  608         dy_dv[3] = -dy_dv[3];
  609     dy_dx = -dcenter;
  610 CALCULATE_DERIV_END(f * k)
  611 
  612 bool FuncVoigtA::get_nonzero_range(double level,
  613                                    realt &left, realt &right) const
  614 {
  615     if (level == 0)
  616         return false;
  617     else if (fabs(level) >= fabs(av_[0]))
  618         left = right = 0;
  619     else {
  620         //TODO estimate Voigt's non-zero range
  621         return false;
  622     }
  623     return true;
  624 }
  625 
  626 bool FuncVoigtA::get_fwhm(realt* a) const
  627 {
  628     *a = voigt_fwhm(av_[2], av_[3]);
  629     return true;
  630 }
  631 
  632 bool FuncVoigtA::get_height(realt* a) const
  633 {
  634     *a = av_[0] / fabs(av_[2] * sqrt(M_PI) * av_[4]);
  635     return true;
  636 }
  637 
  638 
  639 ///////////////////////////////////////////////////////////////////////
  640 
  641 void FuncEMG::more_precomputations()
  642 {
  643 }
  644 
  645 bool FuncEMG::get_nonzero_range(double/*level*/,
  646                                 realt&/*left*/, realt&/*right*/) const
  647 {
  648     return false;
  649 }
  650 
  651 // approximation to erfc(x) * exp(x*x) for |x| > 4;
  652 // if x < -26.something exp(x*x) in the last line returns inf;
  653 // based on "Rational Chebyshev approximations for the error function"
  654 // by W.J. Cody, Math. Comp., 1969, 631-638.
  655 static double erfcexp_x4(double x)
  656 {
  657     double ax = fabs(x);
  658     assert(ax >= 4.);
  659     const double p[4] = { 3.05326634961232344e-1, 3.60344899949804439e-1,
  660                           1.25781726111229246e-1, 1.60837851487422766e-2 };
  661     const double q[4] = { 2.56852019228982242,    1.87295284992346047,
  662                           5.27905102951428412e-1, 6.05183413124413191e-2 };
  663     double rsq = 1 / (ax * ax);
  664     double xnum = 1.63153871373020978e-2 * rsq;
  665     double xden = rsq;
  666     for (int i = 0; i < 4; ++i) {
  667        xnum = (xnum + p[i]) * rsq;
  668        xden = (xden + q[i]) * rsq;
  669     }
  670     double t = rsq * (xnum + 6.58749161529837803e-4)
  671                    / (xden + 2.33520497626869185e-3);
  672     double v = (1/sqrt(M_PI) - t) / ax;
  673     return x >= 0 ? v : 2*exp(x*x) - v;
  674 }
  675 
  676 CALCULATE_VALUE_BEGIN(FuncEMG)
  677     realt a = av_[0];
  678     realt bx = av_[1] - x;
  679     realt c = av_[2];
  680     realt d = av_[3];
  681     realt fact = c*sqrt(M_PI/2)/d;
  682     realt erf_arg = (bx/c + c/d) / M_SQRT2;
  683     // e_arg == bx/d + c*c/(2*d*d)
  684     // erf_arg^2 == bx^2/(2*c^2) + bx/d  + c^2/(2*d^2)
  685     // e_arg == erf_arg^2 - bx^2/(2*c^2)
  686     realt t;
  687     // type double cannot handle erfc(x) for x >= 28
  688     if (fabs(erf_arg) < 20) {
  689         realt e_arg = bx/d + c*c/(2*d*d);
  690         // t = fact * exp(e_arg) * (d >= 0 ? 1-erf(erf_arg) : -1-erf(erf_arg));
  691         t = fact * exp(e_arg) * (d >= 0 ? erfc(erf_arg) : -erfc(-erf_arg));
  692     } else if ((d >= 0 && erf_arg > -26) || (d < 0 && -erf_arg > -26)) {
  693         realt h = exp(-bx*bx/(2*c*c));
  694         realt ee = d >= 0 ? erfcexp_x4(erf_arg) : -erfcexp_x4(-erf_arg);
  695         t = fact * h * ee;
  696     } else
  697         t = 0;
  698 CALCULATE_VALUE_END(a*t)
  699 
  700 CALCULATE_DERIV_BEGIN(FuncEMG)
  701     realt a = av_[0];
  702     realt bx = av_[1] - x;
  703     realt c = av_[2];
  704     realt d = av_[3];
  705     realt fact = c*sqrt(M_PI/2)/d;
  706     realt erf_arg = (bx/c + c/d) / M_SQRT2; //== (c*c+b*d-d*x)/(M_SQRT2*c*d);
  707     realt ee;
  708     if (fabs(erf_arg) < 20)
  709         ee = exp(erf_arg*erf_arg) * (d >= 0 ? erfc(erf_arg) : -erfc(-erf_arg));
  710     else if ((d >= 0 && erf_arg > -26) || (d < 0 && -erf_arg > -26))
  711         ee = d >= 0 ? erfcexp_x4(erf_arg) : -erfcexp_x4(-erf_arg);
  712     else
  713         ee = 0;
  714     realt h = exp(-bx*bx/(2*c*c));
  715     realt t = fact * h * ee;
  716     dy_dv[0] = t;
  717     dy_dv[1] = -a/d * h + a*t/d;
  718     dy_dv[2] = -a/(c*d*d) * (h * (c*c - bx*d) - t * (c*c + d*d));
  719     dy_dv[3] =  a/(d*d*d) * (h * c*c - t * (c*c + d*d + bx*d));
  720     dy_dx = - dy_dv[1];
  721 CALCULATE_DERIV_END(a*t)
  722 
  723 
  724 bool FuncEMG::get_area(realt* a) const
  725 {
  726     *a = av_[0]*av_[2]*sqrt(2*M_PI);
  727     return true;
  728 }
  729 
  730 ///////////////////////////////////////////////////////////////////////
  731 
  732 bool FuncDoniachSunjic::get_nonzero_range(double/*level*/,
  733                                           realt&/*left*/, realt&/*right*/) const
  734 {
  735     return false;
  736 }
  737 
  738 CALCULATE_VALUE_BEGIN(FuncDoniachSunjic)
  739     realt h = av_[0];
  740     realt a = av_[1];
  741     realt F = av_[2];
  742     realt xE = x - av_[3];
  743     realt t = h * cos(M_PI*a/2 + (1-a)*atan(xE/F)) / pow(F*F+xE*xE, (1-a)/2);
  744 CALCULATE_VALUE_END(t)
  745 
  746 CALCULATE_DERIV_BEGIN(FuncDoniachSunjic)
  747     realt h = av_[0];
  748     realt a = av_[1];
  749     realt F = av_[2];
  750     realt xE = x - av_[3];
  751     realt fe2 = F*F+xE*xE;
  752     realt ac = 1-a;
  753     realt p = pow(fe2, -ac/2);
  754     realt at = atan(xE/F);
  755     realt cos_arg = M_PI*a/2 + ac*at;
  756     realt co = cos(cos_arg);
  757     realt si = sin(cos_arg);
  758     realt t = co * p;
  759     dy_dv[0] = t;
  760     dy_dv[1] = h * p * (co/2 * log(fe2) + (at-M_PI/2) * si);
  761     dy_dv[2] = h * ac*p/fe2 * (xE*si - F*co);
  762     dy_dv[3] = h * ac*p/fe2 * (xE*co + si*F);
  763     dy_dx = -dy_dv[3];
  764 CALCULATE_DERIV_END(h*t)
  765 ///////////////////////////////////////////////////////////////////////
  766 
  767 
  768 CALCULATE_VALUE_BEGIN(FuncPielaszekCube)
  769     realt height = av_[0];
  770     realt center = av_[1];
  771     realt R = av_[2];
  772     realt s = av_[3];
  773     realt s2 = s*s;
  774     realt s4 = s2*s2;
  775     realt R2 = R*R;
  776 
  777     realt q = (x-center);
  778     realt q2 = q*q;
  779     realt t = height * (-3*R*(-1 - (R2*(-1 +
  780                           pow(1 + (q2*s4)/R2, 1.5 - R2/(2.*s2))
  781                           * cos(2*(-1.5 + R2/(2.*s2)) * atan((q*s2)/R))))/
  782            (2.*q2*(-1.5 + R2/(2.*s2))* (-1 + R2/(2.*s2))*s4)))/
  783       (sqrt(2*M_PI)*q2*(-0.5 + R2/(2.*s2))* s2);
  784 CALCULATE_VALUE_END(t)
  785 
  786 CALCULATE_DERIV_BEGIN(FuncPielaszekCube)
  787     realt height = av_[0];
  788     realt center = av_[1];
  789     realt R = av_[2];
  790     realt s = av_[3];
  791     realt s2 = s*s;
  792     realt s3 = s*s2;
  793     realt s4 = s2*s2;
  794     realt R2 = R*R;
  795     realt R4 = R2*R2;
  796     realt R3 = R*R2;
  797 
  798     realt q = (x-center);
  799     realt q2 = q*q;
  800     realt t = (-3*R*(-1 - (R2*(-1 +
  801                               pow(1 + (q2*s4)/R2, 1.5 - R2/(2.*s2))
  802                               * cos(2*(-1.5 + R2/(2.*s2)) * atan((q*s2)/R))))/
  803                (2.*q2*(-1.5 + R2/(2.*s2))* (-1 + R2/(2.*s2))*s4)))/
  804           (sqrt(2*M_PI)*q2*(-0.5 + R2/(2.*s2))* s2);
  805 
  806     realt dcenter = height * (
  807             (3*sqrt(2/M_PI)*R*(-1 -
  808                         (R2* (-1 + pow(1 + (q2*s4)/R2, 1.5 - R2/(2.*s2))*
  809              cos(2*(-1.5 + R2/(2.*s2))* atan((q*s2)/R))))/
  810         (2.*q2*(-1.5 + R2/(2.*s2))* (-1 + R2/(2.*s2))*s4)))/
  811    (q*q2*(-0.5 + R2/(2.*s2))*s2) - (3*R*((R2*(-1 +
  812             pow(1 + (q2*s4)/R2, 1.5 - R2/(2.*s2))*
  813              cos(2*(-1.5 + R2/(2.*s2))* atan((q*s2)/R))))/
  814         (q*q2*(-1.5 + R2/(2.*s2))* (-1 + R2/(2.*s2))*s4) -
  815        (R2*((2*q*(1.5 - R2/(2.*s2))* s4*
  816                pow(1 + (q2*s4)/R2, 0.5 - R2/(2.*s2))*
  817                cos(2*(-1.5 + R2/(2.*s2))* atan((q*s2)/R)))/R2 -
  818             (2*(-1.5 + R2/(2.*s2))*s2* pow(1 + (q2*s4)/R2,
  819                 0.5 - R2/(2.*s2))* sin(2*(-1.5 + R2/(2.*s2))*
  820                  atan((q*s2)/R)))/R))/
  821         (2.*q2*(-1.5 + R2/(2.*s2))* (-1 + R2/(2.*s2))*s4)))/
  822    (sqrt(2*M_PI)*q2*(-0.5 + R2/(2.*s2))* s2));
  823 
  824     realt dR = height * (
  825         (3*R2*(-1 - (R2* (-1 + pow(1 + (q2*s4)/R2,
  826               1.5 - R2/(2.*s2))* cos(2*(-1.5 + R2/(2.*s2))*
  827                atan((q*s2)/R))))/ (2.*q2*(-1.5 + R2/(2.*s2))*
  828           (-1 + R2/(2.*s2))*s4)))/ (sqrt(2*M_PI)*q2*pow(-0.5 + R2/(2.*s2),2)*
  829      s4) - (3*(-1 - (R2*(-1 + pow(1 + (q2*s4)/R2,
  830               1.5 - R2/(2.*s2))* cos(2*(-1.5 + R2/(2.*s2))*
  831                atan((q*s2)/R))))/ (2.*q2*(-1.5 + R2/(2.*s2))*
  832           (-1 + R2/(2.*s2))*s4)))/ (sqrt(2*M_PI)*q2*(-0.5 + R2/(2.*s2))*
  833      s2) - (3*R*((R3* (-1 + pow(1 + (q2*s4)/R2,
  834               1.5 - R2/(2.*s2))* cos(2*(-1.5 + R2/(2.*s2))*
  835                atan((q*s2)/R))))/ (2.*q2*(-1.5 + R2/(2.*s2))*
  836           pow(-1 + R2/(2.*s2),2)*s4*s2) + (R3*(-1 + pow(1 + (q2*s4)/R2,
  837               1.5 - R2/(2.*s2))* cos(2*(-1.5 + R2/(2.*s2))* atan((q*s2)/R))))/
  838         (2.*q2*pow(-1.5 + R2/(2.*s2),2)* (-1 + R2/(2.*s2))*(s4*s2)) -
  839        (R*(-1 + pow(1 + (q2*s4)/R2, 1.5 - R2/(2.*s2))*
  840              cos(2*(-1.5 + R2/(2.*s2))* atan((q*s2)/R))))/
  841         (q2*(-1.5 + R2/(2.*s2))* (-1 + R2/(2.*s2))*s4) -
  842        (R2*(pow(1 + (q2*s4)/R2, 1.5 - R2/(2.*s2))*
  843              cos(2*(-1.5 + R2/(2.*s2))* atan((q*s2)/R))*
  844              ((-2*q2*(1.5 - R2/(2.*s2))* s4)/ (R3*
  845                   (1 + (q2*s4)/R2)) - (R*log(1 + (q2*s4)/R2))/ s2) +
  846             pow(1 + (q2*s4)/R2, 1.5 - R2/(2.*s2))* ((2*q*(-1.5 + R2/(2.*s2))*
  847                   s2)/ (R2* (1 + (q2*s4)/R2)) - (2*R*atan((q*s2)/R))/s2)*
  848              sin(2*(-1.5 + R2/(2.*s2))* atan((q*s2)/R))))/
  849         (2.*q2*(-1.5 + R2/(2.*s2))* (-1 + R2/(2.*s2))*s4)))/
  850    (sqrt(2*M_PI)*q2*(-0.5 + R2/(2.*s2))* s2));
  851 
  852     realt ds = height * (
  853             (-3*R3*(-1 - (R2* (-1 + pow(1 + (q2*s4)/R2,
  854               1.5 - R2/(2.*s2))* cos(2*(-1.5 + R2/(2.*s2))*
  855                atan((q*s2)/R))))/ (2.*q2*(-1.5 + R2/(2.*s2))*
  856           (-1 + R2/(2.*s2))*s4)))/
  857    (sqrt(2*M_PI)*q2*pow(-0.5 + R2/(2.*s2),2)* (s4*s)) + (3*sqrt(2/M_PI)*R*
  858      (-1 - (R2*(-1 + pow(1 + (q2*s4)/R2, 1.5 - R2/(2.*s2))*
  859              cos(2*(-1.5 + R2/(2.*s2))* atan((q*s2)/R))))/
  860         (2.*q2*(-1.5 + R2/(2.*s2))* (-1 + R2/(2.*s2))*s4)))/
  861    (q2*(-0.5 + R2/(2.*s2))*s3) - (3*R*(-(R4*(-1 +
  862              pow(1 + (q2*s4)/R2, 1.5 - R2/(2.*s2))*
  863               cos(2*(-1.5 + R2/(2.*s2))* atan((q*s2)/R))))/
  864         (2.*q2*(-1.5 + R2/(2.*s2))* pow(-1 + R2/(2.*s2),2)*(s4*s3)) -
  865        (R4*(-1 + pow(1 + (q2*s4)/R2, 1.5 - R2/(2.*s2))*
  866              cos(2*(-1.5 + R2/(2.*s2))* atan((q*s2)/R))))/
  867         (2.*q2*pow(-1.5 + R2/(2.*s2),2)* (-1 + R2/(2.*s2))*(s4*s3)) +
  868        (2*R2*(-1 + pow(1 + (q2*s4)/R2,
  869               1.5 - R2/(2.*s2))* cos(2*(-1.5 + R2/(2.*s2))*
  870                atan((q*s2)/R))))/ (q2*(-1.5 + R2/(2.*s2))*
  871           (-1 + R2/(2.*s2))*(s4*s)) - (R2*(pow(1 + (q2*s4)/R2,
  872               1.5 - R2/(2.*s2))* cos(2*(-1.5 + R2/(2.*s2))*
  873                atan((q*s2)/R))* ((4*q2*(1.5 - R2/(2.*s2))* s3)/
  874                 (R2* (1 + (q2*s4)/R2)) + (R2*log(1 +
  875                     (q2*s4)/R2))/ s3) +
  876             pow(1 + (q2*s4)/R2, 1.5 - R2/(2.*s2))*
  877              ((-4*q*(-1.5 + R2/(2.*s2))*s)/ (R*(1 + (q2*s4)/R2)) +
  878                (2*R2*atan((q*s2)/R))/ s3)*
  879              sin(2*(-1.5 + R2/(2.*s2))* atan((q*s2)/R))))/
  880         (2.*q2*(-1.5 + R2/(2.*s2))* (-1 + R2/(2.*s2))*s4)))/
  881    (sqrt(2*M_PI)*q2*(-0.5 + R2/(2.*s2))* s2));
  882 
  883     dy_dv[0] = t;
  884     dy_dv[1] = -dcenter;
  885     dy_dv[2] = dR;
  886     dy_dv[3] = ds;
  887     dy_dx = dcenter;
  888 CALCULATE_DERIV_END(height*t)
  889 
  890 ///////////////////////////////////////////////////////////////////////
  891 
  892 // Implemented by Mirko Scholz. The formula is taken from:
  893 // Bingemann, D.; Ernsting, N. P. J. Chem. Phys. 1995, 102, 2691–2700.
  894 
  895 void FuncLogNormal::more_precomputations()
  896 {
  897     if (av_.size() != 4)
  898         av_.resize(4);
  899     if (fabs(av_[2]) < epsilon)
  900         av_[2] = epsilon;
  901     if (fabs(av_[3]) < epsilon)
  902         av_[3] = 0.001;
  903 }
  904 
  905 CALCULATE_VALUE_BEGIN(FuncLogNormal)
  906     realt a = 2.0 * av_[3] * (x - av_[1]) / av_[2];
  907     realt ex = 0.0;
  908     if (a > -1.0) {
  909         realt b = log(1 + a) / av_[3];
  910         ex = av_[0] * exp(-M_LN2 * b * b);
  911     }
  912 CALCULATE_VALUE_END(ex)
  913 
  914 CALCULATE_DERIV_BEGIN(FuncLogNormal)
  915     realt a = 2.0 * av_[3] * (x - av_[1]) / av_[2];
  916     realt ex;
  917     if (a > -1.0) {
  918         realt b = log(1 + a) / av_[3];
  919         ex = exp(-M_LN2 * b * b);
  920         dy_dv[0] = ex;
  921         ex *= av_[0];
  922         dy_dv[1] = 4.0*M_LN2/(av_[2]*(a+1))*ex*b;
  923         dy_dv[2] = 4.0*(x-av_[1])*M_LN2/(av_[2]*av_[2]*(a+1))*ex*b;
  924         dy_dv[3] = ex*(2.0*M_LN2*b*b/av_[3]
  925             -4.0*(x-av_[1])*log(a+1)*M_LN2/(av_[2]*av_[3]*av_[3]*(a+1)));
  926         dy_dx = -4.0*M_LN2/(av_[2]*(a+1))*ex*b;
  927     } else {
  928         ex = 0.0;
  929         dy_dv[0] = 0.0;
  930         dy_dv[1] = 0.0;
  931         dy_dv[2] = 0.0;
  932         dy_dv[3] = 0.0;
  933         dy_dx = 0.0;
  934     }
  935 CALCULATE_DERIV_END(ex)
  936 
  937 bool FuncLogNormal::get_nonzero_range(double level,
  938                                       realt &left, realt &right) const
  939 { /* untested */
  940     if (level == 0)
  941         return false;
  942     else if (fabs(level) >= fabs(av_[0]))
  943         left = right = 0;
  944     else {
  945         //realt w = sqrt (log (fabs(av_[0]/level)) / M_LN2) * av_[2];
  946         realt w1 = (1-exp(sqrt(log(fabs(av_[0]/level))/M_LN2)*av_[3]))*av_[2]
  947             /2.0/av_[3]+av_[1];
  948         realt w0 = (1-exp(-sqrt(log(fabs(av_[0]/level))/M_LN2)*av_[3]))*av_[2]
  949             /2.0/av_[3]+av_[1];
  950         if (w1>w0) {
  951             left = w0;
  952             right = w1;
  953         } else {
  954             left = w1;
  955             right = w0;
  956         }
  957     }
  958     return true;
  959 }
  960 
  961 //cf. eq. 28 of Maroncelli, M.; Fleming, G.R. J. Phys. Chem. 1987, 86, 6221-6239
  962 bool FuncLogNormal::get_fwhm(realt* a) const
  963 {
  964    *a = av_[2]*sinh(av_[3])/av_[3];
  965    return true;
  966 }
  967 
  968 bool FuncLogNormal::get_area(realt* a) const
  969 {
  970     *a = av_[0]/sqrt(M_LN2/M_PI) / (2.0/av_[2]) / exp(-av_[3]*av_[3]/4.0/M_LN2);
  971     return true;
  972 }
  973 
  974 ///////////////////////////////////////////////////////////////////////
  975 
  976 // so far all the va functions have parameters x1,y1,x2,y2,...
  977 std::string VarArgFunction::get_param(int n) const
  978 {
  979     if (is_index(n, av_))
  980         return (n % 2 == 0 ? "x" : "y") + S(n/2 + 1);
  981     else
  982         return "";
  983 }
  984 
  985 void FuncSpline::more_precomputations()
  986 {
  987     q_.resize(nv() / 2);
  988     for (size_t i = 0; i < q_.size(); ++i) {
  989         q_[i].x = av_[2*i];
  990         q_[i].y = av_[2*i+1];
  991     }
  992     prepare_spline_interpolation(q_);
  993 
  994 }
  995 
  996 CALCULATE_VALUE_BEGIN(FuncSpline)
  997     realt t = get_spline_interpolation(q_, x);
  998 CALCULATE_VALUE_END(t)
  999 
 1000 CALCULATE_DERIV_BEGIN(FuncSpline)
 1001     dy_dx = 0; // unused
 1002     //dy_dv[];
 1003     realt t = get_spline_interpolation(q_, x);
 1004 CALCULATE_DERIV_END(t)
 1005 
 1006 ///////////////////////////////////////////////////////////////////////
 1007 
 1008 void FuncPolyline::more_precomputations()
 1009 {
 1010     q_.resize(nv() / 2);
 1011     for (size_t i = 0; i < q_.size(); ++i) {
 1012         q_[i].x = av_[2*i];
 1013         q_[i].y = av_[2*i+1];
 1014     }
 1015 }
 1016 
 1017 CALCULATE_VALUE_BEGIN(FuncPolyline)
 1018     realt t = get_linear_interpolation(q_, x);
 1019 CALCULATE_VALUE_END(t)
 1020 
 1021 CALCULATE_DERIV_BEGIN(FuncPolyline)
 1022     double value;
 1023     if (q_.empty()) {
 1024         dy_dx = 0;
 1025         value = 0.;
 1026     } else if (q_.size() == 1) {
 1027         //dy_dv[0] = 0; // 0 -> p_x
 1028         dy_dv[1] = 1; // 1 -> p_y
 1029         dy_dx = 0;
 1030         value = q_[0].y;
 1031     } else {
 1032         // value = p0.y + (p1.y - p0.y) / (p1.x - p0.x) * (x - p0.x);
 1033         vector<PointD>::iterator pos = get_interpolation_segment(q_, x);
 1034         double lx = (pos + 1)->x - pos->x;
 1035         double ly = (pos + 1)->y - pos->y;
 1036         double d = x - pos->x;
 1037         double a = ly / lx;
 1038         size_t npos = pos - q_.begin();
 1039         dy_dv[2*npos+0] = a*d/lx - a; // p0.x
 1040         dy_dv[2*npos+1] = 1 - d/lx; // p0.y
 1041         dy_dv[2*npos+2] = -a*d/lx; // p1.x
 1042         dy_dv[2*npos+3] = d/lx; // p1.y
 1043         dy_dx = a;
 1044         value = pos->y + a * d;
 1045     }
 1046 CALCULATE_DERIV_END(value)
 1047 
 1048 ///////////////////////////////////////////////////////////////////////
 1049 
 1050 } // namespace fityk