SiStripRecHitConverter/src/VVIObj.cc

0001 //
0002 //  VVIObj.cc  Version 2.0
0003 //
0004 //  Port of CERNLIB G116 Functions vviden/vvidis
0005 //
0006 // Created by Morris Swartz on 1/14/2010.
0007 // 2010 __TheJohnsHopkinsUniversity__.
0008 //
0009 // V1.1 - make dzero call both fcns with a switch
0010 // V1.2 - remove inappriate initializers and add methods to return non-zero/normalized region
0011 // V2.0 - restructuring and speed improvements by V. Innocente
0012 //
0013
0014 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
0015 // put CMSSW location of SimpleHelix.h here
0016 #include "RecoLocalTracker/SiStripRecHitConverter/interface/VVIObj.h"
0017 #else
0018 #include "VVIObj.h"
0019 #endif
0020
0021 #include <cmath>
0022 #include <algorithm>
0023 #include <functional>
0024
0025 namespace sistripvvi {
0026
0027   namespace VVIObjDetails {
0028     void sincosint(double x, double& sint, double& cint);  //! Private version of the cosine and sine integral
0029     double cosint(double x);                               //! Private version of the cosine integral
0030     double sinint(double x);                               //! Private version of the sine integral
0031     double expint(double x);                               //! Private version of the exponential integral
0032
0033     inline double f1(double x, double const* h_) { return h_[0] + h_[1] * std::log(h_[2] * x) - h_[3] * x; }
0034     inline double f2(double x, double const* h_) {
0035       return h_[4] - x + h_[5] * (std::log(std::abs(x)) + expint(x)) - h_[6] * std::exp(-x);
0036     }
0037     template <typename F>
0038     int dzero(double a, double b, double& x0, double& rv, double eps, int mxf, F func);
0039   }  // namespace VVIObjDetails
0040
0041   // ***************************************************************************************************************************************
0042   //! Constructor
0043   //! Set Vavilov parameters kappa and beta2 and define whether to calculate density fcn or distribution fcn
0044   //! \param kappa - (input) Vavilov kappa parameter [0.01 (Landau-like) < kappa < 10. (Gaussian-like)]
0045   //! \param beta2 - (input) Vavilov beta2 parameter (square of particle speed in v/c units)
0046   //! \param mode  - (input) set to 0 to calculate the density function and to 1 to calculate the distribution function
0047   // ***************************************************************************************************************************************
0048
0049   VVIObj::VVIObj(double kappa, double beta2, int mode) : mode_(mode) {
0050     const double xp[9] = {9.29, 2.47, .89, .36, .15, .07, .03, .02, 0.0};
0051     const double xq[7] = {.012, .03, .08, .26, .87, 3.83, 11.0};
0052     double h_[7];
0053     double q, u, x, c1, c2, c3, c4, d1, h4, h5, h6, q2, x1, d, ll, ul, xf1, xf2, rv;
0054     int lp, lq, k, l, n;
0055
0056     // Make sure that the inputs are reasonable
0057
0058     if (kappa < 0.01)
0059       kappa = 0.01;
0060     if (kappa > 10.)
0061       kappa = 10.;
0062     if (beta2 < 0.)
0063       beta2 = 0.;
0064     if (beta2 > 1.)
0065       beta2 = 1.;
0066
0067     h_[4] = 1. - beta2 * 0.42278433999999998 + 7.6 / kappa;
0068     h_[5] = beta2;
0069     h_[6] = 1. - beta2;
0070     h4 = -7.6 / kappa - (beta2 * .57721566 + 1);
0071     h5 = log(kappa);
0072     h6 = 1. / kappa;
0073     t0_ = (h4 - h_[4] * h5 - (h_[4] + beta2) * (log(h_[4]) + VVIObjDetails::expint(h_[4])) + exp(-h_[4])) / h_[4];
0074
0075     // Set up limits for the root search
0076
0077     for (lp = 0; lp < 9; ++lp) {
0078       if (kappa >= xp[lp])
0079         break;
0080     }
0081     ll = -lp - 1.5;
0082     for (lq = 0; lq < 7; ++lq) {
0083       if (kappa <= xq[lq])
0084         break;
0085     }
0086     ul = lq - 6.5;
0087     //  double (*fp2)(double) = reinterpret_cast<double(*)(double)>(&VVIObj::f2);
0088     VVIObjDetails::dzero(ll, ul, u, rv, 1.e-5, 1000, std::bind(&VVIObjDetails::f2, std::placeholders::_1, h_));
0089     q = 1. / u;
0090     t1_ = h4 * q - h5 - (beta2 * q + 1) * (log((fabs(u))) + VVIObjDetails::expint(u)) + exp(-u) * q;
0091     t_ = t1_ - t0_;
0092     omega_ = 6.2831853000000004 / t_;
0093     h_[0] = kappa * (beta2 * .57721566 + 2.) + 9.9166128600000008;
0094     if (kappa >= .07) {
0095       h_[0] += 6.90775527;
0096     }
0097     h_[1] = beta2 * kappa;
0098     h_[2] = h6 * omega_;
0099     h_[3] = omega_ * 1.5707963250000001;
0100     //  double (*fp1)(double) = reinterpret_cast<double(*)(double)>(&VVIObj::f1);
0101     VVIObjDetails::dzero(5., 155., x0_, rv, 1.e-5, 1000, std::bind(&VVIObjDetails::f1, std::placeholders::_1, h_));
0102     n = x0_ + 1.;
0103     d = exp(kappa * (beta2 * (.57721566 - h5) + 1.)) * .31830988654751274;
0104     a_[n - 1] = 0.;
0105     if (mode_ == 0) {
0106       a_[n - 1] = omega_ * .31830988654751274;
0107     }
0108     q = -1.;
0109     q2 = 2.;
0110     for (k = 1; k < n; ++k) {
0111       l = n - k;
0112       x = omega_ * k;
0113       x1 = h6 * x;
0114       VVIObjDetails::sincosint(x1, c2, c1);
0115       c1 = log(x) - c1;
0116       c3 = sin(x1);
0117       c4 = cos(x1);
0118       xf1 = kappa * (beta2 * c1 - c4) - x * c2;
0119       xf2 = x * c1 + kappa * (c3 + beta2 * c2) + t0_ * x;
0120       if (mode_ == 0) {
0121         d1 = q * d * omega_ * exp(xf1);
0122         a_[l - 1] = d1 * cos(xf2);
0123         b_[l - 1] = -d1 * sin(xf2);
0124       } else {
0125         d1 = q * d * exp(xf1) / k;
0126         a_[l - 1] = d1 * sin(xf2);
0127         b_[l - 1] = d1 * cos(xf2);
0128         a_[n - 1] += q2 * a_[l - 1];
0129       }
0130       q = -q;
0131       q2 = -q2;
0132     }
0133
0134   }  // VVIObj
0135
0136   // *************************************************************************************************************************************
0137   //! Vavilov function method
0138   //! Returns density fcn (mode=0) or distribution fcn (mode=1)
0139   //! \param x  - (input) Argument of function [typically defined as (Q-mpv)/sigma]
0140   // *************************************************************************************************************************************
0141
0142   double VVIObj::fcn(double x) const {
0143     // Local variables
0144
0145     double f, u, y, a0, a1;
0146     double a2 = 0.;
0147     double b1, b0, b2, cof;
0148     int k, n, n1;
0149
0150     n = x0_;
0151     if (x < t0_) {
0152       f = 0.;
0153     } else if (x <= t1_) {
0154       y = x - t0_;
0155       u = omega_ * y - 3.141592653589793;
0156       cof = cos(u) * 2.;
0157       a1 = 0.;
0158       a0 = a_[0];
0159       n1 = n + 1;
0160       for (k = 2; k <= n1; ++k) {
0161         a2 = a1;
0162         a1 = a0;
0163         a0 = a_[k - 1] + cof * a1 - a2;
0164       }
0165       b1 = 0.;
0166       b0 = b_[0];
0167       for (k = 2; k <= n; ++k) {
0168         b2 = b1;
0169         b1 = b0;
0170         b0 = b_[k - 1] + cof * b1 - b2;
0171       }
0172       f = (a0 - a2) * .5 + b0 * sin(u);
0173       if (mode_ != 0) {
0174         f += y / t_;
0175       }
0176     } else {
0177       f = 0.;
0178       if (mode_ != 0) {
0179         f = 1.;
0180       }
0181     }
0182     return f;
0183   }  // fcn
0184
0185   // *************************************************************************************************************************************
0186   //! Vavilov limits method
0187   //! \param xl - (output) Smallest value of the argument for the density and the beginning of the normalized region for the distribution
0188   //! \param xu - (output) Largest value of the argument for the density and the end of the normalized region for the distribution
0189   // *************************************************************************************************************************************
0190
0191   void VVIObj::limits(double& xl, double& xu) const {
0192     xl = t0_;
0193     xu = t1_;
0194     return;
0195   }  // limits
0196
0197   namespace VVIObjDetails {
0198     double cosint(double x) {
0199       // Initialized data
0200
0201       const double zero = 0.;
0202       const double one = 1.;
0203       const double two = 2.;
0204       const double eight = 8.;
0205       const double ce = .57721566490153;
0206       const double c__[14] = {1.9405491464836,
0207                               .9413409132865,
0208                               -.579845034293,
0209                               .3091572011159,
0210                               -.0916101792208,
0211                               .0164437407515,
0212                               -.0019713091952,
0213                               1.692538851e-4,
0214                               -1.09393296e-5,
0215                               5.522386e-7,
0216                               -2.23995e-8,
0217                               7.465e-10,
0218                               -2.08e-11,
0219                               5e-13};
0220       const double p[23] = {
0221           .96074783975204, -.0371138962124, .00194143988899, -1.7165988425e-4, 2.112637753e-5, -3.27163257e-6,
0222           6.0069212e-7,    -1.2586794e-7,   2.932563e-8,     -7.45696e-9,      2.04105e-9,     -5.9502e-10,
0223           1.8323e-10,      -5.921e-11,      1.997e-11,       -7e-12,           2.54e-12,       -9.5e-13,
0224           3.7e-13,         -1.4e-13,        6e-14,           -2e-14,           1e-14};
0225       const double q[20] = {.98604065696238, -.0134717382083, 4.5329284117e-4, -3.067288652e-5, 3.13199198e-6,
0226                             -4.2110196e-7,   6.907245e-8,     -1.318321e-8,    2.83697e-9,      -6.7329e-10,
0227                             1.734e-10,       -4.787e-11,      1.403e-11,       -4.33e-12,       1.4e-12,
0228                             -4.7e-13,        1.7e-13,         -6e-14,          2e-14,           -1e-14};
0229
0230       // System generated locals
0231       double d__1;
0232
0233       // Local variables
0234       double h__;
0235       int i__;
0236       double r__, y, b0, b1, b2, pp, qq, alfa;
0237
0238       // If x==0, return same
0239
0240       if (x == zero) {
0241         return zero;
0242       }
0243       if (fabs(x) <= eight) {
0244         y = x / eight;
0245         // Computing 2nd power
0246         d__1 = y;
0247         h__ = two * (d__1 * d__1) - one;
0248         alfa = -two * h__;
0249         b1 = zero;
0250         b2 = zero;
0251         for (i__ = 13; i__ >= 0; --i__) {
0252           b0 = c__[i__] - alfa * b1 - b2;
0253           b2 = b1;
0254           b1 = b0;
0255         }
0256         b1 = ce + log((fabs(x))) - b0 + h__ * b2;
0257       } else {
0258         r__ = one / x;
0259         y = eight * r__;
0260         // Computing 2nd power
0261         d__1 = y;
0262         h__ = two * (d__1 * d__1) - one;
0263         alfa = -two * h__;
0264         b1 = zero;
0265         b2 = zero;
0266         for (i__ = 22; i__ >= 0; --i__) {
0267           b0 = p[i__] - alfa * b1 - b2;
0268           b2 = b1;
0269           b1 = b0;
0270         }
0271         pp = b0 - h__ * b2;
0272         b1 = zero;
0273         b2 = zero;
0274         for (i__ = 19; i__ >= 0; --i__) {
0275           b0 = q[i__] - alfa * b1 - b2;
0276           b2 = b1;
0277           b1 = b0;
0278         }
0279         qq = b0 - h__ * b2;
0280         b1 = r__ * (qq * sin(x) - r__ * pp * cos(x));
0281       }
0282       return b1;
0283     }  // cosint
0284
0285     double sinint(double x) {
0286       // Initialized data
0287
0288       const double zero = 0.;
0289       const double one = 1.;
0290       const double two = 2.;
0291       const double eight = 8.;
0292       const double pih = 1.5707963267949;
0293       const double s[14] = {1.9522209759531,
0294                             -.6884042321257,
0295                             .4551855132256,
0296                             -.1804571236838,
0297                             .0410422133759,
0298                             -.0059586169556,
0299                             6.001427414e-4,
0300                             -4.44708329e-5,
0301                             2.5300782e-6,
0302                             -1.141308e-7,
0303                             4.1858e-9,
0304                             -1.273e-10,
0305                             3.3e-12,
0306                             -1e-13};
0307       const double p[23] = {
0308           .96074783975204, -.0371138962124, .00194143988899, -1.7165988425e-4, 2.112637753e-5, -3.27163257e-6,
0309           6.0069212e-7,    -1.2586794e-7,   2.932563e-8,     -7.45696e-9,      2.04105e-9,     -5.9502e-10,
0310           1.8323e-10,      -5.921e-11,      1.997e-11,       -7e-12,           2.54e-12,       -9.5e-13,
0311           3.7e-13,         -1.4e-13,        6e-14,           -2e-14,           1e-14};
0312       const double q[20] = {.98604065696238, -.0134717382083, 4.5329284117e-4, -3.067288652e-5, 3.13199198e-6,
0313                             -4.2110196e-7,   6.907245e-8,     -1.318321e-8,    2.83697e-9,      -6.7329e-10,
0314                             1.734e-10,       -4.787e-11,      1.403e-11,       -4.33e-12,       1.4e-12,
0315                             -4.7e-13,        1.7e-13,         -6e-14,          2e-14,           -1e-14};
0316
0317       // System generated locals
0318       double d__1;
0319
0320       // Local variables
0321       double h__;
0322       int i__;
0323       double r__, y, b0, b1, b2, pp, qq, alfa;
0324
0325       if (fabs(x) <= eight) {
0326         y = x / eight;
0327         d__1 = y;
0328         h__ = two * (d__1 * d__1) - one;
0329         alfa = -two * h__;
0330         b1 = zero;
0331         b2 = zero;
0332         for (i__ = 13; i__ >= 0; --i__) {
0333           b0 = s[i__] - alfa * b1 - b2;
0334           b2 = b1;
0335           b1 = b0;
0336         }
0337         b1 = y * (b0 - b2);
0338       } else {
0339         r__ = one / x;
0340         y = eight * r__;
0341         d__1 = y;
0342         h__ = two * (d__1 * d__1) - one;
0343         alfa = -two * h__;
0344         b1 = zero;
0345         b2 = zero;
0346         for (i__ = 22; i__ >= 0; --i__) {
0347           b0 = p[i__] - alfa * b1 - b2;
0348           b2 = b1;
0349           b1 = b0;
0350         }
0351         pp = b0 - h__ * b2;
0352         b1 = zero;
0353         b2 = zero;
0354         for (i__ = 19; i__ >= 0; --i__) {
0355           b0 = q[i__] - alfa * b1 - b2;
0356           b2 = b1;
0357           b1 = b0;
0358         }
0359         qq = b0 - h__ * b2;
0360         d__1 = fabs(pih);
0361         if (x < 0.)
0362           d__1 = -d__1;
0363         b1 = d__1 - r__ * (r__ * pp * sin(x) + qq * cos(x));
0364       }
0365
0366       return b1;
0367     }  // sinint
0368
0369     void sincosint(double x, double& sint, double& cint) {
0370       // Initialized data
0371
0372       const double zero = 0.;
0373       const double one = 1.;
0374       const double two = 2.;
0375       const double eight = 8.;
0376       const double ce = .57721566490153;
0377       const double pih = 1.5707963267949;
0378       const double s__[14] = {1.9522209759531,
0379                               -.6884042321257,
0380                               .4551855132256,
0381                               -.1804571236838,
0382                               .0410422133759,
0383                               -.0059586169556,
0384                               6.001427414e-4,
0385                               -4.44708329e-5,
0386                               2.5300782e-6,
0387                               -1.141308e-7,
0388                               4.1858e-9,
0389                               -1.273e-10,
0390                               3.3e-12,
0391                               -1e-13};
0392
0393       const double c__[14] = {1.9405491464836,
0394                               .9413409132865,
0395                               -.579845034293,
0396                               .3091572011159,
0397                               -.0916101792208,
0398                               .0164437407515,
0399                               -.0019713091952,
0400                               1.692538851e-4,
0401                               -1.09393296e-5,
0402                               5.522386e-7,
0403                               -2.23995e-8,
0404                               7.465e-10,
0405                               -2.08e-11,
0406                               5e-13};
0407
0408       const double p[23] = {
0409           .96074783975204, -.0371138962124, .00194143988899, -1.7165988425e-4, 2.112637753e-5, -3.27163257e-6,
0410           6.0069212e-7,    -1.2586794e-7,   2.932563e-8,     -7.45696e-9,      2.04105e-9,     -5.9502e-10,
0411           1.8323e-10,      -5.921e-11,      1.997e-11,       -7e-12,           2.54e-12,       -9.5e-13,
0412           3.7e-13,         -1.4e-13,        6e-14,           -2e-14,           1e-14};
0413       const double q[20] = {.98604065696238, -.0134717382083, 4.5329284117e-4, -3.067288652e-5, 3.13199198e-6,
0414                             -4.2110196e-7,   6.907245e-8,     -1.318321e-8,    2.83697e-9,      -6.7329e-10,
0415                             1.734e-10,       -4.787e-11,      1.403e-11,       -4.33e-12,       1.4e-12,
0416                             -4.7e-13,        1.7e-13,         -6e-14,          2e-14,           -1e-14};
0417
0418       // System generated locals
0419       double d__1;
0420
0421       // Local variables
0422       double h__;
0423       int i__;
0424       double r__, y, b0, b1, b2, pp, qq, alfa;
0425
0426       sint = 0;
0427       cint = 0;
0428
0429       if (fabs(x) <= eight) {
0430         y = x / eight;
0431         // Computing 2nd power
0432         d__1 = y;
0433         h__ = two * (d__1 * d__1) - one;
0434         alfa = -two * h__;
0435
0436         // cos
0437         if (x != 0) {
0438           b1 = zero;
0439           b2 = zero;
0440           for (i__ = 13; i__ >= 0; --i__) {
0441             b0 = c__[i__] - alfa * b1 - b2;
0442             b2 = b1;
0443             b1 = b0;
0444           }
0445           cint = ce + log((fabs(x))) - b0 + h__ * b2;
0446         }
0447         // sin
0448         b1 = zero;
0449         b2 = zero;
0450         for (i__ = 13; i__ >= 0; --i__) {
0451           b0 = s__[i__] - alfa * b1 - b2;
0452           b2 = b1;
0453           b1 = b0;
0454         }
0455         sint = y * (b0 - b2);
0456
0457       } else {
0458         r__ = one / x;
0459         y = eight * r__;
0460         // Computing 2nd power
0461         d__1 = y;
0462         h__ = two * (d__1 * d__1) - one;
0463         alfa = -two * h__;
0464         b1 = zero;
0465         b2 = zero;
0466         for (i__ = 22; i__ >= 0; --i__) {
0467           b0 = p[i__] - alfa * b1 - b2;
0468           b2 = b1;
0469           b1 = b0;
0470         }
0471         pp = b0 - h__ * b2;
0472         b1 = zero;
0473         b2 = zero;
0474         for (i__ = 19; i__ >= 0; --i__) {
0475           b0 = q[i__] - alfa * b1 - b2;
0476           b2 = b1;
0477           b1 = b0;
0478         }
0479         qq = b0 - h__ * b2;
0480         // cos
0481         cint = r__ * (qq * sin(x) - r__ * pp * cos(x));
0482         // sin
0483         d__1 = pih;
0484         if (x < 0.)
0485           d__1 = -d__1;
0486         sint = d__1 - r__ * (r__ * pp * sin(x) + qq * cos(x));
0487       }
0488     }
0489
0490     double expint(double x) {
0491       // Initialized data
0492
0493       const double zero = 0.;
0494       const double q2[7] = {
0495           .10340013040487, 3.319092135933, 20.449478501379, 41.280784189142, 32.426421069514, 10.041164382905, 1.};
0496       const double p3[6] = {
0497           -2.3909964453136, -147.98219500504, -254.3763397689, -119.55761038372, -19.630408535939, -.9999999999036};
0498       const double q3[6] = {177.60070940351, 530.68509610812, 462.23027156148, 156.81843364539, 21.630408494238, 1.};
0499       const double p4[8] = {-8.6693733995107,
0500                             -549.14226552109,
0501                             -4210.0161535707,
0502                             -249301.39345865,
0503                             -119623.66934925,
0504                             -22174462.775885,
0505                             3892804.213112,
0506                             -391546073.8091};
0507       const double q4[8] = {34.171875,
0508                             -1607.0892658722,
0509                             35730.029805851,
0510                             -483547.43616216,
0511                             4285596.2461175,
0512                             -24903337.574054,
0513                             89192576.757561,
0514                             -165254299.72521};
0515       const double a1[8] = {-2.1808638152072,
0516                             -21.901023385488,
0517                             9.3081638566217,
0518                             25.076281129356,
0519                             -33.184253199722,
0520                             60.121799083008,
0521                             -43.253113287813,
0522                             1.0044310922808};
0523       const double b1[8] = {0.,
0524                             3.9370770185272,
0525                             300.89264837292,
0526                             -6.2504116167188,
0527                             1003.6743951673,
0528                             14.325673812194,
0529                             2736.2411988933,
0530                             .52746885196291};
0531       const double a2[8] = {-3.4833465360285,
0532                             -18.65454548834,
0533                             -8.2856199414064,
0534                             -32.34673303054,
0535                             17.960168876925,
0536                             1.7565631546961,
0537                             -1.9502232128966,
0538                             .99999429607471};
0539       const double b2[8] = {0.,
0540                             69.500065588743,
0541                             57.283719383732,
0542                             25.777638423844,
0543                             760.76114800773,
0544                             28.951672792514,
0545                             -3.4394226689987,
0546                             1.0008386740264};
0547       const double a3[6] = {
0548           -27.780928934438, -10.10479081576, -9.1483008216736, -5.0223317461851, -3.0000077799358, 1.0000000000704};
0549       const double one = 1.;
0550       const double b3[6] = {0., 122.39993926823, 2.7276100778779, -7.1897518395045, -2.9990118065262, 1.999999942826};
0551       const double two = 2.;
0552       const double three = 3.;
0553       const double x0 = .37250741078137;
0554       const double xl[6] = {-24., -12., -6., 0., 1., 4.};
0555       const double p1[5] = {4.293125234321, 39.894153870321, 292.52518866921, 425.69682638592, -434.98143832952};
0556       const double q1[5] = {1., 18.899288395003, 150.95038744251, 568.05252718987, 753.58564359843};
0557       const double p2[7] = {.43096783946939,
0558                             6.9052252278444,
0559                             23.019255939133,
0560                             24.378408879132,
0561                             9.0416155694633,
0562                             .99997957705159,
0563                             4.656271079751e-7};
0564
0565       /* Local variables */
0566       double v, y, ap, bp, aq, dp, bq, dq;
0567
0568       if (x <= xl[0]) {
0569         ap = a3[0] - x;
0570         for (int i__ = 2; i__ <= 5; ++i__) {
0571           /* L1: */
0572           ap = a3[i__ - 1] - x + b3[i__ - 1] / ap;
0573         }
0574         y = exp(-x) / x * (one - (a3[5] + b3[5] / ap) / x);
0575       } else if (x <= xl[1]) {
0576         ap = a2[0] - x;
0577         for (int i__ = 2; i__ <= 7; ++i__) {
0578           ap = a2[i__ - 1] - x + b2[i__ - 1] / ap;
0579         }
0580         y = exp(-x) / x * (a2[7] + b2[7] / ap);
0581       } else if (x <= xl[2]) {
0582         ap = a1[0] - x;
0583         for (int i__ = 2; i__ <= 7; ++i__) {
0584           ap = a1[i__ - 1] - x + b1[i__ - 1] / ap;
0585         }
0586         y = exp(-x) / x * (a1[7] + b1[7] / ap);
0587       } else if (x < xl[3]) {
0588         v = -two * (x / three + one);
0589         bp = zero;
0590         dp = p4[0];
0591         for (int i__ = 2; i__ <= 8; ++i__) {
0592           ap = bp;
0593           bp = dp;
0594           dp = p4[i__ - 1] - ap + v * bp;
0595         }
0596         bq = zero;
0597         dq = q4[0];
0598         for (int i__ = 2; i__ <= 8; ++i__) {
0599           aq = bq;
0600           bq = dq;
0601           dq = q4[i__ - 1] - aq + v * bq;
0602         }
0603         y = -log(-x / x0) + (x + x0) * (dp - ap) / (dq - aq);
0604       } else if (x == xl[3]) {
0605         return zero;
0606       } else if (x < xl[4]) {
0607         ap = p1[0];
0608         aq = q1[0];
0609         for (int i__ = 2; i__ <= 5; ++i__) {
0610           ap = p1[i__ - 1] + x * ap;
0611           aq = q1[i__ - 1] + x * aq;
0612         }
0613         y = -log(x) + ap / aq;
0614       } else if (x <= xl[5]) {
0615         y = one / x;
0616         ap = p2[0];
0617         aq = q2[0];
0618         for (int i__ = 2; i__ <= 7; ++i__) {
0619           ap = p2[i__ - 1] + y * ap;
0620           aq = q2[i__ - 1] + y * aq;
0621         }
0622         y = exp(-x) * ap / aq;
0623       } else {
0624         y = one / x;
0625         ap = p3[0];
0626         aq = q3[0];
0627         for (int i__ = 2; i__ <= 6; ++i__) {
0628           ap = p3[i__ - 1] + y * ap;
0629           aq = q3[i__ - 1] + y * aq;
0630         }
0631         y = exp(-x) * y * (one + y * ap / aq);
0632       }
0633       return y;
0634     }  // expint
0635
0636     template <typename F>
0637     int dzero(double a, double b, double& x0, double& rv, double eps, int mxf, F func) {
0638       /* System generated locals */
0639       double d__1, d__2, d__3, d__4;
0640
0641       // Local variables
0642       double f1, f2, f3, u1, u2, x1, x2, u3, u4, x3, ca, cb, cc, fa, fb, ee, ff;
0643       int mc;
0644       double xa, xb, fx, xx, su4;
0645
0646       xa = std::min(a, b);
0647       xb = std::max(a, b);
0648       fa = func(xa);
0649       fb = func(xb);
0650       if (fa * fb > 0.) {
0651         rv = (xb - xa) * -2;
0652         x0 = 0.;
0653         return 1;
0654       }
0655       mc = 0;
0656     L1:
0657       x0 = (xa + xb) * .5;
0658       rv = x0 - xa;
0659       ee = eps * (fabs(x0) + 1);
0660       if (rv <= ee) {
0661         rv = ee;
0662         ff = func(x0);
0663         return 0;
0664       }
0665       f1 = fa;
0666       x1 = xa;
0667       f2 = fb;
0668       x2 = xb;
0669     L2:
0670       fx = func(x0);
0671       ++mc;
0672       if (mc > mxf) {
0673         rv = (d__1 = xb - xa, fabs(d__1)) * -.5;
0674         x0 = 0.;
0675         return 0;
0676       }
0677       if (fx * fa > 0.) {
0678         xa = x0;
0679         fa = fx;
0680       } else {
0681         xb = x0;
0682         fb = fx;
0683       }
0684     L3:
0685       u1 = f1 - f2;
0686       u2 = x1 - x2;
0687       u3 = f2 - fx;
0688       u4 = x2 - x0;
0689       if (u2 == 0. || u4 == 0.) {
0690         goto L1;
0691       }
0692       f3 = fx;
0693       x3 = x0;
0694       u1 /= u2;
0695       u2 = u3 / u4;
0696       ca = u1 - u2;
0697       cb = (x1 + x2) * u2 - (x2 + x0) * u1;
0698       cc = (x1 - x0) * f1 - x1 * (ca * x1 + cb);
0699       if (ca == 0.) {
0700         if (cb == 0.) {
0701           goto L1;
0702         }
0703         x0 = -cc / cb;
0704       } else {
0705         u3 = cb / (ca * 2);
0706         u4 = u3 * u3 - cc / ca;
0707         if (u4 < 0.) {
0708           goto L1;
0709         }
0710         su4 = fabs(u4);
0711         if (x0 + u3 < 0.f) {
0712           su4 = -su4;
0713         }
0714         x0 = -u3 + su4;
0715       }
0716       if (x0 < xa || x0 > xb) {
0717         goto L1;
0718       }
0719       // Computing MIN
0720       d__3 = (d__1 = x0 - x3, fabs(d__1)), d__4 = (d__2 = x0 - x2, fabs(d__2));
0721       rv = std::min(d__3, d__4);
0722       ee = eps * (fabs(x0) + 1);
0723       if (rv > ee) {
0724         f1 = f2;
0725         x1 = x2;
0726         f2 = f3;
0727         x2 = x3;
0728         goto L2;
0729       }
0730       fx = func(x0);
0731       if (fx == 0.) {
0732         rv = ee;
0733         ff = func(x0);
0734         return 0;
0735       }
0736       if (fx * fa < 0.) {
0737         xx = x0 - ee;
0738         if (xx <= xa) {
0739           rv = ee;
0740           ff = func(x0);
0741           return 0;
0742         }
0743         ff = func(xx);
0744         fb = ff;
0745         xb = xx;
0746       } else {
0747         xx = x0 + ee;
0748         if (xx >= xb) {
0749           rv = ee;
0750           ff = func(x0);
0751           return 0;
0752         }
0753         ff = func(xx);
0754         fa = ff;
0755         xa = xx;
0756       }
0757       if (fx * ff > 0.) {
0758         mc += 2;
0759         if (mc > mxf) {
0760           rv = (d__1 = xb - xa, fabs(d__1)) * -.5;
0761           x0 = 0.;
0762           return 0;
0763         }
0764         f1 = f3;
0765         x1 = x3;
0766         f2 = fx;
0767         x2 = x0;
0768         x0 = xx;
0769         fx = ff;
0770         goto L3;
0771       }
0772       /* L4: */
0773       rv = ee;
0774       ff = func(x0);
0775       return 0;
0776     }  // dzero
0777
0778   }  // namespace VVIObjDetails
0779 }  // namespace sistripvvi