Changes in / [d4db147:d8ac2ad] in sasmodels


Ignore:
Location:
sasmodels
Files:
1 deleted
29 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/compare.py

    rff31782 r376b0ee  
    4444from .direct_model import DirectModel, get_mesh 
    4545from .convert import revert_name, revert_pars, constrain_new_to_old 
    46 from .generate import FLOAT_RE, set_integration_size 
     46from .generate import FLOAT_RE 
    4747from .weights import plot_weights 
    4848 
     
    801801    return data, index 
    802802 
    803 def make_engine(model_info, data, dtype, cutoff, ngauss=0): 
     803def make_engine(model_info, data, dtype, cutoff): 
    804804    # type: (ModelInfo, Data, str, float) -> Calculator 
    805805    """ 
     
    809809    than OpenCL. 
    810810    """ 
    811     if ngauss: 
    812         set_integration_size(model_info, ngauss) 
    813  
    814811    if dtype == 'sasview': 
    815812        return eval_sasview(model_info, data) 
     
    10481045    'poly', 'mono', 'cutoff=', 
    10491046    'magnetic', 'nonmagnetic', 
    1050     'accuracy=', 'ngauss=', 
     1047    'accuracy=', 
    10511048    'neval=',  # for timing... 
    10521049 
     
    11821179        'show_weights' : False, 
    11831180        'sphere'    : 0, 
    1184         'ngauss'    : '0', 
    11851181    } 
    11861182    for arg in flags: 
     
    12091205        elif arg.startswith('-engine='):   opts['engine'] = arg[8:] 
    12101206        elif arg.startswith('-neval='):    opts['count'] = arg[7:] 
    1211         elif arg.startswith('-ngauss='):   opts['ngauss'] = arg[8:] 
    12121207        elif arg.startswith('-random='): 
    12131208            opts['seed'] = int(arg[8:]) 
     
    12651260 
    12661261    comparison = any(PAR_SPLIT in v for v in values) 
    1267  
    12681262    if PAR_SPLIT in name: 
    12691263        names = name.split(PAR_SPLIT, 2) 
     
    12781272        return None 
    12791273 
    1280     if PAR_SPLIT in opts['ngauss']: 
    1281         opts['ngauss'] = [int(k) for k in opts['ngauss'].split(PAR_SPLIT, 2)] 
    1282         comparison = True 
    1283     else: 
    1284         opts['ngauss'] = [int(opts['ngauss'])]*2 
    1285  
    12861274    if PAR_SPLIT in opts['engine']: 
    12871275        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2) 
     
    13021290        opts['cutoff'] = [float(opts['cutoff'])]*2 
    13031291 
    1304     base = make_engine(model_info[0], data, opts['engine'][0], 
    1305                        opts['cutoff'][0], opts['ngauss'][0]) 
     1292    base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0]) 
    13061293    if comparison: 
    1307         comp = make_engine(model_info[1], data, opts['engine'][1], 
    1308                            opts['cutoff'][1], opts['ngauss'][1]) 
     1294        comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1]) 
    13091295    else: 
    13101296        comp = None 
  • sasmodels/generate.py

    rff31782 rff10479  
    268268""" 
    269269 
    270  
    271 def set_integration_size(info, n): 
    272     # type: (ModelInfo, int) -> None 
    273     """ 
    274     Update the model definition, replacing the gaussian integration with 
    275     a gaussian integration of a different size. 
    276  
    277     Note: this really ought to be a method in modelinfo, but that leads to 
    278     import loops. 
    279     """ 
    280     if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)): 
    281         import os.path 
    282         from .gengauss import gengauss 
    283         path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n) 
    284         if not os.path.exists(path): 
    285             gengauss(n, path) 
    286         info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss') 
    287                         else lib for lib in info.source] 
    288270 
    289271def 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

    r74768cb rfc0b7aa  
    8080    // outer integral (with gauss points), integration limits = 0, 1 
    8181    double outer_sum = 0; //initialize integral 
    82     for( int i=0; i<GAUSS_N; i++) { 
    83         double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
     82    for( int i=0; i<76; i++) { 
     83        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
    8484        double mu_proj = mu * sqrt(1.0-sigma*sigma); 
    8585 
     
    8888        const double siCt = sas_sinx_x(mu * sigma * tC_over_b); 
    8989        double inner_sum = 0.0; 
    90         for(int j=0; j<GAUSS_N; j++) { 
    91             const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     90        for(int j=0; j<76; j++) { 
     91            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    9292            double sin_uu, cos_uu; 
    9393            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
     
    109109#endif 
    110110 
    111             inner_sum += GAUSS_W[j] * f * f; 
     111            inner_sum += Gauss76Wt[j] * f * f; 
    112112        } 
    113113        inner_sum *= 0.5; 
    114114        // now sum up the outer integral 
    115         outer_sum += GAUSS_W[i] * inner_sum; 
     115        outer_sum += Gauss76Wt[i] * inner_sum; 
    116116    } 
    117117    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

    r74768cb reda8b30  
    4343    P(q) = \text{scale}  <F^2> / V 
    4444 
    45 For 2d data the orientation of the particle is required, described using a different set 
    46 of angles as in the diagrams below, for further details of the calculation and angular 
     45For 2d data the orientation of the particle is required, described using a different set  
     46of angles as in the diagrams below, for further details of the calculation and angular  
    4747dispersions  see :ref:`orientation` . 
    4848 
     
    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

    r74768cb r994d77f  
    77 * 
    88 */ 
    9  #ifdef GAUSS_N 
    10  # undef GAUSS_N 
    11  # undef GAUSS_Z 
    12  # undef GAUSS_W 
    13  #endif 
    14  #define GAUSS_N 150 
    15  #define GAUSS_Z Gauss150Z 
    16  #define GAUSS_W Gauss150Wt 
    17  
    18 // Note: using array size 152 so that it is a multiple of 4 
    199 
    2010// Gaussians 
    21 constant double Gauss150Z[152]={ 
     11constant double Gauss150Z[150]={ 
    2212        -0.9998723404457334, 
    2313        -0.9993274305065947, 
     
    169159        0.9983473449340834, 
    170160        0.9993274305065947, 
    171         0.9998723404457334, 
    172         0., 
    173         0. 
     161        0.9998723404457334 
    174162}; 
    175163 
    176 constant double Gauss150Wt[152]={ 
     164constant double Gauss150Wt[150]={ 
    177165        0.0003276086705538, 
    178166        0.0007624720924706, 
     
    324312        0.0011976474864367, 
    325313        0.0007624720924706, 
    326         0.0003276086705538, 
    327         0., 
    328         0. 
     314        0.0003276086705538 
    329315}; 
  • sasmodels/models/lib/gauss20.c

    r74768cb r994d77f  
    77 * 
    88 */ 
    9  #ifdef GAUSS_N 
    10  # undef GAUSS_N 
    11  # undef GAUSS_Z 
    12  # undef GAUSS_W 
    13  #endif 
    14  #define GAUSS_N 20 
    15  #define GAUSS_Z Gauss20Z 
    16  #define GAUSS_W Gauss20Wt 
    179 
    1810// Gaussians 
  • sasmodels/models/lib/gauss76.c

    r74768cb r66d119f  
    77 * 
    88 */ 
    9  #ifdef GAUSS_N 
    10  # undef GAUSS_N 
    11  # undef GAUSS_Z 
    12  # undef GAUSS_W 
    13  #endif 
    14  #define GAUSS_N 76 
    15  #define GAUSS_Z Gauss76Z 
    16  #define GAUSS_W Gauss76Wt 
     9#define N_POINTS_76 76 
    1710 
    1811// Gaussians 
    19 constant double Gauss76Wt[76]={ 
     12constant double Gauss76Wt[N_POINTS_76]={ 
    2013        .00126779163408536,             //0 
    2114        .00294910295364247, 
     
    9689}; 
    9790 
    98 constant double Gauss76Z[76]={ 
     91constant double Gauss76Z[N_POINTS_76]={ 
    9992        -.999505948362153,              //0 
    10093        -.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 
Note: See TracChangeset for help on using the changeset viewer.