---

glande.cpp

Go to the documentation of this file.

/* glande.f -- translated by f2c (version 19950105).
   You must link the resulting object file with the libraries:
        -lf2c -lm   (in that order)
*/

#ifdef __cplusplus
extern "C" {
#endif
#include "../gheisha_f2c/f2c.h"

doublereal glande_(real *x)
{
    /* Initialized data */

    static real p1[5] = { (float).4259894875,(float)-.124976255,(float)
            .039842437,(float)-.006298287635,(float).001511162253 };
    static real q5[5] = { (float)1.,(float)156.9424537,(float)3745.310488,(
            float)9834.698876,(float)66924.28357 };
    static real p6[5] = { (float)1.000827619,(float)664.9143136,(float)
            62972.92665,(float)475554.6998,(float)-5743609.109 };
    static real q6[5] = { (float)1.,(float)651.4101098,(float)56974.73333,(
            float)165917.4725,(float)-2815759.939 };
    static real a1[3] = { (float).04166666667,(float)-.01996527778,(float)
            .02709538966 };
    static real a2[2] = { (float)-1.84556867,(float)-4.284640743 };
    static real q1[5] = { (float)1.,(float)-.3388260629,(float).09594393323,(
            float)-.01608042283,(float).003778942063 };
    static real p2[5] = { (float).1788541609,(float).1173957403,(float)
            .01488850518,(float)-.001394989411,(float)1.283617211e-4 };
    static real q2[5] = { (float)1.,(float).7428795082,(float).3153932961,(
            float).06694219548,(float).008790609714 };
    static real p3[5] = { (float).1788544503,(float).09359161662,(float)
            .006325387654,(float)6.611667319e-5,(float)-2.031049101e-6 };
    static real q3[5] = { (float)1.,(float).6097809921,(float).2560616665,(
            float).04746722384,(float).006957301675 };
    static real p4[5] = { (float).9874054407,(float)118.6723273,(float)
            849.279436,(float)-743.7792444,(float)427.0262186 };
    static real q4[5] = { (float)1.,(float)106.8615961,(float)337.6496214,(
            float)2016.712389,(float)1597.063511 };
    static real p5[5] = { (float)1.003675074,(float)167.5702434,(float)
            4789.711289,(float)21217.86767,(float)-22324.9491 };

    /* System generated locals */
    real ret_val, r__1;

    /* Builtin functions */
    double exp(doublereal), sqrt(doublereal), log(doublereal);

    /* Local variables */
    static real u, v;

/* . */
/* .    ******************************************************************
 */
/* .    *                                                                *
 */
/* .    *  Copy of the CERN library routine DENLAN (G110)                *
 */
/* .    *                                                                *
 */
/* .    ******************************************************************
 */
/* . */









    v = *x;
    if (v < (float)-5.5) {
        u = exp(v + (float)1.);
        ret_val = exp((float)-1. / u) / sqrt(u) * (float).3989422803 * ((a1[0]
                 + (a1[1] + a1[2] * u) * u) * u + (float)1.);
    } else if (v < (float)-1.) {
        u = exp(-(doublereal)v - (float)1.);
        ret_val = exp(-(doublereal)u) * sqrt(u) * (p1[0] + (p1[1] + (p1[2] + (
                p1[3] + p1[4] * v) * v) * v) * v) / (q1[0] + (q1[1] + (q1[2] 
                + (q1[3] + q1[4] * v) * v) * v) * v);
    } else if (v < (float)1.) {
        ret_val = (p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * 
                v) / (q2[0] + (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v) 
                * v);
    } else if (v < (float)5.) {
        ret_val = (p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * 
                v) / (q3[0] + (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v) 
                * v);
    } else if (v < (float)12.) {
        u = (float)1. / v;
/* Computing 2nd power */
        r__1 = u;
        ret_val = r__1 * r__1 * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u)
                 * u) * u) * u) / (q4[0] + (q4[1] + (q4[2] + (q4[3] + q4[4] * 
                u) * u) * u) * u);
    } else if (v < (float)50.) {
        u = (float)1. / v;
/* Computing 2nd power */
        r__1 = u;
        ret_val = r__1 * r__1 * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u)
                 * u) * u) * u) / (q5[0] + (q5[1] + (q5[2] + (q5[3] + q5[4] * 
                u) * u) * u) * u);
    } else if (v < (float)300.) {
        u = (float)1. / v;
/* Computing 2nd power */
        r__1 = u;
        ret_val = r__1 * r__1 * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u)
                 * u) * u) * u) / (q6[0] + (q6[1] + (q6[2] + (q6[3] + q6[4] * 
                u) * u) * u) * u);
    } else {
        u = (float)1. / (v - v * log(v) / (v + (float)1.));
/* Computing 2nd power */
        r__1 = u;
        ret_val = r__1 * r__1 * ((a2[0] + a2[1] * u) * u + (float)1.);
    }
    return ret_val;
} /* glande_ */

