Index: doc/guide/index.rst
===================================================================
--- doc/guide/index.rst	(revision 19dcb93366741a275baa8ad855ce8a310dca78f8)
+++ doc/guide/index.rst	(revision bb6f0f3f04f718187ac81a56563b8bb7a1d216af)
@@ -1,5 +1,15 @@
-**********************
-SAS Model Organization
-**********************
+**********
+SAS Models
+**********
+
+Small angle X-ray and Neutron (SAXS and SANS) scattering examines the
+scattering patterns produced by a beam travelling through the sample
+and scattering at low angles.  The scattering is computed as a function
+of $q_x$ and $q_y$, which for a given beam wavelength corresponds to
+particular scattering angles. Each pixel on the detector corresponds to
+a different scattering angle. If the sample is unoriented, the scattering
+pattern will appear as rings on the detector.  In this case, a circular
+average can be taken with 1-dimension data at $q = \surd (q_x^2 + q_y^2)$
+compared to the orientationally averaged SAS scattering pattern.
 
 Models have certain features in common.
Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision d15a908796935025260dc12a6d16932a079552d3)
+++ sasmodels/convert.py	(revision 5054e805a7c53aa18f659a0c1613f527f05caab3)
@@ -13,4 +13,5 @@
     'gauss_lorentz_gel',
     'be_polyelectrolyte',
+    'correlation_length',
 ]
 
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision eafc9fa8b200a3b754e7dafb960c8891361eacef)
+++ sasmodels/core.py	(revision d18582e2d66b12afb03e1c3a7c0af6e4449c6bfe)
@@ -73,5 +73,5 @@
     return True
 
-def load_model(model_definition, dtype="single", platform="ocl"):
+def load_model(model_definition, dtype=None, platform="ocl"):
     """
     Prepare the model for the default execution platform.
@@ -87,5 +87,6 @@
     for the calculation. Any valid numpy single or double precision identifier
     is valid, such as 'single', 'f', 'f32', or np.float32 for single, or
-    'double', 'd', 'f64'  and np.float64 for double.
+    'double', 'd', 'f64'  and np.float64 for double.  If *None*, then use
+    'single' unless the model defines single=False.
 
     *platform* should be "dll" to force the dll to be used for C models,
@@ -94,4 +95,6 @@
     if isstr(model_definition):
         model_definition = load_model_definition(model_definition)
+    if dtype is None:
+        dtype = 'single' if getattr(model_definition, 'single', True) else 'double'
     source, info = generate.make(model_definition)
     if callable(info.get('Iq', None)):
Index: sasmodels/data.py
===================================================================
--- sasmodels/data.py	(revision 5c962dff13d0a0a2f6d779be8df775f540958dbc)
+++ sasmodels/data.py	(revision d18582e2d66b12afb03e1c3a7c0af6e4449c6bfe)
@@ -242,5 +242,5 @@
 
 
-def empty_data1D(q, resolution=0.05):
+def empty_data1D(q, resolution=0.0):
     """
     Create empty 1D data using the given *q* as the x value.
@@ -252,4 +252,5 @@
     #dIq = np.sqrt(Iq)
     Iq, dIq = None, None
+    q = np.asarray(q)
     data = Data1D(q, Iq, dx=resolution * q, dy=dIq)
     data.filename = "fake data"
