Index: compare.py
===================================================================
--- compare.py	(revision 87fce00c768607daaeafbbd67ca19eb27459bc2e)
+++ compare.py	(revision 373d1b6113e17053b78f2b86735a454ab02666cf)
@@ -301,4 +301,12 @@
     ]
 
+def get_demo_pars(name):
+    import sasmodels.models
+    __import__('sasmodels.models.'+name)
+    model = getattr(sasmodels.models, name)
+    pars = getattr(model, 'demo', None)
+    if pars is None: pars = dict((p[0],p[2]) for p in model.parameters)
+    return pars
+
 def main():
     opts = [arg for arg in sys.argv[1:] if arg.startswith('-')]
@@ -322,9 +330,5 @@
     # if model does not define demo parameters
     name = args[0]
-    import sasmodels.models
-    __import__('sasmodels.models.'+name)
-    model = getattr(sasmodels.models, name)
-    pars = getattr(model, 'demo', None)
-    if pars is None: pars = dict((p[0],p[2]) for p in model.parameters)
+    pars = get_demo_pars(name)
 
     Nopencl = int(args[1]) if len(args) > 1 else 5
Index: doc/developer/index.rst
===================================================================
--- doc/developer/index.rst	(revision 56fc97aa819cfe686101e6c6b3bac4c2c7cfa01e)
+++ doc/developer/index.rst	(revision 56fc97aa819cfe686101e6c6b3bac4c2c7cfa01e)
@@ -0,0 +1,5 @@
+***************************
+Sasmodels Developers Guide
+***************************
+
+
Index: example/sesansfit.py
===================================================================
--- example/sesansfit.py	(revision e8060779479dfa44414965efc7e91f47d58cc2e8)
+++ example/sesansfit.py	(revision e8060779479dfa44414965efc7e91f47d58cc2e8)
@@ -0,0 +1,44 @@
+import numpy as np
+from bumps.names import *
+
+from sasmodels import bumps_model as sas
+kernel = sas.load_model('sphere', dtype='single')
+
+
+radius = 1000
+if False: # fix when data loader exists
+    from sasmodels.dataloader import load_sesans, load_sans
+    data = load_sesans('mydatfile.pz')
+    sans_data = load_sans('mysansfile.xml')
+else:
+    SElength = np.linspace(0, 2400, 61) # [A]
+    data = np.ones_like(SElength)
+    err_data = np.ones_like(SElength)*0.03
+
+    class SesansData:
+        #q_zmax = 0.23 # [A^-1]
+        q_zmax = 0.1 # [A^-1]
+        SElength = np.linspace(0, 2400, 61) # [A]
+        Rmax = 3*radius # [A]
+        wavelength = 2e-10 # [m]
+        thickness = 0.2 # [cm]
+        SElength = SElength
+        data = data
+        err_data = err_data
+    data = SesansData()
+
+phi = Parameter(0.1, name="phi")
+model = sas.BumpsModel(data, kernel,
+    scale=phi*(1-phi), sld=7.0, solvent_sld=1.0, radius=radius)
+phi.pmp(10)
+model.radius.pmp(40)
+model.sld.pm(2)
+model.background.range(0,5)
+
+
+if False: # have sans data
+    sansmodel = sas.BumpsModel(sans_data, kernel, **model.parameters())
+    problem = FitProblem([model, sansmodel])
+else:
+    problem = FitProblem(model)
+
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision f786ff3daee63fedbada60cec1eb19964c370615)
+++ sasmodels/bumps_model.py	(revision c97724e7592a12c8cc1d1da48b84c88fb8a97093)
@@ -4,4 +4,6 @@
 import sys, os
 import datetime
+
+from sasmodels import sesans
 
 # CRUFT python 2.6
@@ -213,4 +215,16 @@
     #plt.axhline(-1, color='black', ls='--',lw=1, hold=True)
 
+def _plot_sesans(data, theory, view):
+    import matplotlib.pyplot as plt
+    resid = (theory - data.data)/data.err_data
+    plt.subplot(121)
+    plt.errorbar(data.SElength, data.data, yerr=data.err_data)
+    plt.plot(data.SElength, theory, '-', hold=True)
+    plt.xlabel('spin echo length (A)')
+    plt.ylabel('polarization')
+    plt.subplot(122)
+    plt.plot(data.SElength, resid, 'x')
+    plt.xlabel('spin echo length (A)')
+    plt.ylabel('residuals')
 
 def _plot_result2D(data, theory, view):
@@ -230,14 +244,4 @@
     plt.colorbar()
 