/* ** */
/* Subroutine */ int glandg_(real *yran)
{
    /* Initialized data */

    static real xcum1[100] = { (float)-2.57,(float)-2.15412,(float)-1.94167,(
            float)-1.79583,(float)-1.67975,(float)-1.5809,(float)-1.49341,(
            float)-1.41397,(float)-1.34057,(float)-1.27185,(float)-1.20686,(
            float)-1.1449,(float)-1.08545,(float)-1.02809,(float)-.97249,(
            float)-.91839,(float)-.86557,(float)-.81386,(float)-.76308,(float)
            -.71311,(float)-.66384,(float)-.61515,(float)-.56696,(float)
            -.51919,(float)-.47176,(float)-.42462,(float)-.37769,(float)
            -.33092,(float)-.28426,(float)-.23767,(float)-.19109,(float)
            -.14447,(float)-.09779,(float)-.051,(float)-.00405,(float).04308,(
            float).09044,(float).13806,(float).18598,(float).23423,(float)
            .28285,(float).33187,(float).38134,(float).43128,(float).48174,(
            float).53275,(float).58435,(float).63658,(float).68948,(float)
            .74309,(float).79746,(float).85263,(float).90865,(float).96556,(
            float)1.02342,(float)1.08228,(float)1.14219,(float)1.20322,(float)
            1.26542,(float)1.32887,(float)1.39362,(float)1.45976,(float)
            1.52736,(float)1.5965,(float)1.66727,(float)1.73976,(float)
            1.81407,(float)1.89032,(float)1.9686,(float)2.04905,(float)2.1318,
            (float)2.21699,(float)2.30477,(float)2.39531,(float)2.48878,(
            float)2.5854,(float)2.68535,(float)2.78889,(float)2.89626,(float)
            3.00775,(float)3.12364,(float)3.24429,(float)3.37005,(float)
            3.50136,(float)3.63866,(float)3.78246,(float)3.93334,(float)
            4.09194,(float)4.259,(float)4.43533,(float)4.62186,(float)4.8196,(
            float)5.02974,(float)5.25368,(float)5.49312,(float)5.74987,(float)
            6.02605,(float)6.32428,(float)6.64773,(float)7. };
    static real xcum2[100] = { (float)5.5,(float)5.5612,(float)5.62347,(float)
            5.68684,(float)5.75133,(float)5.81699,(float)5.88383,(float)
            5.95191,(float)6.02125,(float)6.0919,(float)6.16391,(float)
            6.23732,(float)6.31219,(float)6.38855,(float)6.46646,(float)
            6.54597,(float)6.62712,(float)6.70998,(float)6.7946,(float)
            6.88103,(float)6.96935,(float)7.05962,(float)7.15192,(float)
            7.24633,(float)7.34294,(float)7.44182,(float)7.54306,(float)
            7.64676,(float)7.753,(float)7.86188,(float)7.97351,(float)8.088,(
            float)8.20548,(float)8.3261,(float)8.44997,(float)8.57725,(float)
            8.70808,(float)8.84262,(float)8.98103,(float)9.12349,(float)
            9.27021,(float)9.42141,(float)9.5773,(float)9.73812,(float)9.9041,
            (float)10.07552,(float)10.25265,(float)10.43584,(float)10.6254,(
            float)10.82169,(float)11.02508,(float)11.23598,(float)11.45487,(
            float)11.68222,(float)11.91855,(float)12.16441,(float)12.42045,(
            float)12.68734,(float)12.9658,(float)13.25663,(float)13.56075,(
            float)13.87912,(float)14.2128,(float)14.56296,(float)14.93088,(
            float)15.31799,(float)15.72593,(float)16.1565,(float)16.6117,(
            float)17.09378,(float)17.60523,(float)18.14884,(float)18.72792,(
            float)19.34622,(float)20.00796,(float)20.71792,(float)21.48176,(
            float)22.30618,(float)23.19862,(float)24.16809,(float)25.22547,(
            float)26.3832,(float)27.65658,(float)29.0649,(float)30.63048,(
            float)32.38211,(float)34.35555,(float)36.59607,(float)39.16212,(
            float)42.13195,(float)45.61204,(float)49.74758,(float)54.74189,(
            float)60.89935,(float)68.6737,(float)78.8148,(float)92.61047,(
            float)112.50807,(float)143.78539,(float)200. };

    /* System generated locals */
    real r__1;

    /* Builtin functions */
    double log(doublereal), sqrt(doublereal);

    /* Local variables */
    extern doublereal rndm_(real *);
    static real a, b;
    static integer j;
    static real p, q, x, q1, tabpo1, tabpo2, qq, yy;
    extern /* Subroutine */ int rannor_(real *, real *);

/* . */
/* .    ******************************************************************
 */
/* .    *                                                                *
 */
/* .    *  Copy of the CERN library routine GENLAN                       *
 */
/* .    *  Generation of LANDAU-distributed random numbers by 4-point    *
 */
/* .    *  interpolation in the previously-tabulated inverse cumulative  *
 */
/* .    *  distribution                                                  *
 */
/* .    *                                                                *
 */
/* .    ******************************************************************
 */
/* . */


/* ====>  1ST TABLE OF INVERSE CUMULATIVE LANDAU POINTS. ( 0.0<P<0.833 ) 
*/

/* ====>  2ND TABLE OF INVERSE CUMULATIVE LANDAU POINTS. ( 0.791<P<0.995 )
 */
/* . */
/* .    ------------------------------------------------------------------
 */
/* . */
    rannor_(&q, &q1);
    x = rndm_(&q);
    if (x < (float).004) {

/* *** Extreme left-hand tail */
        *yran = -(doublereal)sqrt((r__1 = log(x), dabs(r__1)));

    } else if (x <= (float).807069) {

/* *** 4-point interpolation in the first cumulative table */
        tabpo1 = (x + (float)0.) * (float)118.836;
        j = tabpo1 + 1;
        j = max(j,2);
        j = min(j,98);
        p = tabpo1 - (j - 1);
        a = (p + (float)1.) * xcum1[j + 1] - (p - (float)2.) * xcum1[j - 2];
        b = (p - (float)1.) * xcum1[j - 1] - p * xcum1[j];
        *yran = a * p * (p - (float)1.) * (float).16666667 + b * (p + (float)
                1.) * (p - (float)2.) * (float).5;

    } else if (x <= (float).994869) {

/* *** 4-point interpolation in the first cumulative table */
        tabpo2 = (x - (float).79124) * (float)486.178;
        j = tabpo2 + 1;
        j = max(j,2);
        j = min(j,98);
        p = tabpo2 - (j - 1);
        a = (p + (float)1.) * xcum2[j + 1] - (p - (float)2.) * xcum2[j - 2];
        b = (p - (float)1.) * xcum2[j - 1] - p * xcum2[j];
        *yran = a * p * (p - (float)1.) * (float).16666667 + b * (p + (float)
                1.) * (p - (float)2.) * (float).5;

    } else {

/* *** 1/x**2 sampling for extreme Landau tail */
        yy = rndm_(&qq);
        *yran = (float)200. / yy;

    }
    return 0;
} /* glandg_ */

