Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 5753e4ebc8ad092cc172a8a013acd92a6a19e86a)
+++ sasmodels/compare.py	(revision e1ace4d12bc0e49ddee5c117d92a23547094baf4)
@@ -130,4 +130,13 @@
     if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']:
         pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius']
+    if name == 'barbell' and pars['bell_radius'] < pars['radius']:
+        pars['radius'],pars['bell_radius'] = pars['bell_radius'],pars['radius']
+
+    # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff)
+    if name == 'guinier':
+        #q_max = 0.2  # mid q maximum
+        q_max = 1.0  # high q maximum
+        rg_max = np.sqrt(90*np.log(10) + 3*np.log(pars['scale']))/q_max
+        pars['rg'] = min(pars['rg'],rg_max)
 
     # These constraints are only needed for comparison to sasview
@@ -224,5 +233,8 @@
     model_definition = core.load_model_definition(name)
 
-    view = 'linear' if '-linear' in opts else 'log' if '-log' in opts else 'q4' if '-q4' in opts else 'log'
+    view = ('linear' if '-linear' in opts
+            else 'log' if '-log' in opts
+            else 'q4' if '-q4' in opts
+            else 'log')
 
     opt_values = dict(split
@@ -230,5 +242,8 @@
                       if len(split) == 2)
     # Sort out data
-    qmax = 10.0 if '-exq' in opts else 1.0 if '-highq' in opts else 0.2 if '-midq' in opts else 0.05
+    qmax = (10.0 if '-exq' in opts
+            else 1.0 if '-highq' in opts
+            else 0.2 if '-midq' in opts
+            else 0.05)
     Nq = int(opt_values.get('-Nq', '128'))
     res = float(opt_values.get('-res', '0'))
@@ -239,5 +254,7 @@
 
     # modelling accuracy is determined by dtype and cutoff
-    dtype = 'double' if '-double' in opts else 'single'
+    dtype = ('longdouble' if '-quad' in opts
+             else 'double' if '-double' in opts
+             else 'single')
     cutoff = float(opt_values.get('-cutoff','1e-5'))
 
@@ -247,7 +264,7 @@
         seed = int(opt_values['-random']) if '-random' in opt_values else None
         pars, seed = randomize_model(pars, seed=seed)
-        constrain_pars(model_definition, pars)
         print "Randomize using -random=%i"%seed
     pars.update(set_pars)  # set value after random to control value
+    constrain_pars(model_definition, pars)
 
     # parameter selection
@@ -362,5 +379,5 @@
 
     -plot*/-noplot plots or suppress the plot of the model
-    -single*/-double uses double precision for comparison
+    -single*/-double/-quad use single/double/quad precision for comparison
     -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0
     -Nq=128 sets the number of Q points in the data set
@@ -386,5 +403,5 @@
 NAME_OPTIONS = set([
     'plot','noplot',
-    'single','double',
+    'single','double','longdouble',
     'lowq','midq','highq','exq',
     '2d','1d',
Index: sasmodels/compare_many.py
===================================================================
--- sasmodels/compare_many.py	(revision cd3dba0c60af7479ce1bae901d065d3d5786c8ea)
+++ sasmodels/compare_many.py	(revision b514adfe3f715e5ad1efba9c3b18ed3f086b37af)
@@ -57,4 +57,5 @@
     num_good = 0
     first = True
+    max_diff = 0
     for k in range(N):
         print >>sys.stderr, name, k
@@ -73,5 +74,6 @@
             columns.extend(stats)
             labels.append('GPU single')
-            good = good and (stats[0] < 1e-14)
+            max_diff = max(max_diff, stats[0])
+            good = good and (stats[0] < 5e-5)
         if 0 and environment().has_double:
             gpu_double_value = trymodel(eval_opencl, dtype='double', cutoff=cutoff)
@@ -79,5 +81,6 @@
             columns.extend(stats)
             labels.append('GPU double')
-            good = good and (stats[0] < 1e-14)
+            max_diff = max(max_diff, stats[0])
+            good = good and (stats[0] < 1e-12)
         if 1:
             cpu_double_value = trymodel(eval_ctypes, dtype='double', cutoff=cutoff)
@@ -85,10 +88,12 @@
             columns.extend(stats)
             labels.append('CPU double')
-            good = good and (stats[0] < 1e-14)
+            max_diff = max(max_diff, stats[0])
+            good = good and (stats[0] < 1e-12)
         if 0:
             stats = get_stats(cpu_double_value, gpu_single_value, index)
             columns.extend(stats)
             labels.append('single/double')
-            good = good and (stats[0] < 1e-14)
+            max_diff = max(max_diff, stats[0])
+            good = good and (stats[0] < 5e-5)
 
         columns += [v for _,v in sorted(pars_i.items())]
@@ -100,5 +105,5 @@
         else:
             print(("%d,"%seed)+','.join("%g"%v for v in columns))
-    print '"%d/%d good"'%(num_good, N)
+    print '"good","%d/%d","max diff",%g'%(num_good, N, max_diff)
 
 
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision cd3dba0c60af7479ce1bae901d065d3d5786c8ea)
+++ sasmodels/generate.py	(revision e1ace4d12bc0e49ddee5c117d92a23547094baf4)
@@ -197,5 +197,5 @@
 # TODO: identify model files which have changed since loading and reload them.
 
-__all__ = ["make", "doc", "sources", "use_single"]
+__all__ = ["make", "doc", "sources", "use_single", "use_long_double"]
 
 import sys
@@ -206,6 +206,11 @@
 C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c')
 
+F32 = np.dtype('float32')
 F64 = np.dtype('float64')
-F32 = np.dtype('float32')
+try:  # CRUFT: older numpy does not support float128
+    F128 = np.dtype('float128')
+except TypeError:
+    F128 = None
+
 
 # Scale and background, which are parameters common to every form factor
@@ -322,4 +327,19 @@
     source = re.sub(r'[^a-zA-Z_](\d*[.]\d+|\d+[.]\d*)([eE][+-]?\d+)?',
                     r'\g<0>f', source)
+    return source
+
+def use_long_double(source):
+    """
+    Convert code from double precision to long double precision.
+    """
+    # Convert double keyword to float.  Accept an 'n' parameter for vector
+    # values, where n is 2, 4, 8 or 16. Assume complex numbers are represented
+    # as cdouble which is typedef'd to double2.
+    source = re.sub(r'(^|[^a-zA-Z0-9_]c?)double(([248]|16)?($|[^a-zA-Z0-9_]))',
+                    r'\1long double\2', source)
+    # Convert floating point constants to single by adding 'f' to the end.
+    # OS/X driver complains if you don't do this.
+    source = re.sub(r'[^a-zA-Z_](\d*[.]\d+|\d+[.]\d*)([eE][+-]?\d+)?',
+                    r'\g<0>L', source)
     return source
 
Index: sasmodels/kernel_template.c
===================================================================
--- sasmodels/kernel_template.c	(revision 6ee9d39ba499dc3891914588e16e77849c82a434)
+++ sasmodels/kernel_template.c	(revision 062c56df8360e174e5c222181238d564ad661ef2)
@@ -125,5 +125,5 @@
     if (weight > cutoff) {
       const double scattering = Iq(qi, IQ_PARAMETERS);
-      if (scattering >= 0.0) { // scattering cannot be negative
+      //if (scattering >= 0.0) { // scattering cannot be negative
         ret += weight*scattering;
         norm += weight;
@@ -133,5 +133,5 @@
         norm_vol += vol_weight;
       #endif
-      }
+      //}
     //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,scattering); }
     }
@@ -198,5 +198,5 @@
 
       const double scattering = Iqxy(qxi, qyi, IQXY_PARAMETERS);
-      if (scattering >= 0.0) { // scattering cannot be negative
+      //if (scattering >= 0.0) { // scattering cannot be negative
         // TODO: use correct angle for spherical correction
         // Definition of theta and phi are probably reversed relative to the
@@ -219,5 +219,5 @@
       #endif
         norm_vol += vol_weight;
-      }
+      //}
       //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,scattering); }
     }
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision 076300996fb3c63a832e9a4b72e732f2caf60bed)
+++ sasmodels/kernelcl.py	(revision 92da2316c75e24044469af05f4f89967a0422190)
@@ -22,4 +22,29 @@
 devices, where it can be combined with other structure factors and form
 factors and have instrumental resolution effects applied.
+
+In order to use OpenCL for your models, you will need OpenCL drivers for
+your machine.  These should be available from your graphics card vendor.
+Intel provides OpenCL drivers for CPUs as well as their integrated HD
+graphics chipsets.  AMD also provides drivers for Intel CPUs, but as of
+this writing the performance is lacking compared to the Intel drivers.
+NVidia combines drivers for CUDA and OpenCL in one package.  The result
+is a bit messy if you have multiple drivers installed.  You can see which
+drivers are available by starting python and running:
+
+    import pyopencl as cl
+    cl.create_some_context(interactive=True)
+
+Once you have done that, it will show the available drivers which you
+can select.  It will then tell you that you can use these drivers
+automatically by setting the PYOPENCL_CTX environment variable.
+
+Some graphics cards have multiple devices on the same card.  You cannot
+yet use both of them concurrently to evaluate models, but you can run
+the program twice using a different device for each session.
+
+OpenCL kernels are compiled when needed by the device driver.  Some
+drivers produce compiler output even when there is no error.  You
+can see the output by setting PYOPENCL_COMPILER_OUTPUT=1.  It should be
+harmless, albeit annoying.
 """
 import os
Index: sasmodels/kerneldll.py
===================================================================
--- sasmodels/kerneldll.py	(revision a30cdd54ce2ad3ae2e70ee547fed2a33d28a6841)
+++ sasmodels/kerneldll.py	(revision 5edfe12d884117a412d8437e870506cdf8c27006)
@@ -1,3 +1,3 @@
-"""
+r"""
 C types wrapper for sasview models.
 
@@ -5,4 +5,42 @@
 you wish to allow single precision floating point evaluation for the compiled
 models, otherwise it defaults to *False*.
+
+The compiler command line is stored in the attribute *COMPILE*, with string
+substitutions for %(source)s and %(output)s indicating what to compile and
+where to store it.  The actual command is system dependent.
+
+On windows systems, you have a choice of compilers.  *MinGW* is the GNU
+compiler toolchain, available in packages such as anaconda and PythonXY,
+or available stand alone. This toolchain has had difficulties on some
+systems, and may or may not work for you.  In order to build DLLs, *gcc*
+must be on your path.  If the environment variable *SAS_OPENMP* is given
+then -fopenmp is added to the compiler flags.  This requires a version
+of MinGW compiled with OpenMP support.
+
+An alternative toolchain uses the Microsoft Visual C++ compiler, available
+free from microsoft:
+
+    `http://www.microsoft.com/en-us/download/details.aspx?id=44266`_
+
+Again, this requires that the compiler is available on your path.  This is
+done by running vcvarsall.bat in a windows terminal.  Install locations are
+system dependent, such as:
+
+    C:\Program Files (x86)\Common Files\Microsoft\Visual C++ for Python\9.0\vcvarsall.bat
+
+or maybe
+
+    C:\Users\yourname\AppData\Local\Programs\Common\Microsoft\Visual C++ for Python\9.0\vcvarsall.bat
+
+And again, the environment variable *SAS_OPENMP* controls whether OpenMP is
+used to compile the C code.  This requires the Microsoft vcomp90.dll library,
+which doesn't seem to be included with the compiler, nor does there appear
+to be a public download location.  There may be one on your machine already
+in a location such as:
+
+    C:\Windows\winsxs\x86_microsoft.vc90.openmp*\vcomp90.dll
+
+If you copy this onto your path, such as the python directory or the install
+directory for this application, then OpenMP should be supported.
 """
 
@@ -11,5 +49,5 @@
 import tempfile
 import ctypes as ct
-from ctypes import c_void_p, c_int, c_double, c_float
+from ctypes import c_void_p, c_int, c_longdouble, c_double, c_float
 
 import numpy as np
@@ -58,4 +96,8 @@
     if np.dtype(dtype) == generate.F32:
         basename += "32"
+    elif np.dtype(dtype) == generate.F64:
+        basename += "64"
+    else:
+        basename += "128"
     return joinpath(DLL_PATH, basename+'.so')
 
@@ -80,6 +122,7 @@
     models are allowed as DLLs.
     """
-    if not ALLOW_SINGLE_PRECISION_DLLS: dtype = "double"   # Force 64-bit dll
     dtype = np.dtype(dtype)
+    if dtype == generate.F32 and not ALLOW_SINGLE_PRECISION_DLLS:
+        dtype = generate.F64  # Force 64-bit dll
 
     if callable(info.get('Iq',None)):
@@ -89,6 +132,9 @@
         source = generate.use_single(source)
         tempfile_prefix = 'sas_'+info['name']+'32_'
+    elif dtype == generate.F64:
+        tempfile_prefix = 'sas_'+info['name']+'64_'
     else:
-        tempfile_prefix = 'sas_'+info['name']+'_'
+        source = generate.use_long_double(source)
+        tempfile_prefix = 'sas_'+info['name']+'128_'
 
     source_files = generate.sources(info) + [info['filename']]
@@ -158,5 +204,7 @@
             raise
 
-        fp = c_float if self.dtype == generate.F32 else c_double
+        fp = (c_float if self.dtype == generate.F32
+              else c_double if self.dtype == generate.F64
+              else c_longdouble)
         pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else []
         pd_args_2d= [c_void_p, fp] + [c_int]*Npd2d if Npd2d else []
@@ -227,5 +275,7 @@
 
     def __call__(self, fixed_pars, pd_pars, cutoff):
-        real = np.float32 if self.q_input.dtype == generate.F32 else np.float64
+        real = (np.float32 if self.q_input.dtype == generate.F32
+                else np.float64 if self.q_input.dtype == generate.F64
+                else np.float128)
 
         nq = c_int(self.q_input.nq)
Index: sasmodels/kernelpy.py
===================================================================
--- sasmodels/kernelpy.py	(revision 750ffa58d4841615a93ccb846e347c34ab7cb9e9)
+++ sasmodels/kernelpy.py	(revision 062c56df8360e174e5c222181238d564ad661ef2)
@@ -211,5 +211,5 @@
         if w > cutoff:
             I = form(*args)
-            positive = (I >= 0.0)
+            #positive = (I >= 0.0)
 
             # Note: can precompute spherical correction if theta_index is not
@@ -219,6 +219,6 @@
                                     if theta_index >= 0 else 1.0)
             #spherical_correction = 1.0
-            ret += w * I * spherical_correction * positive
-            norm += w * positive
+            ret += w * I * spherical_correction #* positive
+            norm += w #* positive
 
             # Volume normalization.
@@ -230,6 +230,6 @@
                 vol_args = [args[index] for index in vol_index]
                 vol_weight = np.prod(weight[vol_weight_index])
-                vol += vol_weight * form_volume(*vol_args) * positive
-                vol_norm += vol_weight * positive
+                vol += vol_weight * form_volume(*vol_args) #* positive
+                vol_norm += vol_weight #* positive
 
     positive = (vol * vol_norm != 0.0)
Index: sasmodels/resolution.py
===================================================================
--- sasmodels/resolution.py	(revision eb588ef084c46a7548c1de175696f106522cf6bb)
+++ sasmodels/resolution.py	(revision 7f7f99fcf117177d5281d8a7c0ef25462f0d7f04)
@@ -1074,4 +1074,4 @@
 
 if __name__ == "__main__":
-    demo()
-    #main()
+    #demo()
+    main()
