Index: .gitignore
===================================================================
--- .gitignore	(revision 8a5f0216970a47f3e897a508054840132bd87909)
+++ .gitignore	(revision 9248bf738cf2cd8f3c0ef73ee96e3bae5099f6a1)
@@ -22,2 +22,3 @@
 /example/Fit_*/
 /example/batch_fit.csv
+/sasmodels/models/lib/gauss*.c
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/compare.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
@@ -42,5 +42,5 @@
 from .data import plot_theory, empty_data1D, empty_data2D, load_data
 from .direct_model import DirectModel, get_mesh
-from .generate import FLOAT_RE
+from .generate import FLOAT_RE, set_integration_size
 from .weights import plot_weights
 
@@ -706,5 +706,5 @@
     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
     """
@@ -714,4 +714,7 @@
     than OpenCL.
     """
+    if ngauss:
+        set_integration_size(model_info, ngauss)
+
     if dtype is None or not dtype.endswith('!'):
         return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff)
@@ -954,5 +957,5 @@
     'poly', 'mono', 'cutoff=',
     'magnetic', 'nonmagnetic',
-    'accuracy=',
+    'accuracy=', 'ngauss=',
     'neval=',  # for timing...
 
@@ -1089,4 +1092,5 @@
         'show_weights' : False,
         'sphere'    : 0,
+        'ngauss'    : '0',
     }
     for arg in flags:
@@ -1115,4 +1119,5 @@
         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:])
@@ -1169,4 +1174,5 @@
 
     comparison = any(PAR_SPLIT in v for v in values)
+
     if PAR_SPLIT in name:
         names = name.split(PAR_SPLIT, 2)
@@ -1181,4 +1187,10 @@
         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)
