Index: sasmodels/__init__.py
===================================================================
--- sasmodels/__init__.py	(revision 19dcb93366741a275baa8ad855ce8a310dca78f8)
+++ sasmodels/__init__.py	(revision 989005365096670e54afec254875979f35bd395a)
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision 5134b2c3f75004ca3a032a5e4ca3c5dcc3fe9156)
+++ sasmodels/bumps_model.py	(revision 989005365096670e54afec254875979f35bd395a)
@@ -437,5 +437,5 @@
     data = load_data('DEC07086.DAT')
     set_beam_stop(data, 0.004)
-    plot_data(data)
+    plot_data(data, data.data)
     import matplotlib.pyplot as plt; plt.show()
 
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 84e01d2e0e3a15d7bf739f532e0ba213dd222062)
+++ sasmodels/core.py	(revision 989005365096670e54afec254875979f35bd395a)
@@ -53,4 +53,19 @@
     return v,w/np.sum(w)
 
+def dispersion_mesh(pars):
+    """
+    Create a mesh grid of dispersion parameters and weights.
+
+    Returns [p1,p2,...],w where pj is a vector of values for parameter j
+    and w is a vector containing the products for weights for each
+    parameter set in the vector.
+    """
+    values, weights = zip(*pars)
+    if len(values) > 1:
+        values = [v.flatten() for v in np.meshgrid(*values)]
+        weights = np.vstack([v.flatten() for v in np.meshgrid(*weights)])
+        weights = np.prod(weights, axis=0)
+    return values, weights
+
 def call_kernel(kernel, pars, cutoff=1e-5):
     fixed_pars = [pars.get(name, kernel.info['defaults'][name])
@@ -59,2 +74,25 @@
     return kernel(fixed_pars, pd_pars, cutoff=cutoff)
 
+def call_ER(kernel, pars):
+    ER = kernel.info.get('ER', None)
+    if ER is None:
+        return 1.0
+    else:
+        vol_pars = [get_weights(kernel, pars, name)
+                    for name in kernel.info['partype']['volume']]
+        values, weights = dispersion_mesh(vol_pars)
+        fv = ER(*values)
+        #print values[0].shape, weights.shape, fv.shape
+        return np.sum(weights*fv) / np.sum(weights)
+
+def call_VR(kernel, pars):
+    VR = kernel.info.get('VR', None)
+    if VR is None:
+        return 1.0
+    else:
+        vol_pars = [get_weights(kernel, pars, name)
+                    for name in kernel.info['partype']['volume']]
+        values, weights = dispersion_mesh(vol_pars)
+        whole,part = VR(*values)
+        return np.sum(weights*part)/np.sum(weights*whole)
+
Index: sasmodels/model_test.py
===================================================================
--- sasmodels/model_test.py	(revision 84e01d2e0e3a15d7bf739f532e0ba213dd222062)
+++ sasmodels/model_test.py	(revision 989005365096670e54afec254875979f35bd395a)
@@ -1,7 +1,44 @@
 # -*- coding: utf-8 -*-
 """
-Created on Tue Feb 17 11:43:56 2015
-
-@author: David
+Run model unit tests.
+
+Usage::
+
+     python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ...
+
+     if model1 is 'all', then all except the remaining models will be tested
+
+Each model is tested using the default parameters at q=0.1, (qx,qy)=(0.1,0.1),
+and the ER and VR are computed.  The return values at these points are not
+considered.  The test is only to verify that the models run to completion,
+and do not produce inf or NaN.
+
+Tests are defined with the *tests* attribute in the model.py file.  *tests*
+is a list of individual tests to run, where each test consists of the
+parameter values for the test, the q-values and the expected results.  For
+the effective radius test, the q-value should be 'ER'.  For the VR test,
+the q-value should be 'VR'.  For 1-D tests, either specify the q value or
+a list of q-values, and the corresponding I(q) value, or list of I(q) values.
+
+That is::
+
+    tests = [
+        [ {parameters}, q, I(q)],
+        [ {parameters}, [q], [I(q)] ],
+        [ {parameters}, [q1, q2, ...], [I(q1), I(q2), ...]],
+
+        [ {parameters}, (qx, qy), I(qx, Iqy)],
+        [ {parameters}, [(qx1, qy1), (qx2, qy2), ...], [I(qx1,qy1), I(qx2,qy2), ...],
+
+        [ {parameters}, 'ER', ER(pars) ],
+        [ {parameters}, 'VR', VR(pars) ],
+        ...
+    ]
+
+Parameters are *key:value* pairs, where key is one of the parameters of the
+model and value is the value to use for the test.  Any parameters not given
+in the parameter list will take on the default parameter value.
+
+Precision defaults to 5 digits (relative).
 """
 
@@ -13,5 +50,5 @@
 from .core import list_models, load_model_definition
 from .core import load_model_cl, load_model_dll
-from .core import make_kernel, call_kernel
+from .core import make_kernel, call_kernel, call_ER, call_VR
 
 def annotate_exception(exc, msg):
@@ -54,5 +91,10 @@
         model_definition = load_model_definition(model_name)
 
-        smoke_tests = [[{},0.1,None],[{},(0.1,0.1),None]]
+        smoke_tests = [
+            [{},0.1,None],
+            [{},(0.1,0.1),None],
+            [{},'ER',None],
+            [{},'VR',None],
+            ]
         tests = smoke_tests + getattr(model_definition, 'tests', [])
         
@@ -102,19 +144,24 @@
                 pars, Q, I = test
 
+                if not isinstance(I, list):
+                    I = [I]
                 if not isinstance(Q, list):
                     Q = [Q]
-                if not isinstance(I, list):
-                    I = [I]
-                    
-                if isinstance(Q[0], tuple):
+
+                self.assertEqual(len(I), len(Q))
+
+                if Q[0] == 'ER':
+                    Iq = [call_ER(kernel, pars)]
+                elif Q[0] == 'VR':
+                    Iq = [call_VR(kernel, pars)]
+                elif isinstance(Q[0], tuple):
                     Qx,Qy = zip(*Q)
                     Q_vectors = [np.array(Qx), np.array(Qy)]
+                    kernel = make_kernel(model, Q_vectors)
+                    Iq = call_kernel(kernel, pars)
                 else:
                     Q_vectors = [np.array(Q)]
-
-                self.assertEqual(len(I), len(Q))
-
-                kernel = make_kernel(model, Q_vectors)
-                Iq = call_kernel(kernel, pars)
+                    kernel = make_kernel(model, Q_vectors)
+                    Iq = call_kernel(kernel, pars)
             
                 self.assertGreater(len(Iq), 0)    
@@ -122,9 +169,11 @@
                 
                 for q, i, iq in zip(Q, I, Iq):
-                    if i is None: continue # smoke test --- make sure it runs
-                    err = abs(i - iq)
-                    nrm = abs(i)
-            
-                    self.assertLess(err * 10**5, nrm, 'q:%s; expected:%s; actual:%s' % (q, i, iq))
+                    if i is None:
+                        # smoke test --- make sure it runs and produces a value
+                        self.assertTrue(np.isfinite(iq), 'q:%s; not finite; actual:%s' % (q, iq))
+                    else:
+                        err = abs(i - iq)
+                        nrm = abs(i)
+                        self.assertLess(err * 10**5, nrm, 'q:%s; expected:%s; actual:%s' % (q, i, iq))
                     
         except Exception,exc: 
Index: sasmodels/models/barbell.py
===================================================================
--- sasmodels/models/barbell.py	(revision a5d0d008778409dd90a86aa3db6037187ed3ff57)
+++ sasmodels/models/barbell.py	(revision 989005365096670e54afec254875979f35bd395a)
@@ -127,7 +127,4 @@
 source = [ "lib/J1.c", "lib/gauss76.c", "barbell.c" ]
 
-def ER(radius, length):
-    return 1.0
-
 # parameters for demo
 demo = dict(
Index: sasmodels/models/bcc.py
===================================================================
--- sasmodels/models/bcc.py	(revision 6272968de11aabcab9d51c7e02bcc3dfbb821018)
+++ sasmodels/models/bcc.py	(revision 989005365096670e54afec254875979f35bd395a)
@@ -128,7 +128,4 @@
 source = [ "lib/J1.c", "lib/gauss150.c", "bcc.c" ]
 
-def ER(radius, length):
-    return 0
-
 # parameters for demo
 demo = dict(
Index: sasmodels/models/cylinder.py
===================================================================
--- sasmodels/models/cylinder.py	(revision a5d0d008778409dd90a86aa3db6037187ed3ff57)
+++ sasmodels/models/cylinder.py	(revision 989005365096670e54afec254875979f35bd395a)
@@ -172,16 +172,4 @@
 
 
-# test values:
-# [
-#   [ {parameters}, q, I(q)],
-#   [ {parameters}, [q], [I(q)] ],
-#   [ {parameters}, [q1, q2, ...], [I(q1), I(q2), ...]],
-
-#   [ {parameters}, (qx, qy), I(qx, Iqy)],
-#   [ {parameters}, [(qx1, qy1), (qx2, qy2), ...], [I(qx1,qy1), I(qx2,qy2), ...],
-#   ...
-# ]
-# Precision defaults to 7 digits (relative) for single, 14 for double
-# May want a limited precision version that checks for 8-n or 15-n digits respectively
 qx,qy = 0.2*np.cos(2.5), 0.2*np.sin(2.5)
 tests = [
Index: sasmodels/models/fcc.py
===================================================================
--- sasmodels/models/fcc.py	(revision a5d0d008778409dd90a86aa3db6037187ed3ff57)
+++ sasmodels/models/fcc.py	(revision 989005365096670e54afec254875979f35bd395a)
@@ -103,7 +103,4 @@
 source = [ "lib/J1.c", "lib/gauss150.c", "fcc.c" ]
 
-def ER(radius, length):
-    return 0
-
 # parameters for demo
 demo = dict(
Index: sasmodels/models/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision a5d0d008778409dd90a86aa3db6037187ed3ff57)
+++ sasmodels/models/parallelepiped.py	(revision 989005365096670e54afec254875979f35bd395a)
@@ -101,5 +101,5 @@
 """
 
-from numpy import pi, inf
+from numpy import pi, inf, sqrt
 
 name = "parallelepiped"
