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