@@ -1199,7 +1211,9 @@
         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
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/generate.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
@@ -270,4 +270,22 @@
 """
 
+
+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):
Index: sasmodels/gengauss.py
===================================================================
--- sasmodels/gengauss.py	(revision ff31782cc25cfbfea355acdd6d0204c959b41dcc)
+++ sasmodels/gengauss.py	(revision ff31782cc25cfbfea355acdd6d0204c959b41dcc)
@@ -0,0 +1,43 @@
+#!/usr/bin/env python
+from __future__ import division, print_function
+
+import sys
+
+import numpy as np
+from numpy.polynomial.legendre import leggauss
+
+def gengauss(n, path):
+    z, w = leggauss(n)
+
+    # Make sure array size is a multiple of 4
+    if n%4:
+        array_size = n + (4 - n%4)
+        z, w = [np.hstack((v, [0.]*(4-n%4))) for v in (z, w)]
+    else:
+        array_size = n
+
+    with open(path, "w") as fid:
+        fid.write("""\
+// Generated by sasmodels.gengauss.gengauss(%d)
+
+#ifdef GAUSS_N
+# undef GAUSS_N
+# undef GAUSS_Z
+# undef GAUSS_W
+#endif
+#define GAUSS_N %d
+#define GAUSS_Z Gauss%dZ
+#define GAUSS_W Gauss%dWt
+
+"""%(n, n, n, n))
+
+        if array_size != n:
+            fid.write("// Note: using array size %d so that it is a multiple of 4\n\n"%array_size)
+
+        fid.write("constant double Gauss%dWt[%d]={\n"%(n, array_size))
+        fid.write(",\n".join("\t% .15e"%v for v in w))
+        fid.write("\n};\n")
+
+        fid.write("constant double Gauss%dZ[%d]={\n"%(n, array_size))
+        fid.write(",\n".join("\t% .15e"%v for v in z))
+        fid.write("\n};")
Index: sasmodels/models/barbell.c
===================================================================
--- sasmodels/models/barbell.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/barbell.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -23,10 +23,10 @@
     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,10 +73,10 @@
     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]
Index: sasmodels/models/bcc_paracrystal.c
===================================================================
--- sasmodels/models/bcc_paracrystal.c	(revision ea60e084e3e40eeb3abfc5f728359135acaed020)
+++ sasmodels/models/bcc_paracrystal.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -81,13 +81,13 @@
 
     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,8 +95,8 @@
             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;
Index: sasmodels/models/capped_cylinder.c
===================================================================
--- sasmodels/models/capped_cylinder.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/capped_cylinder.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -30,10 +30,10 @@
     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,6 +95,6 @@
     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,5 +103,5 @@
         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]
Index: sasmodels/models/core_shell_bicelle.c
===================================================================
--- sasmodels/models/core_shell_bicelle.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/core_shell_bicelle.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -52,11 +52,11 @@
 
     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;
     }
 
Index: sasmodels/models/core_shell_bicelle_elliptical.c
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical.c	(revision 82592da3721b7e2098cf6346cf224814189fadb3)
+++ sasmodels/models/core_shell_bicelle_elliptical.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -37,10 +37,10 @@
     //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,11 +49,11 @@
         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,8 +61,8 @@
             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;
     }
 
Index: sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c
===================================================================
--- sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c	(revision 82592da3721b7e2098cf6346cf224814189fadb3)
+++ sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -7,5 +7,5 @@
         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,9 +47,9 @@
     //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,8 +58,8 @@
         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,10 +67,10 @@
             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;
     }
 
Index: sasmodels/models/core_shell_cylinder.c
===================================================================
--- sasmodels/models/core_shell_cylinder.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/core_shell_cylinder.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -30,9 +30,9 @@
     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,5 +40,5 @@
         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]
Index: sasmodels/models/core_shell_ellipsoid.c
===================================================================
--- sasmodels/models/core_shell_ellipsoid.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/core_shell_ellipsoid.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -59,6 +59,6 @@
     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,5 +66,5 @@
             equat_shell, polar_shell,
             sld_core_shell, sld_shell_solvent);
-        total += Gauss76Wt[i] * fq * fq;
+        total += GAUSS_W[i] * fq * fq;
     }
     total *= m;
Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision 4493288f50c1e858772d5058d6ca027bfa72f50e)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision 4493288f50c1e858772d5058d6ca027bfa72f50e)
@@ -60,6 +60,6 @@
     // outer integral (with gauss points), integration limits = 0, 1
     double outer_sum = 0; //initialize integral
-    for( int i=0; i<76; i++) {
-        const double cos_alpha = 0.5 * ( Gauss76Z[i] + 1.0 );
+    for( int i=0; i<GAUSS_N; i++) {
+        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 );
         const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha);
 
@@ -68,6 +68,6 @@
         const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q);
         double inner_sum = 0.0;
-        for(int j=0; j<76; j++) {
-            const double beta = 0.5 * ( Gauss76Z[j] + 1.0 );
+        for(int j=0; j<GAUSS_N; j++) {
+            const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 );
             double sin_beta, cos_beta;
             SINCOS(M_PI_2*beta, sin_beta, cos_beta);
@@ -89,9 +89,9 @@
 #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;
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/cylinder.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -21,11 +21,11 @@
 
     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]
Index: sasmodels/models/ellipsoid.c
===================================================================
--- sasmodels/models/ellipsoid.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/ellipsoid.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -22,13 +22,13 @@
 
     // 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]
Index: sasmodels/models/elliptical_cylinder.c
===================================================================
--- sasmodels/models/elliptical_cylinder.c	(revision 82592da3721b7e2098cf6346cf224814189fadb3)
+++ sasmodels/models/elliptical_cylinder.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -22,16 +22,15 @@
     //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,5 +39,5 @@
         //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);
Index: sasmodels/models/elliptical_cylinder.py
===================================================================
--- sasmodels/models/elliptical_cylinder.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
+++ sasmodels/models/elliptical_cylinder.py	(revision 2d81cfe3be42579d563407fb9275b7b516fee234)
@@ -121,6 +121,5 @@
 # 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,
Index: sasmodels/models/fcc_paracrystal.c
===================================================================
--- sasmodels/models/fcc_paracrystal.c	(revision f728001e783d4f2cad1f52bfa5374c457648c268)
+++ sasmodels/models/fcc_paracrystal.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -53,13 +53,13 @@
 
     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,8 +67,8 @@
             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;
Index: sasmodels/models/flexible_cylinder_elliptical.c
===================================================================
--- sasmodels/models/flexible_cylinder_elliptical.c	(revision 592343f7d47ec3aad7ffa1ca242acb1346553ede)
+++ sasmodels/models/flexible_cylinder_elliptical.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -17,11 +17,11 @@
     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;
Index: sasmodels/models/hollow_cylinder.c
===================================================================
--- sasmodels/models/hollow_cylinder.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/hollow_cylinder.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -38,10 +38,10 @@
 
     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;
     }
 
Index: sasmodels/models/hollow_rectangular_prism.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism.c	(revision 8de147789a271b990dd7ffc4244215e5def052c8)
+++ sasmodels/models/hollow_rectangular_prism.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -39,7 +39,7 @@
 
     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,7 +49,7 @@
 
         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,9 +66,9 @@
             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);
Index: sasmodels/models/hollow_rectangular_prism_thin_walls.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision ab2aea81f7fb392e5f24f1c76dd989914d8f7022)
+++ sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -1,4 +1,4 @@
 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,8 +29,8 @@
     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,6 +44,6 @@
 
         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,9 +62,9 @@
                 * ( 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;
     }
 
Index: sasmodels/models/lib/gauss150.c
===================================================================
--- sasmodels/models/lib/gauss150.c	(revision 994d77f9b713691e41106d1f9af9661dfb2dc419)
+++ sasmodels/models/lib/gauss150.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -7,7 +7,17 @@
  *
  */
+ #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,8 +169,10 @@
   	0.9983473449340834,
   	0.9993274305065947,
-  	0.9998723404457334
+  	0.9998723404457334,
+  	0.,
+  	0.
 };
 
-constant double Gauss150Wt[150]={
+constant double Gauss150Wt[152]={
   	0.0003276086705538,
   	0.0007624720924706,
@@ -312,4 +324,6 @@
   	0.0011976474864367,
   	0.0007624720924706,
-  	0.0003276086705538
+  	0.0003276086705538,
+  	0.,
+  	0.
 };
Index: sasmodels/models/lib/gauss20.c
===================================================================
--- sasmodels/models/lib/gauss20.c	(revision 994d77f9b713691e41106d1f9af9661dfb2dc419)
+++ sasmodels/models/lib/gauss20.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -7,4 +7,12 @@
  *
  */
+ #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
Index: sasmodels/models/lib/gauss76.c
===================================================================
--- sasmodels/models/lib/gauss76.c	(revision 66d119fbbe121774089e17e3086b4809c9cf2a80)
+++ sasmodels/models/lib/gauss76.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -7,8 +7,15 @@
  *
  */
-#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,5 +96,5 @@
 };
 
-constant double Gauss76Z[N_POINTS_76]={
+constant double Gauss76Z[76]={
 	-.999505948362153,		//0
 	-.997397786355355,
Index: sasmodels/models/parallelepiped.c
===================================================================
--- sasmodels/models/parallelepiped.c	(revision 9b7b23fffe93c3cd39b094e9e17d0bf62c39b2bc)
+++ sasmodels/models/parallelepiped.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -23,6 +23,6 @@
     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,16 +30,16 @@
         // 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;
Index: sasmodels/models/pringle.c
===================================================================
--- sasmodels/models/pringle.c	(revision 1e7b0db0b7b4315c62b59e1d659332bcbc962b8f)
+++ sasmodels/models/pringle.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -29,11 +29,11 @@
     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,6 +86,6 @@
 
     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,5 +93,5 @@
         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;
     }
 
Index: sasmodels/models/rectangular_prism.c
===================================================================
--- sasmodels/models/rectangular_prism.c	(revision 8de147789a271b990dd7ffc4244215e5def052c8)
+++ sasmodels/models/rectangular_prism.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -28,6 +28,6 @@
 
     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,6 +36,6 @@
 
         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,8 +45,8 @@
             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;
     }
 
Index: sasmodels/models/sc_paracrystal.c
===================================================================
--- sasmodels/models/sc_paracrystal.c	(revision f728001e783d4f2cad1f52bfa5374c457648c268)
+++ sasmodels/models/sc_paracrystal.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -54,13 +54,13 @@
 
     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,8 +68,8 @@
             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;
Index: sasmodels/models/stacked_disks.c
===================================================================
--- sasmodels/models/stacked_disks.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/stacked_disks.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -81,6 +81,6 @@
     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,5 +95,5 @@
                            solvent_sld,
                            d);
-        summ += Gauss76Wt[i] * yyy * sin_alpha;
+        summ += GAUSS_W[i] * yyy * sin_alpha;
     }
 
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 74768cb764a04223ec9e992e316cb8c211adfe34)
@@ -21,7 +21,7 @@
     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,12 +29,12 @@
         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