@@ -257,5 +258,5 @@
 
 
-def empty_data2D(qx, qy=None, resolution=0.05):
+def empty_data2D(qx, qy=None, resolution=0.0):
     """
     Create empty 2D data using the given mesh.
@@ -267,4 +268,5 @@
     if qy is None:
         qy = qx
+    qx, qy = np.asarray(qx), np.asarray(qy)
     # 5% dQ/Q resolution
     Qx, Qy = np.meshgrid(qx, qy)
Index: sasmodels/direct_model.py
===================================================================
--- sasmodels/direct_model.py	(revision eafc9fa8b200a3b754e7dafb960c8891361eacef)
+++ sasmodels/direct_model.py	(revision d18582e2d66b12afb03e1c3a7c0af6e4449c6bfe)
@@ -234,5 +234,5 @@
 
     model_definition = load_model_definition(model_name)
-    model = load_model(model_definition, dtype='single')
+    model = load_model(model_definition)
     calculator = DirectModel(data, model)
     pars = dict((k, float(v))
Index: sasmodels/kernel_template.c
===================================================================
--- sasmodels/kernel_template.c	(revision 9c79c32da6740e2bbac1a67d4d4794a71533de93)
+++ sasmodels/kernel_template.c	(revision caf768d5b2c191061248bcdc8e74f6e85385e393)
@@ -16,10 +16,12 @@
      using namespace std;
      #if defined(_MSC_VER)
-     #   define kernel extern "C" __declspec( dllexport )
+         #include <float.h>
+         #define kernel extern "C" __declspec( dllexport )
          inline double trunc(double x) { return x>=0?floor(x):-floor(-x); }
-	 inline double fmin(double x, double y) { return x>y ? y : x; }
-	 inline double fmax(double x, double y) { return x<y ? y : x; }
+	     inline double fmin(double x, double y) { return x>y ? y : x; }
+	     inline double fmax(double x, double y) { return x<y ? y : x; }
+	     inline double isnan(double x) { return _isnan(x); }
      #else
-     #   define kernel extern "C"
+         #define kernel extern "C"
      #endif
      inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); }
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision eafc9fa8b200a3b754e7dafb960c8891361eacef)
+++ sasmodels/kernelcl.py	(revision e6a5556cf35ef5c00c42dc56786102efc2b7d5c1)
@@ -172,6 +172,6 @@
         #self.data_boundary = max(d.min_data_type_align_size
         #                         for d in self.context.devices)
-        self.queues = [cl.CommandQueue(self.context, d)
-                       for d in self.context.devices]
+        self.queues = [cl.CommandQueue(context, context.devices[0])
+                       for context in self.context]
         self.compiled = {}
 
@@ -181,5 +181,23 @@
         """
         dtype = generate.F32 if dtype == 'fast' else np.dtype(dtype)
-        return all(has_type(d, dtype) for d in self.context.devices)
+        return any(has_type(d, dtype)
+                   for context in self.context
+                   for d in context.devices)
+
+    def get_queue(self, dtype):
+        """
+        Return a command queue for the kernels of type dtype.
+        """
+        for context, queue in zip(self.context, self.queues):
+            if all(has_type(d, dtype) for d in context.devices):
+                return queue
+
+    def get_context(self, dtype):
+        """
+        Return a OpenCL context for the kernels of type dtype.
+        """
+        for context, queue in zip(self.context, self.queues):
+            if all(has_type(d, dtype) for d in context.devices):
+                return context
 
     def _create_some_context(self):
@@ -190,5 +208,5 @@
         """
         try:
-            self.context = cl.create_some_context(interactive=False)
+            self.context = [cl.create_some_context(interactive=False)]
         except Exception as exc:
             warnings.warn(str(exc))
@@ -204,5 +222,5 @@
             #print("compiling",name)
             dtype = np.dtype(dtype)
-            program = compile_model(self.context, source, dtype, fast)
+            program = compile_model(self.get_context(dtype), source, dtype, fast)
             self.compiled[key] = program
         return self.compiled[key]
@@ -218,18 +236,51 @@
 def _get_default_context():
     """
-    Get an OpenCL context, preferring GPU over CPU.
-    """
-    default = None
+    Get an OpenCL context, preferring GPU over CPU, and preferring Intel
+    drivers over AMD drivers.
+    """
+    # Note: on mobile devices there is automatic clock scaling if either the
+    # CPU or the GPU is underutilized; probably doesn't affect us, but we if
+    # it did, it would mean that putting a busy loop on the CPU while the GPU
+    # is running may increase throughput.
+    #
+    # Macbook pro, base install:
+    #     {'Apple': [Intel CPU, NVIDIA GPU]}
+    # Macbook pro, base install:
+    #     {'Apple': [Intel CPU, Intel GPU]}
+    # 2 x nvidia 295 with Intel and NVIDIA opencl drivers installed
+    #     {'Intel': [CPU], 'NVIDIA': [GPU, GPU, GPU, GPU]}
+    gpu, cpu = None, None
     for platform in cl.get_platforms():
+        # AMD provides a much weaker CPU driver than Intel/Apple, so avoid it.
+        # If someone has bothered to install the AMD/NVIDIA drivers, prefer them over the integrated
+        # graphics driver that may have been supplied with the CPU chipset.
+        preferred_cpu = platform.vendor.startswith('Intel') or platform.vendor.startswith('Apple')
+        preferred_gpu = platform.vendor.startswith('Advanced') or platform.vendor.startswith('NVIDIA')
         for device in platform.get_devices():
             if device.type == cl.device_type.GPU:
