Changes in / [c1bccff:925b3b5] in sasmodels


Ignore:
Files:
1 deleted
32 edited

Legend:

Unmodified
Added
Removed

  • .gitignore

    re9ed2de r8a5f021  
    2222/example/Fit_*/ 
    2323/example/batch_fit.csv 
    24 # Note: changes to gauss20, gauss76 and gauss150 are still tracked since they 
    25 # are part of the repo and required by some models. 
    26 /sasmodels/models/lib/gauss*.c 

  • doc/guide/plugin.rst

    rd0dc9a3 r3048ec6  
    543543    M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E: 
    544544        $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$ 
    545     exp, log, pow(x,y), expm1, log1p, sqrt, cbrt: 
    546         Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\ln 1 + x$, 
    547         $\sqrt{x}$, $\sqrt[3]{x}$. The functions expm1(x) and log1p(x) 
    548         are accurate across all $x$, including $x$ very close to zero. 
     545    exp, log, pow(x,y), expm1, sqrt: 
     546        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\sqrt{x}$. 
     547        The function expm1(x) is accurate across all $x$, including $x$ 
     548        very close to zero. 
    549549    sin, cos, tan, asin, acos, atan: 
    550550        Trigonometry functions and inverses, operating on radians. 
    … …  
    557557        atan(y/x) would return a value in quadrant I. Similarly for 
    558558        quadrants II and IV when $x$ and $y$ have opposite sign. 
    559     fabs(x), fmin(x,y), fmax(x,y), trunc, rint: 
     559    fmin(x,y), fmax(x,y), trunc, rint: 
    560560        Floating point functions.  rint(x) returns the nearest integer. 
    561561    NAN: 
    562562        NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that 
    563563        you cannot use :code:`x == NAN` to test for NaN values since that 
    564         will always return false.  NAN does not equal NAN!  The alternative, 
    565         :code:`x != x` may fail if the compiler optimizes the test away. 
     564        will always return false.  NAN does not equal NAN! 
    566565    INFINITY: 
    567566        $\infty, 1/0$.  Use isinf(x) to test for infinity, or isfinite(x) 
    … …  
    735734        Similar arrays are available in :code:`gauss20.c` for 20-point 
    736735        quadrature and in :code:`gauss150.c` for 150-point quadrature. 
    737         The macros :code:`GAUSS_N`, :code:`GAUSS_Z` and :code:`GAUSS_W` are 
    738         defined so that you can change the order of the integration by 
    739         selecting an different source without touching the C code. 
    740736 
    741737        :code:`source = ["lib/gauss76.c", ...]` 

  • sasmodels/compare.py

    r2d81cfe r2d81cfe  
    4242from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    4343from .direct_model import DirectModel, get_mesh 
    44 from .generate import FLOAT_RE, set_integration_size 
     44from .generate import FLOAT_RE 
    4545from .weights import plot_weights 
    4646 
    … …  
    706706    return data, index 
    707707 
    708 def make_engine(model_info, data, dtype, cutoff, ngauss=0): 
     708def make_engine(model_info, data, dtype, cutoff): 
    709709    # type: (ModelInfo, Data, str, float) -> Calculator 
    710710    """ 
    … …  
    714714    than OpenCL. 
    715715    """ 
    716     if ngauss: 
    717         set_integration_size(model_info, ngauss) 
    718  
    719716    if dtype is None or not dtype.endswith('!'): 
    720717        return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff) 
    … …  
    957954    'poly', 'mono', 'cutoff=', 
    958955    'magnetic', 'nonmagnetic', 
    959     'accuracy=', 'ngauss=', 
     956    'accuracy=', 
    960957    'neval=',  # for timing... 
    961958 
    … …  
    10921089        'show_weights' : False, 
    10931090        'sphere'    : 0, 
    1094         'ngauss'    : '0', 
    10951091    } 
    10961092    for arg in flags: 
    … …  
    11191115        elif arg.startswith('-engine='):   opts['engine'] = arg[8:] 
    11201116        elif arg.startswith('-neval='):    opts['count'] = arg[7:] 
    1121         elif arg.startswith('-ngauss='):   opts['ngauss'] = arg[8:] 
    11221117        elif arg.startswith('-random='): 
    11231118            opts['seed'] = int(arg[8:]) 
    … …  
    11741169 
    11751170    comparison = any(PAR_SPLIT in v for v in values) 
    1176  
    11771171    if PAR_SPLIT in name: 
    11781172        names = name.split(PAR_SPLIT, 2) 
    … …  
    11871181        return None 
    11881182 
    1189     if PAR_SPLIT in opts['ngauss']: 
    1190         opts['ngauss'] = [int(k) for k in opts['ngauss'].split(PAR_SPLIT, 2)] 
    1191         comparison = True 
    1192     else: 
    1193         opts['ngauss'] = [int(opts['ngauss'])]*2 
    1194  
    11951183    if PAR_SPLIT in opts['engine']: 
    11961184        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2) 
    … …  
    12111199        opts['cutoff'] = [float(opts['cutoff'])]*2 
    12121200 
    1213     base = make_engine(model_info[0], data, opts['engine'][0], 
    1214                        opts['cutoff'][0], opts['ngauss'][0]) 
     1201    base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0]) 
    12151202    if comparison: 
    1216         comp = make_engine(model_info[1], data, opts['engine'][1], 
    1217                            opts['cutoff'][1], opts['ngauss'][1]) 
     1203        comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1]) 
    12181204    else: 
    12191205        comp = None 

  • sasmodels/generate.py

    r2d81cfe r2d81cfe  
    270270""" 
    271271 
    272  
    273 def set_integration_size(info, n): 
    274     # type: (ModelInfo, int) -> None 
    275     """ 
    276     Update the model definition, replacing the gaussian integration with 
    277     a gaussian integration of a different size. 
    278  
    279     Note: this really ought to be a method in modelinfo, but that leads to 
    280     import loops. 
    281     """ 
    282     if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)): 
    283         import os.path 
    284         from .gengauss import gengauss 
    285         path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n) 
    286         if not os.path.exists(path): 
    287             gengauss(n, path) 
    288         info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss') 
    289                         else lib for lib in info.source] 
    290272 
    291273def format_units(units): 

  • sasmodels/models/barbell.c

    r74768cb rbecded3  
    2323    const double qab_r = radius_bell*qab; // Q*R*sin(theta) 
    2424    double total = 0.0; 
    25     for (int i = 0; i < GAUSS_N; i++){ 
    26         const double t = GAUSS_Z[i]*zm + zb; 
     25    for (int i = 0; i < 76; i++){ 
     26        const double t = Gauss76Z[i]*zm + zb; 
    2727        const double radical = 1.0 - t*t; 
    2828        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    2929        const double Fq = cos(m*t + b) * radical * bj; 
    30         total += GAUSS_W[i] * Fq; 
     30        total += Gauss76Wt[i] * Fq; 
    3131    } 
    3232    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    7373    const double zb = M_PI_4; 
    7474    double total = 0.0; 
    75     for (int i = 0; i < GAUSS_N; i++){ 
    76         const double alpha= GAUSS_Z[i]*zm + zb; 
     75    for (int i = 0; i < 76; i++){ 
     76        const double alpha= Gauss76Z[i]*zm + zb; 
    7777        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    7878        SINCOS(alpha, sin_alpha, cos_alpha); 
    7979        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    80         total += GAUSS_W[i] * Aq * Aq * sin_alpha; 
     80        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8181    } 
    8282    // translate dx in [-1,1] to dx in [lower,upper] 

  • sasmodels/models/bcc_paracrystal.c

    r74768cb rea60e08  
    8181 
    8282    double outer_sum = 0.0; 
    83     for(int i=0; i<GAUSS_N; i++) { 
     83    for(int i=0; i<150; i++) { 
    8484        double inner_sum = 0.0; 
    85         const double theta = GAUSS_Z[i]*theta_m + theta_b; 
     85        const double theta = Gauss150Z[i]*theta_m + theta_b; 
    8686        double sin_theta, cos_theta; 
    8787        SINCOS(theta, sin_theta, cos_theta); 
    8888        const double qc = q*cos_theta; 
    8989        const double qab = q*sin_theta; 
    90         for(int j=0;j<GAUSS_N;j++) { 
    91             const double phi = GAUSS_Z[j]*phi_m + phi_b; 
     90        for(int j=0;j<150;j++) { 
     91            const double phi = Gauss150Z[j]*phi_m + phi_b; 
    9292            double sin_phi, cos_phi; 
    9393            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    9595            const double qb = qab*sin_phi; 
    9696            const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 
    97             inner_sum += GAUSS_W[j] * form; 
     97            inner_sum += Gauss150Wt[j] * form; 
    9898        } 
    9999        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    100         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     100        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
    101101    } 
    102102    outer_sum *= theta_m; 

  • sasmodels/models/capped_cylinder.c

    r74768cb rbecded3  
    3030    const double qab_r = radius_cap*qab; // Q*R*sin(theta) 
    3131    double total = 0.0; 
    32     for (int i=0; i<GAUSS_N; i++) { 
    33         const double t = GAUSS_Z[i]*zm + zb; 
     32    for (int i=0; i<76 ;i++) { 
     33        const double t = Gauss76Z[i]*zm + zb; 
    3434        const double radical = 1.0 - t*t; 
    3535        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    3636        const double Fq = cos(m*t + b) * radical * bj; 
    37         total += GAUSS_W[i] * Fq; 
     37        total += Gauss76Wt[i] * Fq; 
    3838    } 
    3939    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    9595    const double zb = M_PI_4; 
    9696    double total = 0.0; 
    97     for (int i=0; i<GAUSS_N ;i++) { 
    98         const double theta = GAUSS_Z[i]*zm + zb; 
     97    for (int i=0; i<76 ;i++) { 
     98        const double theta = Gauss76Z[i]*zm + zb; 
    9999        double sin_theta, cos_theta; // slots to hold sincos function output 
    100100        SINCOS(theta, sin_theta, cos_theta); 
    … …  
    103103        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
    104104        // scale by sin_theta for spherical coord integration 
    105         total += GAUSS_W[i] * Aq * Aq * sin_theta; 
     105        total += Gauss76Wt[i] * Aq * Aq * sin_theta; 
    106106    } 
    107107    // translate dx in [-1,1] to dx in [lower,upper] 

  • sasmodels/models/core_shell_bicelle.c

    r74768cb rbecded3  
    5252 
    5353    double total = 0.0; 
    54     for(int i=0;i<GAUSS_N;i++) { 
    55         double theta = (GAUSS_Z[i] + 1.0)*uplim; 
     54    for(int i=0;i<N_POINTS_76;i++) { 
     55        double theta = (Gauss76Z[i] + 1.0)*uplim; 
    5656        double sin_theta, cos_theta; // slots to hold sincos function output 
    5757        SINCOS(theta, sin_theta, cos_theta); 
    5858        double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
    5959                                   halflength, sld_core, sld_face, sld_rim, sld_solvent); 
    60         total += GAUSS_W[i]*fq*fq*sin_theta; 
     60        total += Gauss76Wt[i]*fq*fq*sin_theta; 
    6161    } 
    6262 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r74768cb r82592da  
    3737    //initialize integral 
    3838    double outer_sum = 0.0; 
    39     for(int i=0;i<GAUSS_N;i++) { 
     39    for(int i=0;i<76;i++) { 
    4040        //setup inner integral over the ellipsoidal cross-section 
    4141        //const double va = 0.0; 
    4242        //const double vb = 1.0; 
    43         //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
    44         const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0; 
     43        //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     44        const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0; 
    4545        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
    4646        const double qab = q*sin_theta; 
    … …  
    4949        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
    5050        double inner_sum=0.0; 
    51         for(int j=0;j<GAUSS_N;j++) { 
     51        for(int j=0;j<76;j++) { 
    5252            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    5353            // inner integral limits 
    5454            //const double vaj=0.0; 
    5555            //const double vbj=M_PI; 
    56             //const double phi = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    57             const double phi = ( GAUSS_Z[j] +1.0)*M_PI_2; 
     56            //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     57            const double phi = ( Gauss76Z[j] +1.0)*M_PI_2; 
    5858            const double rr = sqrt(r2A - r2B*cos(phi)); 
    5959            const double be1 = sas_2J1x_x(rr*qab); 
    … …  
    6161            const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
    6262 
    63             inner_sum += GAUSS_W[j] * fq * fq; 
     63            inner_sum += Gauss76Wt[j] * fq * fq; 
    6464        } 
    6565        //now calculate outer integral 
    66         outer_sum += GAUSS_W[i] * inner_sum; 
     66        outer_sum += Gauss76Wt[i] * inner_sum; 
    6767    } 
    6868 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r74768cb r82592da  
    77        double length) 
    88{ 
    9     return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length + 
     9    return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length +  
    1010                 square(r_minor)*x_core*2.0*thick_face  ); 
    1111} 
    … …  
    4747    //initialize integral 
    4848    double outer_sum = 0.0; 
    49     for(int i=0;i<GAUSS_N;i++) { 
     49    for(int i=0;i<76;i++) { 
    5050        //setup inner integral over the ellipsoidal cross-section 
    5151        // since we generate these lots of times, why not store them somewhere? 
    52         //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
    53         const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0; 
     52        //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     53        const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0; 
    5454        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    5555        double inner_sum=0; 
    … …  
    5858        si1 = sas_sinx_x(sinarg1); 
    5959        si2 = sas_sinx_x(sinarg2); 
    60         for(int j=0;j<GAUSS_N;j++) { 
     60        for(int j=0;j<76;j++) { 
    6161            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    62             //const double beta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    63             const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2; 
     62            //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     63            const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
    6464            const double rr = sqrt(r2A - r2B*cos(beta)); 
    6565            double besarg1 = q*rr*sin_alpha; 
    … …  
    6767            be1 = sas_2J1x_x(besarg1); 
    6868            be2 = sas_2J1x_x(besarg2); 
    69             inner_sum += GAUSS_W[j] *square(dr1*si1*be1 + 
     69            inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 + 
    7070                                              dr2*si1*be2 + 
    7171                                              dr3*si2*be1); 
    7272        } 
    7373        //now calculate outer integral 
    74         outer_sum += GAUSS_W[i] * inner_sum; 
     74        outer_sum += Gauss76Wt[i] * inner_sum; 
    7575    } 
    7676 

  • sasmodels/models/core_shell_cylinder.c

    r74768cb rbecded3  
    3030    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    3131    double total = 0.0; 
    32     for (int i=0; i<GAUSS_N ;i++) { 
     32    for (int i=0; i<76 ;i++) { 
    3333        // translate a point in [-1,1] to a point in [0, pi/2] 
    34         //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0; 
     34        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    3535        double sin_theta, cos_theta; 
    36         const double theta = GAUSS_Z[i]*M_PI_4 + M_PI_4; 
     36        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 
    3737        SINCOS(theta, sin_theta,  cos_theta); 
    3838        const double qab = q*sin_theta; 
    … …  
    4040        const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
    4141            + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
    42         total += GAUSS_W[i] * fq * fq * sin_theta; 
     42        total += Gauss76Wt[i] * fq * fq * sin_theta; 
    4343    } 
    4444    // translate dx in [-1,1] to dx in [lower,upper] 

  • sasmodels/models/core_shell_ellipsoid.c

    r74768cb rbecded3  
    5959    const double b = 0.5; 
    6060    double total = 0.0;     //initialize intergral 
    61     for(int i=0;i<GAUSS_N;i++) { 
    62         const double cos_theta = GAUSS_Z[i]*m + b; 
     61    for(int i=0;i<76;i++) { 
     62        const double cos_theta = Gauss76Z[i]*m + b; 
    6363        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
    6464        double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta, 
    … …  
    6666            equat_shell, polar_shell, 
    6767            sld_core_shell, sld_shell_solvent); 
    68         total += GAUSS_W[i] * fq * fq; 
     68        total += Gauss76Wt[i] * fq * fq; 
    6969    } 
    7070    total *= m; 

  • sasmodels/models/core_shell_parallelepiped.c

    r4493288 r4493288  
    6060    // outer integral (with gauss points), integration limits = 0, 1 
    6161    double outer_sum = 0; //initialize integral 
    62     for( int i=0; i<GAUSS_N; i++) { 
    63         const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
     62    for( int i=0; i<76; i++) { 
     63        const double cos_alpha = 0.5 * ( Gauss76Z[i] + 1.0 ); 
    6464        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
    6565 
    … …  
    6868        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 
    6969        double inner_sum = 0.0; 
    70         for(int j=0; j<GAUSS_N; j++) { 
    71             const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     70        for(int j=0; j<76; j++) { 
     71            const double beta = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    7272            double sin_beta, cos_beta; 
    7373            SINCOS(M_PI_2*beta, sin_beta, cos_beta); 
    … …  
    8989#endif 
    9090 
    91             inner_sum += GAUSS_W[j] * f * f; 
     91            inner_sum += Gauss76Wt[j] * f * f; 
    9292        } 
    9393        inner_sum *= 0.5; 
    9494        // now sum up the outer integral 
    95         outer_sum += GAUSS_W[i] * inner_sum; 
     95        outer_sum += Gauss76Wt[i] * inner_sum; 
    9696    } 
    9797    outer_sum *= 0.5; 

  • sasmodels/models/cylinder.c

    r74768cb rbecded3  
    2121 
    2222    double total = 0.0; 
    23     for (int i=0; i<GAUSS_N ;i++) { 
    24         const double theta = GAUSS_Z[i]*zm + zb; 
     23    for (int i=0; i<76 ;i++) { 
     24        const double theta = Gauss76Z[i]*zm + zb; 
    2525        double sin_theta, cos_theta; // slots to hold sincos function output 
    2626        // theta (theta,phi) the projection of the cylinder on the detector plane 
    2727        SINCOS(theta , sin_theta, cos_theta); 
    2828        const double form = fq(q*sin_theta, q*cos_theta, radius, length); 
    29         total += GAUSS_W[i] * form * form * sin_theta; 
     29        total += Gauss76Wt[i] * form * form * sin_theta; 
    3030    } 
    3131    // translate dx in [-1,1] to dx in [lower,upper] 

  • sasmodels/models/ellipsoid.c

    r74768cb rbecded3  
    2222 
    2323    // translate a point in [-1,1] to a point in [0, 1] 
    24     // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
     24    // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2; 
    2525    const double zm = 0.5; 
    2626    const double zb = 0.5; 
    2727    double total = 0.0; 
    28     for (int i=0;i<GAUSS_N;i++) { 
    29         const double u = GAUSS_Z[i]*zm + zb; 
     28    for (int i=0;i<76;i++) { 
     29        const double u = Gauss76Z[i]*zm + zb; 
    3030        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 
    3131        const double f = sas_3j1x_x(q*r); 
    32         total += GAUSS_W[i] * f * f; 
     32        total += Gauss76Wt[i] * f * f; 
    3333    } 
    3434    // translate dx in [-1,1] to dx in [lower,upper] 

  • sasmodels/models/elliptical_cylinder.c

    r74768cb r82592da  
    2222    //initialize integral 
    2323    double outer_sum = 0.0; 
    24     for(int i=0;i<GAUSS_N;i++) { 
     24    for(int i=0;i<76;i++) { 
    2525        //setup inner integral over the ellipsoidal cross-section 
    26         const double cos_val = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
     26        const double cos_val = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    2727        const double sin_val = sqrt(1.0 - cos_val*cos_val); 
    2828        //const double arg = radius_minor*sin_val; 
    2929        double inner_sum=0; 
    30         for(int j=0;j<GAUSS_N;j++) { 
    31             const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     30        for(int j=0;j<76;j++) { 
     31            //20 gauss points for the inner integral, increase to 76, RKH 6Nov2017 
     32            const double theta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    3233            const double r = sin_val*sqrt(rA - rB*cos(theta)); 
    3334            const double be = sas_2J1x_x(q*r); 
    34             inner_sum += GAUSS_W[j] * be * be; 
     35            inner_sum += Gauss76Wt[j] * be * be; 
    3536        } 
    3637        //now calculate the value of the inner integral 
    … …  
    3940        //now calculate outer integral 
    4041        const double si = sas_sinx_x(q*0.5*length*cos_val); 
    41         outer_sum += GAUSS_W[i] * inner_sum * si * si; 
     42        outer_sum += Gauss76Wt[i] * inner_sum * si * si; 
    4243    } 
    4344    outer_sum *= 0.5*(vb-va); 

  • sasmodels/models/elliptical_cylinder.py

    r2d81cfe r2d81cfe  
    121121# pylint: enable=bad-whitespace, line-too-long 
    122122 
    123 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
     123source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "lib/gauss20.c", 
     124          "elliptical_cylinder.c"] 
    124125 
    125126demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 

  • sasmodels/models/fcc_paracrystal.c

    r74768cb rf728001  
    5353 
    5454    double outer_sum = 0.0; 
    55     for(int i=0; i<GAUSS_N; i++) { 
     55    for(int i=0; i<150; i++) { 
    5656        double inner_sum = 0.0; 
    57         const double theta = GAUSS_Z[i]*theta_m + theta_b; 
     57        const double theta = Gauss150Z[i]*theta_m + theta_b; 
    5858        double sin_theta, cos_theta; 
    5959        SINCOS(theta, sin_theta, cos_theta); 
    6060        const double qc = q*cos_theta; 
    6161        const double qab = q*sin_theta; 
    62         for(int j=0;j<GAUSS_N;j++) { 
    63             const double phi = GAUSS_Z[j]*phi_m + phi_b; 
     62        for(int j=0;j<150;j++) { 
     63            const double phi = Gauss150Z[j]*phi_m + phi_b; 
    6464            double sin_phi, cos_phi; 
    6565            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    6767            const double qb = qab*sin_phi; 
    6868            const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 
    69             inner_sum += GAUSS_W[j] * form; 
     69            inner_sum += Gauss150Wt[j] * form; 
    7070        } 
    7171        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    72         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     72        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
    7373    } 
    7474    outer_sum *= theta_m; 

  • sasmodels/models/flexible_cylinder_elliptical.c

    r74768cb r592343f  
    1717    double sum=0.0; 
    1818 
    19     for(int i=0;i<GAUSS_N;i++) { 
    20         const double zi = ( GAUSS_Z[i] + 1.0 )*M_PI_4; 
     19    for(int i=0;i<N_POINTS_76;i++) { 
     20        const double zi = ( Gauss76Z[i] + 1.0 )*M_PI_4; 
    2121        double sn, cn; 
    2222        SINCOS(zi, sn, cn); 
    2323        const double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn); 
    2424        const double yyy = sas_2J1x_x(arg); 
    25         sum += GAUSS_W[i] * yyy * yyy; 
     25        sum += Gauss76Wt[i] * yyy * yyy; 
    2626    } 
    2727    sum *= 0.5; 

  • sasmodels/models/hollow_cylinder.c

    r74768cb rbecded3  
    3838 
    3939    double summ = 0.0;            //initialize intergral 
    40     for (int i=0;i<GAUSS_N;i++) { 
    41         const double cos_theta = 0.5*( GAUSS_Z[i] * (upper-lower) + lower + upper ); 
     40    for (int i=0;i<76;i++) { 
     41        const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
    4242        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
    4343        const double form = _fq(q*sin_theta, q*cos_theta, 
    4444                                radius, thickness, length); 
    45         summ += GAUSS_W[i] * form * form; 
     45        summ += Gauss76Wt[i] * form * form; 
    4646    } 
    4747 

  • sasmodels/models/hollow_rectangular_prism.c

    r74768cb r8de1477  
    3939 
    4040    double outer_sum = 0.0; 
    41     for(int i=0; i<GAUSS_N; i++) { 
     41    for(int i=0; i<76; i++) { 
    4242 
    43         const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
     43        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
    4444        double sin_theta, cos_theta; 
    4545        SINCOS(theta, sin_theta, cos_theta); 
    … …  
    4949 
    5050        double inner_sum = 0.0; 
    51         for(int j=0; j<GAUSS_N; j++) { 
     51        for(int j=0; j<76; j++) { 
    5252 
    53             const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
     53            const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
    5454            double sin_phi, cos_phi; 
    5555            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    6666            const double AP2 = vol_core * termA2 * termB2 * termC2; 
    6767 
    68             inner_sum += GAUSS_W[j] * square(AP1-AP2); 
     68            inner_sum += Gauss76Wt[j] * square(AP1-AP2); 
    6969        } 
    7070        inner_sum *= 0.5 * (v2b-v2a); 
    7171 
    72         outer_sum += GAUSS_W[i] * inner_sum * sin(theta); 
     72        outer_sum += Gauss76Wt[i] * inner_sum * sin(theta); 
    7373    } 
    7474    outer_sum *= 0.5*(v1b-v1a); 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    r74768cb rab2aea8  
    11double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, 
     2double Iq(double q, double sld, double solvent_sld, double length_a,  
    33          double b2a_ratio, double c2a_ratio); 
    44 
    … …  
    2929    const double v2a = 0.0; 
    3030    const double v2b = M_PI_2;  //phi integration limits 
    31  
     31     
    3232    double outer_sum = 0.0; 
    33     for(int i=0; i<GAUSS_N; i++) { 
    34         const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
     33    for(int i=0; i<76; i++) { 
     34        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
    3535 
    3636        double sin_theta, cos_theta; 
    … …  
    4444 
    4545        double inner_sum = 0.0; 
    46         for(int j=0; j<GAUSS_N; j++) { 
    47             const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
     46        for(int j=0; j<76; j++) { 
     47            const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
    4848 
    4949            double sin_phi, cos_phi; 
    … …  
    6262                * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi ); 
    6363 
    64             inner_sum += GAUSS_W[j] * square(AL+AT); 
     64            inner_sum += Gauss76Wt[j] * square(AL+AT); 
    6565        } 
    6666 
    6767        inner_sum *= 0.5 * (v2b-v2a); 
    68         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     68        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
    6969    } 
    7070 

  • sasmodels/models/lib/gauss150.c

    r99b84ec r994d77f  
    1 // Created by Andrew Jackson on 4/23/07 
     1/* 
     2 *  GaussWeights.c 
     3 *  SANSAnalysis 
     4 * 
     5 *  Created by Andrew Jackson on 4/23/07. 
     6 *  Copyright 2007 __MyCompanyName__. All rights reserved. 
     7 * 
     8 */ 
    29 
    3  #ifdef GAUSS_N 
    4  # undef GAUSS_N 
    5  # undef GAUSS_Z 
    6  # undef GAUSS_W 
    7  #endif 
    8  #define GAUSS_N 150 
    9  #define GAUSS_Z Gauss150Z 
    10  #define GAUSS_W Gauss150Wt 
    11  
    12  
    13 // Note: using array size 152 rather than 150 so that it is a multiple of 4. 
    14 // Some OpenCL devices prefer that vectors start and end on nice boundaries. 
    15 constant double Gauss150Z[152]={ 
     10// Gaussians 
     11constant double Gauss150Z[150]={ 
    1612        -0.9998723404457334, 
    1713        -0.9993274305065947, 
    … …  
    163159        0.9983473449340834, 
    164160        0.9993274305065947, 
    165         0.9998723404457334, 
    166         0., // zero padding is ignored 
    167         0.  // zero padding is ignored 
     161        0.9998723404457334 
    168162}; 
    169163 
    170 constant double Gauss150Wt[152]={ 
     164constant double Gauss150Wt[150]={ 
    171165        0.0003276086705538, 
    172166        0.0007624720924706, 
    … …  
    318312        0.0011976474864367, 
    319313        0.0007624720924706, 
    320         0.0003276086705538, 
    321         0., // zero padding is ignored 
    322         0.  // zero padding is ignored 
     314        0.0003276086705538 
    323315}; 

  • sasmodels/models/lib/gauss20.c

    r99b84ec r994d77f  
    1 // Created by Andrew Jackson on 4/23/07 
    2  
    3  #ifdef GAUSS_N 
    4  # undef GAUSS_N 
    5  # undef GAUSS_Z 
    6  # undef GAUSS_W 
    7  #endif 
    8  #define GAUSS_N 20 
    9  #define GAUSS_Z Gauss20Z 
    10  #define GAUSS_W Gauss20Wt 
     1/* 
     2 *  GaussWeights.c 
     3 *  SANSAnalysis 
     4 * 
     5 *  Created by Andrew Jackson on 4/23/07. 
     6 *  Copyright 2007 __MyCompanyName__. All rights reserved. 
     7 * 
     8 */ 
    119 
    1210// Gaussians 

  • sasmodels/models/lib/gauss76.c

    r99b84ec r66d119f  
    1 // Created by Andrew Jackson on 4/23/07 
    2  
    3  #ifdef GAUSS_N 
    4  # undef GAUSS_N 
    5  # undef GAUSS_Z 
    6  # undef GAUSS_W 
    7  #endif 
    8  #define GAUSS_N 76 
    9  #define GAUSS_Z Gauss76Z 
    10  #define GAUSS_W Gauss76Wt 
     1/* 
     2 *  GaussWeights.c 
     3 *  SANSAnalysis 
     4 * 
     5 *  Created by Andrew Jackson on 4/23/07. 
     6 *  Copyright 2007 __MyCompanyName__. All rights reserved. 
     7 * 
     8 */ 
     9#define N_POINTS_76 76 
    1110 
    1211// Gaussians 
    13 constant double Gauss76Wt[76]={ 
     12constant double Gauss76Wt[N_POINTS_76]={ 
    1413        .00126779163408536,             //0 
    1514        .00294910295364247, 
    … …  
    9089}; 
    9190 
    92 constant double Gauss76Z[76]={ 
     91constant double Gauss76Z[N_POINTS_76]={ 
    9392        -.999505948362153,              //0 
    9493        -.997397786355355, 

  • sasmodels/models/parallelepiped.c

    r74768cb r9b7b23f  
    2323    double outer_total = 0; //initialize integral 
    2424 
    25     for( int i=0; i<GAUSS_N; i++) { 
    26         const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
     25    for( int i=0; i<76; i++) { 
     26        const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
    2727        const double mu_proj = mu * sqrt(1.0-sigma*sigma); 
    2828 
    … …  
    3030        // corresponding to angles from 0 to pi/2. 
    3131        double inner_total = 0.0; 
    32         for(int j=0; j<GAUSS_N; j++) { 
    33             const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     32        for(int j=0; j<76; j++) { 
     33            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    3434            double sin_uu, cos_uu; 
    3535            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
    3636            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
    3737            const double si2 = sas_sinx_x(mu_proj * cos_uu); 
    38             inner_total += GAUSS_W[j] * square(si1 * si2); 
     38            inner_total += Gauss76Wt[j] * square(si1 * si2); 
    3939        } 
    4040        inner_total *= 0.5; 
    4141 
    4242        const double si = sas_sinx_x(mu * c_scaled * sigma); 
    43         outer_total += GAUSS_W[i] * inner_total * si * si; 
     43        outer_total += Gauss76Wt[i] * inner_total * si * si; 
    4444    } 
    4545    outer_total *= 0.5; 

  • sasmodels/models/pringle.c

    r74768cb r1e7b0db0  
    2929    double sumC = 0.0;          // initialize integral 
    3030    double r; 
    31     for (int i=0; i < GAUSS_N; i++) { 
    32         r = GAUSS_Z[i]*zm + zb; 
     31    for (int i=0; i < 76; i++) { 
     32        r = Gauss76Z[i]*zm + zb; 
    3333 
    3434        const double qrs = r*q_sin_psi; 
    3535        const double qrrc = r*r*q_cos_psi; 
    3636 
    37         double y = GAUSS_W[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 
     37        double y = Gauss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 
    3838        double S, C; 
    3939        SINCOS(alpha*qrrc, S, C); 
    … …  
    8686 
    8787    double sum = 0.0; 
    88     for (int i = 0; i < GAUSS_N; i++) { 
    89         double psi = GAUSS_Z[i]*zm + zb; 
     88    for (int i = 0; i < 76; i++) { 
     89        double psi = Gauss76Z[i]*zm + zb; 
    9090        double sin_psi, cos_psi; 
    9191        SINCOS(psi, sin_psi, cos_psi); 
    … …  
    9393        double sinc_term = square(sas_sinx_x(q * thickness * cos_psi / 2.0)); 
    9494        double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term; 
    95         sum += GAUSS_W[i] * pringle_kernel; 
     95        sum += Gauss76Wt[i] * pringle_kernel; 
    9696    } 
    9797 

  • sasmodels/models/rectangular_prism.c

    r74768cb r8de1477  
    2828 
    2929    double outer_sum = 0.0; 
    30     for(int i=0; i<GAUSS_N; i++) { 
    31         const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
     30    for(int i=0; i<76; i++) { 
     31        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
    3232        double sin_theta, cos_theta; 
    3333        SINCOS(theta, sin_theta, cos_theta); 
    … …  
    3636 
    3737        double inner_sum = 0.0; 
    38         for(int j=0; j<GAUSS_N; j++) { 
    39             double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
     38        for(int j=0; j<76; j++) { 
     39            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
    4040            double sin_phi, cos_phi; 
    4141            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    4545            const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi); 
    4646            const double AP = termA * termB * termC; 
    47             inner_sum += GAUSS_W[j] * AP * AP; 
     47            inner_sum += Gauss76Wt[j] * AP * AP; 
    4848        } 
    4949        inner_sum = 0.5 * (v2b-v2a) * inner_sum; 
    50         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     50        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
    5151    } 
    5252 

  • sasmodels/models/sc_paracrystal.c

    r74768cb rf728001  
    5454 
    5555    double outer_sum = 0.0; 
    56     for(int i=0; i<GAUSS_N; i++) { 
     56    for(int i=0; i<150; i++) { 
    5757        double inner_sum = 0.0; 
    58         const double theta = GAUSS_Z[i]*theta_m + theta_b; 
     58        const double theta = Gauss150Z[i]*theta_m + theta_b; 
    5959        double sin_theta, cos_theta; 
    6060        SINCOS(theta, sin_theta, cos_theta); 
    6161        const double qc = q*cos_theta; 
    6262        const double qab = q*sin_theta; 
    63         for(int j=0;j<GAUSS_N;j++) { 
    64             const double phi = GAUSS_Z[j]*phi_m + phi_b; 
     63        for(int j=0;j<150;j++) { 
     64            const double phi = Gauss150Z[j]*phi_m + phi_b; 
    6565            double sin_phi, cos_phi; 
    6666            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    6868            const double qb = qab*sin_phi; 
    6969            const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 
    70             inner_sum += GAUSS_W[j] * form; 
     70            inner_sum += Gauss150Wt[j] * form; 
    7171        } 
    7272        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    73         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     73        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
    7474    } 
    7575    outer_sum *= theta_m; 

  • sasmodels/models/stacked_disks.c

    r74768cb rbecded3  
    8181    double halfheight = 0.5*thick_core; 
    8282 
    83     for(int i=0; i<GAUSS_N; i++) { 
    84         double zi = (GAUSS_Z[i] + 1.0)*M_PI_4; 
     83    for(int i=0; i<N_POINTS_76; i++) { 
     84        double zi = (Gauss76Z[i] + 1.0)*M_PI_4; 
    8585        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8686        SINCOS(zi, sin_alpha, cos_alpha); 
    … …  
    9595                           solvent_sld, 
    9696                           d); 
    97         summ += GAUSS_W[i] * yyy * sin_alpha; 
     97        summ += Gauss76Wt[i] * yyy * sin_alpha; 
    9898    } 
    9999 

  • sasmodels/models/triaxial_ellipsoid.c

    r74768cb rbecded3  
    2121    const double zb = M_PI_4; 
    2222    double outer = 0.0; 
    23     for (int i=0;i<GAUSS_N;i++) { 
    24         //const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper + lower)/2; 
    25         const double phi = GAUSS_Z[i]*zm + zb; 
     23    for (int i=0;i<76;i++) { 
     24        //const double u = Gauss76Z[i]*(upper-lower)/2 + (upper + lower)/2; 
     25        const double phi = Gauss76Z[i]*zm + zb; 
    2626        const double pa_sinsq_phi = pa*square(sin(phi)); 
    2727 
    … …  
    2929        const double um = 0.5; 
    3030        const double ub = 0.5; 
    31         for (int j=0;j<GAUSS_N;j++) { 
     31        for (int j=0;j<76;j++) { 
    3232            // translate a point in [-1,1] to a point in [0, 1] 
    33             const double usq = square(GAUSS_Z[j]*um + ub); 
     33            const double usq = square(Gauss76Z[j]*um + ub); 
    3434            const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq); 
    3535            const double fq = sas_3j1x_x(q*r); 
    36             inner += GAUSS_W[j] * fq * fq; 
     36            inner += Gauss76Wt[j] * fq * fq; 
    3737        } 
    38         outer += GAUSS_W[i] * inner;  // correcting for dx later 
     38        outer += Gauss76Wt[i] * inner;  // correcting for dx later 
    3939    } 
    4040    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 

  • sasmodels/special.py

    rdf69efa re65c3ba  
    33................. 
    44 
    5 This following standard C99 math functions are available: 
     5The C code follows the C99 standard, with the usual math functions, 
     6as defined in 
     7`OpenCL <https://www.khronos.org/registry/cl/sdk/1.1/docs/man/xhtml/mathFunctions.html>`_. 
     8This includes the following: 
    69 
    710    M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E: 
    811        $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$ 
    912 
    10     exp, log, pow(x,y), expm1, log1p, sqrt, cbrt: 
    11         Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\ln 1 + x$, 
    12         $\sqrt{x}$, $\sqrt[3]{x}$. The functions expm1(x) and log1p(x) 
    13         are accurate across all $x$, including $x$ very close to zero. 
     13    exp, log, pow(x,y), expm1, sqrt: 
     14        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\sqrt{x}$. 
     15        The function expm1(x) is accurate across all $x$, including $x$ 
     16        very close to zero. 
    1417 
    1518    sin, cos, tan, asin, acos, atan: 
    … …  
    2629        quadrants II and IV when $x$ and $y$ have opposite sign. 
    2730 
    28     fabs(x), fmin(x,y), fmax(x,y), trunc, rint: 
     31    fmin(x,y), fmax(x,y), trunc, rint: 
    2932        Floating point functions.  rint(x) returns the nearest integer. 
    3033 
    … …  
    3235        NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that 
    3336        you cannot use :code:`x == NAN` to test for NaN values since that 
    34         will always return false.  NAN does not equal NAN!  The alternative, 
    35         :code:`x != x` may fail if the compiler optimizes the test away. 
     37        will always return false.  NAN does not equal NAN! 
    3638 
    3739    INFINITY: 
    … …  
    8789        for forcing a constant to stay double precision. 
    8890 
    89 The following special functions and scattering calculations are defined. 
     91The following special functions and scattering calculations are defined in 
     92`sasmodels/models/lib <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib>`_. 
    9093These functions have been tuned to be fast and numerically stable down 
    9194to $q=0$ even in single precision.  In some cases they work around bugs 
    … …  
    181184 
    182185 
    183     gauss76.n, gauss76.z[i], gauss76.w[i]: 
     186    Gauss76Z[i], Gauss76Wt[i]: 
    184187        Points $z_i$ and weights $w_i$ for 76-point Gaussian quadrature, respectively, 
    185188        computing $\int_{-1}^1 f(z)\,dz \approx \sum_{i=1}^{76} w_i\,f(z_i)$. 
    186         When translating the model to C, include 'lib/gauss76.c' in the source 
    187         and use :code:`GAUSS_N`, :code:`GAUSS_Z`, and :code:`GAUSS_W`. 
    188  
    189         Similar arrays are available in :code:`gauss20` for 20-point quadrature 
    190         and :code:`gauss150.c` for 150-point quadrature. By using 
    191         :code:`import gauss76 as gauss` it is easy to change the number of 
    192         points in the integration. 
     189 
     190        Similar arrays are available in :code:`gauss20.c` for 20-point 
     191        quadrature and in :code:`gauss150.c` for 150-point quadrature. 
     192 
    193193""" 
    194194# pylint: disable=unused-import 
    … …  
    200200 
    201201# C99 standard math library functions 
    202 from numpy import exp, log, power as pow, expm1, log1p, sqrt, cbrt 
     202from numpy import exp, log, power as pow, expm1, sqrt 
    203203from numpy import sin, cos, tan, arcsin as asin, arccos as acos, arctan as atan 
    204204from numpy import sinh, cosh, tanh, arcsinh as asinh, arccosh as acosh, arctanh as atanh 
    205205from numpy import arctan2 as atan2 
    206 from numpy import fabs, fmin, fmax, trunc, rint 
    207 from numpy import pi, nan, inf 
     206from numpy import fmin, fmax, trunc, rint 
     207from numpy import NAN, inf as INFINITY 
     208 
    208209from scipy.special import gamma as sas_gamma 
    209210from scipy.special import erf as sas_erf 
    … …  
    217218# C99 standard math constants 
    218219M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E = np.pi, np.pi/2, np.pi/4, np.sqrt(0.5), np.e 
    219 NAN = nan 
    220 INFINITY = inf 
    221220 
    222221# non-standard constants 
    … …  
    227226    """return sin(x), cos(x)""" 
    228227    return sin(x), cos(x) 
    229 sincos = SINCOS 
    230228 
    231229def square(x): 
    … …  
    296294 
    297295# Gaussians 
    298 class Gauss: 
    299     def __init__(self, w, z): 
    300         self.n = len(w) 
    301         self.w = w 
    302         self.z = z 
    303  
    304 gauss20 = Gauss( 
    305     w=np.array([ 
    306         .0176140071391521, 
    307         .0406014298003869, 
    308         .0626720483341091, 
    309         .0832767415767047, 
    310         .10193011981724, 
    311         .118194531961518, 
    312         .131688638449177, 
    313         .142096109318382, 
    314         .149172986472604, 
    315         .152753387130726, 
    316         .152753387130726, 
    317         .149172986472604, 
    318         .142096109318382, 
    319         .131688638449177, 
    320         .118194531961518, 
    321         .10193011981724, 
    322         .0832767415767047, 
    323         .0626720483341091, 
    324         .0406014298003869, 
    325         .0176140071391521 
    326     ]), 
    327     z=np.array([ 
    328         -.993128599185095, 
    329         -.963971927277914, 
    330         -.912234428251326, 
    331         -.839116971822219, 
    332         -.746331906460151, 
    333         -.636053680726515, 
    334         -.510867001950827, 
    335         -.37370608871542, 
    336         -.227785851141645, 
    337         -.076526521133497, 
    338         .0765265211334973, 
    339         .227785851141645, 
    340         .37370608871542, 
    341         .510867001950827, 
    342         .636053680726515, 
    343         .746331906460151, 
    344         .839116971822219, 
    345         .912234428251326, 
    346         .963971927277914, 
    347         .993128599185095 
    348     ]) 
    349 ) 
    350  
    351 gauss76 = Gauss( 
    352     w=np.array([ 
    353         .00126779163408536,             #0 
    354         .00294910295364247, 
    355         .00462793522803742, 
    356         .00629918049732845, 
    357         .00795984747723973, 
    358         .00960710541471375, 
    359         .0112381685696677, 
    360         .0128502838475101, 
    361         .0144407317482767, 
    362         .0160068299122486, 
    363         .0175459372914742,              #10 
    364         .0190554584671906, 
    365         .020532847967908, 
    366         .0219756145344162, 
    367         .0233813253070112, 
    368         .0247476099206597, 
    369         .026072164497986, 
    370         .0273527555318275, 
    371         .028587223650054, 
    372         .029773487255905, 
    373         .0309095460374916,              #20 
    374         .0319934843404216, 
    375         .0330234743977917, 
    376         .0339977794120564, 
    377         .0349147564835508, 
    378         .0357728593807139, 
    379         .0365706411473296, 
    380         .0373067565423816, 
    381         .0379799643084053, 
    382         .0385891292645067, 
    383         .0391332242205184,              #30 
    384         .0396113317090621, 
    385         .0400226455325968, 
    386         .040366472122844, 
    387         .0406422317102947, 
    388         .0408494593018285, 
    389         .040987805464794, 
    390         .0410570369162294, 
    391         .0410570369162294, 
    392         .040987805464794, 
    393         .0408494593018285,              #40 
    394         .0406422317102947, 
    395         .040366472122844, 
    396         .0400226455325968, 
    397         .0396113317090621, 
    398         .0391332242205184, 
    399         .0385891292645067, 
    400         .0379799643084053, 
    401         .0373067565423816, 
    402         .0365706411473296, 
    403         .0357728593807139,              #50 
    404         .0349147564835508, 
    405         .0339977794120564, 
    406         .0330234743977917, 
    407         .0319934843404216, 
    408         .0309095460374916, 
    409         .029773487255905, 
    410         .028587223650054, 
    411         .0273527555318275, 
    412         .026072164497986, 
    413         .0247476099206597,              #60 
    414         .0233813253070112, 
    415         .0219756145344162, 
    416         .020532847967908, 
    417         .0190554584671906, 
    418         .0175459372914742, 
    419         .0160068299122486, 
    420         .0144407317482767, 
    421         .0128502838475101, 
    422         .0112381685696677, 
    423         .00960710541471375,             #70 
    424         .00795984747723973, 
    425         .00629918049732845, 
    426         .00462793522803742, 
    427         .00294910295364247, 
    428         .00126779163408536              #75 (indexed from 0) 
    429     ]), 
    430     z=np.array([ 
    431         -.999505948362153,              #0 
    432         -.997397786355355, 
    433         -.993608772723527, 
    434         -.988144453359837, 
    435         -.981013938975656, 
    436         -.972229228520377, 
    437         -.961805126758768, 
    438         -.949759207710896, 
    439         -.936111781934811, 
    440         -.92088586125215, 
    441         -.904107119545567,              #10 
    442         -.885803849292083, 
    443         -.866006913771982, 
    444         -.844749694983342, 
    445         -.822068037328975, 
    446         -.7980001871612, 
    447         -.77258672828181, 
    448         -.74587051350361, 
    449         -.717896592387704, 
    450         -.688712135277641, 
    451         -.658366353758143,              #20 
    452         -.626910417672267, 
    453         -.594397368836793, 
    454         -.560882031601237, 
    455         -.526420920401243, 
    456         -.491072144462194, 
    457         -.454895309813726, 
    458         -.417951418780327, 
    459         -.380302767117504, 
    460         -.342012838966962, 
    461         -.303146199807908,              #30 
    462         -.263768387584994, 
    463         -.223945802196474, 
    464         -.183745593528914, 
    465         -.143235548227268, 
    466         -.102483975391227, 
    467         -.0615595913906112, 
    468         -.0205314039939986, 
    469         .0205314039939986, 
    470         .0615595913906112, 
    471         .102483975391227,                       #40 
    472         .143235548227268, 
    473         .183745593528914, 
    474         .223945802196474, 
    475         .263768387584994, 
    476         .303146199807908, 
    477         .342012838966962, 
    478         .380302767117504, 
    479         .417951418780327, 
    480         .454895309813726, 
    481         .491072144462194,               #50 
    482         .526420920401243, 
    483         .560882031601237, 
    484         .594397368836793, 
    485         .626910417672267, 
    486         .658366353758143, 
    487         .688712135277641, 
    488         .717896592387704, 
    489         .74587051350361, 
    490         .77258672828181, 
    491         .7980001871612, #60 
    492         .822068037328975, 
    493         .844749694983342, 
    494         .866006913771982, 
    495         .885803849292083, 
    496         .904107119545567, 
    497         .92088586125215, 
    498         .936111781934811, 
    499         .949759207710896, 
    500         .961805126758768, 
    501         .972229228520377,               #70 
    502         .981013938975656, 
    503         .988144453359837, 
    504         .993608772723527, 
    505         .997397786355355, 
    506         .999505948362153                #75 
    507     ]) 
    508 ) 
    509  
    510 gauss150 = Gauss( 
    511     z=np.array([ 
    512         -0.9998723404457334, 
    513         -0.9993274305065947, 
    514         -0.9983473449340834, 
    515         -0.9969322929775997, 
    516         -0.9950828645255290, 
    517         -0.9927998590434373, 
    518         -0.9900842691660192, 
    519         -0.9869372772712794, 
    520         -0.9833602541697529, 
    521         -0.9793547582425894, 
    522         -0.9749225346595943, 
    523         -0.9700655145738374, 
    524         -0.9647858142586956, 
    525         -0.9590857341746905, 
    526         -0.9529677579610971, 
    527         -0.9464345513503147, 
    528         -0.9394889610042837, 
    529         -0.9321340132728527, 
    530         -0.9243729128743136, 
    531         -0.9162090414984952, 
    532         -0.9076459563329236, 
    533         -0.8986873885126239, 
    534         -0.8893372414942055, 
    535         -0.8795995893549102, 
    536         -0.8694786750173527, 
    537         -0.8589789084007133, 
    538         -0.8481048644991847, 
    539         -0.8368612813885015, 
    540         -0.8252530581614230, 
    541         -0.8132852527930605, 
    542         -0.8009630799369827, 
    543         -0.7882919086530552, 
    544         -0.7752772600680049, 
    545         -0.7619248049697269, 
    546         -0.7482403613363824, 
    547         -0.7342298918013638, 
    548         -0.7198995010552305, 
    549         -0.7052554331857488, 
    550         -0.6903040689571928, 
    551         -0.6750519230300931, 
    552         -0.6595056411226444, 
    553         -0.6436719971150083, 
    554         -0.6275578900977726, 
    555         -0.6111703413658551, 
    556         -0.5945164913591590, 
    557         -0.5776035965513142, 
    558         -0.5604390262878617, 
    559         -0.5430302595752546, 
    560         -0.5253848818220803, 
    561         -0.5075105815339176, 
    562         -0.4894151469632753, 
    563         -0.4711064627160663, 
    564         -0.4525925063160997, 
    565         -0.4338813447290861, 
    566         -0.4149811308476706, 
    567         -0.3959000999390257, 
    568         -0.3766465660565522, 
    569         -0.3572289184172501, 
    570         -0.3376556177463400, 
    571         -0.3179351925907259, 
    572         -0.2980762356029071, 
    573         -0.2780873997969574, 
    574         -0.2579773947782034, 
    575         -0.2377549829482451, 
    576         -0.2174289756869712, 
    577         -0.1970082295132342, 
    578         -0.1765016422258567, 
    579         -0.1559181490266516, 
    580         -0.1352667186271445, 
    581         -0.1145563493406956, 
    582         -0.0937960651617229, 
    583         -0.0729949118337358, 
    584         -0.0521619529078925, 
    585         -0.0313062657937972, 
    586         -0.0104369378042598, 
    587         0.0104369378042598, 
    588         0.0313062657937972, 
    589         0.0521619529078925, 
    590         0.0729949118337358, 
    591         0.0937960651617229, 
    592         0.1145563493406956, 
    593         0.1352667186271445, 
    594         0.1559181490266516, 
    595         0.1765016422258567, 
    596         0.1970082295132342, 
    597         0.2174289756869712, 
    598         0.2377549829482451, 
    599         0.2579773947782034, 
    600         0.2780873997969574, 
    601         0.2980762356029071, 
    602         0.3179351925907259, 
    603         0.3376556177463400, 
    604         0.3572289184172501, 
    605         0.3766465660565522, 
    606         0.3959000999390257, 
    607         0.4149811308476706, 
    608         0.4338813447290861, 
    609         0.4525925063160997, 
    610         0.4711064627160663, 
    611         0.4894151469632753, 
    612         0.5075105815339176, 
    613         0.5253848818220803, 
    614         0.5430302595752546, 
    615         0.5604390262878617, 
    616         0.5776035965513142, 
    617         0.5945164913591590, 
    618         0.6111703413658551, 
    619         0.6275578900977726, 
    620         0.6436719971150083, 
    621         0.6595056411226444, 
    622         0.6750519230300931, 
    623         0.6903040689571928, 
    624         0.7052554331857488, 
    625         0.7198995010552305, 
    626         0.7342298918013638, 
    627         0.7482403613363824, 
    628         0.7619248049697269, 
    629         0.7752772600680049, 
    630         0.7882919086530552, 
    631         0.8009630799369827, 
    632         0.8132852527930605, 
    633         0.8252530581614230, 
    634         0.8368612813885015, 
    635         0.8481048644991847, 
    636         0.8589789084007133, 
    637         0.8694786750173527, 
    638         0.8795995893549102, 
    639         0.8893372414942055, 
    640         0.8986873885126239, 
    641         0.9076459563329236, 
    642         0.9162090414984952, 
    643         0.9243729128743136, 
    644         0.9321340132728527, 
    645         0.9394889610042837, 
    646         0.9464345513503147, 
    647         0.9529677579610971, 
    648         0.9590857341746905, 
    649         0.9647858142586956, 
    650         0.9700655145738374, 
    651         0.9749225346595943, 
    652         0.9793547582425894, 
    653         0.9833602541697529, 
    654         0.9869372772712794, 
    655         0.9900842691660192, 
    656         0.9927998590434373, 
    657         0.9950828645255290, 
    658         0.9969322929775997, 
    659         0.9983473449340834, 
    660         0.9993274305065947, 
    661         0.9998723404457334 
    662     ]), 
    663     w=np.array([ 
    664         0.0003276086705538, 
    665         0.0007624720924706, 
    666         0.0011976474864367, 
    667         0.0016323569986067, 
    668         0.0020663664924131, 
    669         0.0024994789888943, 
    670         0.0029315036836558, 
    671         0.0033622516236779, 
    672         0.0037915348363451, 
    673         0.0042191661429919, 
    674         0.0046449591497966, 
    675         0.0050687282939456, 
    676         0.0054902889094487, 
    677         0.0059094573005900, 
    678         0.0063260508184704, 
    679         0.0067398879387430, 
    680         0.0071507883396855, 
    681         0.0075585729801782, 
    682         0.0079630641773633, 
    683         0.0083640856838475, 
    684         0.0087614627643580, 
    685         0.0091550222717888, 
    686         0.0095445927225849, 
    687         0.0099300043714212, 
    688         0.0103110892851360, 
    689         0.0106876814158841, 
    690         0.0110596166734735, 
    691         0.0114267329968529, 
    692         0.0117888704247183, 
    693         0.0121458711652067, 
    694         0.0124975796646449, 
    695         0.0128438426753249, 
    696         0.0131845093222756, 
    697         0.0135194311690004, 
    698         0.0138484622795371, 
    699         0.0141714592928592, 
    700         0.0144882814685445, 
    701         0.0147987907597169, 
    702         0.0151028518701744, 
    703         0.0154003323133401, 
    704         0.0156911024699895, 
    705         0.0159750356447283, 
    706         0.0162520081211971, 
    707         0.0165218992159766, 
    708         0.0167845913311726, 
    709         0.0170399700056559, 
    710         0.0172879239649355, 
    711         0.0175283451696437, 
    712         0.0177611288626114, 
    713         0.0179861736145128, 
    714         0.0182033813680609, 
    715         0.0184126574807331, 
    716         0.0186139107660094, 
    717         0.0188070535331042, 
    718         0.0189920016251754, 
    719         0.0191686744559934, 
    720         0.0193369950450545, 
    721         0.0194968900511231, 
    722         0.0196482898041878, 
    723         0.0197911283358190, 
    724         0.0199253434079123, 
    725         0.0200508765398072, 
    726         0.0201676730337687, 
    727         0.0202756819988200, 
    728         0.0203748563729175, 
    729         0.0204651529434560, 
    730         0.0205465323660984, 
    731         0.0206189591819181, 
    732         0.0206824018328499, 
    733         0.0207368326754401, 
    734         0.0207822279928917, 
    735         0.0208185680053983, 
    736         0.0208458368787627, 
    737         0.0208640227312962, 
    738         0.0208731176389954, 
    739         0.0208731176389954, 
    740         0.0208640227312962, 
    741         0.0208458368787627, 
    742         0.0208185680053983, 
    743         0.0207822279928917, 
    744         0.0207368326754401, 
    745         0.0206824018328499, 
    746         0.0206189591819181, 
    747         0.0205465323660984, 
    748         0.0204651529434560, 
    749         0.0203748563729175, 
    750         0.0202756819988200, 
    751         0.0201676730337687, 
    752         0.0200508765398072, 
    753         0.0199253434079123, 
    754         0.0197911283358190, 
    755         0.0196482898041878, 
    756         0.0194968900511231, 
    757         0.0193369950450545, 
    758         0.0191686744559934, 
    759         0.0189920016251754, 
    760         0.0188070535331042, 
    761         0.0186139107660094, 
    762         0.0184126574807331, 
    763         0.0182033813680609, 
    764         0.0179861736145128, 
    765         0.0177611288626114, 
    766         0.0175283451696437, 
    767         0.0172879239649355, 
    768         0.0170399700056559, 
    769         0.0167845913311726, 
    770         0.0165218992159766, 
    771         0.0162520081211971, 
    772         0.0159750356447283, 
    773         0.0156911024699895, 
    774         0.0154003323133401, 
    775         0.0151028518701744, 
    776         0.0147987907597169, 
    777         0.0144882814685445, 
    778         0.0141714592928592, 
    779         0.0138484622795371, 
    780         0.0135194311690004, 
    781         0.0131845093222756, 
    782         0.0128438426753249, 
    783         0.0124975796646449, 
    784         0.0121458711652067, 
    785         0.0117888704247183, 
    786         0.0114267329968529, 
    787         0.0110596166734735, 
    788         0.0106876814158841, 
    789         0.0103110892851360, 
    790         0.0099300043714212, 
    791         0.0095445927225849, 
    792         0.0091550222717888, 
    793         0.0087614627643580, 
    794         0.0083640856838475, 
    795         0.0079630641773633, 
    796         0.0075585729801782, 
    797         0.0071507883396855, 
    798         0.0067398879387430, 
    799         0.0063260508184704, 
    800         0.0059094573005900, 
    801         0.0054902889094487, 
    802         0.0050687282939456, 
    803         0.0046449591497966, 
    804         0.0042191661429919, 
    805         0.0037915348363451, 
    806         0.0033622516236779, 
    807         0.0029315036836558, 
    808         0.0024994789888943, 
    809         0.0020663664924131, 
    810         0.0016323569986067, 
    811         0.0011976474864367, 
    812         0.0007624720924706, 
    813         0.0003276086705538 
    814     ]) 
    815 ) 
     296 
     297Gauss20Wt = np.array([ 
     298    .0176140071391521, 
     299    .0406014298003869, 
     300    .0626720483341091, 
     301    .0832767415767047, 
     302    .10193011981724, 
     303    .118194531961518, 
     304    .131688638449177, 
     305    .142096109318382, 
     306    .149172986472604, 
     307    .152753387130726, 
     308    .152753387130726, 
     309    .149172986472604, 
     310    .142096109318382, 
     311    .131688638449177, 
     312    .118194531961518, 
     313    .10193011981724, 
     314    .0832767415767047, 
     315    .0626720483341091, 
     316    .0406014298003869, 
     317    .0176140071391521 
     318]) 
     319 
     320Gauss20Z = np.array([ 
     321    -.993128599185095, 
     322    -.963971927277914, 
     323    -.912234428251326, 
     324    -.839116971822219, 
     325    -.746331906460151, 
     326    -.636053680726515, 
     327    -.510867001950827, 
     328    -.37370608871542, 
     329    -.227785851141645, 
     330    -.076526521133497, 
     331    .0765265211334973, 
     332    .227785851141645, 
     333    .37370608871542, 
     334    .510867001950827, 
     335    .636053680726515, 
     336    .746331906460151, 
     337    .839116971822219, 
     338    .912234428251326, 
     339    .963971927277914, 
     340    .993128599185095 
     341]) 
     342 
     343Gauss76Wt = np.array([ 
     344    .00126779163408536,         #0 
     345    .00294910295364247, 
     346    .00462793522803742, 
     347    .00629918049732845, 
     348    .00795984747723973, 
     349    .00960710541471375, 
     350    .0112381685696677, 
     351    .0128502838475101, 
     352    .0144407317482767, 
     353    .0160068299122486, 
     354    .0175459372914742,          #10 
     355    .0190554584671906, 
     356    .020532847967908, 
     357    .0219756145344162, 
     358    .0233813253070112, 
     359    .0247476099206597, 
     360    .026072164497986, 
     361    .0273527555318275, 
     362    .028587223650054, 
     363    .029773487255905, 
     364    .0309095460374916,          #20 
     365    .0319934843404216, 
     366    .0330234743977917, 
     367    .0339977794120564, 
     368    .0349147564835508, 
     369    .0357728593807139, 
     370    .0365706411473296, 
     371    .0373067565423816, 
     372    .0379799643084053, 
     373    .0385891292645067, 
     374    .0391332242205184,          #30 
     375    .0396113317090621, 
     376    .0400226455325968, 
     377    .040366472122844, 
     378    .0406422317102947, 
     379    .0408494593018285, 
     380    .040987805464794, 
     381    .0410570369162294, 
     382    .0410570369162294, 
     383    .040987805464794, 
     384    .0408494593018285,          #40 
     385    .0406422317102947, 
     386    .040366472122844, 
     387    .0400226455325968, 
     388    .0396113317090621, 
     389    .0391332242205184, 
     390    .0385891292645067, 
     391    .0379799643084053, 
     392    .0373067565423816, 
     393    .0365706411473296, 
     394    .0357728593807139,          #50 
     395    .0349147564835508, 
     396    .0339977794120564, 
     397    .0330234743977917, 
     398    .0319934843404216, 
     399    .0309095460374916, 
     400    .029773487255905, 
     401    .028587223650054, 
     402    .0273527555318275, 
     403    .026072164497986, 
     404    .0247476099206597,          #60 
     405    .0233813253070112, 
     406    .0219756145344162, 
     407    .020532847967908, 
     408    .0190554584671906, 
     409    .0175459372914742, 
     410    .0160068299122486, 
     411    .0144407317482767, 
     412    .0128502838475101, 
     413    .0112381685696677, 
     414    .00960710541471375,         #70 
     415    .00795984747723973, 
     416    .00629918049732845, 
     417    .00462793522803742, 
     418    .00294910295364247, 
     419    .00126779163408536          #75 (indexed from 0) 
     420]) 
     421 
     422Gauss76Z = np.array([ 
     423    -.999505948362153,          #0 
     424    -.997397786355355, 
     425    -.993608772723527, 
     426    -.988144453359837, 
     427    -.981013938975656, 
     428    -.972229228520377, 
     429    -.961805126758768, 
     430    -.949759207710896, 
     431    -.936111781934811, 
     432    -.92088586125215, 
     433    -.904107119545567,          #10 
     434    -.885803849292083, 
     435    -.866006913771982, 
     436    -.844749694983342, 
     437    -.822068037328975, 
     438    -.7980001871612, 
     439    -.77258672828181, 
     440    -.74587051350361, 
     441    -.717896592387704, 
     442    -.688712135277641, 
     443    -.658366353758143,          #20 
     444    -.626910417672267, 
     445    -.594397368836793, 
     446    -.560882031601237, 
     447    -.526420920401243, 
     448    -.491072144462194, 
     449    -.454895309813726, 
     450    -.417951418780327, 
     451    -.380302767117504, 
     452    -.342012838966962, 
     453    -.303146199807908,          #30 
     454    -.263768387584994, 
     455    -.223945802196474, 
     456    -.183745593528914, 
     457    -.143235548227268, 
     458    -.102483975391227, 
     459    -.0615595913906112, 
     460    -.0205314039939986, 
     461    .0205314039939986, 
     462    .0615595913906112, 
     463    .102483975391227,                   #40 
     464    .143235548227268, 
     465    .183745593528914, 
     466    .223945802196474, 
     467    .263768387584994, 
     468    .303146199807908, 
     469    .342012838966962, 
     470    .380302767117504, 
     471    .417951418780327, 
     472    .454895309813726, 
     473    .491072144462194,           #50 
     474    .526420920401243, 
     475    .560882031601237, 
     476    .594397368836793, 
     477    .626910417672267, 
     478    .658366353758143, 
     479    .688712135277641, 
     480    .717896592387704, 
     481    .74587051350361, 
     482    .77258672828181, 
     483    .7980001871612,     #60 
     484    .822068037328975, 
     485    .844749694983342, 
     486    .866006913771982, 
     487    .885803849292083, 
     488    .904107119545567, 
     489    .92088586125215, 
     490    .936111781934811, 
     491    .949759207710896, 
     492    .961805126758768, 
     493    .972229228520377,           #70 
     494    .981013938975656, 
     495    .988144453359837, 
     496    .993608772723527, 
     497    .997397786355355, 
     498    .999505948362153            #75 
     499]) 
     500 
     501Gauss150Z = np.array([ 
     502    -0.9998723404457334, 
     503    -0.9993274305065947, 
     504    -0.9983473449340834, 
     505    -0.9969322929775997, 
     506    -0.9950828645255290, 
     507    -0.9927998590434373, 
     508    -0.9900842691660192, 
     509    -0.9869372772712794, 
     510    -0.9833602541697529, 
     511    -0.9793547582425894, 
     512    -0.9749225346595943, 
     513    -0.9700655145738374, 
     514    -0.9647858142586956, 
     515    -0.9590857341746905, 
     516    -0.9529677579610971, 
     517    -0.9464345513503147, 
     518    -0.9394889610042837, 
     519    -0.9321340132728527, 
     520    -0.9243729128743136, 
     521    -0.9162090414984952, 
     522    -0.9076459563329236, 
     523    -0.8986873885126239, 
     524    -0.8893372414942055, 
     525    -0.8795995893549102, 
     526    -0.8694786750173527, 
     527    -0.8589789084007133, 
     528    -0.8481048644991847, 
     529    -0.8368612813885015, 
     530    -0.8252530581614230, 
     531    -0.8132852527930605, 
     532    -0.8009630799369827, 
     533    -0.7882919086530552, 
     534    -0.7752772600680049, 
     535    -0.7619248049697269, 
     536    -0.7482403613363824, 
     537    -0.7342298918013638, 
     538    -0.7198995010552305, 
     539    -0.7052554331857488, 
     540    -0.6903040689571928, 
     541    -0.6750519230300931, 
     542    -0.6595056411226444, 
     543    -0.6436719971150083, 
     544    -0.6275578900977726, 
     545    -0.6111703413658551, 
     546    -0.5945164913591590, 
     547    -0.5776035965513142, 
     548    -0.5604390262878617, 
     549    -0.5430302595752546, 
     550    -0.5253848818220803, 
     551    -0.5075105815339176, 
     552    -0.4894151469632753, 
     553    -0.4711064627160663, 
     554    -0.4525925063160997, 
     555    -0.4338813447290861, 
     556    -0.4149811308476706, 
     557    -0.3959000999390257, 
     558    -0.3766465660565522, 
     559    -0.3572289184172501, 
     560    -0.3376556177463400, 
     561    -0.3179351925907259, 
     562    -0.2980762356029071, 
     563    -0.2780873997969574, 
     564    -0.2579773947782034, 
     565    -0.2377549829482451, 
     566    -0.2174289756869712, 
     567    -0.1970082295132342, 
     568    -0.1765016422258567, 
     569    -0.1559181490266516, 
     570    -0.1352667186271445, 
     571    -0.1145563493406956, 
     572    -0.0937960651617229, 
     573    -0.0729949118337358, 
     574    -0.0521619529078925, 
     575    -0.0313062657937972, 
     576    -0.0104369378042598, 
     577    0.0104369378042598, 
     578    0.0313062657937972, 
     579    0.0521619529078925, 
     580    0.0729949118337358, 
     581    0.0937960651617229, 
     582    0.1145563493406956, 
     583    0.1352667186271445, 
     584    0.1559181490266516, 
     585    0.1765016422258567, 
     586    0.1970082295132342, 
     587    0.2174289756869712, 
     588    0.2377549829482451, 
     589    0.2579773947782034, 
     590    0.2780873997969574, 
     591    0.2980762356029071, 
     592    0.3179351925907259, 
     593    0.3376556177463400, 
     594    0.3572289184172501, 
     595    0.3766465660565522, 
     596    0.3959000999390257, 
     597    0.4149811308476706, 
     598    0.4338813447290861, 
     599    0.4525925063160997, 
     600    0.4711064627160663, 
     601    0.4894151469632753, 
     602    0.5075105815339176, 
     603    0.5253848818220803, 
     604    0.5430302595752546, 
     605    0.5604390262878617, 
     606    0.5776035965513142, 
     607    0.5945164913591590, 
     608    0.6111703413658551, 
     609    0.6275578900977726, 
     610    0.6436719971150083, 
     611    0.6595056411226444, 
     612    0.6750519230300931, 
     613    0.6903040689571928, 
     614    0.7052554331857488, 
     615    0.7198995010552305, 
     616    0.7342298918013638, 
     617    0.7482403613363824, 
     618    0.7619248049697269, 
     619    0.7752772600680049, 
     620    0.7882919086530552, 
     621    0.8009630799369827, 
     622    0.8132852527930605, 
     623    0.8252530581614230, 
     624    0.8368612813885015, 
     625    0.8481048644991847, 
     626    0.8589789084007133, 
     627    0.8694786750173527, 
     628    0.8795995893549102, 
     629    0.8893372414942055, 
     630    0.8986873885126239, 
     631    0.9076459563329236, 
     632    0.9162090414984952, 
     633    0.9243729128743136, 
     634    0.9321340132728527, 
     635    0.9394889610042837, 
     636    0.9464345513503147, 
     637    0.9529677579610971, 
     638    0.9590857341746905, 
     639    0.9647858142586956, 
     640    0.9700655145738374, 
     641    0.9749225346595943, 
     642    0.9793547582425894, 
     643    0.9833602541697529, 
     644    0.9869372772712794, 
     645    0.9900842691660192, 
     646    0.9927998590434373, 
     647    0.9950828645255290, 
     648    0.9969322929775997, 
     649    0.9983473449340834, 
     650    0.9993274305065947, 
     651    0.9998723404457334 
     652]) 
     653 
     654Gauss150Wt = np.array([ 
     655    0.0003276086705538, 
     656    0.0007624720924706, 
     657    0.0011976474864367, 
     658    0.0016323569986067, 
     659    0.0020663664924131, 
     660    0.0024994789888943, 
     661    0.0029315036836558, 
     662    0.0033622516236779, 
     663    0.0037915348363451, 
     664    0.0042191661429919, 
     665    0.0046449591497966, 
     666    0.0050687282939456, 
     667    0.0054902889094487, 
     668    0.0059094573005900, 
     669    0.0063260508184704, 
     670    0.0067398879387430, 
     671    0.0071507883396855, 
     672    0.0075585729801782, 
     673    0.0079630641773633, 
     674    0.0083640856838475, 
     675    0.0087614627643580, 
     676    0.0091550222717888, 
     677    0.0095445927225849, 
     678    0.0099300043714212, 
     679    0.0103110892851360, 
     680    0.0106876814158841, 
     681    0.0110596166734735, 
     682    0.0114267329968529, 
     683    0.0117888704247183, 
     684    0.0121458711652067, 
     685    0.0124975796646449, 
     686    0.0128438426753249, 
     687    0.0131845093222756, 
     688    0.0135194311690004, 
     689    0.0138484622795371, 
     690    0.0141714592928592, 
     691    0.0144882814685445, 
     692    0.0147987907597169, 
     693    0.0151028518701744, 
     694    0.0154003323133401, 
     695    0.0156911024699895, 
     696    0.0159750356447283, 
     697    0.0162520081211971, 
     698    0.0165218992159766, 
     699    0.0167845913311726, 
     700    0.0170399700056559, 
     701    0.0172879239649355, 
     702    0.0175283451696437, 
     703    0.0177611288626114, 
     704    0.0179861736145128, 
     705    0.0182033813680609, 
     706    0.0184126574807331, 
     707    0.0186139107660094, 
     708    0.0188070535331042, 
     709    0.0189920016251754, 
     710    0.0191686744559934, 
     711    0.0193369950450545, 
     712    0.0194968900511231, 
     713    0.0196482898041878, 
     714    0.0197911283358190, 
     715    0.0199253434079123, 
     716    0.0200508765398072, 
     717    0.0201676730337687, 
     718    0.0202756819988200, 
     719    0.0203748563729175, 
     720    0.0204651529434560, 
     721    0.0205465323660984, 
     722    0.0206189591819181, 
     723    0.0206824018328499, 
     724    0.0207368326754401, 
     725    0.0207822279928917, 
     726    0.0208185680053983, 
     727    0.0208458368787627, 
     728    0.0208640227312962, 
     729    0.0208731176389954, 
     730    0.0208731176389954, 
     731    0.0208640227312962, 
     732    0.0208458368787627, 
     733    0.0208185680053983, 
     734    0.0207822279928917, 
     735    0.0207368326754401, 
     736    0.0206824018328499, 
     737    0.0206189591819181, 
     738    0.0205465323660984, 
     739    0.0204651529434560, 
     740    0.0203748563729175, 
     741    0.0202756819988200, 
     742    0.0201676730337687, 
     743    0.0200508765398072, 
     744    0.0199253434079123, 
     745    0.0197911283358190, 
     746    0.0196482898041878, 
     747    0.0194968900511231, 
     748    0.0193369950450545, 
     749    0.0191686744559934, 
     750    0.0189920016251754, 
     751    0.0188070535331042, 
     752    0.0186139107660094, 
     753    0.0184126574807331, 
     754    0.0182033813680609, 
     755    0.0179861736145128, 
     756    0.0177611288626114, 
     757    0.0175283451696437, 
     758    0.0172879239649355, 
     759    0.0170399700056559, 
     760    0.0167845913311726, 
     761    0.0165218992159766, 
     762    0.0162520081211971, 
     763    0.0159750356447283, 
     764    0.0156911024699895, 
     765    0.0154003323133401, 
     766    0.0151028518701744, 
     767    0.0147987907597169, 
     768    0.0144882814685445, 
     769    0.0141714592928592, 
     770    0.0138484622795371, 
     771    0.0135194311690004, 
     772    0.0131845093222756, 
     773    0.0128438426753249, 
     774    0.0124975796646449, 
     775    0.0121458711652067, 
     776    0.0117888704247183, 
     777    0.0114267329968529, 
     778    0.0110596166734735, 
     779    0.0106876814158841, 
     780    0.0103110892851360, 
     781    0.0099300043714212, 
     782    0.0095445927225849, 
     783    0.0091550222717888, 
     784    0.0087614627643580, 
     785    0.0083640856838475, 
     786    0.0079630641773633, 
     787    0.0075585729801782, 
     788    0.0071507883396855, 
     789    0.0067398879387430, 
     790    0.0063260508184704, 
     791    0.0059094573005900, 
     792    0.0054902889094487, 
     793    0.0050687282939456, 
     794    0.0046449591497966, 
     795    0.0042191661429919, 
     796    0.0037915348363451, 
     797    0.0033622516236779, 
     798    0.0029315036836558, 
     799    0.0024994789888943, 
     800    0.0020663664924131, 
     801    0.0016323569986067, 
     802    0.0011976474864367, 
     803    0.0007624720924706, 
     804    0.0003276086705538 
     805]) 
Note: See TracChangeset for help on using the changeset viewer.