Changes in / [d8ac2ad:d4db147] in sasmodels


Ignore:
Location:
sasmodels
Files:
1 added
29 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/compare.py

    r376b0ee rff31782  
    4444from .direct_model import DirectModel, get_mesh 
    4545from .convert import revert_name, revert_pars, constrain_new_to_old 
    46 from .generate import FLOAT_RE 
     46from .generate import FLOAT_RE, set_integration_size 
    4747from .weights import plot_weights 
    4848 
    … …  
    801801    return data, index 
    802802 
    803 def make_engine(model_info, data, dtype, cutoff): 
     803def make_engine(model_info, data, dtype, cutoff, ngauss=0): 
    804804    # type: (ModelInfo, Data, str, float) -> Calculator 
    805805    """ 
    … …  
    809809    than OpenCL. 
    810810    """ 
     811    if ngauss: 
     812        set_integration_size(model_info, ngauss) 
     813 
    811814    if dtype == 'sasview': 
    812815        return eval_sasview(model_info, data) 
    … …  
    10451048    'poly', 'mono', 'cutoff=', 
    10461049    'magnetic', 'nonmagnetic', 
    1047     'accuracy=', 
     1050    'accuracy=', 'ngauss=', 
    10481051    'neval=',  # for timing... 
    10491052 
    … …  
    11791182        'show_weights' : False, 
    11801183        'sphere'    : 0, 
     1184        'ngauss'    : '0', 
    11811185    } 
    11821186    for arg in flags: 
    … …  
    12051209        elif arg.startswith('-engine='):   opts['engine'] = arg[8:] 
    12061210        elif arg.startswith('-neval='):    opts['count'] = arg[7:] 
     1211        elif arg.startswith('-ngauss='):   opts['ngauss'] = arg[8:] 
    12071212        elif arg.startswith('-random='): 
    12081213            opts['seed'] = int(arg[8:]) 
    … …  
    12601265 
    12611266    comparison = any(PAR_SPLIT in v for v in values) 
     1267 
    12621268    if PAR_SPLIT in name: 
    12631269        names = name.split(PAR_SPLIT, 2) 
    … …  
    12721278        return None 
    12731279 
     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 
    12741286    if PAR_SPLIT in opts['engine']: 
    12751287        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2) 
    … …  
    12901302        opts['cutoff'] = [float(opts['cutoff'])]*2 
    12911303 
    1292     base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0]) 
     1304    base = make_engine(model_info[0], data, opts['engine'][0], 
     1305                       opts['cutoff'][0], opts['ngauss'][0]) 
    12931306    if comparison: 
    1294         comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1]) 
     1307        comp = make_engine(model_info[1], data, opts['engine'][1], 
     1308                           opts['cutoff'][1], opts['ngauss'][1]) 
    12951309    else: 
    12961310        comp = None 

  • sasmodels/generate.py

    rff10479 rff31782  
    268268""" 
    269269 
     270 
     271def 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] 
    270288 
    271289def format_units(units): 

  • sasmodels/models/barbell.c

    rbecded3 r74768cb  
    2323    const double qab_r = radius_bell*qab; // Q*R*sin(theta) 
    2424    double total = 0.0; 
    25     for (int i = 0; i < 76; i++){ 
    26         const double t = Gauss76Z[i]*zm + zb; 
     25    for (int i = 0; i < GAUSS_N; i++){ 
     26        const double t = GAUSS_Z[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 += Gauss76Wt[i] * Fq; 
     30        total += GAUSS_W[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 < 76; i++){ 
    76         const double alpha= Gauss76Z[i]*zm + zb; 
     75    for (int i = 0; i < GAUSS_N; i++){ 
     76        const double alpha= GAUSS_Z[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 += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
     80        total += GAUSS_W[i] * Aq * Aq * sin_alpha; 
    8181    } 
    8282    // translate dx in [-1,1] to dx in [lower,upper] 

  • sasmodels/models/bcc_paracrystal.c

    rea60e08 r74768cb  
    8181 
    8282    double outer_sum = 0.0; 
    83     for(int i=0; i<150; i++) { 
     83    for(int i=0; i<GAUSS_N; i++) { 
    8484        double inner_sum = 0.0; 
    85         const double theta = Gauss150Z[i]*theta_m + theta_b; 
     85        const double theta = GAUSS_Z[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<150;j++) { 
    91             const double phi = Gauss150Z[j]*phi_m + phi_b; 
     90        for(int j=0;j<GAUSS_N;j++) { 
     91            const double phi = GAUSS_Z[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 += Gauss150Wt[j] * form; 
     97            inner_sum += GAUSS_W[j] * form; 
    9898        } 
    9999        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    100         outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     100        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    101101    } 
    102102    outer_sum *= theta_m; 

  • sasmodels/models/capped_cylinder.c

    rbecded3 r74768cb  
    3030    const double qab_r = radius_cap*qab; // Q*R*sin(theta) 
    3131    double total = 0.0; 
    32     for (int i=0; i<76 ;i++) { 
    33         const double t = Gauss76Z[i]*zm + zb; 
     32    for (int i=0; i<GAUSS_N; i++) { 
     33        const double t = GAUSS_Z[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 += Gauss76Wt[i] * Fq; 
     37        total += GAUSS_W[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<76 ;i++) { 
    98         const double theta = Gauss76Z[i]*zm + zb; 
     97    for (int i=0; i<GAUSS_N ;i++) { 
     98        const double theta = GAUSS_Z[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 += Gauss76Wt[i] * Aq * Aq * sin_theta; 
     105        total += GAUSS_W[i] * Aq * Aq * sin_theta; 
    106106    } 
    107107    // translate dx in [-1,1] to dx in [lower,upper] 

  • sasmodels/models/core_shell_bicelle.c

    rbecded3 r74768cb  
    5252 
    5353    double total = 0.0; 
    54     for(int i=0;i<N_POINTS_76;i++) { 
    55         double theta = (Gauss76Z[i] + 1.0)*uplim; 
     54    for(int i=0;i<GAUSS_N;i++) { 
     55        double theta = (GAUSS_Z[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 += Gauss76Wt[i]*fq*fq*sin_theta; 
     60        total += GAUSS_W[i]*fq*fq*sin_theta; 
    6161    } 
    6262 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r82592da r74768cb  
    3737    //initialize integral 
    3838    double outer_sum = 0.0; 
    39     for(int i=0;i<76;i++) { 
     39    for(int i=0;i<GAUSS_N;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 = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    44         const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.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; 
    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<76;j++) { 
     51        for(int j=0;j<GAUSS_N;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 = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    57             const double phi = ( Gauss76Z[j] +1.0)*M_PI_2; 
     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; 
    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 += Gauss76Wt[j] * fq * fq; 
     63            inner_sum += GAUSS_W[j] * fq * fq; 
    6464        } 
    6565        //now calculate outer integral 
    66         outer_sum += Gauss76Wt[i] * inner_sum; 
     66        outer_sum += GAUSS_W[i] * inner_sum; 
    6767    } 
    6868 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r82592da r74768cb  
    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<76;i++) { 
     49    for(int i=0;i<GAUSS_N;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 = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    53         const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0; 
     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; 
    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<76;j++) { 
     60        for(int j=0;j<GAUSS_N;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 = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    63             const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
     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; 
    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 += Gauss76Wt[j] *square(dr1*si1*be1 + 
     69            inner_sum += GAUSS_W[j] *square(dr1*si1*be1 + 
    7070                                              dr2*si1*be2 + 
    7171                                              dr3*si2*be1); 
    7272        } 
    7373        //now calculate outer integral 
    74         outer_sum += Gauss76Wt[i] * inner_sum; 
     74        outer_sum += GAUSS_W[i] * inner_sum; 
    7575    } 
    7676 

  • sasmodels/models/core_shell_cylinder.c

    rbecded3 r74768cb  
    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<76 ;i++) { 
     32    for (int i=0; i<GAUSS_N ;i++) { 
    3333        // translate a point in [-1,1] to a point in [0, pi/2] 
    34         //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     34        //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0; 
    3535        double sin_theta, cos_theta; 
    36         const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 
     36        const double theta = GAUSS_Z[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 += Gauss76Wt[i] * fq * fq * sin_theta; 
     42        total += GAUSS_W[i] * fq * fq * sin_theta; 
    4343    } 
    4444    // translate dx in [-1,1] to dx in [lower,upper] 

  • sasmodels/models/core_shell_ellipsoid.c

    rbecded3 r74768cb  
    5959    const double b = 0.5; 
    6060    double total = 0.0;     //initialize intergral 
    61     for(int i=0;i<76;i++) { 
    62         const double cos_theta = Gauss76Z[i]*m + b; 
     61    for(int i=0;i<GAUSS_N;i++) { 
     62        const double cos_theta = GAUSS_Z[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 += Gauss76Wt[i] * fq * fq; 
     68        total += GAUSS_W[i] * fq * fq; 
    6969    } 
    7070    total *= m; 

  • sasmodels/models/core_shell_parallelepiped.c

    rfc0b7aa r74768cb  
    8080    // outer integral (with gauss points), integration limits = 0, 1 
    8181    double outer_sum = 0; //initialize integral 
    82     for( int i=0; i<76; i++) { 
    83         double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     82    for( int i=0; i<GAUSS_N; i++) { 
     83        double sigma = 0.5 * ( GAUSS_Z[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<76; j++) { 
    91             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     90        for(int j=0; j<GAUSS_N; j++) { 
     91            const double uu = 0.5 * ( GAUSS_Z[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 += Gauss76Wt[j] * f * f; 
     111            inner_sum += GAUSS_W[j] * f * f; 
    112112        } 
    113113        inner_sum *= 0.5; 
    114114        // now sum up the outer integral 
    115         outer_sum += Gauss76Wt[i] * inner_sum; 
     115        outer_sum += GAUSS_W[i] * inner_sum; 
    116116    } 
    117117    outer_sum *= 0.5; 

  • sasmodels/models/cylinder.c

    rbecded3 r74768cb  
    2121 
    2222    double total = 0.0; 
    23     for (int i=0; i<76 ;i++) { 
    24         const double theta = Gauss76Z[i]*zm + zb; 
     23    for (int i=0; i<GAUSS_N ;i++) { 
     24        const double theta = GAUSS_Z[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 += Gauss76Wt[i] * form * form * sin_theta; 
     29        total += GAUSS_W[i] * form * form * sin_theta; 
    3030    } 
    3131    // translate dx in [-1,1] to dx in [lower,upper] 

  • sasmodels/models/ellipsoid.c

    rbecded3 r74768cb  
    2222 
    2323    // translate a point in [-1,1] to a point in [0, 1] 
    24     // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2; 
     24    // const double u = GAUSS_Z[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<76;i++) { 
    29         const double u = Gauss76Z[i]*zm + zb; 
     28    for (int i=0;i<GAUSS_N;i++) { 
     29        const double u = GAUSS_Z[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 += Gauss76Wt[i] * f * f; 
     32        total += GAUSS_W[i] * f * f; 
    3333    } 
    3434    // translate dx in [-1,1] to dx in [lower,upper] 

  • sasmodels/models/elliptical_cylinder.c

    r82592da r74768cb  
    2222    //initialize integral 
    2323    double outer_sum = 0.0; 
    24     for(int i=0;i<76;i++) { 
     24    for(int i=0;i<GAUSS_N;i++) { 
    2525        //setup inner integral over the ellipsoidal cross-section 
    26         const double cos_val = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     26        const double cos_val = ( GAUSS_Z[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<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; 
     30        for(int j=0;j<GAUSS_N;j++) { 
     31            const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    3332            const double r = sin_val*sqrt(rA - rB*cos(theta)); 
    3433            const double be = sas_2J1x_x(q*r); 
    35             inner_sum += Gauss76Wt[j] * be * be; 
     34            inner_sum += GAUSS_W[j] * be * be; 
    3635        } 
    3736        //now calculate the value of the inner integral 
    … …  
    4039        //now calculate outer integral 
    4140        const double si = sas_sinx_x(q*0.5*length*cos_val); 
    42         outer_sum += Gauss76Wt[i] * inner_sum * si * si; 
     41        outer_sum += GAUSS_W[i] * inner_sum * si * si; 
    4342    } 
    4443    outer_sum *= 0.5*(vb-va); 

  • sasmodels/models/elliptical_cylinder.py

    reda8b30 r74768cb  
    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", "lib/gauss20.c", 
    124           "elliptical_cylinder.c"] 
     123source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
    125124 
    126125demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 

  • sasmodels/models/fcc_paracrystal.c

    rf728001 r74768cb  
    5353 
    5454    double outer_sum = 0.0; 
    55     for(int i=0; i<150; i++) { 
     55    for(int i=0; i<GAUSS_N; i++) { 
    5656        double inner_sum = 0.0; 
    57         const double theta = Gauss150Z[i]*theta_m + theta_b; 
     57        const double theta = GAUSS_Z[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<150;j++) { 
    63             const double phi = Gauss150Z[j]*phi_m + phi_b; 
     62        for(int j=0;j<GAUSS_N;j++) { 
     63            const double phi = GAUSS_Z[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 += Gauss150Wt[j] * form; 
     69            inner_sum += GAUSS_W[j] * form; 
    7070        } 
    7171        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    72         outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     72        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    7373    } 
    7474    outer_sum *= theta_m; 

  • sasmodels/models/flexible_cylinder_elliptical.c

    r592343f r74768cb  
    1717    double sum=0.0; 
    1818 
    19     for(int i=0;i<N_POINTS_76;i++) { 
    20         const double zi = ( Gauss76Z[i] + 1.0 )*M_PI_4; 
     19    for(int i=0;i<GAUSS_N;i++) { 
     20        const double zi = ( GAUSS_Z[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 += Gauss76Wt[i] * yyy * yyy; 
     25        sum += GAUSS_W[i] * yyy * yyy; 
    2626    } 
    2727    sum *= 0.5; 

  • sasmodels/models/hollow_cylinder.c

    rbecded3 r74768cb  
    3838 
    3939    double summ = 0.0;            //initialize intergral 
    40     for (int i=0;i<76;i++) { 
    41         const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
     40    for (int i=0;i<GAUSS_N;i++) { 
     41        const double cos_theta = 0.5*( GAUSS_Z[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 += Gauss76Wt[i] * form * form; 
     45        summ += GAUSS_W[i] * form * form; 
    4646    } 
    4747 

  • sasmodels/models/hollow_rectangular_prism.c

    r8de1477 r74768cb  
    3939 
    4040    double outer_sum = 0.0; 
    41     for(int i=0; i<76; i++) { 
     41    for(int i=0; i<GAUSS_N; i++) { 
    4242 
    43         const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     43        const double theta = 0.5 * ( GAUSS_Z[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<76; j++) { 
     51        for(int j=0; j<GAUSS_N; j++) { 
    5252 
    53             const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     53            const double phi = 0.5 * ( GAUSS_Z[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 += Gauss76Wt[j] * square(AP1-AP2); 
     68            inner_sum += GAUSS_W[j] * square(AP1-AP2); 
    6969        } 
    7070        inner_sum *= 0.5 * (v2b-v2a); 
    7171 
    72         outer_sum += Gauss76Wt[i] * inner_sum * sin(theta); 
     72        outer_sum += GAUSS_W[i] * inner_sum * sin(theta); 
    7373    } 
    7474    outer_sum *= 0.5*(v1b-v1a); 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    rab2aea8 r74768cb  
    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<76; i++) { 
    34         const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     33    for(int i=0; i<GAUSS_N; i++) { 
     34        const double theta = 0.5 * ( GAUSS_Z[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<76; j++) { 
    47             const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     46        for(int j=0; j<GAUSS_N; j++) { 
     47            const double phi = 0.5 * ( GAUSS_Z[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 += Gauss76Wt[j] * square(AL+AT); 
     64            inner_sum += GAUSS_W[j] * square(AL+AT); 
    6565        } 
    6666 
    6767        inner_sum *= 0.5 * (v2b-v2a); 
    68         outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
     68        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    6969    } 
    7070 

  • sasmodels/models/lib/gauss150.c

    r994d77f r74768cb  
    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 
    919 
    1020// Gaussians 
    11 constant double Gauss150Z[150]={ 
     21constant double Gauss150Z[152]={ 
    1222        -0.9998723404457334, 
    1323        -0.9993274305065947, 
    … …  
    159169        0.9983473449340834, 
    160170        0.9993274305065947, 
    161         0.9998723404457334 
     171        0.9998723404457334, 
     172        0., 
     173        0. 
    162174}; 
    163175 
    164 constant double Gauss150Wt[150]={ 
     176constant double Gauss150Wt[152]={ 
    165177        0.0003276086705538, 
    166178        0.0007624720924706, 
    … …  
    312324        0.0011976474864367, 
    313325        0.0007624720924706, 
    314         0.0003276086705538 
     326        0.0003276086705538, 
     327        0., 
     328        0. 
    315329}; 

  • sasmodels/models/lib/gauss20.c

    r994d77f r74768cb  
    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 
    917 
    1018// Gaussians 

  • sasmodels/models/lib/gauss76.c

    r66d119f r74768cb  
    77 * 
    88 */ 
    9 #define N_POINTS_76 76 
     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 
    1017 
    1118// Gaussians 
    12 constant double Gauss76Wt[N_POINTS_76]={ 
     19constant double Gauss76Wt[76]={ 
    1320        .00126779163408536,             //0 
    1421        .00294910295364247, 
    … …  
    8996}; 
    9097 
    91 constant double Gauss76Z[N_POINTS_76]={ 
     98constant double Gauss76Z[76]={ 
    9299        -.999505948362153,              //0 
    93100        -.997397786355355, 

  • sasmodels/models/parallelepiped.c

    r9b7b23f r74768cb  
    2323    double outer_total = 0; //initialize integral 
    2424 
    25     for( int i=0; i<76; i++) { 
    26         const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     25    for( int i=0; i<GAUSS_N; i++) { 
     26        const double sigma = 0.5 * ( GAUSS_Z[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<76; j++) { 
    33             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     32        for(int j=0; j<GAUSS_N; j++) { 
     33            const double uu = 0.5 * ( GAUSS_Z[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 += Gauss76Wt[j] * square(si1 * si2); 
     38            inner_total += GAUSS_W[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 += Gauss76Wt[i] * inner_total * si * si; 
     43        outer_total += GAUSS_W[i] * inner_total * si * si; 
    4444    } 
    4545    outer_total *= 0.5; 

  • sasmodels/models/pringle.c

    r1e7b0db0 r74768cb  
    2929    double sumC = 0.0;          // initialize integral 
    3030    double r; 
    31     for (int i=0; i < 76; i++) { 
    32         r = Gauss76Z[i]*zm + zb; 
     31    for (int i=0; i < GAUSS_N; i++) { 
     32        r = GAUSS_Z[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 = Gauss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 
     37        double y = GAUSS_W[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 < 76; i++) { 
    89         double psi = Gauss76Z[i]*zm + zb; 
     88    for (int i = 0; i < GAUSS_N; i++) { 
     89        double psi = GAUSS_Z[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 += Gauss76Wt[i] * pringle_kernel; 
     95        sum += GAUSS_W[i] * pringle_kernel; 
    9696    } 
    9797 

  • sasmodels/models/rectangular_prism.c

    r8de1477 r74768cb  
    2828 
    2929    double outer_sum = 0.0; 
    30     for(int i=0; i<76; i++) { 
    31         const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     30    for(int i=0; i<GAUSS_N; i++) { 
     31        const double theta = 0.5 * ( GAUSS_Z[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<76; j++) { 
    39             double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     38        for(int j=0; j<GAUSS_N; j++) { 
     39            double phi = 0.5 * ( GAUSS_Z[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 += Gauss76Wt[j] * AP * AP; 
     47            inner_sum += GAUSS_W[j] * AP * AP; 
    4848        } 
    4949        inner_sum = 0.5 * (v2b-v2a) * inner_sum; 
    50         outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
     50        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    5151    } 
    5252 

  • sasmodels/models/sc_paracrystal.c

    rf728001 r74768cb  
    5454 
    5555    double outer_sum = 0.0; 
    56     for(int i=0; i<150; i++) { 
     56    for(int i=0; i<GAUSS_N; i++) { 
    5757        double inner_sum = 0.0; 
    58         const double theta = Gauss150Z[i]*theta_m + theta_b; 
     58        const double theta = GAUSS_Z[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<150;j++) { 
    64             const double phi = Gauss150Z[j]*phi_m + phi_b; 
     63        for(int j=0;j<GAUSS_N;j++) { 
     64            const double phi = GAUSS_Z[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 += Gauss150Wt[j] * form; 
     70            inner_sum += GAUSS_W[j] * form; 
    7171        } 
    7272        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    73         outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     73        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    7474    } 
    7575    outer_sum *= theta_m; 

  • sasmodels/models/stacked_disks.c

    rbecded3 r74768cb  
    8181    double halfheight = 0.5*thick_core; 
    8282 
    83     for(int i=0; i<N_POINTS_76; i++) { 
    84         double zi = (Gauss76Z[i] + 1.0)*M_PI_4; 
     83    for(int i=0; i<GAUSS_N; i++) { 
     84        double zi = (GAUSS_Z[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 += Gauss76Wt[i] * yyy * sin_alpha; 
     97        summ += GAUSS_W[i] * yyy * sin_alpha; 
    9898    } 
    9999 

  • sasmodels/models/triaxial_ellipsoid.c

    rbecded3 r74768cb  
    2121    const double zb = M_PI_4; 
    2222    double outer = 0.0; 
    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; 
     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; 
    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<76;j++) { 
     31        for (int j=0;j<GAUSS_N;j++) { 
    3232            // translate a point in [-1,1] to a point in [0, 1] 
    33             const double usq = square(Gauss76Z[j]*um + ub); 
     33            const double usq = square(GAUSS_Z[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 += Gauss76Wt[j] * fq * fq; 
     36            inner += GAUSS_W[j] * fq * fq; 
    3737        } 
    38         outer += Gauss76Wt[i] * inner;  // correcting for dx later 
     38        outer += GAUSS_W[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.