/* ** */
doublereal glands_(real *x)
{
    /* Initialized data */

    static real p1[5] = { (float).2514091491,(float)-.06250580444,(float)
            .0145838123,(float)-.002108817737,(float)7.41124729e-4 };
    static real q1[5] = { (float)1.,(float)-.005571175625,(float).06225310236,
            (float)-.003137378427,(float).001931496439 };
    static real p2[4] = { (float).2868328584,(float).3564363231,(float)
            .1523518695,(float).02251304883 };
    static real q2[4] = { (float)1.,(float).6191136137,(float).1720721448,(
            float).02278594771 };
    static real p3[4] = { (float).2868329066,(float).3003828436,(float)
            .09950951941,(float).008733827185 };
    static real q3[4] = { (float)1.,(float).4237190502,(float).1095631512,(
            float).008693851567 };
    static real p4[4] = { (float)1.00035163,(float)4.503592498,(float)
            10.8588388,(float)7.536052269 };
    static real q4[4] = { (float)1.,(float)5.539969678,(float)19.33581111,(
            float)27.21321508 };
    static real p5[4] = { (float)1.000006517,(float)49.09414111,(float)
            85.05544753,(float)153.2153455 };
    static real q5[4] = { (float)1.,(float)50.09928881,(float)139.9819104,(
            float)420.0002909 };
    static real p6[4] = { (float)1.000000983,(float)132.9868456,(float)
            916.2149244,(float)-960.5054274 };
    static real q6[4] = { (float)1.,(float)133.9887843,(float)1055.990413,(
            float)553.2224619 };
    static real a1[3] = { (float)-.4583333333,(float).6675347222,(float)
            -1.641741416 };
    static real a2[3] = { (float)1.,(float)-.4227843351,(float)-2.043403138 };

    /* System generated locals */
    real ret_val;

    /* Builtin functions */
    double exp(doublereal), sqrt(doublereal), log(doublereal);

    /* Local variables */
    static real u, v;

/* . */
/* .    ******************************************************************
 */
/* .    *                                                                *
 */
/* .    *  Copy of the CERN library routine DSTLAN (G110)                *
 */
/* .    *                                                                *
 */
/* .    ******************************************************************
 */
/* . */









    v = *x;
    if (v < (float)-5.5) {
        u = exp(v + (float)1.);
        ret_val = exp((float)-1. / u) * (float).3989422803 * sqrt(u) * ((a1[0]
                 + (a1[1] + a1[2] * u) * u) * u + (float)1.);
    } else if (v < (float)-1.) {
        u = exp(-(doublereal)v - (float)1.);
        ret_val = exp(-(doublereal)u) / sqrt(u) * (p1[0] + (p1[1] + (p1[2] + (
                p1[3] + p1[4] * v) * v) * v) * v) / (q1[0] + (q1[1] + (q1[2] 
                + (q1[3] + q1[4] * v) * v) * v) * v);
    } else if (v < (float)1.) {
        ret_val = (p2[0] + (p2[1] + (p2[2] + p2[3] * v) * v) * v) / (q2[0] + (
                q2[1] + (q2[2] + q2[3] * v) * v) * v);
    } else if (v < (float)4.) {
        ret_val = (p3[0] + (p3[1] + (p3[2] + p3[3] * v) * v) * v) / (q3[0] + (
                q3[1] + (q3[2] + q3[3] * v) * v) * v);
    } else if (v < (float)12.) {
        u = (float)1. / v;
        ret_val = (p4[0] + (p4[1] + (p4[2] + p4[3] * u) * u) * u) / (q4[0] + (
                q4[1] + (q4[2] + q4[3] * u) * u) * u);
    } else if (v < (float)50.) {
        u = (float)1. / v;
        ret_val = (p5[0] + (p5[1] + (p5[2] + p5[3] * u) * u) * u) / (q5[0] + (
                q5[1] + (q5[2] + q5[3] * u) * u) * u);
    } else if (v < (float)300.) {
        u = (float)1. / v;
        ret_val = (p6[0] + (p6[1] + (p6[2] + p6[3] * u) * u) * u) / (q6[0] + (
                q6[1] + (q6[2] + q6[3] * u) * u) * u);
    } else {
        u = (float)1. / (v - v * log(v) / (v + (float)1.));
        ret_val = (float)1. - (a2[0] + (a2[1] + a2[2] * u) * u) * u;
    }
    return ret_val;
} /* glands_ */

