Changeset 0b07b91 in sasmodels

Timestamp:

Nov 20, 2017 9:45:59 AM (6 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:

a261a83

Parents:

435cc89 (diff), 688d315 (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:

• sasmodels/models/core_shell_parallelepiped.py (modified) (1 diff)
• .gitignore (modified) (1 diff)
• sasmodels/compare.py (modified) (9 diffs)
• sasmodels/generate.py (modified) (1 diff)
• sasmodels/gengauss.py (added)
• sasmodels/models/barbell.c (modified) (2 diffs)
• sasmodels/models/bcc_paracrystal.c (modified) (2 diffs)
• sasmodels/models/capped_cylinder.c (modified) (3 diffs)
• sasmodels/models/core_shell_bicelle.c (modified) (1 diff)
• sasmodels/models/core_shell_bicelle_elliptical.c (modified) (3 diffs)
• sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c (modified) (4 diffs)
• sasmodels/models/core_shell_cylinder.c (modified) (2 diffs)
• sasmodels/models/core_shell_ellipsoid.c (modified) (2 diffs)
• sasmodels/models/core_shell_parallelepiped.c (modified) (3 diffs)
• sasmodels/models/cylinder.c (modified) (1 diff)
• sasmodels/models/ellipsoid.c (modified) (1 diff)
• sasmodels/models/elliptical_cylinder.c (modified) (2 diffs)
• sasmodels/models/elliptical_cylinder.py (modified) (2 diffs)
• sasmodels/models/fcc_paracrystal.c (modified) (2 diffs)
• sasmodels/models/flexible_cylinder_elliptical.c (modified) (1 diff)
• sasmodels/models/hollow_cylinder.c (modified) (1 diff)
• sasmodels/models/hollow_rectangular_prism.c (modified) (3 diffs)
• sasmodels/models/hollow_rectangular_prism_thin_walls.c (modified) (4 diffs)
• sasmodels/models/lib/gauss150.c (modified) (3 diffs)
• sasmodels/models/lib/gauss20.c (modified) (1 diff)
• sasmodels/models/lib/gauss76.c (modified) (2 diffs)
• sasmodels/models/parallelepiped.c (modified) (2 diffs)
• sasmodels/models/pringle.c (modified) (3 diffs)
• sasmodels/models/rectangular_prism.c (modified) (3 diffs)
• sasmodels/models/sc_paracrystal.c (modified) (2 diffs)
• sasmodels/models/stacked_disks.c (modified) (2 diffs)
• sasmodels/models/triaxial_ellipsoid.c (modified) (2 diffs)

sasmodels/models/core_shell_parallelepiped.py

@@ -215 +215 @@
 
 # rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models
-qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
-tests = [[{}, 0.2, 0.533149288477],
-         [{}, [0.2], [0.533149288477]],
-         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
-         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
-        ]
-del tests  # TODO: fix the tests
-del qx, qy  # not necessary to delete, but cleaner
+if 0:  # pak: model rewrite; need to update tests
+    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
+    tests = [[{}, 0.2, 0.533149288477],
+            [{}, [0.2], [0.533149288477]],
+            [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
+            [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
+            ]
+    del qx, qy  # not necessary to delete, but cleaner

.gitignore

@@ -22 +22 @@
 /example/Fit_*/
 /example/batch_fit.csv
+/sasmodels/models/lib/gauss*.c

sasmodels/compare.py

@@ -44 +44 @@
 from .direct_model import DirectModel, get_mesh
 from .convert import revert_name, revert_pars, constrain_new_to_old
-from .generate import FLOAT_RE
+from .generate import FLOAT_RE, set_integration_size
 from .weights import plot_weights
 
@@ -801 +801 @@
     return data, index
 
-def make_engine(model_info, data, dtype, cutoff):
+def make_engine(model_info, data, dtype, cutoff, ngauss=0):
     # type: (ModelInfo, Data, str, float) -> Calculator
     """
@@ -809 +809 @@
     than OpenCL.
     """
+    if ngauss:
+        set_integration_size(model_info, ngauss)
+
     if dtype == 'sasview':
         return eval_sasview(model_info, data)
@@ -1045 +1048 @@
     'poly', 'mono', 'cutoff=',
     'magnetic', 'nonmagnetic',
-    'accuracy=',
+    'accuracy=', 'ngauss=',
     'neval=',  # for timing...
 
@@ -1179 +1182 @@
         'show_weights' : False,
         'sphere'    : 0,
+        'ngauss'    : '0',
     }
     for arg in flags:
@@ -1205 +1209 @@
         elif arg.startswith('-engine='):   opts['engine'] = arg[8:]
         elif arg.startswith('-neval='):    opts['count'] = arg[7:]
+        elif arg.startswith('-ngauss='):   opts['ngauss'] = arg[8:]
         elif arg.startswith('-random='):
             opts['seed'] = int(arg[8:])
@@ -1260 +1265 @@
 
     comparison = any(PAR_SPLIT in v for v in values)
+
     if PAR_SPLIT in name:
         names = name.split(PAR_SPLIT, 2)
@@ -1272 +1278 @@
         return None
 
+    if PAR_SPLIT in opts['ngauss']:
+        opts['ngauss'] = [int(k) for k in opts['ngauss'].split(PAR_SPLIT, 2)]
+        comparison = True
+    else:
+        opts['ngauss'] = [int(opts['ngauss'])]*2
+
     if PAR_SPLIT in opts['engine']:
         opts['engine'] = opts['engine'].split(PAR_SPLIT, 2)
@@ -1290 +1302 @@
         opts['cutoff'] = [float(opts['cutoff'])]*2
 
-    base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0])
+    base = make_engine(model_info[0], data, opts['engine'][0],
+                       opts['cutoff'][0], opts['ngauss'][0])
     if comparison:
-        comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1])
+        comp = make_engine(model_info[1], data, opts['engine'][1],
+                           opts['cutoff'][1], opts['ngauss'][1])
     else:
         comp = None