-                return cl.Context([device])
-            if default is None:
-                default = device
-
-    if not default:
-        raise RuntimeError("OpenCL device not found")
-
-    return cl.Context([default])
+                # If the existing type is not GPU then it will be CUSTOM or ACCELERATOR,
+                # so don't override it.
+                if gpu is None or (preferred_gpu and gpu.type == cl.device_type.GPU):
+                    gpu = device
+            elif device.type == cl.device_type.CPU:
+                if cpu is None or preferred_cpu:
+                    cpu = device
+            else:
+                # System has cl.device_type.ACCELERATOR or cl.device_type.CUSTOM
+                # Intel Phi for example registers as an accelerator
+                # Since the user installed a custom device on their system and went through the
+                # pain of sorting out OpenCL drivers for it, lets assume they really do want to
+                # use it as their primary compute device.
+                gpu = device
+
+    # order the devices by gpu then by cpu; when searching for an available device by data type they
+    # will be checked in this order, which means that if the gpu supports double then the cpu will never
+    # be used (though we may make it possible to explicitly request the cpu at some point).
+    devices = []
+    if gpu is not None:
+        devices.append(gpu)
+    if cpu is not None:
+        devices.append(cpu)
+    return [cl.Context([d]) for d in devices]
 
 
@@ -314,6 +365,7 @@
         # architectures tested so far.
         self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors]
+        context = env.get_context(self.dtype)
         self.q_buffers = [
-            cl.Buffer(env.context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)
+            cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)
             for q in self.q_vectors
         ]
@@ -363,11 +415,12 @@
         # Note: res may be shorter than res_b if global_size != nq
         env = environment()
-        self.loops_b = [cl.Buffer(env.context, mf.READ_WRITE,
-                                  2 * MAX_LOOPS * q_input.dtype.itemsize)
-                        for _ in env.queues]
-        self.res_b = [cl.Buffer(env.context, mf.READ_WRITE,
-                                q_input.global_size[0] * q_input.dtype.itemsize)
-                      for _ in env.queues]
+        self.queue = env.get_queue(dtype)
+        self.loops_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
+                                 2 * MAX_LOOPS * q_input.dtype.itemsize)
+        self.res_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
+                               q_input.global_size[0] * q_input.dtype.itemsize)
         self.q_input = q_input
+
+        self._need_release = [self.loops_b, self.res_b, self.q_input]
 
     def __call__(self, fixed_pars, pd_pars, cutoff=1e-5):
@@ -377,7 +430,5 @@
                 else np.float32)  # will never get here, so use np.float32
 
-        device_num = 0
-        queuei = environment().queues[device_num]
-        res_bi = self.res_b[device_num]
+        res_bi = self.res_b
         nq = np.uint32(self.q_input.nq)
         if pd_pars:
@@ -394,6 +445,6 @@
                 raise ValueError("too many polydispersity points")
 
-            loops_bi = self.loops_b[device_num]
-            cl.enqueue_copy(queuei, loops_bi, loops)
+            loops_bi = self.loops_b
+            cl.enqueue_copy(self.queue, loops_bi, loops)
             loops_l = cl.LocalMemory(len(loops.data))
             #ctx = environment().context
@@ -404,6 +455,6 @@
         fixed = [real(p) for p in fixed_pars]
         args = self.q_input.q_buffers + [res_bi, nq] + dispersed + fixed
-        self.kernel(queuei, self.q_input.global_size, None, *args)
-        cl.enqueue_copy(queuei, self.res, res_bi)
+        self.kernel(self.queue, self.q_input.global_size, None, *args)
+        cl.enqueue_copy(self.queue, self.res, res_bi)
 
         return self.res
@@ -413,11 +464,7 @@
         Release resources associated with the kernel.
         """
-        for b in self.loops_b:
-            b.release()
-        self.loops_b = []
-        for b in self.res_b:
-            b.release()
-        self.res_b = []
-        self.q_input.release()
+        for v in self._need_release:
+            v.release()
+        self._need_release = []
 
     def __del__(self):
Index: sasmodels/model_test.py
===================================================================
--- sasmodels/model_test.py	(revision 5c962dff13d0a0a2f6d779be8df775f540958dbc)
+++ sasmodels/model_test.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -100,7 +100,11 @@
                 test_name = "Model: %s, Kernel: OpenCL"%model_name
                 test_method_name = "test_%s_opencl" % model_name
+                # Using dtype=None so that the models that are only
+                # correct for double precision are not tested using
+                # single precision.  The choice is determined by the
+                # presence of *single=False* in the model file.
                 test = ModelTestCase(test_name, model_definition,
                                      test_method_name,
-                                     platform="ocl", dtype='single')
+                                     platform="ocl", dtype=None)
                 #print("defining", test_name)
                 suite.addTest(test)
@@ -157,5 +161,5 @@
                     ## values an error.  Only do so for the "dll" tests
                     ## to reduce noise from both opencl and dll, and because
-                    ## python kernels us
+                    ## python kernels use platform="dll".
                     #raise Exception("No test cases provided")
                     pass
@@ -198,7 +202,5 @@
                                     'invalid f(%s): %s' % (xi, actual_yi))
                 else:
-                    err = abs(yi - actual_yi)
-                    nrm = abs(yi)
-                    self.assertLess(err * 10**5, nrm,
+                    self.assertTrue(is_near(yi, actual_yi, 5),
                                     'f(%s); expected:%s; actual:%s'
                                     % (xi, yi, actual_yi))
@@ -206,5 +208,11 @@
     return ModelTestCase
 
-
+def is_near(target, actual, digits=5):
+    """
+    Returns true if *actual* is within *digits* significant digits of *target*.
+    """
+    import math
+    shift = 10**math.ceil(math.log10(abs(target)))
+    return abs(target-actual)/shift < 1.5*10**-digits
 
 def main():
@@ -217,4 +225,9 @@
 
     models = sys.argv[1:]
+    if models and models[0] == '-v':
+        verbosity = 2
+        models = models[1:]
+    else:
+        verbosity = 1
     if models and models[0] == 'opencl':
         if not HAVE_OPENCL:
@@ -235,14 +248,18 @@
         print("""\
 usage:
-  python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ...
+  python -m sasmodels.model_test [-v] [opencl|dll] model1 model2 ...
+
+If -v is included on the
+If neither opencl nor dll is specified, then models will be tested with
+both opencl and dll; the compute target is ignored for pure python models.
 
 If model1 is 'all', then all except the remaining models will be tested.
-If no compute target is specified, then models will be tested with both opencl
-and dll; the compute target is ignored for pure python models.""")
+
+""")
 
         return 1
 
     #runner = unittest.TextTestRunner()
-    runner = xmlrunner.XMLTestRunner(output='logs')
+    runner = xmlrunner.XMLTestRunner(output='logs', verbosity=verbosity)
     result = runner.run(make_suite(loaders, models))
     return 1 if result.failures or result.errors else 0
Index: sasmodels/models/HayterMSAsq.py
===================================================================
--- sasmodels/models/HayterMSAsq.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
+++ sasmodels/models/HayterMSAsq.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -56,4 +56,5 @@
         parameters used in P(Q).
 """
+single = False  # double precision only for now
 #             [ "name", "units", default, [lower, upper], "type", "description" ],
 parameters = [["effect_radius", "Ang", 20.75, [0, inf], "volume",
Index: sasmodels/models/bcc.py
===================================================================
--- sasmodels/models/bcc.py	(revision dcdf29d606cfaa2f3bb9d578ac394c6f167fa360)
+++ sasmodels/models/bcc.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -116,4 +116,7 @@
     """
 category = "shape:paracrystal"
+
+single = False
+
 # pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type","description" ],
Index: sasmodels/models/core_shell_ellipsoid.py
===================================================================
--- sasmodels/models/core_shell_ellipsoid.py	(revision 81dd61963c319ce8b1c38d14a6f3056a517fccdd)
+++ sasmodels/models/core_shell_ellipsoid.py	(revision 177c1a138481a9833b4929af687e86583b5c16d5)
@@ -96,4 +96,5 @@
 category = "shape:ellipsoid"
 
+single = False  # TODO: maybe using sph_j1c inside gfn would help?
 # pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type", "description"],
@@ -111,5 +112,5 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/gfn.c", "lib/gauss76.c", "core_shell_ellipsoid.c"]
+source = ["lib/sph_j1c.c", "lib/gfn.c", "lib/gauss76.c", "core_shell_ellipsoid.c"]
 
 demo = dict(scale=1, background=0.001,
Index: sasmodels/models/correlation_length.py
===================================================================
--- sasmodels/models/correlation_length.py	(revision 5054e805a7c53aa18f659a0c1613f527f05caab3)
+++ sasmodels/models/correlation_length.py	(revision 5054e805a7c53aa18f659a0c1613f527f05caab3)
@@ -0,0 +1,92 @@
+#correlation length model
+# Note: model title and parameter table are inserted automatically
+r"""
+Definition
+----------
+
+The scattering intensity I(q) is calculated as
+
+.. math::
+    I(Q) = \frac{A}{Q^n} + \frac{C}{1 + (Q\xi)^m} + B
+
+The first term describes Porod scattering from clusters (exponent = n) and the
+second term is a Lorentzian function describing scattering from polymer chains
+(exponent = m). This second term characterizes the polymer/solvent interactions
+and therefore the thermodynamics. The two multiplicative factors A and C, the
+incoherent background B and the two exponents n and m are used as fitting
+parameters. The final parameter &xi; is a correlation length for the polymer
+chains. Note that when m=2 this functional form becomes the familiar Lorentzian
+function.
+
+For 2D data: The 2D scattering intensity is calculated in the same way as 1D,
+where the q vector is defined as
+
+.. math::
+    q = \sqrt{q_x^2 + q_y^2}
+
+.. figure:: img/correlation_length_1d.jpg
+
+    1D plot using the default values (w/500 data points).
+
+REFERENCE
+B Hammouda, D L Ho and S R Kline, Insight into Clustering in
+Poly(ethylene oxide) Solutions, Macromolecules, 37 (2004) 6932-6937
+"""
+
+from numpy import inf, sqrt
+
+name = "correlation_length"
+title = """Calculates an empirical functional form for SAS data characterized
+by a low-Q signal and a high-Q signal."""
+description = """
+"""
+category = "shape-independent"
+# pylint: disable=bad-continuation, line-too-long
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [
+              ["lorentz_scale", "", 10.0, [0, inf], "", "Lorentzian Scaling Factor"],
+              ["porod_scale", "", 1e-06, [0, inf], "", "Porod Scaling Factor"],
+              ["cor_length", "Ang", 50.0, [0, inf], "", "Correlation length"],
+              ["exponent_p", "", 3.0, [0, inf], "", "Porod Exponent"],
+              ["exponent_l", "1/Ang^2", 2.0, [0, inf], "", "Lorentzian Exponent"],
+             ]
+# pylint: enable=bad-continuation, line-too-long
+
+def Iq(q, lorentz_scale, porod_scale, cor_length, exponent_p, exponent_l):
+    """
+    1D calculation of the Correlation length model
+    """
+    porod = porod_scale / pow(q, exponent_p)
+    lorentz = lorentz_scale / (1.0 + pow(q * cor_length, exponent_l))
+    inten = porod + lorentz
+    return inten
+
+def Iqxy(qx, qy, lorentz_scale, porod_scale, cor_length, exponent_p, exponent_l):
+    """
+    2D calculation of the Correlation length model
+    There is no orientation contribution.
+    """
+    q = sqrt(qx ** 2 + qy ** 2)
+    return Iq(q, lorentz_scale, porod_scale, cor_length, exponent_p, exponent_l)
+
+# parameters for demo
+demo = dict(lorentz_scale=10.0, porod_scale=1.0e-06, cor_length=50.0,
+            exponent_p=3.0, exponent_l=2.0, background=0.1,
+           )
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname = 'CorrLengthModel'
+# pylint: disable=bad-continuation
+oldpars = dict(
+               lorentz_scale='scale_l', porod_scale='scale_p',
+               cor_length='length_l', exponent_p='exponent_p',
+               exponent_l='exponent_l'
+               )
+
+tests = [
+         [{}, 0.001, 1009.98],
+         [{}, 0.150141, 0.174645],
+         [{}, 0.442528, 0.0203957]
+         ]
+# pylint: enable=bad-continuation
Index: sasmodels/models/fcc.c
===================================================================
--- sasmodels/models/fcc.c	(revision 82d239a5cb4d8005a4b04bf42b1425792a28383d)
+++ sasmodels/models/fcc.c	(revision eeb8bacf06af69341f80d5e76d7b10a7239c4e67)
@@ -101,5 +101,5 @@
 
   double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
-  double q_z;
+  // double q_z;
   double cos_val_b3, cos_val_b2, cos_val_b1;
   double a1_dot_q, a2_dot_q,a3_dot_q;
@@ -124,5 +124,5 @@
   const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
   // q vector
-  q_z = 0.0; // for SANS; assuming qz is negligible
+  // q_z = 0.0; // for SANS; assuming qz is negligible
   /// Angles here are respect to detector coordinate
   ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
Index: sasmodels/models/fcc.py
===================================================================
--- sasmodels/models/fcc.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/fcc.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -112,4 +112,6 @@
 category = "shape:paracrystal"
 
+single = False
+
 #             ["name", "units", default, [lower, upper], "type","description"],
 parameters = [["dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"],
Index: sasmodels/models/flexible_cylinder.py
===================================================================
--- sasmodels/models/flexible_cylinder.py	(revision 168052c5791506d7dfad00173f5636fc5701ad11)
+++ sasmodels/models/flexible_cylinder.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -78,4 +78,5 @@
 
 category = "shape:cylinder"
+single = False
 
 # pylint: disable=bad-whitespace, line-too-long
Index: sasmodels/models/flexible_cylinder_ex.py
===================================================================
--- sasmodels/models/flexible_cylinder_ex.py	(revision 504abeee8bc5d67a14d7598227542961d3d748ed)
+++ sasmodels/models/flexible_cylinder_ex.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -98,4 +98,5 @@
         during model fitting.
         """
+single = False
 
 category = "shape:cylinder"
Index: sasmodels/models/gaussian_peak.py
===================================================================
--- sasmodels/models/gaussian_peak.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/gaussian_peak.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -42,4 +42,5 @@
 category = "shape-independent"
 
+single = False
 #             ["name", "units", default, [lower, upper], "type","description"],
 parameters = [["q0", "1/Ang", 0.05, [-inf, inf], "", "Peak position"],
Index: sasmodels/models/hardsphere.py
===================================================================
--- sasmodels/models/hardsphere.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
+++ sasmodels/models/hardsphere.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -55,4 +55,5 @@
                "volume fraction of hard spheres"],
              ]
+single = False
 
 # No volume normalization despite having a volume parameter
Index: sasmodels/models/lamellarCaille.py
===================================================================
--- sasmodels/models/lamellarCaille.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
+++ sasmodels/models/lamellarCaille.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -87,4 +87,6 @@
 category = "shape:lamellae"
 
+single = False
+
 #             ["name", "units", default, [lower, upper], "type","description"],
 parameters = [["thickness", "Ang",  30.0, [0, inf], "volume", "sheet thickness"],
Index: sasmodels/models/lamellarCailleHG.py
===================================================================
--- sasmodels/models/lamellarCailleHG.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
+++ sasmodels/models/lamellarCailleHG.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -91,4 +91,5 @@
 category = "shape:lamellae"
 
+single = False
 parameters = [
     #   [ "name", "units", default, [lower, upper], "type",
Index: sasmodels/models/lamellarPC.py
===================================================================
--- sasmodels/models/lamellarPC.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
+++ sasmodels/models/lamellarPC.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -111,4 +111,6 @@
 category = "shape:lamellae"
 
+single = False
+
 #             ["name", "units", default, [lower, upper], "type","description"],
 parameters = [["thickness", "Ang", 33.0, [0, inf], "volume",
Index: sasmodels/models/lib/gfn.c
===================================================================
--- sasmodels/models/lib/gfn.c	(revision 81dd61963c319ce8b1c38d14a6f3056a517fccdd)
+++ sasmodels/models/lib/gfn.c	(revision 177c1a138481a9833b4929af687e86583b5c16d5)
@@ -7,34 +7,30 @@
 // function gfn4 for oblate ellipsoids
 double
+gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq);
+double
 gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
 {
-	// local variables
-	double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi;
+    // local variables
+    const double pi43=4.0/3.0*M_PI;
+    const double aa = crmaj;
+    const double bb = crmin;
+    const double u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx));
+    const double uq = sqrt(u2)*qq;
+    // changing to more accurate sph_j1c since the following inexplicably fails on Radeon Nano.
+    //const double siq = (uq == 0.0 ? 1.0 : 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq);
+    const double siq = sph_j1c(uq);
+    const double vc = pi43*aa*aa*bb;
+    const double gfnc = siq*vc*delpc;
 
-	Pi = 4.0*atan(1.0);
-	pi43=4.0/3.0*Pi;
-  	aa = crmaj;
- 	bb = crmin;
- 	u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx));
- 	ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx));
-   	uq = sqrt(u2)*qq;
- 	ut= sqrt(ut2)*qq;
-	vc = pi43*aa*aa*bb;
-   	vt = pi43*trmaj*trmaj*trmin;
-   	if (uq == 0.0){
-   		siq = 1.0/3.0;
-   	}else{
-   		siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
-   	}
-   	if (ut == 0.0){
-   		sit = 1.0/3.0;
-   	}else{
-   		sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
-   	}
-   	gfnc = 3.0*siq*vc*delpc;
-  	gfnt = 3.0*sit*vt*delps;
-  	tgfn = gfnc+gfnt;
-  	gfn4 = tgfn*tgfn;
+    const double ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx));
+    const double ut= sqrt(ut2)*qq;
+    const double vt = pi43*trmaj*trmaj*trmin;
+    //const double sit = (ut == 0.0 ? 1.0 : 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut);
+    const double sit = sph_j1c(ut);
+    const double gfnt = sit*vt*delps;
 
-  	return (gfn4);
+    const double tgfn = gfnc + gfnt;
+    const double result = tgfn*tgfn;
+
+    return (result);
 }
Index: sasmodels/models/lib/wrc_cyl.c
===================================================================
--- sasmodels/models/lib/wrc_cyl.c	(revision 504abeee8bc5d67a14d7598227542961d3d748ed)
+++ sasmodels/models/lib/wrc_cyl.c	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -379,4 +379,5 @@
 }
 
+double Sk_WR(double q, double L, double b);
 double Sk_WR(double q, double L, double b)
 {
Index: sasmodels/models/pearl_necklace.py
===================================================================
--- sasmodels/models/pearl_necklace.py	(revision f12357fd54eaa1b861fbe7a32191f62f1f1d5afd)
+++ sasmodels/models/pearl_necklace.py	(revision d18582e2d66b12afb03e1c3a7c0af6e4449c6bfe)
@@ -96,6 +96,5 @@
 
 source = ["lib/Si.c", "pearl_necklace.c"]
-# new flag to let the compiler know to never use single precision
-single = False
+single = False  # use double precision unless told otherwise
 
 def volume(radius, edge_separation, string_thickness, number_of_pearls):
Index: sasmodels/models/rpa.c
===================================================================
--- sasmodels/models/rpa.c	(revision 82c299fa51ab091ee8ead0a235d4177ccfeecb99)
+++ sasmodels/models/rpa.c	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -205,5 +205,5 @@
   const double Kbb = 0.0;
   const double Kcc = 0.0;
-  const double Kdd = 0.0;
+  //const double Kdd = 0.0;
   const double Zaa = Kaa - Kad - Kad;
   const double Zab = Kab - Kad - Kbd;
@@ -278,9 +278,9 @@
   const double Q12 = (-Mab*Mcc + Mac*Mcb)/DenQ;
   const double Q13 = ( Mab*Mbc - Mac*Mbb)/DenQ;
-  const double Q21 = (-Mba*Mcc + Mbc*Mca)/DenQ;
+  //const double Q21 = (-Mba*Mcc + Mbc*Mca)/DenQ;
   const double Q22 = ( Maa*Mcc - Mac*Mca)/DenQ;
   const double Q23 = (-Maa*Mbc + Mac*Mba)/DenQ;
-  const double Q31 = ( Mba*Mcb - Mbb*Mca)/DenQ;
-  const double Q32 = (-Maa*Mcb + Mab*Mca)/DenQ;
+  //const double Q31 = ( Mba*Mcb - Mbb*Mca)/DenQ;
+  //const double Q32 = (-Maa*Mcb + Mab*Mca)/DenQ;
   const double Q33 = ( Maa*Mbb - Mab*Mba)/DenQ;
 
Index: sasmodels/models/star_polymer.py
===================================================================
--- sasmodels/models/star_polymer.py	(revision 168052c5791506d7dfad00173f5636fc5701ad11)
+++ sasmodels/models/star_polymer.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -55,4 +55,5 @@
         """
 category = "shape-independent"
+single = False
 # pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type","description"],
Index: sasmodels/models/stickyhardsphere.py
===================================================================
--- sasmodels/models/stickyhardsphere.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
+++ sasmodels/models/stickyhardsphere.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
@@ -85,4 +85,5 @@
 category = "structure-factor"
 
+single = False
 #             ["name", "units", default, [lower, upper], "type","description"],
 parameters = [