-def plot_result(data, theory, view='log'):
-    """
-    Plot the data and residuals.
-    """
-    if hasattr(data, 'qx_data'):
-        _plot_result2D(data, theory, view)
-    else:
-        _plot_result1D(data, theory, view)
-
-
 class BumpsModel(object):
     """
@@ -263,9 +267,22 @@
         self.model = model
         self.cutoff = cutoff
+        if hasattr(data, 'SElength'):
+            self.data_type = 'sesans'
+        elif hasattr(data, 'qx_data'):
+            self.data_type = 'Iqxy'
+        else:
+            self.data_type = 'Iq'
 
         partype = model.info['partype']
 
         # interpret data
-        if hasattr(data, 'qx_data'):
+        if self.data_type == 'sesans':
+            q = sesans.make_q(data.q_zmax, data.Rmax)
+            self.index = slice(None,None)
+            self.iq = data.data
+            self.diq = data.err_data
+            self._theory = np.zeros_like(q)
+            q_vectors = [q]
+        elif self.data_type == 'Iqxy':
             self.index = (data.mask==0) & (~np.isnan(data.data))
             self.iq = data.data[self.index]
@@ -276,5 +293,5 @@
             else:
                 q_vectors = [data.qx_data, data.qy_data]
-        else:
+        elif self.data_type == 'Iq':
             self.index = (data.x>=data.qmin) & (data.x<=data.qmax) & ~np.isnan(data.y)
             self.iq = data.y[self.index]
@@ -282,4 +299,6 @@
             self._theory = np.zeros_like(data.y)
             q_vectors = [data.x]
+        else:
+            raise ValueError("Unknown data type") # never gets here
 
         # Remember function inputs so we can delay loading the function and
@@ -333,5 +352,11 @@
             self._theory[self.index] = self._fn(pars, pd_pars, self.cutoff)
             #self._theory[:] = self._fn.eval(pars, pd_pars)
-            self._cache['theory'] = self._theory
+            if self.data_type == 'sesans':
+                P = sesans.hankel(self.data.SElength, self.data.wavelength,
+                                  self.data.thickness, self._fn_inputs[0],
+                                  self._theory)
+                self._cache['theory'] = P
+            else:
+                self._cache['theory'] = self._theory
         return self._cache['theory']
 
@@ -349,5 +374,34 @@
 
     def plot(self, view='log'):
-        plot_result(self.data, self.theory(), view=view)
+        """
+        Plot the data and residuals.
+        """
+        data, theory = self.data, self.theory()
+        if self.data_type == 'Iq':
+            _plot_result1D(data, theory, view)
+        elif self.data_type == 'Iqxy':
+            _plot_result2D(data, theory, view)
+        elif self.data_type == 'sesans':
+            _plot_sesans(data, theory, view)
+        else:
+            raise ValueError("Unknown data type")
+
+    def simulate_data(self, noise=None):
+        print "noise", noise
+        if noise is None:
+            noise = self.diq[self.index]
+        else:
+            noise = 0.01*noise
+            self.diq[self.index] = noise
+        y = self.theory()
+        y += y*np.random.randn(*y.shape)*noise
+        if self.data_type == 'Iq':
+            self.data.y[self.index] = y
+        elif self.data_type == 'Iqxy':
+            self.data.data[self.index] = y
+        elif self.data_type == 'sesans':
+            self.data.data[self.index] = y
+        else:
+            raise ValueError("Unknown model")
 
     def save(self, basename):
Index: smodels/hankel.py
===================================================================
--- sasmodels/hankel.py	(revision 10576d11d250be03daca38ed2d4cdce37228346c)
+++ 	(revision )
@@ -1,46 +1,0 @@
-"""
-Conversion of scattering cross section from SANS in absolute
-units into SESANS using a Hankel transformation
-
-Everything is in units of metres except specified otherwise
-
-Wim Bouwman (w.g.bouwman@tudelft.nl), June 2013
-"""
-
-from __future__ import division
-
-import numpy as np
-from numpy import pi, exp
-from scipy.special import jv as besselj
-
-def hankel(zz, wavelength, th, q, I):
-    """
-    Compute the expected SESANS polarization for a given SANS pattern.
-
-    Uses the hankel transform followed by the exponential.  The values
-    for zz (or spin echo length, or delta), wavelength and sample thickness
-    information should come from the dataset.  *q* should be chosen such
-    that the oscillations in *I(q)* are well sampled (e.g., 5*2*pi/d_max).
-
-    *zz* [nm] is the set of z points at which to compute the hankel transform
-
-    *wavelength* [m]  is the wavelength of each individual point *zz*
-
-    *th* [m] is the sample thickness.
-
-    *q* [nm^{-1}] is the set of q points at which the model has been computed.
-    These should be equally spaced.
-
-    *qwidth* is the width of the integration
-
-    *I* [m^{-1}] is the value of the SANS model at *q*
-    """
-    dq=(q[1]-q[0])*1e9   # [m^-1] step size in q, needed for integration
-    G = np.zeros(len(zz), 'd')
-    for i in range(len(zz)):
-        integr = besselj(0,q*zz[i])*I*q
-        G[i] = np.sum(integr)
-    G *= dq*1e9*2*pi # integr step, conver q into [m**-1] and 2 pi circle integr
-    PP = exp(th*wavelength**2/(4*pi**2)*(G-G[0]))
-
-    return PP
Index: sasmodels/sesans.py
===================================================================
--- sasmodels/sesans.py	(revision c97724e7592a12c8cc1d1da48b84c88fb8a97093)
+++ sasmodels/sesans.py	(revision c97724e7592a12c8cc1d1da48b84c88fb8a97093)
@@ -0,0 +1,68 @@
+"""
+Conversion of scattering cross section from SANS in absolute
+units into SESANS using a Hankel transformation
+
+Everything is in units of metres except specified otherwise
+
+Wim Bouwman (w.g.bouwman@tudelft.nl), June 2013
+"""
+
+from __future__ import division
+
+import numpy as np
+from numpy import pi, exp
+
+from scipy.special import jv as besselj
+
+def make_q(q_zmax, Rmax):
+    q_min = dq = 0.1 * 2*pi / Rmax
+    #q_min = 0.00003
+    return np.arange(q_min, q_zmax, dq)
+
+# TODO: dead code; for now the call to the hankel transform happens in BumpsModel
+class SesansCalculator:
+    def __init__(self, sans_kernel, q_zmax, Rmax, SElength, wavelength, thickness):
+        self._set_kernel(sans_kernel, q_zmax, Rmax)
+        self.SElength = SElength
+        self.wavelength = wavelength
+        self.thickness = thickness
+
+    def _set_kernel(self, sans_kernel, q_zmax, Rmax):
+        input = sans_kernel.make_input([make_q(q_zmax, Rmax)])
+        self.sans_calculator = sans_kernel(input)
+
+    def __call__(self, pars, pd_pars, cutoff=1e-5):
+        Iq = self.sans_calculator(pars, pd_pars, cutoff)
+        P = hankel(self.SElength, self.wavelength, self.thickness, self.q, Iq)
+        self.Iq = Iq
+        return P
+
+def hankel(SElength, wavelength, thickness, q, Iq):
+    """
+    Compute the expected SESANS polarization for a given SANS pattern.
+
+    Uses the hankel transform followed by the exponential.  The values
+    for zz (or spin echo length, or delta), wavelength and sample thickness
+    information should come from the dataset.  *q* should be chosen such
+    that the oscillations in *I(q)* are well sampled (e.g., 5*2*pi/d_max).
+
+    *SElength* [A] is the set of z points at which to compute the hankel transform
+
+    *wavelength* [m]  is the wavelength of each individual point *zz*
+
+    *thickness* [cm] is the sample thickness.
+
+    *q* [A^{-1}] is the set of q points at which the model has been computed.
+    These should be equally spaced.
+
+    *I* [cm^{-1}] is the value of the SANS model at *q*
+    """
+    G = np.zeros(len(SElength), 'd')
+    for i in range(len(SElength)):
+        integr = besselj(0,q*SElength[i])*Iq*q
+        G[i] = np.sum(integr)
+    dq=(q[1]-q[0])*1e10   # [m^-1] step size in q, needed for integration
+    G *= dq*1e10*2*pi # integr step, conver q into [m**-1] and 2 pi circle integr
+    P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0]))
+
+    return P
Index: sesansdemo.py
===================================================================
--- sesansdemo.py	(revision 10576d11d250be03daca38ed2d4cdce37228346c)
+++ sesansdemo.py	(revision c97724e7592a12c8cc1d1da48b84c88fb8a97093)
@@ -10,4 +10,5 @@
 
 # q-range parameters
+
 q = arange(0.0003, 1.0, 0.0003);    # [nm^-1] range wide enough for  Hankel transform
 dq=(q[1]-q[0])*1e9;   # [m^-1] step size in q, needed for integration
@@ -34,5 +35,5 @@
 
 clf()
-subplot(1,2,1)  # plot the SANS calculation
+subplot(211)  # plot the SANS calculation
 plot(q,I,'k')
 loglog(q,I)
@@ -55,5 +56,5 @@
 PP=exp(th*Lambda**2/4/pi/pi*(G-G[0]));
 
-subplot(1,2,2)
+subplot(212)
 plot(zz,PP,'k',label="Hankel transform") # Hankel transform 1D
 xlabel('spin-echo length [nm]')
@@ -62,5 +63,5 @@
 
 # Cosine transformation of 2D scattering patern
-if True:
+if False:
     qy,qz = meshgrid(q,q)
     qr=R*sqrt(qy**2 + qz**2); # reuse variable names Hankel transform, but now 2D