sasmodels/generate.py

@@ -268 +268 @@
 """
 
+
+def set_integration_size(info, n):
+    # type: (ModelInfo, int) -> None
+    """
+    Update the model definition, replacing the gaussian integration with
+    a gaussian integration of a different size.
+
+    Note: this really ought to be a method in modelinfo, but that leads to
+    import loops.
+    """
+    if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)):
+        import os.path
+        from .gengauss import gengauss
+        path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n)
+        if not os.path.exists(path):
+            gengauss(n, path)
+        info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss')
+                        else lib for lib in info.source]
 
 def format_units(units):

sasmodels/models/barbell.c

@@ -23 +23 @@
     const double qab_r = radius_bell*qab; // Q*R*sin(theta)
     double total = 0.0;
-    for (int i = 0; i < 76; i++){
-        const double t = Gauss76Z[i]*zm + zb;
+    for (int i = 0; i < GAUSS_N; i++){
+        const double t = GAUSS_Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
         const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
-        total += Gauss76Wt[i] * Fq;
+        total += GAUSS_W[i] * Fq;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -73 +73 @@
     const double zb = M_PI_4;
     double total = 0.0;
-    for (int i = 0; i < 76; i++){
-        const double alpha= Gauss76Z[i]*zm + zb;
+    for (int i = 0; i < GAUSS_N; i++){
+        const double alpha= GAUSS_Z[i]*zm + zb;
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(alpha, sin_alpha, cos_alpha);
         const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length);
-        total += Gauss76Wt[i] * Aq * Aq * sin_alpha;
+        total += GAUSS_W[i] * Aq * Aq * sin_alpha;
     }
     // translate dx in [-1,1] to dx in [lower,upper]

sasmodels/models/bcc_paracrystal.c

@@ -81 +81 @@
 
     double outer_sum = 0.0;
-    for(int i=0; i<150; i++) {
+    for(int i=0; i<GAUSS_N; i++) {
         double inner_sum = 0.0;
-        const double theta = Gauss150Z[i]*theta_m + theta_b;
+        const double theta = GAUSS_Z[i]*theta_m + theta_b;
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
         const double qc = q*cos_theta;
         const double qab = q*sin_theta;
-        for(int j=0;j<150;j++) {
-            const double phi = Gauss150Z[j]*phi_m + phi_b;
+        for(int j=0;j<GAUSS_N;j++) {
+            const double phi = GAUSS_Z[j]*phi_m + phi_b;
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -95 +95 @@
             const double qb = qab*sin_phi;
             const double form = bcc_Zq(qa, qb, qc, dnn, d_factor);
-            inner_sum += Gauss150Wt[j] * form;
+            inner_sum += GAUSS_W[j] * form;
         }
         inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
-        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
     outer_sum *= theta_m;

sasmodels/models/capped_cylinder.c

@@ -30 +30 @@
     const double qab_r = radius_cap*qab; // Q*R*sin(theta)
     double total = 0.0;
-    for (int i=0; i<76 ;i++) {
-        const double t = Gauss76Z[i]*zm + zb;
+    for (int i=0; i<GAUSS_N; i++) {
+        const double t = GAUSS_Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
         const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
-        total += Gauss76Wt[i] * Fq;
+        total += GAUSS_W[i] * Fq;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -95 +95 @@
     const double zb = M_PI_4;
     double total = 0.0;
-    for (int i=0; i<76 ;i++) {
-        const double theta = Gauss76Z[i]*zm + zb;
+    for (int i=0; i<GAUSS_N ;i++) {
+        const double theta = GAUSS_Z[i]*zm + zb;
         double sin_theta, cos_theta; // slots to hold sincos function output
         SINCOS(theta, sin_theta, cos_theta);
@@ -103 +103 @@
         const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length);
         // scale by sin_theta for spherical coord integration
-        total += Gauss76Wt[i] * Aq * Aq * sin_theta;
+        total += GAUSS_W[i] * Aq * Aq * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]

sasmodels/models/core_shell_bicelle.c

@@ -52 +52 @@
 
     double total = 0.0;
-    for(int i=0;i<N_POINTS_76;i++) {
-        double theta = (Gauss76Z[i] + 1.0)*uplim;
+    for(int i=0;i<GAUSS_N;i++) {
+        double theta = (GAUSS_Z[i] + 1.0)*uplim;
         double sin_theta, cos_theta; // slots to hold sincos function output
         SINCOS(theta, sin_theta, cos_theta);
         double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face,
                                    halflength, sld_core, sld_face, sld_rim, sld_solvent);
-        total += Gauss76Wt[i]*fq*fq*sin_theta;
+        total += GAUSS_W[i]*fq*fq*sin_theta;
     }
 

sasmodels/models/core_shell_bicelle_elliptical.c

@@ -37 +37 @@
     //initialize integral
     double outer_sum = 0.0;
-    for(int i=0;i<76;i++) {
+    for(int i=0;i<GAUSS_N;i++) {
         //setup inner integral over the ellipsoidal cross-section
         //const double va = 0.0;
         //const double vb = 1.0;
-        //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
-        const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0;
+        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0;
         const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
         const double qab = q*sin_theta;
@@ -49 +49 @@
         const double si2 = sas_sinx_x((halfheight+thick_face)*qc);
         double inner_sum=0.0;
-        for(int j=0;j<76;j++) {
+        for(int j=0;j<GAUSS_N;j++) {
             //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
             // inner integral limits
             //const double vaj=0.0;
             //const double vbj=M_PI;
-            //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
-            const double phi = ( Gauss76Z[j] +1.0)*M_PI_2;
+            //const double phi = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+            const double phi = ( GAUSS_Z[j] +1.0)*M_PI_2;
             const double rr = sqrt(r2A - r2B*cos(phi));
             const double be1 = sas_2J1x_x(rr*qab);
@@ -61 +61 @@
             const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1;
 
-            inner_sum += Gauss76Wt[j] * fq * fq;
+            inner_sum += GAUSS_W[j] * fq * fq;
         }
         //now calculate outer integral
-        outer_sum += Gauss76Wt[i] * inner_sum;
+        outer_sum += GAUSS_W[i] * inner_sum;
     }
 

sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

@@ -7 +7 @@
         double length)
 {
-    return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length + 
+    return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length +
                  square(r_minor)*x_core*2.0*thick_face  );
 }
@@ -47 +47 @@
     //initialize integral
     double outer_sum = 0.0;
-    for(int i=0;i<76;i++) {
+    for(int i=0;i<GAUSS_N;i++) {
         //setup inner integral over the ellipsoidal cross-section
         // since we generate these lots of times, why not store them somewhere?
-        //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
-        const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0;
+        //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0;
         const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
         double inner_sum=0;
@@ -58 +58 @@
         si1 = sas_sinx_x(sinarg1);
         si2 = sas_sinx_x(sinarg2);
-        for(int j=0;j<76;j++) {
+        for(int j=0;j<GAUSS_N;j++) {
             //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
-            //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
-            const double beta = ( Gauss76Z[j] +1.0)*M_PI_2;
+            //const double beta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2;
             const double rr = sqrt(r2A - r2B*cos(beta));
             double besarg1 = q*rr*sin_alpha;
@@ -67 +67 @@
             be1 = sas_2J1x_x(besarg1);
             be2 = sas_2J1x_x(besarg2);
-            inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 +
+            inner_sum += GAUSS_W[j] *square(dr1*si1*be1 +
                                               dr2*si1*be2 +
                                               dr3*si2*be1);
         }
         //now calculate outer integral
-        outer_sum += Gauss76Wt[i] * inner_sum;
+        outer_sum += GAUSS_W[i] * inner_sum;
     }
 

sasmodels/models/core_shell_cylinder.c

@@ -30 +30 @@
     const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
     double total = 0.0;
-    for (int i=0; i<76 ;i++) {
+    for (int i=0; i<GAUSS_N ;i++) {
         // translate a point in [-1,1] to a point in [0, pi/2]
-        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
+        //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0;
         double sin_theta, cos_theta;
-        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4;
+        const double theta = GAUSS_Z[i]*M_PI_4 + M_PI_4;
         SINCOS(theta, sin_theta,  cos_theta);
         const double qab = q*sin_theta;
@@ -40 +40 @@
         const double fq = _cyl(core_vd, core_r*qab, core_h*qc)
             + _cyl(shell_vd, shell_r*qab, shell_h*qc);
-        total += Gauss76Wt[i] * fq * fq * sin_theta;
+        total += GAUSS_W[i] * fq * fq * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]

sasmodels/models/core_shell_ellipsoid.c

@@ -59 +59 @@
     const double b = 0.5;
     double total = 0.0;     //initialize intergral
-    for(int i=0;i<76;i++) {
-        const double cos_theta = Gauss76Z[i]*m + b;
+    for(int i=0;i<GAUSS_N;i++) {
+        const double cos_theta = GAUSS_Z[i]*m + b;
         const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
         double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta,
@@ -66 +66 @@
             equat_shell, polar_shell,
             sld_core_shell, sld_shell_solvent);
-        total += Gauss76Wt[i] * fq * fq;
+        total += GAUSS_W[i] * fq * fq;
     }
     total *= m;

sasmodels/models/core_shell_parallelepiped.c

@@ -80 +80 @@
     // outer integral (with gauss points), integration limits = 0, 1
     double outer_sum = 0; //initialize integral
-    for( int i=0; i<76; i++) {
-        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
+    for( int i=0; i<GAUSS_N; i++) {
+        double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 );
         double mu_proj = mu * sqrt(1.0-sigma*sigma);
 
@@ -88 +88 @@
         const double siCt = sas_sinx_x(mu * sigma * tC_over_b);
         double inner_sum = 0.0;
-        for(int j=0; j<76; j++) {
-            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
+        for(int j=0; j<GAUSS_N; j++) {
+            const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 );
             double sin_uu, cos_uu;
             SINCOS(M_PI_2*uu, sin_uu, cos_uu);
@@ -109 +109 @@
 #endif
 
-            inner_sum += Gauss76Wt[j] * f * f;
+            inner_sum += GAUSS_W[j] * f * f;
         }
         inner_sum *= 0.5;
         // now sum up the outer integral
-        outer_sum += Gauss76Wt[i] * inner_sum;
+        outer_sum += GAUSS_W[i] * inner_sum;
     }
     outer_sum *= 0.5;

sasmodels/models/cylinder.c

@@ -21 +21 @@
 
     double total = 0.0;
-    for (int i=0; i<76 ;i++) {
-        const double theta = Gauss76Z[i]*zm + zb;
+    for (int i=0; i<GAUSS_N ;i++) {
+        const double theta = GAUSS_Z[i]*zm + zb;
         double sin_theta, cos_theta; // slots to hold sincos function output
         // theta (theta,phi) the projection of the cylinder on the detector plane
         SINCOS(theta , sin_theta, cos_theta);
         const double form = fq(q*sin_theta, q*cos_theta, radius, length);
-        total += Gauss76Wt[i] * form * form * sin_theta;
+        total += GAUSS_W[i] * form * form * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]

sasmodels/models/ellipsoid.c

@@ -22 +22 @@
 
     // translate a point in [-1,1] to a point in [0, 1]
-    // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2;
+    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2;
     const double zm = 0.5;
     const double zb = 0.5;
     double total = 0.0;
-    for (int i=0;i<76;i++) {
-        const double u = Gauss76Z[i]*zm + zb;
+    for (int i=0;i<GAUSS_N;i++) {
+        const double u = GAUSS_Z[i]*zm + zb;
         const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one);
         const double f = sas_3j1x_x(q*r);
-        total += Gauss76Wt[i] * f * f;
+        total += GAUSS_W[i] * f * f;
     }
     // translate dx in [-1,1] to dx in [lower,upper]

sasmodels/models/elliptical_cylinder.c

@@ -22 +22 @@
     //initialize integral
     double outer_sum = 0.0;
-    for(int i=0;i<76;i++) {
+    for(int i=0;i<GAUSS_N;i++) {
         //setup inner integral over the ellipsoidal cross-section
-        const double cos_val = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_val = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
         const double sin_val = sqrt(1.0 - cos_val*cos_val);
         //const double arg = radius_minor*sin_val;
         double inner_sum=0;
-        for(int j=0;j<76;j++) {
-            //20 gauss points for the inner integral, increase to 76, RKH 6Nov2017
-            const double theta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+        for(int j=0;j<GAUSS_N;j++) {
+            const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
             const double r = sin_val*sqrt(rA - rB*cos(theta));
             const double be = sas_2J1x_x(q*r);
-            inner_sum += Gauss76Wt[j] * be * be;
+            inner_sum += GAUSS_W[j] * be * be;
         }
         //now calculate the value of the inner integral
@@ -40 +39 @@
         //now calculate outer integral
         const double si = sas_sinx_x(q*0.5*length*cos_val);
-        outer_sum += Gauss76Wt[i] * inner_sum * si * si;
+        outer_sum += GAUSS_W[i] * inner_sum * si * si;
     }
     outer_sum *= 0.5*(vb-va);

sasmodels/models/elliptical_cylinder.py

@@ -43 +43 @@
     P(q) = \text{scale}  <F^2> / V
 
-For 2d data the orientation of the particle is required, described using a different set 
-of angles as in the diagrams below, for further details of the calculation and angular 
+For 2d data the orientation of the particle is required, described using a different set
+of angles as in the diagrams below, for further details of the calculation and angular
 dispersions  see :ref:`orientation` .
 
@@ -121 +121 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "lib/gauss20.c",
-          "elliptical_cylinder.c"]
+source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"]
 
 demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0,

sasmodels/models/fcc_paracrystal.c

@@ -53 +53 @@
 
     double outer_sum = 0.0;
-    for(int i=0; i<150; i++) {
+    for(int i=0; i<GAUSS_N; i++) {
         double inner_sum = 0.0;
-        const double theta = Gauss150Z[i]*theta_m + theta_b;
+        const double theta = GAUSS_Z[i]*theta_m + theta_b;
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
         const double qc = q*cos_theta;
         const double qab = q*sin_theta;
-        for(int j=0;j<150;j++) {
-            const double phi = Gauss150Z[j]*phi_m + phi_b;
+        for(int j=0;j<GAUSS_N;j++) {
+            const double phi = GAUSS_Z[j]*phi_m + phi_b;
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -67 +67 @@
             const double qb = qab*sin_phi;
             const double form = fcc_Zq(qa, qb, qc, dnn, d_factor);
-            inner_sum += Gauss150Wt[j] * form;
+            inner_sum += GAUSS_W[j] * form;
         }
         inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
-        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
     outer_sum *= theta_m;

sasmodels/models/flexible_cylinder_elliptical.c

@@ -17 +17 @@
     double sum=0.0;
 
-    for(int i=0;i<N_POINTS_76;i++) {
-        const double zi = ( Gauss76Z[i] + 1.0 )*M_PI_4;
+    for(int i=0;i<GAUSS_N;i++) {
+        const double zi = ( GAUSS_Z[i] + 1.0 )*M_PI_4;
         double sn, cn;
         SINCOS(zi, sn, cn);
         const double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn);
         const double yyy = sas_2J1x_x(arg);
-        sum += Gauss76Wt[i] * yyy * yyy;
+        sum += GAUSS_W[i] * yyy * yyy;
     }
     sum *= 0.5;

sasmodels/models/hollow_cylinder.c

@@ -38 +38 @@
 
     double summ = 0.0;            //initialize intergral
-    for (int i=0;i<76;i++) {
-        const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper );
+    for (int i=0;i<GAUSS_N;i++) {
+        const double cos_theta = 0.5*( GAUSS_Z[i] * (upper-lower) + lower + upper );
         const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
         const double form = _fq(q*sin_theta, q*cos_theta,
                                 radius, thickness, length);
-        summ += Gauss76Wt[i] * form * form;
+        summ += GAUSS_W[i] * form * form;
     }
 

sasmodels/models/hollow_rectangular_prism.c

@@ -39 +39 @@
 
     double outer_sum = 0.0;
-    for(int i=0; i<76; i++) {
+    for(int i=0; i<GAUSS_N; i++) {
 
-        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );
+        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b );
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
@@ -49 +49 @@
 
         double inner_sum = 0.0;
-        for(int j=0; j<76; j++) {
+        for(int j=0; j<GAUSS_N; j++) {
 
-            const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );
+            const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b );
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -66 +66 @@
             const double AP2 = vol_core * termA2 * termB2 * termC2;
 
-            inner_sum += Gauss76Wt[j] * square(AP1-AP2);
+            inner_sum += GAUSS_W[j] * square(AP1-AP2);
         }
         inner_sum *= 0.5 * (v2b-v2a);
 
-        outer_sum += Gauss76Wt[i] * inner_sum * sin(theta);
+        outer_sum += GAUSS_W[i] * inner_sum * sin(theta);
     }
     outer_sum *= 0.5*(v1b-v1a);

sasmodels/models/hollow_rectangular_prism_thin_walls.c

@@ -1 +1 @@
 double form_volume(double length_a, double b2a_ratio, double c2a_ratio);
-double Iq(double q, double sld, double solvent_sld, double length_a, 
+double Iq(double q, double sld, double solvent_sld, double length_a,
           double b2a_ratio, double c2a_ratio);
 
@@ -29 +29 @@
     const double v2a = 0.0;
     const double v2b = M_PI_2;  //phi integration limits
-    
+
     double outer_sum = 0.0;
-    for(int i=0; i<76; i++) {
-        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );
+    for(int i=0; i<GAUSS_N; i++) {
+        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b );
 
         double sin_theta, cos_theta;
@@ -44 +44 @@
 
         double inner_sum = 0.0;
-        for(int j=0; j<76; j++) {
-            const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );
+        for(int j=0; j<GAUSS_N; j++) {
+            const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b );
 
             double sin_phi, cos_phi;
@@ -62 +62 @@
                 * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi );
 
-            inner_sum += Gauss76Wt[j] * square(AL+AT);
+            inner_sum += GAUSS_W[j] * square(AL+AT);
         }
 
         inner_sum *= 0.5 * (v2b-v2a);
-        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
 

sasmodels/models/lib/gauss150.c

@@ -7 +7 @@
  *
  */
+ #ifdef GAUSS_N
+ # undef GAUSS_N
+ # undef GAUSS_Z
+ # undef GAUSS_W
+ #endif
+ #define GAUSS_N 150
+ #define GAUSS_Z Gauss150Z
+ #define GAUSS_W Gauss150Wt
+
+// Note: using array size 152 so that it is a multiple of 4
 
 // Gaussians
-constant double Gauss150Z[150]={
+constant double Gauss150Z[152]={
         -0.9998723404457334,
         -0.9993274305065947,
@@ -159 +169 @@
         0.9983473449340834,
         0.9993274305065947,
-        0.9998723404457334
+        0.9998723404457334,
+        0.,
+        0.
 };
 
-constant double Gauss150Wt[150]={
+constant double Gauss150Wt[152]={
         0.0003276086705538,
         0.0007624720924706,
@@ -312 +324 @@
         0.0011976474864367,
         0.0007624720924706,
-        0.0003276086705538
+        0.0003276086705538,
+        0.,
+        0.
 };

sasmodels/models/lib/gauss20.c

@@ -7 +7 @@
  *
  */
+ #ifdef GAUSS_N
+ # undef GAUSS_N
+ # undef GAUSS_Z
+ # undef GAUSS_W
+ #endif
+ #define GAUSS_N 20
+ #define GAUSS_Z Gauss20Z
+ #define GAUSS_W Gauss20Wt
 
 // Gaussians

sasmodels/models/lib/gauss76.c

@@ -7 +7 @@
  *
  */
-#define N_POINTS_76 76
+ #ifdef GAUSS_N
+ # undef GAUSS_N
+ # undef GAUSS_Z
+ # undef GAUSS_W
+ #endif
+ #define GAUSS_N 76
+ #define GAUSS_Z Gauss76Z
+ #define GAUSS_W Gauss76Wt
 
 // Gaussians
-constant double Gauss76Wt[N_POINTS_76]={
+constant double Gauss76Wt[76]={
         .00126779163408536,             //0
         .00294910295364247,
@@ -89 +96 @@
 };
 
-constant double Gauss76Z[N_POINTS_76]={
+constant double Gauss76Z[76]={
         -.999505948362153,              //0
         -.997397786355355,

sasmodels/models/parallelepiped.c

@@ -23 +23 @@
     double outer_total = 0; //initialize integral
 
-    for( int i=0; i<76; i++) {
-        const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
+    for( int i=0; i<GAUSS_N; i++) {
+        const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 );
         const double mu_proj = mu * sqrt(1.0-sigma*sigma);
 
@@ -30 +30 @@
         // corresponding to angles from 0 to pi/2.
         double inner_total = 0.0;
-        for(int j=0; j<76; j++) {
-            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
+        for(int j=0; j<GAUSS_N; j++) {
+            const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 );
             double sin_uu, cos_uu;
             SINCOS(M_PI_2*uu, sin_uu, cos_uu);
             const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
             const double si2 = sas_sinx_x(mu_proj * cos_uu);
-            inner_total += Gauss76Wt[j] * square(si1 * si2);
+            inner_total += GAUSS_W[j] * square(si1 * si2);
         }
         inner_total *= 0.5;
 
         const double si = sas_sinx_x(mu * c_scaled * sigma);
-        outer_total += Gauss76Wt[i] * inner_total * si * si;
+        outer_total += GAUSS_W[i] * inner_total * si * si;
     }
     outer_total *= 0.5;

sasmodels/models/pringle.c

@@ -29 +29 @@
     double sumC = 0.0;          // initialize integral
     double r;
-    for (int i=0; i < 76; i++) {
-        r = Gauss76Z[i]*zm + zb;
+    for (int i=0; i < GAUSS_N; i++) {
+        r = GAUSS_Z[i]*zm + zb;
 
         const double qrs = r*q_sin_psi;
         const double qrrc = r*r*q_cos_psi;
 
-        double y = Gauss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs);
+        double y = GAUSS_W[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs);
         double S, C;
         SINCOS(alpha*qrrc, S, C);
@@ -86 +86 @@
 
     double sum = 0.0;
-    for (int i = 0; i < 76; i++) {
-        double psi = Gauss76Z[i]*zm + zb;
+    for (int i = 0; i < GAUSS_N; i++) {
+        double psi = GAUSS_Z[i]*zm + zb;
         double sin_psi, cos_psi;
         SINCOS(psi, sin_psi, cos_psi);
@@ -93 +93 @@
         double sinc_term = square(sas_sinx_x(q * thickness * cos_psi / 2.0));
         double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term;
-        sum += Gauss76Wt[i] * pringle_kernel;
+        sum += GAUSS_W[i] * pringle_kernel;
     }
 

sasmodels/models/rectangular_prism.c

@@ -28 +28 @@
 
     double outer_sum = 0.0;
-    for(int i=0; i<76; i++) {
-        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );
+    for(int i=0; i<GAUSS_N; i++) {
+        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b );
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
@@ -36 +36 @@
 
         double inner_sum = 0.0;
-        for(int j=0; j<76; j++) {
-            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );
+        for(int j=0; j<GAUSS_N; j++) {
+            double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b );
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -45 +45 @@
             const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi);
             const double AP = termA * termB * termC;
-            inner_sum += Gauss76Wt[j] * AP * AP;
+            inner_sum += GAUSS_W[j] * AP * AP;
         }
         inner_sum = 0.5 * (v2b-v2a) * inner_sum;
-        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
 

sasmodels/models/sc_paracrystal.c

@@ -54 +54 @@
 
     double outer_sum = 0.0;
-    for(int i=0; i<150; i++) {
+    for(int i=0; i<GAUSS_N; i++) {
         double inner_sum = 0.0;
-        const double theta = Gauss150Z[i]*theta_m + theta_b;
+        const double theta = GAUSS_Z[i]*theta_m + theta_b;
         double sin_theta, cos_theta;
         SINCOS(theta, sin_theta, cos_theta);
         const double qc = q*cos_theta;
         const double qab = q*sin_theta;
-        for(int j=0;j<150;j++) {
-            const double phi = Gauss150Z[j]*phi_m + phi_b;
+        for(int j=0;j<GAUSS_N;j++) {
+            const double phi = GAUSS_Z[j]*phi_m + phi_b;
             double sin_phi, cos_phi;
             SINCOS(phi, sin_phi, cos_phi);
@@ -68 +68 @@
             const double qb = qab*sin_phi;
             const double form = sc_Zq(qa, qb, qc, dnn, d_factor);
-            inner_sum += Gauss150Wt[j] * form;
+            inner_sum += GAUSS_W[j] * form;
         }
         inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
-        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta;
+        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
     }
     outer_sum *= theta_m;

sasmodels/models/stacked_disks.c

@@ -81 +81 @@
     double halfheight = 0.5*thick_core;
 
-    for(int i=0; i<N_POINTS_76; i++) {
-        double zi = (Gauss76Z[i] + 1.0)*M_PI_4;
+    for(int i=0; i<GAUSS_N; i++) {
+        double zi = (GAUSS_Z[i] + 1.0)*M_PI_4;
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(zi, sin_alpha, cos_alpha);
@@ -95 +95 @@
                            solvent_sld,
                            d);
-        summ += Gauss76Wt[i] * yyy * sin_alpha;
+        summ += GAUSS_W[i] * yyy * sin_alpha;
     }
 

sasmodels/models/triaxial_ellipsoid.c

@@ -21 +21 @@
     const double zb = M_PI_4;
     double outer = 0.0;
-    for (int i=0;i<76;i++) {
-        //const double u = Gauss76Z[i]*(upper-lower)/2 + (upper + lower)/2;
-        const double phi = Gauss76Z[i]*zm + zb;
+    for (int i=0;i<GAUSS_N;i++) {
+        //const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper + lower)/2;
+        const double phi = GAUSS_Z[i]*zm + zb;
         const double pa_sinsq_phi = pa*square(sin(phi));
 
@@ -29 +29 @@
         const double um = 0.5;
         const double ub = 0.5;
-        for (int j=0;j<76;j++) {
+        for (int j=0;j<GAUSS_N;j++) {
             // translate a point in [-1,1] to a point in [0, 1]
-            const double usq = square(Gauss76Z[j]*um + ub);
+            const double usq = square(GAUSS_Z[j]*um + ub);
             const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq);
             const double fq = sas_3j1x_x(q*r);
-            inner += Gauss76Wt[j] * fq * fq;
+            inner += GAUSS_W[j] * fq * fq;
         }
-        outer += Gauss76Wt[i] * inner;  // correcting for dx later
+        outer += GAUSS_W[i] * inner;  // correcting for dx later
     }
     // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi