Index: doc/genmodel.py
===================================================================
--- doc/genmodel.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ doc/genmodel.py	(revision c094758ad0df0ba9f8a4aa5b788a7c4f1bdf249f)
@@ -10,83 +10,115 @@
 # Convert ../sasmodels/models/name.py to name
 model_name = os.path.basename(sys.argv[1])[:-3]
+model_info = core.load_model_info(model_name)
+model = core.build_model(model_info)
 
 # Load the doc string from the module definition file and store it in rst
-docstr = generate.make_doc(core.load_model_info(model_name))
-    
-# Generate automatically plot of the model and add it to rst documentation
+docstr = generate.make_doc(model_info)
 
-info = core.load_model_info(model_name)
 
 # Calculate 1D curve for default parameters
-pars = dict((p[0], p[2]) for p in info['parameters'])
+pars = dict((p.name, p.default) for p in model_info['parameters'])
 
 # Plotting ranges and options
 opts = {
-        'xscale'    : 'log',
-        'yscale'    : 'log' if not info['structure_factor'] else 'linear',
-        'qmin'      : 0.001,
-        'qmax'      : 1.0,
-        'nq'        : 1000,
-        'nq2d'      : 100,
+    'xscale'    : 'log',
+    'yscale'    : 'log' if not model_info['structure_factor'] else 'linear',
+    'zscale'    : 'log' if not model_info['structure_factor'] else 'linear',
+    'q_min'     : 0.001,
+    'q_max'     : 1.0,
+    'nq'        : 1000,
+    'nq2d'      : 400,
+    'vmin'      : 1e-3,  # floor for the 2D data results
+    'qx_max'    : 0.5,
 }
 
-qmin, qmax, nq = opts['qmin'], opts['qmax'], opts['nq']
-qmin = math.log10(qmin)
-qmax = math.log10(qmax)
-q = np.logspace(qmin, qmax, nq)
-data = empty_data1D(q)
-model = core.load_model(model_name)
-calculator = DirectModel(data, model)
-Iq1D = calculator()
 
-# TO DO: Generation of 2D plots 
-# Problem in sasmodels.direct_model._calc_theory
-# There self._kernel.q_input.nq  gets a value of 0 in the 2D case
-# and returns a 0 numpy array (it does not call the C code)
+def plot_1d(model, opts, ax):
+    q_min, q_max, nq = opts['q_min'], opts['q_max'], opts['nq']
+    q_min = math.log10(q_min)
+    q_max = math.log10(q_max)
+    q = np.logspace(q_min, q_max, nq)
+    data = empty_data1D(q)
+    calculator = DirectModel(data, model)
+    Iq1D = calculator()
 
-# If 2D model, compute 2D image
-#if info['has_2d'] != []:
-#    qmax, nq2d = opts['qmax'], opts['nq2d']
-#    data2d = empty_data2D(np.linspace(-qmax, qmax, nq2d), resolution=0.0) 
-#    #model = core.load_model(model_name)
-#    calculator = DirectModel(data2d, model)
-#    Iq2D = calculator()
+    ax.plot(q, Iq1D, color='blue', lw=2, label=model_info['name'])
+    ax.set_xlabel(r'$Q \/(\AA^{-1})$')
+    ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$')
+    ax.set_xscale(opts['xscale'])
+    ax.set_yscale(opts['yscale'])
+    #ax.legend(loc='best')
+
+def plot_2d(model, opts, ax):
+    qx_max, nq2d = opts['qx_max'], opts['nq2d']
+    q = np.linspace(-qx_max, qx_max, nq2d)
+    data2d = empty_data2D(q, resolution=0.0)
+    calculator = DirectModel(data2d, model)
+    Iq2D = calculator() #background=0)
+    Iq2D = Iq2D.reshape(nq2d, nq2d)
+    if opts['zscale'] == 'log':
+        Iq2D = np.log(np.clip(Iq2D, opts['vmin'], np.inf))
+    h = ax.imshow(Iq2D, interpolation='nearest', aspect=1, origin='upper',
+           extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=ice_cm())
+    # , vmin=vmin, vmax=vmax)
+    ax.set_xlabel(r'$Q_x \/(\AA^{-1})$')
+    ax.set_ylabel(r'$Q_y \/(\AA^{-1})$')
+
+def ice_cm():
+    from matplotlib._cm import _Blues_data
+    from matplotlib import colors
+    from matplotlib import rcParams
+    def from_white(segments):
+        scale = 1.0/segments[0][1]
+        return [(k, v*scale, w*scale) for k, v, w in segments]
+    ice_data = dict((k,from_white(v)) for k,v in _Blues_data.items())
+    ice = colors.LinearSegmentedColormap("ice", ice_data, rcParams['image.lut'])
+    return ice
+
 
 # Generate image (comment IF for 1D/2D for the moment) and generate only 1D
-#if info['has_2d'] == []:
-#    fig = plt.figure()
-#    ax = fig.add_subplot(1,1,1)
-#    ax.plot(q, Iq1D, color='blue', lw=2, label=model_name)
-#    ax.set_xlabel(r'$Q \/(\AA^{-1})$')
-#    ax.set_xscale(opts['xscale'])   
-#    ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$')
-#    ax.set_yscale(opts['yscale'])  
-#    ax.legend()
-#else:
-#    # need figure with 1D + 2D
-#    pass
-fig = plt.figure()
-ax = fig.add_subplot(1,1,1)
-ax.plot(q, Iq1D, color='blue', lw=2, label=info['name'])
-ax.set_xlabel(r'$Q \/(\AA^{-1})$')
-ax.set_xscale(opts['xscale'])   
-ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$')
-ax.set_yscale(opts['yscale'])  
-ax.legend()
- 
+fig_height = 3.0 # in
+fig_left = 0.6 # in
+fig_right = 0.5 # in
+fig_top = 0.6*0.25 # in
+fig_bottom = 0.6*0.75
+if model_info['has_2d']:
+    plot_height = fig_height - (fig_top+fig_bottom)
+    plot_width = plot_height
+    fig_width = 2*(plot_width + fig_left + fig_right)
+    aspect = (fig_width, fig_height)
+    ratio = aspect[0]/aspect[1]
+    ax_left = fig_left/fig_width
+    ax_bottom = fig_bottom/fig_height
+    ax_height = plot_height/fig_height
+    ax_width = ax_height/ratio # square axes
+    fig = plt.figure(figsize=aspect)
+    ax2d = fig.add_axes([0.5+ax_left, ax_bottom, ax_width, ax_height])
+    plot_2d(model, opts, ax2d)
+    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height])
+    #ax.set_aspect('square')
+else:
+    plot_height = fig_height - (fig_top+fig_bottom)
+    plot_width = (1+np.sqrt(5))/2*fig_height
+    fig_width = plot_width + fig_left + fig_right
+    ax_left = fig_left/fig_width
+    ax_bottom = fig_bottom/fig_height
+    ax_width = plot_width/fig_width
+    ax_height = plot_height/fig_height
+    aspect = (fig_width, fig_height)
+    fig = plt.figure(figsize=aspect)
+    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height])
+plot_1d(model, opts, ax1d)
 
 # Save image in model/img
 figname = model_name + '_autogenfig.png'
 filename = os.path.join('model', 'img', figname)
-plt.savefig(filename)
+plt.savefig(filename, bbox_inches='tight')
+#print "figure saved in",filename
 
 # Auto caption for figure
 captionstr = '\n'
-captionstr += '.. figure:: img/' + model_name + '_autogenfig.png\n'
+captionstr += '.. figure:: img/' + model_info['id'] + '_autogenfig.png\n'
 captionstr += '\n'
-#if info['has_2d'] == []:
-#    captionstr += '    1D plot corresponding to the default parameters of the model.\n'
-#else:
-#    captionstr += '    1D and 2D plots corresponding to the default parameters of the model.\n'
 captionstr += '    1D plot corresponding to the default parameters of the model.\n'
 captionstr += '\n'
@@ -102,5 +134,5 @@
 else:
     print '------------------------------------------------------------------'
-    print 'References NOT FOUND for model: ', model_name
+    print 'References NOT FOUND for model: ', model_info['id']
     print '------------------------------------------------------------------'
     docstr = docstr + captionstr
Index: explore/J1c.py
===================================================================
--- explore/J1c.py	(revision 0a6da3c2133a98c9475f95364de124e647838b54)
+++ explore/J1c.py	(revision cbd37a7e9b2cd2085f3e1d42c1ab4a48c4db56a7)
@@ -9,21 +9,24 @@
 
 
-SHOW_DIFF = True  # True if show diff rather than function value
+SHOW_DIFF = True # True if show diff rather than function value
+#SHOW_DIFF = False # True if show diff rather than function value
 LINEAR_X = False  # True if q is linearly spaced instead of log spaced
+#LINEAR_X = True # True if q is linearly spaced instead of log spaced
+FUNCTION = "2*J1(x)/x"
 
-def mp_J1c(vec, bits=500):
+def mp_fn(vec, bits=500):
     """
     Direct calculation using sympy multiprecision library.
     """
     with mp.workprec(bits):
-        return [_mp_J1c(mp.mpf(x)) for x in vec]
+        return [_mp_fn(mp.mpf(x)) for x in vec]
 
-def _mp_J1c(x):
+def _mp_fn(x):
     """
-    Helper funciton for mp_j1c
+    Actual function that gets evaluated.  The caller just vectorizes.
     """
     return mp.mpf(2)*mp.j1(x)/x
 
-def np_J1c(x, dtype):
+def np_fn(x, dtype):
     """
     Direct calculation using scipy.
@@ -33,76 +36,12 @@
     return np.asarray(2, dtype)*J1(x)/x
 
-def cephes_J1c(x, dtype, n):
+def sasmodels_fn(x, dtype, platform='ocl'):
     """
     Calculation using pade approximant.
     """
-    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32
-    x = np.asarray(x, dtype)
-    ans = np.empty_like(x)
-    ax = abs(x)
-
-    # Branch a
-    a_idx = ax < f(8.0)
-    a_xsq = x[a_idx]**2
-    a_coeff1 = list(reversed((72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606)))
-    a_coeff2 = list(reversed((144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0)))
-    a_ans1 = np.polyval(np.asarray(a_coeff1[n:], dtype), a_xsq)
-    a_ans2 = np.polyval(np.asarray(a_coeff2[n:], dtype), a_xsq)
-    ans[a_idx] = f(2.0)*a_ans1/a_ans2
-
-    # Branch b
-    b_idx = ~a_idx
-    b_ax = ax[b_idx]
-    b_x = x[b_idx]
-
-    b_y = f(64.0)/(b_ax**2)
-    b_xx = b_ax - f(2.356194491)
-    b_coeff1 = list(reversed((1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6)))
-    b_coeff2 = list(reversed((0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6)))
-    b_ans1 = np.polyval(np.asarray(b_coeff1[n:], dtype),b_y)
-    b_ans2 = np.polyval(np.asarray(b_coeff2[n:], dtype), b_y)
-    b_sn, b_cn = np.sin(b_xx), np.cos(b_xx)
-    ans[b_idx] = np.sign(b_x)*np.sqrt(f(0.636619772)/b_ax) * (b_cn*b_ans1 - (f(8.0)/b_ax)*b_sn*b_ans2)*f(2.0)/b_x
-
-    return ans
-
-def div_J1c(x, dtype):
-    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32
-    x = np.asarray(x, dtype)
-    return f(2.0)*np.asarray([_J1(xi, f)/xi for xi in x], dtype)
-
-def _J1(x, f):
-    ax = abs(x)
-    if ax < f(8.0):
-        y = x*x
-        ans1 = x*(f(72362614232.0)
-                  + y*(f(-7895059235.0)
-                  + y*(f(242396853.1)
-                  + y*(f(-2972611.439)
-                  + y*(f(15704.48260)
-                  + y*(f(-30.16036606)))))))
-        ans2 = (f(144725228442.0)
-                  + y*(f(2300535178.0)
-                  + y*(f(18583304.74)
-                  + y*(f(99447.43394)
-                  + y*(f(376.9991397)
-                  + y)))))
-        return ans1/ans2
-    else:
-        y = f(64.0)/(ax*ax)
-        xx = ax - f(2.356194491)
-        ans1 = (f(1.0)
-                  + y*(f(0.183105e-2)
-                  + y*(f(-0.3516396496e-4)
-                  + y*(f(0.2457520174e-5)
-                  + y*f(-0.240337019e-6)))))
-        ans2 = (f(0.04687499995)
-                  + y*(f(-0.2002690873e-3)
-                  + y*(f(0.8449199096e-5)
-                  + y*(f(-0.88228987e-6)
-                  + y*f(0.105787412e-6)))))
-        sn, cn = np.sin(xx), np.cos(xx)
-        ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1 - (f(8.0)/ax)*sn*ans2)
-        return -ans if (x < f(0.0)) else ans
+    from sasmodels import core, data, direct_model
+    model = core.load_model('bessel', dtype=dtype)
+    calculator = direct_model.DirectModel(data.empty_data1D(x), model)
+    return calculator(background=0)
 
 def plotdiff(x, target, actual, label):
@@ -113,5 +52,6 @@
     """
     if SHOW_DIFF:
-        err = np.clip(abs((target-actual)/target), 0, 1)
+        err = abs((target-actual)/target)
+        #err = np.clip(err, 0, 1)
         pylab.loglog(x, err, '-', label=label)
     else:
@@ -119,19 +59,16 @@
         pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
 
-def compare(x, precision):
+def compare(x, precision, target):
     r"""
     Compare the different computation methods using the given precision.
     """
-    target = np.asarray(mp_J1c(x), 'double')
-    #plotdiff(x, target, mp_J1c(x, 11), 'mp 11 bits')
-    plotdiff(x, target, np_J1c(x, precision), 'direct '+precision)
-    plotdiff(x, target, cephes_J1c(x, precision, 0), 'cephes '+precision)
-    #plotdiff(x, target, cephes_J1c(x, precision, 1), 'cephes '+precision)
-    #plotdiff(x, target, div_J1c(x, precision), 'cephes 2 J1(x)/x '+precision)
+    #plotdiff(x, target, mp_fn(x, 11), 'mp 11 bits')
+    plotdiff(x, target, np_fn(x, precision), 'numpy '+precision)
+    plotdiff(x, target, sasmodels_fn(x, precision, 0), 'sasmodels '+precision)
     pylab.xlabel("qr (1/Ang)")
     if SHOW_DIFF:
         pylab.ylabel("relative error")
     else:
-        pylab.ylabel("2 J1(x)/x")
+        pylab.ylabel(FUNCTION)
         pylab.semilogx(x, target,  '-', label="true value")
     if LINEAR_X:
@@ -147,11 +84,12 @@
     else:
         qr = np.logspace(-3,5,400)
+    target = np.asarray(mp_fn(qr), 'double')
     pylab.subplot(121)
-    compare(qr, 'single')
+    compare(qr, 'single', target)
     pylab.legend(loc='best')
     pylab.subplot(122)
-    compare(qr, 'double')
+    compare(qr, 'double', target)
     pylab.legend(loc='best')
-    pylab.suptitle('2 J1(x)/x')
+    pylab.suptitle(FUNCTION)
 
 if __name__ == "__main__":
Index: sasmodels/data.py
===================================================================
--- sasmodels/data.py	(revision 84db7a50ae40a38a16b4f9c6430e31c26f8407f4)
+++ sasmodels/data.py	(revision c094758ad0df0ba9f8a4aa5b788a7c4f1bdf249f)
@@ -178,6 +178,6 @@
         self.qy_data, self.dqy_data = y, dy
         self.data, self.err_data = z, dz
-        self.mask = (~np.isnan(z) if z is not None
-                     else np.ones_like(x) if x is not None
+        self.mask = (np.isnan(z) if z is not None
+                     else np.zeros_like(x, dtype='bool') if x is not None
                      else None)
         self.q_data = np.sqrt(x**2 + y**2)
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision 17bbadde9de73711781f748816b06e6956dc8cfa)
+++ sasmodels/kernelcl.py	(revision c094758ad0df0ba9f8a4aa5b788a7c4f1bdf249f)
@@ -367,9 +367,10 @@
         self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors]
         context = env.get_context(self.dtype)
+        self.global_size = [self.q_vectors[0].size]
+        #print("creating inputs of size", self.global_size)
         self.q_buffers = [
             cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)
             for q in self.q_vectors
         ]
-        self.global_size = [self.q_vectors[0].size]
 
     def release(self):
Index: sasmodels/models/bessel.py
===================================================================
--- sasmodels/models/bessel.py	(revision 07142f3db170174ebf3c5f6b28931ae502435967)
+++ sasmodels/models/bessel.py	(revision cbd37a7e9b2cd2085f3e1d42c1ab4a48c4db56a7)
@@ -67,7 +67,8 @@
 #Bessel
 parameters = [
+    ["ignored", "", 0.0, [-inf, inf], "", "no parameterless functions"],
              ]
 
-source = ["lib/polevl.c", "lib/j1d.c"]
+source = ["lib/polevl.c", "lib/j1_cephes.c"]
 
 # No volume normalization despite having a volume parameter
@@ -77,5 +78,5 @@
 
 Iq = """
-    return j1(q);
+    return 2.0*j1(q)/q;
     """
 
Index: sasmodels/models/lib/j0_cephes.c
===================================================================
--- sasmodels/models/lib/j0_cephes.c	(revision bfef528b72ff0998652ebfa1704aa03b745dc7af)
+++ sasmodels/models/lib/j0_cephes.c	(revision 094e320030609bfc152d3be6fc35776ebb90d361)
@@ -44,45 +44,4 @@
  */
 
-/*							y0.c
- *
- *	Bessel function of the second kind, order zero
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, y0();
- *
- * y = y0( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns Bessel function of the second kind, of order
- * zero, of the argument.
- *
- * The domain is divided into the intervals [0, 5] and
- * (5, infinity). In the first interval a rational approximation
- * R(x) is employed to compute
- *   y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.
- * Thus a call to j0() is required.
- *
- * In the second interval, the Hankel asymptotic expansion
- * is employed with two rational functions of degree 6/6
- * and 7/7.
- *
- *
- *
- * ACCURACY:
- *
- *  Absolute error, when y0(x) < 1; else relative error:
- *
- * arithmetic   domain     # trials      peak         rms
- *    DEC       0, 30        9400       7.0e-17     7.9e-18
- *    IEEE      0, 30       30000       1.3e-15     1.6e-16
- *
- */
-
 
 /*
@@ -95,5 +54,4 @@
 
 double j0( double );
-
 double j0(double x) {
 
@@ -291,10 +249,10 @@
 
     q = 1.0/x;
-    w = sqrtf(q);
+    w = sqrt(q);
 
     p = w * polevl( q, MO, 7);
     w = q*q;
     xn = q * polevl( w, PH, 7) - PIO4F;
-    p = p * cosf(xn + xx);
+    p = p * cos(xn + xx);
     return(p);
 #endif
Index: sasmodels/models/lib/j1_cephes.c
===================================================================
--- sasmodels/models/lib/j1_cephes.c	(revision bfef528b72ff0998652ebfa1704aa03b745dc7af)
+++ sasmodels/models/lib/j1_cephes.c	(revision e2af2a953458fc0c3f5a1c0c5754c8455115eac8)
@@ -32,41 +32,4 @@
  *    IEEE      0, 30       30000       2.6e-16     1.1e-16
  *
- *
- */
-/*							y1.c
- *
- *	Bessel function of second kind of order one
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, y1();
- *
- * y = y1( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns Bessel function of the second kind of order one
- * of the argument.
- *
- * The domain is divided into the intervals [0, 8] and
- * (8, infinity). In the first interval a 25 term Chebyshev
- * expansion is used, and a call to j1() is required.
- * In the second, the asymptotic trigonometric representation
- * is employed using two rational functions of degree 5/5.
- *
- *
- *
- * ACCURACY:
- *
- *                      Absolute error:
- * arithmetic   domain      # trials      peak         rms
- *    DEC       0, 30       10000       8.6e-17     1.3e-17
- *    IEEE      0, 30       30000       1.0e-15     1.3e-16
- *
- * (error criterion relative when |y1| > 1).
  *
  */
Index: sasmodels/models/lib/polevl.c
===================================================================
--- sasmodels/models/lib/polevl.c	(revision 3936ad3134459cadbcc9696bac14e0ffc9da869b)
+++ sasmodels/models/lib/polevl.c	(revision e2af2a953458fc0c3f5a1c0c5754c8455115eac8)
@@ -50,4 +50,5 @@
 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 */
+
 double polevl( double x, double coef[8], int N );
 double p1evl( double x, double coef[8], int N );
