Changeset d4db147 in sasmodels


Ignore:
Timestamp:
Nov 6, 2017 4:11:48 PM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
9248bf7
Parents:
ff31782 (diff), d8ac2ad (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'ticket-786' into generic_integration_loop

Files:
1 added
1 deleted
30 edited

Legend:

Unmodified
Added
Removed

  • doc/guide/orientation/orientation.rst

    r3d40839 r82592da  
    6262care for large ranges of angle. 
    6363 
    64 Note that the form factors for asymmetric particles are also performing 
    65 numerical integrations over one or more variables, so care should be taken, 
    66 especially with very large particles or more extreme aspect ratios. Users can 
    67 experiment with the values of *Npts* and *Nsigs*, the number of steps used in the 
    68 integration and the range spanned in number of standard deviations. The 
    69 standard deviation is entered in units of degrees. For a rectangular 
    70 (uniform) distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard 
    71 deviations (this may be changed soon). 
     64.. note:: 
     65    Note that the form factors for oriented particles are also performing 
     66    numerical integrations over one or more variables, so care should be taken, 
     67    especially with very large particles or more extreme aspect ratios. In such  
     68    cases results may not be accurate, particularly at very high Q, unless the model 
     69    has been specifically coded to use limiting forms of the scattering equations. 
     70     
     71    For best numerical results keep the $\theta$ distribution narrower than the $\phi$  
     72    distribution. Thus for asymmetric particles, such as elliptical_cylinder, you may  
     73    need to reorder the sizes of the three axes to acheive the desired result.  
     74    This is due to the issues of mapping a rectangular distribution onto the  
     75    surface of a sphere. 
    7276 
    73 Where appropriate, for best numerical results, keep $a < b < c$ and the 
    74 $\theta$ distribution narrower than the $\phi$ distribution. 
     77Users can experiment with the values of *Npts* and *Nsigs*, the number of steps  
     78used in the integration and the range spanned in number of standard deviations.  
     79The standard deviation is entered in units of degrees. For a "rectangular"  
     80distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations.  
     81The new "uniform" distribution avoids this by letting you directly specify the  
     82half width. 
     83 
     84The angular distributions will be truncated outside of the range -180 to +180  
     85degrees, so beware of using saying a broad Gaussian distribution with large value 
     86of *Nsigs*, as the array of *Npts* may be truncated to many fewer points than would  
     87give a good integration,as well as becoming rather meaningless. (At some point  
     88in the future the actual polydispersity arrays may be made available to the user  
     89for inspection.) 
    7590 
    7691Some more detailed technical notes are provided in the developer section of 
    … …  
    7994*Document History* 
    8095 
    81 | 2017-10-27 Richard Heenan 
     96| 2017-11-06 Richard Heenan 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r74768cb rd4db147  
    9292 
    9393    // Compute effective radius in rotated coordinates 
    94     const double qr_hat = sqrt(square(r_major*qa) + square(r_minor*qb)); 
    95     const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa) 
    96                                    + square((r_minor+thick_rim)*qb)); 
     94    const double qr_hat = sqrt(square(r_major*qb) + square(r_minor*qa)); 
     95    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qb) 
     96                                   + square((r_minor+thick_rim)*qa)); 
    9797    const double be1 = sas_2J1x_x( qr_hat ); 
    9898    const double be2 = sas_2J1x_x( qrshell_hat ); 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r74768cb rd4db147  
    9191          double sigma) 
    9292{ 
    93     // THIS NEEDS TESTING 
     93    // integrated 2d seems to match 1d reasonably well, except perhaps at very high Q 
    9494    // Vol1,2,3 and dr1,2,3 are now for Vcore, Vcore+rim, Vcore+face, 
    9595    const double dr1 = -rhor - rhoh + rhoc + rhosolv; 
    … …  
    103103 
    104104    // Compute effective radius in rotated coordinates 
    105     const double qr_hat = sqrt(square(r_major*qa) + square(r_minor*qb)); 
     105    const double qr_hat = sqrt(square(r_major*qb) + square(r_minor*qa)); 
    106106    // does this need to be changed for the "missing corners" where there there is no "belt" ? 
    107     const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa) 
    108                                    + square((r_minor+thick_rim)*qb)); 
     107    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qb) 
     108                                   + square((r_minor+thick_rim)*qa)); 
    109109    const double be1 = sas_2J1x_x( qr_hat ); 
    110110    const double be2 = sas_2J1x_x( qrshell_hat ); 

  • sasmodels/models/elliptical_cylinder.c

    r74768cb rd4db147  
    6060    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
    6161    // Given:    radius_major = r_ratio * radius_minor 
    62     const double qr = radius_minor*sqrt(square(r_ratio*qa) + square(qb)); 
     62    const double qr = radius_minor*sqrt(square(r_ratio*qb) + square(qa)); 
    6363    const double be = sas_2J1x_x(qr); 
    6464    const double si = sas_sinx_x(qc*0.5*length); 

  • 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

    r4991048 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_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.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.