/* ** */
doublereal gvaviv_(real *rkappa, real *beta2, real *ran)
{
    /* Initialized data */

    static real fninv[5] = { (float)1.,(float).5,(float).33333333,(float).25,(
            float).2 };
    static real edgec[6] = { (float).16666667,(float).041666667,(float)
            .0083333333,(float).013888889,(float).0069444444,(float)
            7.7160493e-4 };
    static real u1[13] = { (float).25850868,(float).032477982,(float)
            -.0059020496,(float)0.,(float).024880692,(float).0047404356,(
            float)-7.444513e-4,(float).0073225731,(float)0.,(float)
            .0011668284,(float)0.,(float)-.0015727318,(float)-.0011210142 };
    static real u2[13] = { (float).43142611,(float).040797543,(float)
            -.0091490215,(float)0.,(float).042127077,(float).0073167928,(
            float)-.0014026047,(float).016195241,(float).0024714789,(float)
            .0020751278,(float)0.,(float)-.0025141668,(float)-.0014064022 };
    static real u3[13] = { (float).25225955,(float).064820468,(float)
            -.023615759,(float)0.,(float).023834176,(float).0021624675,(float)
            -.0026865597,(float)-.0054891384,(float).0039800522,(float)
            .0048447456,(float)-.0089439554,(float)-.0062756944,(float)
            -.0024655436 };
    static real u4[12] = { (float)1.2593231,(float)-.20374501,(float)
            .095055662,(float)-.020771531,(float)-.04686518,(float)
            -.0077222986,(float).0032241039,(float).008988292,(float)
            -.0067167236,(float)-.013049241,(float).018786468,(float)
            .014484097 };
    static real u5[13] = { (float)-.024864376,(float)-.0010368495,(float)
            .0014330117,(float)2.005273e-4,(float).0018751903,(float)
            .0012668869,(float)4.8736023e-4,(float).0034850854,(float)0.,(
            float)-3.6597173e-4,(float).0019372124,(float)7.0761825e-4,(float)
            4.6898375e-4 };
    static real u6[13] = { (float).035855696,(float)-.027542114,(float)
            .012631023,(float)-.0030188807,(float)-8.4479939e-4,(float)0.,(
            float)4.5675843e-4,(float)-.0069836141,(float).0039876546,(float)
            -.0036055679,(float)0.,(float).0015298434,(float).0019247256 };
    static real u7[8] = { (float)10.234691,(float)-3.5619655,(float).69387764,
            (float)-.14047599,(float)-1.995239,(float)-.45679694,(float)0.,(
            float).50505298 };
    static real u8[13] = { (float)21.487518,(float)-11.825253,(float)
            4.3133087,(float)-1.4500543,(float)-3.4343169,(float)-1.1063164,(
            float)-.21000819,(float)1.7891643,(float)-.89601916,(float)
            .39120793,(float).73410606,(float)0.,(float)-.32454506 };
    static real v1[12] = { (float).27827257,(float)-.0014227603,(float)
            .0024848327,(float)0.,(float).045091424,(float).0080559636,(float)
            -.0038974523,(float)0.,(float)-.0030634124,(float)7.5633702e-4,(
            float).0054730726,(float).0019792507 };
    static real v2[12] = { (float).41421789,(float)-.030061649,(float)
            .0052249697,(float)0.,(float).12693873,(float).022999801,(float)
            -.0086792801,(float).031875584,(float)-.0061757928,(float)0.,(
            float).019716857,(float).0032596742 };
    static real v3[13] = { (float).20191056,(float)-.046831422,(float)
            .0096777473,(float)-.0017995317,(float).053921588,(float)
            .003506874,(float)-.012621494,(float)-.0054996531,(float)
            -.0090029985,(float).0034958743,(float).018513506,(float)
            .0068332334,(float)-.0012940502 };
    static real v4[12] = { (float)1.3206081,(float).10036618,(float)
            -.022015201,(float).0061667091,(float)-.14986093,(float)
            -.012720568,(float).024972042,(float)-.0097751962,(float)
            .026087455,(float)-.011399062,(float)-.048282515,(float)
            -.0098552378 };
    static real v5[13] = { (float).016435243,(float).0360514,(float)
            .002303652,(float)-6.1666343e-4,(float)-.010775802,(float)
            .0051476061,(float).0056856517,(float)-.013438433,(float)0.,(
            float)0.,(float)-.0025421507,(float).0020169108,(float)
            -.0015144931 };
    static real v6[13] = { (float).033432405,(float).0060583916,(float)
            -.0023381379,(float)8.3846081e-4,(float)-.013346861,(float)
            -.0017402116,(float).0021052496,(float).0015528195,(float)
            .002190067,(float)-.0013202847,(float)-.0045124157,(float)
            -.0015629454,(float)2.2499176e-4 };
    static real v7[11] = { (float)5.4529572,(float)-.90906096,(float)
            .086122438,(float)0.,(float)-1.2218009,(float)-.3232412,(float)
            -.027373591,(float).12173464,(float)0.,(float)0.,(float)
            .040917471 };
    static real v8[11] = { (float)9.3841352,(float)-1.6276904,(float)
            .16571423,(float)0.,(float)-1.8160479,(float)-.50919193,(float)
            -.051384654,(float).21413992,(float)0.,(float)0.,(float)
            .066596366 };
    static real w1[13] = { (float).29712951,(float).0097572934,(float)0.,(
            float)-.0015291686,(float).035707399,(float).0096221631,(float)
            -.0018402821,(float)-.0049821585,(float).0018831112,(float)
            .0043541673,(float).0020301312,(float)-.0018723311,(float)
            -7.3403108e-4 };
    static real w2[11] = { (float).40882635,(float).014474912,(float)
            .0025023704,(float)-.0037707379,(float).18719727,(float)
            .056954987,(float)0.,(float).023020158,(float).0050574313,(float)
            .009455014,(float).019300232 };
    static real w3[13] = { (float).16861629,(float)0.,(float).0036317285,(
            float)-.0043657818,(float).030144338,(float).013891826,(float)
            -.0058030495,(float)-.0038717547,(float).0085359607,(float)
            .014507659,(float).0082387775,(float)-.010116105,(float)
            -.005513567 };
    static real w4[13] = { (float)1.3493891,(float)-.0026863185,(float)
            -.003521604,(float).024434909,(float)-.083447911,(float)
            -.04806136,(float).0076473951,(float).02449443,(float)-.0162092,(
            float)-.037768479,(float)-.047890063,(float).017778596,(float)
            .013179324 };
    static real w5[13] = { (float).10264945,(float).032738857,(float)0.,(
            float).0043608779,(float)-.043097757,(float)-.0022647176,(float)
            .009453129,(float)-.012442571,(float)-.0032283517,(float)
            -.0075640352,(float)-.0088293329,(float).0052537299,(float)
            .0013340546 };
    static real w6[13] = { (float).029568177,(float)-.001630006,(float)
            -2.1119745e-4,(float).0023599053,(float)-.0048515387,(float)
            -.0040797531,(float)4.0403265e-4,(float).0018200105,(float)
            -.0014346306,(float)-.0039165276,(float)-.0037432073,(float)
            .001995038,(float).0012222675 };
    static real w8[8] = { (float)6.6184645,(float)-.73866379,(float)
            .044693973,(float)0.,(float)-1.4540925,(float)-.39529833,(float)
            -.044293243,(float).088741049 };

    /* System generated locals */
    integer i__1;
    real ret_val, r__1;

    /* Builtin functions */
    double sqrt(doublereal), log(doublereal), exp(doublereal);

    /* Local variables */
    static real alfa[4], rlam, h[9];
    static integer j, k, n;
    static real s, t, x, y, dsigm[5];
    static integer itype;
    static real p2, p3, q2, q3, x2, x3, y2, ac[13], y3, s0, hc[9], fl, fn;
    extern doublereal glande_(real *);
    static real fu, pq, wk;
    extern doublereal glands_(real *);
    static real xx, xy, yy, drk[5];
    static integer npt;

/* . */
/* .    ******************************************************************
 */
/* .    *                                                                *
 */
/* .    *                                                                *
 */
/* .    *  GVAVIV                                                        *
 */
/* .    *  ------                                                        *
 */
/* .    *                                                                *
 */
/* .    *  Initializes the parameters for a Vavilov distribution         *
 */
/* .    *                                                                *
 */
/* .    *  BETA2 is the particle velocity squared                        *
 */
/* .    *  RKAPPA the usual straggling parameter                         *
 */
/* .    *                                                                *
 */
/* .    *  This routine has been extracted from a proposed CERNLIB       *
 */
/* .    *  set of routines (VAVCOE,VAVFSM).                              *
 */
/* .    *  The authors of these routines have submitted an article       *
 */
/* .    *  to NIM:                                                       *
 */
/* .    *                                                                *
 */
/* .    *    Alberto Rotondi, Paolo Montagna                             *
 */
/* .    *    Fast calculation of Vavilov distribution NIMB 1990          *
 */
/* .    *                                                                *
 */
/* .    *       Author    K.S.Koelbig     ********                       *
 */
/* .    *                                                                *
 */
/* .    ******************************************************************
 */
/* . */
    ret_val = (float)0.;
    if (*rkappa < (float).01 || *rkappa > (float)12.) {
        return ret_val;
    }
    if (*rkappa >= (float).29) {
        itype = 1;
        npt = 100;
        wk = 1 / sqrt(*rkappa);
        ac[0] = (*beta2 * (float)-.032227 - (float).074275) * *rkappa + (*
                beta2 * (float).24533 + (float).070152) * wk + (*beta2 * (
                float)-.5561 - (float)3.1579);
        ac[8] = (*beta2 * (float)-.013483 - (float).048801) * *rkappa + (*
                beta2 * (float)-1.6921 + (float)8.3656) * wk + (*beta2 * (
                float)-.73275 - (float)3.5226);
/* Computing 2nd power */
        r__1 = wk;
        drk[0] = r__1 * r__1;
        dsigm[0] = sqrt(*rkappa / (1 - *beta2 * (float).5));
        for (j = 1; j <= 4; ++j) {
            drk[j] = drk[0] * drk[j - 1];
            dsigm[j] = dsigm[0] * dsigm[j - 1];
/* L1: */
            alfa[j - 1] = (fninv[j - 1] - *beta2 * fninv[j]) * drk[j - 1];
        }
        hc[0] = log(*rkappa) + *beta2 + (float).42278434;
        hc[1] = dsigm[0];
        hc[2] = alfa[1] * dsigm[2];
/* Computing 2nd power */
        r__1 = alfa[0];
        hc[3] = (r__1 * r__1 * 3 + alfa[2]) * dsigm[3] - 3;
        hc[4] = (alfa[0] * 10 * alfa[1] + alfa[3]) * dsigm[4] - hc[2] * 10;
/* Computing 2nd power */
        r__1 = hc[2];
        hc[5] = r__1 * r__1;
        hc[6] = hc[2] * hc[3];
        hc[7] = hc[2] * hc[5];
        for (j = 2; j <= 7; ++j) {
/* L2: */
            hc[j] = edgec[j - 2] * hc[j];
        }
        hc[8] = hc[1] * (float).39894228;
    } else if (*rkappa >= (float).22) {
        itype = 2;
        npt = 150;
        x = (*rkappa - (float).3) * (float)25.000000000000004 + 1;
        y = (sqrt(*beta2) - (float)1.) * (float)2.1052631578947367 + 1;
        xx = x * 2;
        yy = y * 2;
        x2 = xx * x - 1;
        x3 = xx * x2 - x;
        y2 = yy * y - 1;
        y3 = yy * y2 - y;
        xy = x * y;
        p2 = x2 * y;
        p3 = x3 * y;
        q2 = y2 * x;
        q3 = y3 * x;
        pq = x2 * y2;
        ac[1] = w1[0] + w1[1] * x + w1[3] * x3 + w1[4] * y + w1[5] * y2 + w1[
                6] * y3 + w1[7] * xy + w1[8] * p2 + w1[9] * p3 + w1[10] * q2 
                + w1[11] * q3 + w1[12] * pq;
        ac[2] = w2[0] + w2[1] * x + w2[2] * x2 + w2[3] * x3 + w2[4] * y + w2[
                5] * y2 + w2[7] * xy + w2[8] * p2 + w2[9] * p3 + w2[10] * q2;
        ac[3] = w3[0] + w3[2] * x2 + w3[3] * x3 + w3[4] * y + w3[5] * y2 + w3[
                6] * y3 + w3[7] * xy + w3[8] * p2 + w3[9] * p3 + w3[10] * q2 
                + w3[11] * q3 + w3[12] * pq;
        ac[4] = w4[0] + w4[1] * x + w4[2] * x2 + w4[3] * x3 + w4[4] * y + w4[
                5] * y2 + w4[6] * y3 + w4[7] * xy + w4[8] * p2 + w4[9] * p3 + 
                w4[10] * q2 + w4[11] * q3 + w4[12] * pq;
        ac[5] = w5[0] + w5[1] * x + w5[3] * x3 + w5[4] * y + w5[5] * y2 + w5[
                6] * y3 + w5[7] * xy + w5[8] * p2 + w5[9] * p3 + w5[10] * q2 
                + w5[11] * q3 + w5[12] * pq;
        ac[6] = w6[0] + w6[1] * x + w6[2] * x2 + w6[3] * x3 + w6[4] * y + w6[
                5] * y2 + w6[6] * y3 + w6[7] * xy + w6[8] * p2 + w6[9] * p3 + 
                w6[10] * q2 + w6[11] * q3 + w6[12] * pq;
        ac[8] = w8[0] + w8[1] * x + w8[2] * x2 + w8[4] * y + w8[5] * y2 + w8[
                6] * y3 + w8[7] * xy;
        ac[0] = (float)-3.05;
    } else if (*rkappa >= (float).12) {
        itype = 3;
        npt = 200;
        x = (*rkappa - (float).2) * (float)24.999999999999996 + 1;
        y = (sqrt(*beta2) - (float)1.) * (float)2.1052631578947367 + 1;
        xx = x * 2;
        yy = y * 2;
        x2 = xx * x - 1;
        x3 = xx * x2 - x;
        y2 = yy * y - 1;
        y3 = yy * y2 - y;
        xy = x * y;
        p2 = x2 * y;
        p3 = x3 * y;
        q2 = y2 * x;
        q3 = y3 * x;
        pq = x2 * y2;
        ac[1] = v1[0] + v1[1] * x + v1[2] * x2 + v1[4] * y + v1[5] * y2 + v1[
                6] * y3 + v1[8] * p2 + v1[9] * p3 + v1[10] * q2 + v1[11] * q3;
        ac[2] = v2[0] + v2[1] * x + v2[2] * x2 + v2[4] * y + v2[5] * y2 + v2[
                6] * y3 + v2[7] * xy + v2[8] * p2 + v2[10] * q2 + v2[11] * q3;
        ac[3] = v3[0] + v3[1] * x + v3[2] * x2 + v3[3] * x3 + v3[4] * y + v3[
                5] * y2 + v3[6] * y3 + v3[7] * xy + v3[8] * p2 + v3[9] * p3 + 
                v3[10] * q2 + v3[11] * q3 + v3[12] * pq;
        ac[4] = v4[0] + v4[1] * x + v4[2] * x2 + v4[3] * x3 + v4[4] * y + v4[
                5] * y2 + v4[6] * y3 + v4[7] * xy + v4[8] * p2 + v4[9] * p3 + 
                v4[10] * q2 + v4[11] * q3;
        ac[5] = v5[0] + v5[1] * x + v5[2] * x2 + v5[3] * x3 + v5[4] * y + v5[
                5] * y2 + v5[6] * y3 + v5[7] * xy + v5[10] * q2 + v5[11] * q3 
                + v5[12] * pq;
        ac[6] = v6[0] + v6[1] * x + v6[2] * x2 + v6[3] * x3 + v6[4] * y + v6[
                5] * y2 + v6[6] * y3 + v6[7] * xy + v6[8] * p2 + v6[9] * p3 + 
                v6[10] * q2 + v6[11] * q3 + v6[12] * pq;
        ac[7] = v7[0] + v7[1] * x + v7[2] * x2 + v7[4] * y + v7[5] * y2 + v7[
                6] * y3 + v7[7] * xy + v7[10] * q2;
        ac[8] = v8[0] + v8[1] * x + v8[2] * x2 + v8[4] * y + v8[5] * y2 + v8[
                6] * y3 + v8[7] * xy + v8[10] * q2;
        ac[0] = (float)-3.04;
    } else {
        itype = 4;
        if (*rkappa >= (float).02) {
            itype = 3;
        }
        npt = 200;
        x = (*rkappa - (float).1) * (float)25. + 1;
        y = (sqrt(*beta2) - (float)1.) * (float)2.1052631578947367 + 1;
        xx = x * 2;
        yy = y * 2;
        x2 = xx * x - 1;
        x3 = xx * x2 - x;
        y2 = yy * y - 1;
        y3 = yy * y2 - y;
        xy = x * y;
        p2 = x2 * y;
        p3 = x3 * y;
        q2 = y2 * x;
        q3 = y3 * x;
        pq = x2 * y2;
        if (itype == 4) {
            goto L4;
        }
        ac[1] = u1[0] + u1[1] * x + u1[2] * x2 + u1[4] * y + u1[5] * y2 + u1[
                6] * y3 + u1[7] * xy + u1[9] * p3 + u1[11] * q3 + u1[12] * pq;
        ac[2] = u2[0] + u2[1] * x + u2[2] * x2 + u2[4] * y + u2[5] * y2 + u2[
                6] * y3 + u2[7] * xy + u2[8] * p2 + u2[9] * p3 + u2[11] * q3 
                + u2[12] * pq;
        ac[3] = u3[0] + u3[1] * x + u3[2] * x2 + u3[4] * y + u3[5] * y2 + u3[
                6] * y3 + u3[7] * xy + u3[8] * p2 + u3[9] * p3 + u3[10] * q2 
                + u3[11] * q3 + u3[12] * pq;
        ac[4] = u4[0] + u4[1] * x + u4[2] * x2 + u4[3] * x3 + u4[4] * y + u4[
                5] * y2 + u4[6] * y3 + u4[7] * xy + u4[8] * p2 + u4[9] * p3 + 
                u4[10] * q2 + u4[11] * q3;
        ac[5] = u5[0] + u5[1] * x + u5[2] * x2 + u5[3] * x3 + u5[4] * y + u5[
                5] * y2 + u5[6] * y3 + u5[7] * xy + u5[9] * p3 + u5[10] * q2 
                + u5[11] * q3 + u5[12] * pq;
        ac[6] = u6[0] + u6[1] * x + u6[2] * x2 + u6[3] * x3 + u6[4] * y + u6[
                6] * y3 + u6[7] * xy + u6[8] * p2 + u6[9] * p3 + u6[11] * q3 
                + u6[12] * pq;
L4:
        ac[7] = u7[0] + u7[1] * x + u7[2] * x2 + u7[3] * x3 + u7[4] * y + u7[
                5] * y2 + u7[7] * xy;
        ac[8] = u8[0] + u8[1] * x + u8[2] * x2 + u8[3] * x3 + u8[4] * y + u8[
                5] * y2 + u8[6] * y3 + u8[7] * xy + u8[8] * p2 + u8[9] * p3 + 
                u8[10] * q2 + u8[12] * pq;
        ac[0] = (float)-3.03;
    }
    ac[9] = (ac[8] - ac[0]) / npt;
    if (itype == 3) {
        x = (ac[7] - ac[8]) / (ac[7] * ac[8]);
        y = 1 / log(ac[8] / ac[7]);
/* Computing 2nd power */
        r__1 = ac[7];
        p2 = r__1 * r__1;
        ac[11] = p2 * (ac[1] * exp(-(doublereal)ac[2] * (ac[7] + ac[5] * p2) 
                - ac[3] * exp(-(doublereal)ac[4] * (ac[7] + ac[6] * p2))) - y 
                * (float).045 / ac[7]) / (x * y * ac[7] + 1);
        ac[12] = (x * ac[11] + (float).045) * y;
    }
    if (itype == 4) {
        ac[10] = (float).995 / glands_(&ac[8]);
    }
    t = *ran * 2 / ac[9];
    rlam = ac[0];
    fl = (float)0.;
    s = (float)0.;
    i__1 = npt;
    for (n = 1; n <= i__1; ++n) {
        rlam += ac[9];
        if (itype == 1) {
            fn = (float)1.;
            x = (rlam + hc[0]) * hc[1];
            h[0] = x;
/* Computing 2nd power */
            r__1 = x;
            h[1] = r__1 * r__1 - 1;
            for (k = 2; k <= 8; ++k) {
                fn += 1;
/* L31: */
                h[k] = x * h[k - 1] - fn * h[k - 2];
            }
            y = hc[7] * h[8] + 1;
            for (k = 2; k <= 6; ++k) {
/* L32: */
                y += hc[k] * h[k];
            }
/* Computing 2nd power */
            r__1 = x;
            fu = hc[8] * exp(r__1 * r__1 * (float)-.5) * dmax(y,(float)0.);
        } else if (itype == 2) {
/* Computing 2nd power */
            r__1 = rlam;
            x = r__1 * r__1;
            fu = ac[1] * exp(-(doublereal)ac[2] * (rlam + ac[5] * x) - ac[3] *
                     exp(-(doublereal)ac[4] * (rlam + ac[6] * x)));
        } else if (itype == 3) {
            if (rlam < ac[7]) {
/* Computing 2nd power */
                r__1 = rlam;
                x = r__1 * r__1;
                fu = ac[1] * exp(-(doublereal)ac[2] * (rlam + ac[5] * x) - ac[
                        3] * exp(-(doublereal)ac[4] * (rlam + ac[6] * x)));
            } else {
                x = 1 / rlam;
                fu = (ac[11] * x + ac[12]) * x;
            }
        } else {
            fu = ac[10] * glande_(&rlam);
        }
        s = s + fl + fu;
        if (s > t) {
            goto L22;
        }
/* L21: */
        fl = fu;
    }
L22:
    s0 = s - fl - fu;
    ret_val = rlam - ac[9];
    if (s > s0) {
        ret_val += ac[9] * (t - s0) / (s - s0);
    }
    return ret_val;
} /* gvaviv_ */

#ifdef __cplusplus
        }
#endif

---