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: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision a503bfddd742ec35fe2b801f1d75a0a5f11eccfa)
+++ sasmodels/generate.py	(revision 3271e207340eb2fd5f01b1d8ba963ab4fe9ede8c)
@@ -722,4 +722,6 @@
         parameters = COMMON_PARAMETERS + kernel_module.parameters,
         source = getattr(kernel_module, 'source', []),
+        oldname = kernel_module.oldname,
+        oldpars = kernel_module.oldpars,
         )
     # Fill in attributes which default to None
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/models/fcc.c
===================================================================
--- sasmodels/models/fcc.c	(revision 3271e207340eb2fd5f01b1d8ba963ab4fe9ede8c)
+++ sasmodels/models/fcc.c	(revision 3271e207340eb2fd5f01b1d8ba963ab4fe9ede8c)
@@ -0,0 +1,175 @@
+double form_volume(double radius);
+double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld);
+double Iqxy(double qx, double qy, double dnn,
+    double d_factor, double radius,double sld, double solvent_sld,
+    double theta, double phi, double psi);
+
+double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
+double _BCCeval(double Theta, double Phi, double temp1, double temp3);
+double _sphereform(double q, double radius, double sld, double solvent_sld);
+
+
+double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
+
+	const double Da = d_factor*dnn;
+	const double temp1 = q*q*Da*Da;
+	const double temp3 = q*dnn;
+
+	double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
+	return(retVal);
+}
+
+double _BCCeval(double Theta, double Phi, double temp1, double temp3) {
+
+	double temp6,temp7,temp8,temp9,temp10;
+	double result;
+
+	temp6 = sin(Theta);
+	temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta);
+	temp8 = -1.0*sin(Theta)*cos(Phi)-sin(Theta)*sin(Phi)+cos(Theta);
+	temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)-cos(Theta);
+	temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
+	result = pow(1.0-(temp10*temp10),3)*temp6/((1.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10)));
+
+	return (result);
+}
+
+double _sphereform(double q, double radius, double sld, double solvent_sld){
+    const double qr = q*radius;
+    double sn, cn;
+    SINCOS(qr, sn, cn);
+    const double bes = (qr == 0.0 ? 1.0 : 3.0*(sn-qr*cn)/(qr*qr*qr));
+    const double fq = bes * (sld - solvent_sld)*form_volume(radius);
+    return 1.0e-4*fq*fq;
+}
+
+double form_volume(double radius){
+    return 1.333333333333333*M_PI*radius*radius*radius;
+}
+
+
+double Iq(double q, double dnn,
+  double d_factor, double radius,
+  double sld, double solvent_sld){
+
+	//Volume fraction calculated from lattice symmetry and sphere radius
+	const double s1 = dnn/sqrt(0.75);
+	const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
+
+    const double va = 0.0;
+    const double vb = 2.0*M_PI;
+    const double vaj = 0.0;
+    const double vbj = M_PI;
+
+    double summ = 0.0;
+    double answer = 0.0;
+	for(int i=0; i<150; i++) {
+		//setup inner integral over the ellipsoidal cross-section
+		double summj=0.0;
+		const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;		//the outer dummy is phi
+		for(int j=0;j<150;j++) {
+			//20 gauss points for the inner integral
+			double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;		//the inner dummy is theta
+			double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi);
+			summj += yyy;
+		}
+		//now calculate the value of the inner integral
+		double answer = (vbj-vaj)/2.0*summj;
+
+		//now calculate outer integral
+		summ = summ+(Gauss150Wt[i] * answer);
+	}		//final scaling is done at the end of the function, after the NT_FP64 case
+
+	answer = (vb-va)/2.0*summ;
+	answer = answer*_sphereform(q,radius,sld,solvent_sld)*latticescale;
+
+    return answer;
+
+
+}
+
+
+double Iqxy(double qx, double qy, double dnn,
+    double d_factor, double radius,double sld, double solvent_sld,
+    double theta, double phi, double psi){
+
+  double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
+  double q_z;
+  double cos_val_b3, cos_val_b2, cos_val_b1;
+  double a1_dot_q, a2_dot_q,a3_dot_q;
+  double answer;
+  double Zq, Fkq, Fkq_2;
+
+  //convert to q and make scaled values
+  double q = sqrt(qx*qx+qy*qy);
+  double q_x = qx/q;
+  double q_y = qy/q;
+
+  //convert angle degree to radian
+  theta = theta * M_PI_180;
+  phi = phi * M_PI_180;
+  psi = psi * M_PI_180;
+
+  const double Da = d_factor*dnn;
+  const double qDa_2 = q*q*Da*Da;
+  const double s1 = dnn/sqrt(0.75);
+
+
+  //the occupied volume of the lattice
+  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
+  /// Angles here are respect to detector coordinate
+  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
+    // b3 axis orientation
+    b3_x = cos(theta) * cos(phi);
+    b3_y = sin(theta);
+    //b3_z = -cos(theta) * sin(phi);
+    cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
+
+    //alpha = acos(cos_val_b3);
+    // b1 axis orientation
+    b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+    b1_y = sin(psi)*cos(theta);
+    cos_val_b1 = b1_x*q_x + b1_y*q_y;
+    // b2 axis orientation
+    b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+  	b2_y = cos(theta)*cos(psi);
+    cos_val_b2 = b2_x*q_x + b2_y*q_y;
+
+    // The following test should always pass
+    if (fabs(cos_val_b3)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b3 = 1.0;
+    }
+    if (fabs(cos_val_b2)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b2 = 1.0;
+    }
+    if (fabs(cos_val_b1)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b1 = 1.0;
+    }
+    // Compute the angle btw vector q and the a3 axis
+    a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
+
+    // a1 axis
+    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
+
+    // a2 axis
+    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
+
+
+    // Get Fkq and Fkq_2
+    Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
+    Fkq_2 = Fkq*Fkq;
+    // Call Zq=Z1*Z2*Z3
+    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
+
+  // Use SphereForm directly from libigor
+  answer = _sphereform(q,radius,sld,solvent_sld)*Zq*latticescale;
+
+  return answer;
+ }
Index: sasmodels/models/fcc.py
===================================================================
--- sasmodels/models/fcc.py	(revision 3271e207340eb2fd5f01b1d8ba963ab4fe9ede8c)
+++ sasmodels/models/fcc.py	(revision 3271e207340eb2fd5f01b1d8ba963ab4fe9ede8c)
@@ -0,0 +1,123 @@
+#bcc paracrystal model
+#note model title and parameter table are automatically inserted
+#note - calculation requires double precision
+"""
+Calculates the scattering from a **body-centered cubic lattice** with paracrystalline distortion. Thermal vibrations
+are considered to be negligible, and the size of the paracrystal is infinitely large. Paracrystalline distortion is
+assumed to be isotropic and characterized by a Gaussian distribution.
+
+The returned value is scaled to units of |cm^-1|\ |sr^-1|, absolute scale.
+
+Definition
+----------
+
+The scattering intensity *I(q)* is calculated as
+
+.. image:: img/image167.jpg
+
+where *scale* is the volume fraction of spheres, *Vp* is the volume of the primary particle, *V(lattice)* is a volume
+correction for the crystal structure, *P(q)* is the form factor of the sphere (normalized), and *Z(q)* is the
+paracrystalline structure factor for a body-centered cubic structure.
+
+Equation (1) of the 1990 reference is used to calculate *Z(q)*, using equations (29)-(31) from the 1987 paper for
+*Z1*\ , *Z2*\ , and *Z3*\ .
+
+The lattice correction (the occupied volume of the lattice) for a body-centered cubic structure of particles of radius
+*R* and nearest neighbor separation *D* is
+
+.. image:: img/image159.jpg
+
+The distortion factor (one standard deviation) of the paracrystal is included in the calculation of *Z(q)*
+
+.. image:: img/image160.jpg
+
+where *g* is a fractional distortion based on the nearest neighbor distance.
+
+The body-centered cubic lattice is
+
+.. image:: img/image168.jpg
+
+For a crystal, diffraction peaks appear at reduced q-values given by
+
+.. image:: img/image162.jpg
+
+where for a body-centered cubic lattice, only reflections where (\ *h* + *k* + *l*\ ) = even are allowed and
+reflections where (\ *h* + *k* + *l*\ ) = odd are forbidden. Thus the peak positions correspond to (just the first 5)
+
+.. image:: img/image169.jpg
+
+**NB: The calculation of** *Z(q)* **is a double numerical integral that must be carried out with a high density of**
+**points to properly capture the sharp peaks of the paracrystalline scattering.** So be warned that the calculation is
+SLOW. Go get some coffee. Fitting of any experimental data must be resolution smeared for any meaningful fit. This
+makes a triple integral. Very, very slow. Go get lunch!
+
+This example dataset is produced using 200 data points, *qmin* = 0.001 |Ang^-1|, *qmax* = 0.1 |Ang^-1| and the above
+default values.
+
+.. image:: img/image170.jpg
+
+*Figure. 1D plot in the linear scale using the default values (w/200 data point).*
+
+The 2D (Anisotropic model) is based on the reference below where *I(q)* is approximated for 1d scattering. Thus the
+scattering pattern for 2D may not be accurate. Note that we are not responsible for any incorrectness of the 2D model
+computation.
+
+.. image:: img/image165.gif
+
+.. image:: img/image171.jpg
+
+*Figure. 2D plot using the default values (w/200X200 pixels).*
+
+REFERENCE
+
+Hideki Matsuoka et. al. *Physical Review B*, 36 (1987) 1754-1765
+(Original Paper)
+
+Hideki Matsuoka et. al. *Physical Review B*, 41 (1990) 3854 -3856
+(Corrections to FCC and BCC lattice structure calculation)
+"""
+
+from numpy import pi, inf
+
+name = "BCCparacrystal"
+title = "Body-centred cubic lattic with paracrystalline distortion"
+description = """
+    Calculates the scattering from a **body-centered cubic lattice** with paracrystalline distortion. Thermal vibrations
+    are considered to be negligible, and the size of the paracrystal is infinitely large. Paracrystalline distortion is
+    assumed to be isotropic and characterized by a Gaussian distribution.
+    """
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type","description" ],
+    [ "dnn", "Ang", 220, [-inf,inf],"","Nearest neighbour distance"],
+    [ "d_factor", "", 0.06,[-inf,inf],"","Paracrystal distortion factor" ],
+    [ "radius", "Ang",  40, [0, inf], "volume","Particle radius" ],
+    [ "sld", "1e-6/Ang^2", 4, [-inf,inf], "", "Particle scattering length density" ],
+    [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "","Solvent scattering length density" ],
+    [ "theta", "degrees", 60, [-inf, inf], "orientation","In plane angle" ],
+    [ "phi", "degrees", 60, [-inf, inf], "orientation","Out of plane angle" ],
+    [ "psi", "degrees", 60, [-inf,inf], "orientation","Out of plane angle"]
+    ]
+
+source = [ "lib/J1.c", "lib/gauss150.c", "bcc.c" ]
+
+def ER(radius, length):
+    return 0
+
+# parameters for demo
+demo = dict(
+    scale=1, background=0,
+    dnn=220, d_factor=0.06, sld=4, solvent_sld=1,
+    radius=40,
+    theta=60, phi=60, psi=60,
+    radius_pd=.2, radius_pd_n=0.2,
+    theta_pd=15, theta_pd_n=0,
+    phi_pd=15, phi_pd_n=0,
+    psi_pd=15, psi_pd_n=0,
+    )
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname='BCCrystalModel'
+oldpars=dict(sld='sldSph',
+             solvent_sld='sldSolv')
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision 250fa25ad675ea9a85cddbc91dd32cd7a1255fb0)
+++ sasmodels/sasview_model.py	(revision 0a82216fe41763d3d9ab2277f58228a0921edcce)
@@ -32,9 +32,12 @@
 
 
-def make_class(kernel_module, dtype='single'):
+def make_class(kernel_module, dtype='single', namestyle='name'):
     """
     Load the sasview model defined in *kernel_module*.
 
     Returns a class that can be used directly as a sasview model.
+
+    Defaults to using the new name for a model. Setting namestyle='name'
+    will produce a class with a name compatible with SasView
     """
     model =  load_model(kernel_module, dtype=dtype)
@@ -42,5 +45,5 @@
         SasviewModel.__init__(self, model)
     attrs = dict(__init__=__init__)
-    ConstructedModel = type(model.info['name'],  (SasviewModel,), attrs)
+    ConstructedModel = type(model.info[namestyle],  (SasviewModel,), attrs)
     return ConstructedModel
 
@@ -54,4 +57,5 @@
 
         self.name = model.info['name']
+        self.oldname = model.info['oldname']
         self.description = model.info['description']
         self.category = None
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
