Index: compare
===================================================================
--- compare	(revision 5852c866b0b68e7ed626252d76561ef3fe4baace)
+++ compare	(revision 5852c866b0b68e7ed626252d76561ef3fe4baace)
@@ -0,0 +1,20 @@
+#!/usr/bin/env python
+
+import sys
+import os
+import glob
+
+def main():
+    sasmodels = os.path.dirname(os.path.abspath(__file__))
+    root = os.path.dirname(sasmodels)
+    sasview=glob.glob(os.path.join(root, 'sasview', 'build', 'lib.*' ))[0]
+    sys.path.insert(0, os.path.join(root, 'bumps'))
+    sys.path.insert(0, os.path.join(root, 'periodictable'))
+    sys.path.insert(0, sasview)
+    sys.path.insert(0, sasmodels)
+
+    import sasmodels.compare
+    sasmodels.compare.main()
+
+if __name__ == "__main__":
+    main()
Index: doc/genmodel.py
===================================================================
--- doc/genmodel.py	(revision c094758ad0df0ba9f8a4aa5b788a7c4f1bdf249f)
+++ doc/genmodel.py	(revision 044a9c12e31225b960e033a80e969530db1641d4)
@@ -2,4 +2,5 @@
 import numpy as np
 import matplotlib.pyplot as plt
+import pylab
 sys.path.insert(0, os.path.abspath('..'))
 from sasmodels import generate, core
@@ -28,5 +29,5 @@
     'q_max'     : 1.0,
     'nq'        : 1000,
-    'nq2d'      : 400,
+    'nq2d'      : 100,
     'vmin'      : 1e-3,  # floor for the 2D data results
     'qx_max'    : 0.5,
@@ -59,9 +60,15 @@
     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.imshow(Iq2D, interpolation='nearest', aspect=1, origin='lower',
+        extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=pylab.cm.jet)   
     ax.set_xlabel(r'$Q_x \/(\AA^{-1})$')
     ax.set_ylabel(r'$Q_y \/(\AA^{-1})$')
+
+#            im = self.subplot.imshow(output, interpolation='nearest',
+#                                     origin='lower',
+#                                     vmin=zmin_temp, vmax=self.zmax_2D,
+#                                     cmap=self.cmap,
+#                                     extent=(self.xmin_2D, self.xmax_2D,
+#                                             self.ymin_2D, self.ymax_2D))
 
 def ice_cm():
@@ -77,5 +84,5 @@
 
 
-# Generate image (comment IF for 1D/2D for the moment) and generate only 1D
+# Generate image 
 fig_height = 3.0 # in
 fig_left = 0.6 # in
@@ -97,4 +104,5 @@
     plot_2d(model, opts, ax2d)
     ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height])
+    plot_1d(model, opts, ax1d)
     #ax.set_aspect('square')
 else:
@@ -109,5 +117,5 @@
     fig = plt.figure(figsize=aspect)
     ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height])
-plot_1d(model, opts, ax1d)
+    plot_1d(model, opts, ax1d)
 
 # Save image in model/img
@@ -121,5 +129,8 @@
 captionstr += '.. figure:: img/' + model_info['id'] + '_autogenfig.png\n'
 captionstr += '\n'
-captionstr += '    1D plot corresponding to the default parameters of the model.\n'
+if model_info['has_2d']:
+    captionstr += '    1D and 2D plots corresponding to the default parameters of the model.\n'
+else:
+    captionstr += '    1D plot corresponding to the default parameters of the model.\n'
 captionstr += '\n'
 
Index: example/sesansfit.py
===================================================================
--- example/sesansfit.py	(revision 84db7a50ae40a38a16b4f9c6430e31c26f8407f4)
+++ example/sesansfit.py	(revision 170ea69cd9ce6a059d8adf8d340e81b5064042b9)
@@ -1,4 +1,5 @@
 from bumps.names import *
 from sasmodels import core, bumps_model, sesans
+from sas.sascalc.dataloader.loader import Loader
 
 HAS_CONVERTER = True
@@ -8,4 +9,5 @@
     HAS_CONVERTER = False
 
+
 def get_bumps_model(model_name):
     kernel = core.load_model(model_name)
@@ -13,6 +15,7 @@
     return model
 
-def sesans_fit(file, model, initial_vals={}, custom_params={}, param_range=[]):
+def sesans_fit(file, model, initial_vals={}, custom_params={}, param_range=[], acceptance_angle=None):
     """
+
     @param file: SESANS file location
     @param model: Bumps model object or model name - can be model, model_1 * model_2, and/or model_1 + model_2
@@ -24,5 +27,4 @@
     """
     try:
-        from sas.sascalc.dataloader.loader import Loader
         loader = Loader()
         data = loader.load(file)
@@ -55,6 +57,9 @@
             dy = err_data
             sample = Sample()
+            needs_all_q = acceptance_angle is not None
         data = SESANSData1D()
+        data.acceptance_angle = acceptance_angle
 
+    data.needs_all_q = acceptance_angle is not None
     if "radius" in initial_vals:
         radius = initial_vals.get("radius")
Index: example/sesansfit.sh
===================================================================
--- example/sesansfit.sh	(revision bdb653cc842ebbb2483def7e99946c8db6fb3346)
+++ example/sesansfit.sh	(revision 170ea69cd9ce6a059d8adf8d340e81b5064042b9)
@@ -1,6 +1,15 @@
 #!/bin/bash
 
-SASVIEW=$PWD/../../sasview/src
-PYTHONPATH=$PWD/../:$PWD/../../bumps:$PWD/../../periodictable:$SASVIEW
+# Need to fix the paths to sasmodels and sasview if no eggs present
+echo $PWD
+ONEUP="$(dirname "$PWD")"
+PROJECTS="$(dirname "$ONEUP")"
+CCOLON="C:/"
+CSLASH="/c/"
+SASMODELSBASE=$PROJECTS/sasmodels/
+SASMODELS="${SASMODELSBASE/$CSLASH/$CCOLON}"
+SASVIEWBASE=$PROJECTS/sasview/src/
+SASVIEW="${SASVIEWBASE/$CSLASH/$CCOLON}"
+PYTHONPATH="$SASVIEW;$SASMODELS"
 export PYOPENCL_CTX PYTHONPATH
 
Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision 78d3341c5d095cf5c713ac48cb7e918807415d33)
+++ sasmodels/convert.py	(revision 74f7238cb6c48971beea1fac687b23814387b471)
@@ -15,5 +15,5 @@
     'be_polyelectrolyte',
     'correlation_length',
-    'fractal_core_shell'
+    'fractal_core_shell',
     'binary_hard_sphere',
 ]
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 7b3e62cb529ae72b9e299f503716e786281f6f83)
+++ sasmodels/core.py	(revision d5ba8419b0a797faa3fe41fde9334c5a19282992)
@@ -170,11 +170,10 @@
     """
     value, weight = zip(*pars)
-    if len(value) > 1:
-        value = [v.flatten() for v in np.meshgrid(*value)]
-        weight = np.vstack([v.flatten() for v in np.meshgrid(*weight)])
-        weight = np.prod(weight, axis=0)
+    value = [v.flatten() for v in np.meshgrid(*value)]
+    weight = np.vstack([v.flatten() for v in np.meshgrid(*weight)])
+    weight = np.prod(weight, axis=0)
     return value, weight
 
-def call_kernel(kernel, pars, cutoff=0):
+def call_kernel(kernel, pars, cutoff=0, mono=False):
     """
     Call *kernel* returned from :func:`make_kernel` with parameters *pars*.
@@ -189,5 +188,9 @@
     fixed_pars = [pars.get(name, kernel.info['defaults'][name])
                   for name in kernel.fixed_pars]
-    pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars]
+    if mono:
+        pd_pars = [( np.array([pars[name]]), np.array([1.0]) )
+                   for name in kernel.pd_pars]
+    else:
+        pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars]
     return kernel(fixed_pars, pd_pars, cutoff=cutoff)
 
Index: sasmodels/direct_model.py
===================================================================
--- sasmodels/direct_model.py	(revision 17bbadde9de73711781f748816b06e6956dc8cfa)
+++ sasmodels/direct_model.py	(revision 02e70ffbd2f37e9a91d5d027972b30069bd05c6b)
@@ -69,4 +69,5 @@
 
         if self.data_type == 'sesans':
+            
             q = sesans.make_q(data.sample.zacceptance, data.Rmax)
             index = slice(None, None)
@@ -77,5 +78,6 @@
                 Iq, dIq = None, None
             #self._theory = np.zeros_like(q)
-            q_vectors = [q]
+            q_vectors = [q]            
+            q_mono = sesans.make_all_q(data)
         elif self.data_type == 'Iqxy':
             if not partype['orientation'] and not partype['magnetic']:
@@ -96,4 +98,5 @@
             #self._theory = np.zeros_like(self.Iq)
             q_vectors = res.q_calc
+            q_mono = []
         elif self.data_type == 'Iq':
             index = (data.x >= data.qmin) & (data.x <= data.qmax)
@@ -120,4 +123,5 @@
             #self._theory = np.zeros_like(self.Iq)
             q_vectors = [res.q_calc]
+            q_mono = []
         else:
             raise ValueError("Unknown data type") # never gets here
@@ -125,5 +129,6 @@
         # Remember function inputs so we can delay loading the function and
         # so we can save/restore state
-        self._kernel_inputs = [v for v in q_vectors]
+        self._kernel_inputs = q_vectors
+        self._kernel_mono_inputs = q_mono
         self._kernel = None
         self.Iq, self.dIq, self.index = Iq, dIq, index
@@ -149,14 +154,16 @@
     def _calc_theory(self, pars, cutoff=0.0):
         if self._kernel is None:
-            self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-defined-outside-init
+            self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-dedata_type
+            self._kernel_mono = make_kernel(self._model, self._kernel_mono_inputs) if self._kernel_mono_inputs else None
 
         Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff)
+        Iq_mono = call_kernel(self._kernel_mono, pars, mono=True) if self._kernel_mono_inputs else None
         if self.data_type == 'sesans':
-            result = sesans.hankel(self._data.x, self._data.lam * 1e-9,
-                                   self._data.sample.thickness / 10,
-                                   self._kernel_inputs[0], Iq_calc)
+            result = sesans.transform(self._data,
+                                   self._kernel_inputs[0], Iq_calc, 
+                                   self._kernel_mono_inputs, Iq_mono)
         else:
             result = self.resolution.apply(Iq_calc)
-        return result
+        return result        
 
 
Index: sasmodels/models/adsorbed_layer.py
===================================================================
--- sasmodels/models/adsorbed_layer.py	(revision f10bc52ab6508ac7b449d2990c9f805461b8996c)
+++ sasmodels/models/adsorbed_layer.py	(revision c079f5048e0520234276a4b55ef0538f7b83852c)
@@ -6,5 +6,5 @@
 
 r"""
-This model describes the scattering from a layer of surfactant or polymer adsorbed on spherical particles under the conditions that (i) the particles (cores) are contrast-matched to the dispersion medium, (ii) *S(Q)* ~ 1 (ie, the particle volume fraction is dilute), (iii) the particle radius is >> layer thickness (ie, the interface is locally flat), and (iv) scattering from excess unadsorbed adsorbate in the bulk medium is absent or has been corrected for.
+This model describes the scattering from a layer of surfactant or polymer adsorbed on large, smooth, notionally spherical particles under the conditions that (i) the particles (cores) are contrast-matched to the dispersion medium, (ii) *S(Q)* ~ 1 (ie, the particle volume fraction is dilute), (iii) the particle radius is >> layer thickness (ie, the interface is locally flat), and (iv) scattering from excess unadsorbed adsorbate in the bulk medium is absent or has been corrected for.
 
 Unlike many other core-shell models, this model does not assume any form for the density distribution of the adsorbed species normal to the interface (cf, a core-shell model normally assumes the density distribution to be a homogeneous step-function). For comparison, if the thickness of a (traditional core-shell like) step function distribution is *t*, the second moment about the mean of the density distribution (ie, the distance of the centre-of-mass of the distribution from the interface), |sigma| =  sqrt((*t* :sup:`2` )/12).
@@ -34,5 +34,5 @@
 
 description =  """
-    Evaluates the scattering from particles
+    Evaluates the scattering from large particles
     with an adsorbed layer of surfactant or
     polymer, independent of the form of the
@@ -42,22 +42,26 @@
 
 #             ["name", "units", default, [lower, upper], "type", "description"],
-parameters =  [["second_moment", "Ang", 23.0, [0.0, inf], "", "Second moment"],
-              ["adsorbed_amount", "mg/m2", 1.9, [0.0, inf], "", "Adsorbed amount"],
-              ["density_poly", "g/cm3", 0.7, [0.0, inf], "", "Polymer density"],
-              ["radius", "Ang", 500.0, [0.0, inf], "", "Particle radius"],
-              ["volfraction", "None", 0.14, [0.0, inf], "", "Particle vol fraction"],
-              ["polymer_sld", "1/Ang^2", 1.5e-06, [-inf, inf], "", "Polymer SLD"],
-              ["solvent_sld", "1/Ang^2", 6.3e-06, [-inf, inf], "", "Solvent SLD"]]
+parameters =  [["second_moment", "Ang", 23.0, [0.0, inf], "", "Second moment of polymer distribution"],
+              ["adsorbed_amount", "mg/m2", 1.9, [0.0, inf], "", "Adsorbed amount of polymer"],
+              ["density_shell", "g/cm3", 0.7, [0.0, inf], "", "Bulk density of polymer in the shell"],
+              ["radius", "Ang", 500.0, [0.0, inf], "", "Core particle radius"],
+              ["volfraction", "None", 0.14, [0.0, inf], "", "Core particle volume fraction"],
+              ["sld_shell", "1e-6/Ang^2", 1.5, [-inf, inf], "sld", "Polymer shell SLD"],
+              ["sld_solvent", "1e-6/Ang^2", 6.3, [-inf, inf], "sld", "Solvent SLD"]]
 
 # NB: Scale and Background are implicit parameters on every model
-def Iq(q, second_moment, adsorbed_amount, density_poly, radius, 
-        volfraction, polymer_sld, solvent_sld):
+def Iq(q, second_moment, adsorbed_amount, density_shell, radius, 
+        volfraction, sld_shell, sld_solvent):
     # pylint: disable = missing-docstring
-    deltarhosqrd =  (polymer_sld - solvent_sld) * (polymer_sld - solvent_sld)
-    numerator =  6.0 * pi * volfraction * (adsorbed_amount * adsorbed_amount)
-    denominator =  (q * q) * (density_poly * density_poly) * radius
-    eterm =  exp(-1.0 * (q * q) * (second_moment * second_moment))
-    #scale by 10^10 for units conversion to cm^-1
-    inten =  1.0e+10 * deltarhosqrd * ((numerator / denominator) * eterm)
+#    deltarhosqrd =  (sld_shell - sld_solvent) * (sld_shell - sld_solvent)
+#    numerator =  6.0 * pi * volfraction * (adsorbed_amount * adsorbed_amount)
+#    denominator =  (q * q) * (density_shell * density_shell) * radius
+#    eterm =  exp(-1.0 * (q * q) * (second_moment * second_moment))
+#    #scale by 10^-2 for units conversion to cm^-1
+#    inten =  1.0e-02 * deltarhosqrd * ((numerator / denominator) * eterm)
+    aa =  (sld_shell - sld_solvent) * adsorbed_amount / q / density_shell 
+    bb = q * second_moment
+    #scale by 10^-2 for units conversion to cm^-1
+    inten =  6.0e-02 * pi * volfraction * aa * aa * exp(-bb * bb) / radius
     return inten
 Iq.vectorized =  True  # Iq accepts an array of q values
@@ -71,9 +75,9 @@
             second_moment = 23.0,
             adsorbed_amount = 1.9,
-            density_poly = 0.7,
+            density_shell = 0.7,
             radius = 500.0,
             volfraction = 0.14,
-            polymer_sld = 1.5e-06,
-            solvent_sld = 6.3e-06,
+            sld_shell = 1.5,
+            sld_solvent = 6.3,
             background = 0.0)
 
@@ -82,9 +86,9 @@
                second_moment = 'second_moment',
                adsorbed_amount = 'ads_amount',
-               density_poly = 'density_poly',
+               density_shell = 'density_poly',
                radius = 'radius_core',
                volfraction = 'volf_cores',
-               polymer_sld = 'sld_poly',
-               solvent_sld = 'sld_solv',
+               sld_shell = 'sld_poly',
+               sld_solvent = 'sld_solv',
                background = 'background')
 
@@ -92,6 +96,7 @@
 tests =  [
     [{'scale': 1.0, 'second_moment': 23.0, 'adsorbed_amount': 1.9, 
-     'density_poly': 0.7, 'radius': 500.0, 'volfraction': 0.14, 
-     'polymer_sld': 1.5e-06, 'solvent_sld': 6.3e-06, 'background': 0.0},
+     'density_shell': 0.7, 'radius': 500.0, 'volfraction': 0.14, 
+     'sld_shell': 1.5, 'sld_solvent': 6.3, 'background': 0.0},
      [0.0106939, 0.469418], [73.741, 9.65391e-53]],
     ]
+# ADDED by: SMK  ON: 16Mar2016  convert from sasview, check vs SANDRA, 18Mar2016 RKH some edits & renaming
Index: sasmodels/models/binary_hard_sphere.c
===================================================================
--- sasmodels/models/binary_hard_sphere.c	(revision 115d0f015d008c59c827203a96ee8894b7ec49b4)
+++ sasmodels/models/binary_hard_sphere.c	(revision 115d0f015d008c59c827203a96ee8894b7ec49b4)
@@ -0,0 +1,193 @@
+double form_volume(void);
+
+double Iq(double q, 
+    double lg_radius, double sm_radius,
+    double lg_vol_frac, double sm_vol_frac,
+    double lg_sld, double sm_sld, double solvent_sld
+    );
+    
+double Iqxy(double qx, double qy,
+    double lg_radius, double sm_radius,
+    double lg_vol_frac, double sm_vol_frac,
+    double lg_sld, double sm_sld, double solvent_sld
+    );
+
+int calculate_psfs(double qval,
+    double r2, double nf2,
+    double aa, double phi,
+    double *s11, double *s22, double *s12
+    );
+
+double form_volume(void)
+{
+    return 1.0;
+}
+
+double Iq(double q,
+    double lg_radius, double sm_radius,
+    double lg_vol_frac, double sm_vol_frac,
+    double lg_sld, double sm_sld, double solvent_sld)
+    
+{
+    double r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten;       //my local names
+    double psf11,psf12,psf22;
+    double phi1,phi2,phr,a3;
+    double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2;
+    int err;
+    
+    r2 = lg_radius;
+    r1 = sm_radius;
+    phi2 = lg_vol_frac;
+    phi1 = sm_vol_frac;
+    rho2 = lg_sld;
+    rho1 = sm_sld;
+    rhos = solvent_sld;
+    
+    
+    phi = phi1 + phi2;
+    aa = r1/r2;
+    //calculate the number fraction of larger spheres (eqn 2 in reference)
+    a3=aa*aa*aa;
+    phr=phi2/phi;
+    nf2 = phr*a3/(1.0-phr+phr*a3);
+    // calculate the PSF's here
+    err = calculate_psfs(q,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
+    
+    // /* do form factor calculations  */
+    
+    v1 = 4.0*M_PI/3.0*r1*r1*r1;
+    v2 = 4.0*M_PI/3.0*r2*r2*r2;
+    
+    n1 = phi1/v1;
+    n2 = phi2/v2;
+    
+    qr1 = r1*q;
+    qr2 = r2*q;
+
+    //if (qr1 == 0){
+        //sc1 = 1.0/3.0;
+    //}else{
+        //sc1 = (sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1;
+    //}
+    //if (qr2 == 0){
+        //sc2 = 1.0/3.0;
+    //}else{
+        //sc2 = (sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2;
+    //}
+    sc1 = sph_j1c(qr1);
+    sc2 = sph_j1c(qr2);
+    b1 = r1*r1*r1*(rho1-rhos)*4.0/3.0*M_PI*sc1;
+    b2 = r2*r2*r2*(rho2-rhos)*4.0/3.0*M_PI*sc2;
+    inten = n1*b1*b1*psf11;
+    inten += sqrt(n1*n2)*2.0*b1*b2*psf12;
+    inten += n2*b2*b2*psf22;
+    ///* convert I(1/A) to (1/cm)  */
+    inten *= 1.0e8;
+    ///*convert rho^2 in 10^-6A to A*/
+    inten *= 1.0e-12;    
+    return(inten);
+}
+
+
+double Iqxy(double qx, double qy,
+    double lg_radius, double sm_radius,
+    double lg_vol_frac, double sm_vol_frac,
+    double lg_sld, double sm_sld, double solvent_sld)
+    
+{
+    double q = sqrt(qx*qx + qy*qy);
+    return Iq(q,
+        lg_radius, sm_radius,
+        lg_vol_frac, sm_vol_frac,
+        lg_sld, sm_sld, solvent_sld);
+}
+
+int calculate_psfs(double qval,
+    double r2, double nf2,
+    double aa, double phi,
+    double *s11, double *s22, double *s12)
+
+{
+    //  variable qval,r2,nf2,aa,phi,&s11,&s22,&s12
+    
+    //   calculate constant terms
+    double s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4;
+    double a1,a2i,a2,b1,b2,b12,gm1,gm12;
+    double err=0.0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3;
+    double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13;
+    double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1;
+    
+    s2 = 2.0*r2;
+    s1 = aa*s2;
+    v = phi;
+    a3 = aa*aa*aa;
+    v1=((1.-nf2)*a3/(nf2+(1.-nf2)*a3))*v;
+    v2=(nf2/(nf2+(1.-nf2)*a3))*v;
+    g11=((1.+.5*v)+1.5*v2*(aa-1.))/(1.-v)/(1.-v);
+    g22=((1.+.5*v)+1.5*v1*(1./aa-1.))/(1.-v)/(1.-v);
+    g12=((1.+.5*v)+1.5*(1.-aa)*(v1-v2)/(1.+aa))/(1.-v)/(1.-v);
+    wmv = 1/(1.-v);
+    wmv3 = wmv*wmv*wmv;
+    wmv4 = wmv*wmv3;
+    a1=3.*wmv4*((v1+a3*v2)*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2))) + ((v1+a3*v2)*(1.+2.*v)+(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)-3.*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*wmv3;
+    a2i=((v1+a3*v2)*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*3*wmv4 + ((v1+a3*v2)*(1.+2.*v)+a3*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*aa-3.*v1*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*wmv3;
+    a2=a2i/a3;
+    b1=-6.*(v1*g11*g11+.25*v2*(1.+aa)*(1.+aa)*aa*g12*g12);
+    b2=-6.*(v2*g22*g22+.25*v1/a3*(1.+aa)*(1.+aa)*g12*g12);
+    b12=-3.*aa*(1.+aa)*(v1*g11/aa/aa+v2*g22)*g12;
+    gm1=(v1*a1+a3*v2*a2)*.5;
+    gm12=2.*gm1*(1.-aa)/aa;
+    //c  
+    //c   calculate the direct correlation functions and print results
+    //c
+    //  do 20 j=1,npts
+    
+    yy=qval*s2;
+    //c   calculate direct correlation functions
+    //c   ----c11
+    ay=aa*yy;
+    ay2 = ay*ay;
+    ay3 = ay*ay*ay;
+    t1=a1*(sin(ay)-ay*cos(ay));
+    t2=b1*(2.*ay*sin(ay)-(ay2-2.)*cos(ay)-2.)/ay;
+    t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3;
+    f11=24.*v1*(t1+t2+t3)/ay3;
+    
+    //c ------c22
+    y2=yy*yy;
+    y3=yy*y2;
+    tt1=a2*(sin(yy)-yy*cos(yy));
+    tt2=b2*(2.*yy*sin(yy)-(y2-2.)*cos(yy)-2.)/yy;
+    tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3;
+    f22=24.*v2*(tt1+tt2+tt3)/y3;
+    
+    //c   -----c12
+    yl=.5*yy*(1.-aa);
+    yl3=yl*yl*yl;
+    wma3 = (1.-aa)*(1.-aa)*(1.-aa);
+    y1=aa*yy;
+    y13 = y1*y1*y1;
+    ttt1=3.*wma3*v*sqrt(nf2)*sqrt(1.-nf2)*a1*(sin(yl)-yl*cos(yl))/((nf2+(1.-nf2)*a3)*yl3);
+    t21=b12*(2.*y1*cos(y1)+(y1*y1-2.)*sin(y1));
+    t22=gm12*((3.*y1*y1-6.)*cos(y1)+(y1*y1*y1-6.*y1)*sin(y1)+6.)/y1;
+    t23=gm1*((4.*y13-24.*y1)*cos(y1)+(y13*y1-12.*y1*y1+24.)*sin(y1))/(y1*y1);
+    t31=b12*(2.*y1*sin(y1)-(y1*y1-2.)*cos(y1)-2.);
+    t32=gm12*((3.*y1*y1-6.)*sin(y1)-(y1*y1*y1-6.*y1)*cos(y1))/y1;
+    t33=gm1*((4.*y13-24.*y1)*sin(y1)-(y13*y1-12.*y1*y1+24.)*cos(y1)+24.)/(y1*y1);
+    t41=cos(yl)*((sin(y1)-y1*cos(y1))/(y1*y1) + (1.-aa)/(2.*aa)*(1.-cos(y1))/y1);
+    t42=sin(yl)*((cos(y1)+y1*sin(y1)-1.)/(y1*y1) + (1.-aa)/(2.*aa)*sin(y1)/y1);
+    ttt2=sin(yl)*(t21+t22+t23)/(y13*y1);
+    ttt3=cos(yl)*(t31+t32+t33)/(y13*y1);
+    ttt4=a1*(t41+t42)/y1;
+    f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3);
+    
+    c11=f11;
+    c22=f22;
+    c12=f12;
+    *s11=1./(1.+c11-(c12)*c12/(1.+c22));
+    *s22=1./(1.+c22-(c12)*c12/(1.+c11)); 
+    *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));   
+    
+    return(err);
+}
+
Index: sasmodels/models/binary_hard_sphere.py
===================================================================
--- sasmodels/models/binary_hard_sphere.py	(revision 115d0f015d008c59c827203a96ee8894b7ec49b4)
+++ sasmodels/models/binary_hard_sphere.py	(revision 115d0f015d008c59c827203a96ee8894b7ec49b4)
@@ -0,0 +1,128 @@
+r"""
+
+The binary hard sphere model provides the scattering intensity, for binary
+mixture of hard spheres including hard sphere interaction between those
+particles, using rhw Percus-Yevick closure. The calculation is an exact
+multi-component solution that properly accounts for the 3 partial structure
+factors.
+
+Definition
+----------
+
+.. math::
+
+    \begin{eqnarray}
+    I(q) = (1-x)f_1^2(q) S_{11}(q) + 2[x(1-x)]^{1/2} f_1(q)f_2(q)S_{12}(q) + \\
+    x\,f_2^2(q)S_{22}(q) \\
+    \end{eqnarray}
+
+where $S_{ij}$ are the partial structure factors and $f_i$ are the scattering
+amplitudes of the particles. The subscript 1 is for the smaller particle and 2
+is for the larger. The number fraction of the larger particle,
+($x = n2/(n1+n2)$, where $n$ = the number density) is internally calculated
+based on the diameter ratio and the volume fractions.
+
+.. math::
+
+    \begin{eqnarray}
+    x &=& \frac{(\phi_2 / \phi)\alpha^3}{(1-(\phi_2/\phi) + (\phi_2/\phi) \\
+    \alpha^3)} \phi &=& \phi_1 + \phi_2 = \text{total volume fraction} \\
+    \alpha &=& R_1/R_2 = \text{size ratio}
+    \end{eqnarray}
+
+The 2D scattering intensity is the same as 1D, regardless of the orientation of
+the *q* vector which is defined as
+
+.. math::
+
+    q = \sqrt{q_x^2 + q_y^2}
+
+
+**NOTE 1:** The volume fractions and the scattering contrasts are loosely
+correlated, so holding as many parameters fixed to known values during fitting
+will improve the robustness of the fit.
+
+**NOTE 2:** Since the calculation uses the Percus-Yevick closure, all of the
+limitations of that closure relation apply here. Specifically, one should be
+wary of results for (total) volume fractions greater than approximately 40%.
+Depending on the size ratios or number fractions, the limit on total volume
+fraction may be lower.
+
+**NOTE 3:** The heavy arithmatic operations also mean that at present the
+function is poorly behaved at very low qr.  In some cases very large qr may
+also be poorly behaved.  These should however be outside any useful region of
+qr.
+
+The code for this model is based originally on a c-library implementation by the
+NIST Center for Neutron Research (Kline, 2006).
+
+See the references for details.
+
+References
+----------
+
+N W Ashcroft and D C Langreth, *Physical Review*, 156 (1967) 685-692
+[Errata found in *Phys. Rev.* 166 (1968) 934]
+
+S R Kline, *J Appl. Cryst.*, 39 (2006) 895
+
+**Author:** N/A **on:**
+
+**Modified by:** Paul Butler **on:** March 18, 2016
+
+**Reviewed by:** Paul Butler **on:** March 18, 2016
+"""
+
+import numpy as np
+from numpy import pi, inf
+
+name = "binary_hard_sphere"
+title = "binary mixture of hard spheres with hard sphere interactions."
+description = """Describes the scattering from a mixture of two distinct
+monodisperse, hard sphere particles.
+        [Parameters];
+        radius_lg: large radius of binary hard sphere,
+        radius_sm: small radius of binary hard sphere,
+        volfraction_lg: volume fraction of large spheres,
+        volfraction_sm: volume fraction of small spheres,
+        sld_lg: large sphere  scattering length density,
+        sld_sm: small sphere scattering length density,
+        sld_solvent: solvent scattering length density.
+"""
+category = "shape:sphere"
+
+#             ["name", "units", default, [lower, upper], "type", "description"],
+parameters = [["radius_lg", "Ang", 100, [0, inf], "",
+               "radius of large particle"],
+              ["radius_sm", "Ang", 25, [0, inf], "",
+               "radius of small particle"],
+              ["volfraction_lg", "", 0.1, [0, 1], "",
+               "volume fraction of large particle"],
+              ["volfraction_sm", "", 0.2, [0, 1], "",
+               "volume fraction of small particle"],
+              ["sld_lg", "1e-6/Ang^2", 3.5, [-inf, inf], "",
+               "scattering length density of large particle"],
+              ["sld_sm", "1e-6/Ang^2", 0.5, [-inf, inf], "",
+               "scattering length density of small particle"],
+              ["sld_solvent", "1e-6/Ang^2", 6.36, [-inf, inf], "",
+               "Solvent scattering length density"],
+             ]
+
+source = ["lib/sph_j1c.c", "binary_hard_sphere.c"]
+
+# parameters for demo and documentation
+demo = dict(scale=1, background=0,
+            sld_lg=3.5, sld_sm=0.5, sld_solvent=6.36,
+            radius_lg=100, radius_sm=20,
+            volfraction_lg=0.1, volfraction_sm=0.2)
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname = 'BinaryHSModel'
+oldpars = dict(volfraction_lg='vol_frac_ls', volfraction_sm='vol_frac_ss',
+               sld_lg='ls_sld', sld_sm='ss_sld', sld_solvent='solvent_sld',
+               radius_lg='l_radius', radius_sm='s_radius')
+
+# NOTE: test results taken from values returned by SasView 3.1.2
+tests = [[{'scale':1.0}, 0.001, 25.8927262013]]
+
Index: sasmodels/models/core_shell_ellipsoid.c
===================================================================
--- sasmodels/models/core_shell_ellipsoid.c	(revision e7678b2f993854334aca4c45e1c7022b2fdbb509)
+++ sasmodels/models/core_shell_ellipsoid.c	(revision 29172aa52f4fdcefe534e792b14ec43221514fe1)
@@ -8,7 +8,7 @@
           double equat_shell,
           double polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld);
+          double sld_core,
+          double sld_shell,
+          double sld_solvent);
 
 
@@ -18,7 +18,7 @@
           double equat_shell,
           double polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
+          double sld_core,
+          double sld_shell,
+          double sld_solvent,
           double theta,
           double phi);
@@ -40,7 +40,7 @@
           double equat_shell,
           double polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld)
+          double sld_core,
+          double sld_shell,
+          double sld_solvent)
 {
 
@@ -51,6 +51,6 @@
     double summ = 0.0;	 //initialize intergral
 
-    const double delpc = core_sld - shell_sld;    //core - shell
-    const double delps = shell_sld - solvent_sld; //shell - solvent
+    const double delpc = sld_core - sld_shell;    //core - shell
+    const double delps = sld_shell - sld_solvent; //shell - solvent
 
     for(int i=0;i<N_POINTS_76;i++) {
@@ -81,7 +81,7 @@
           double equat_shell,
           double polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
+          double sld_core,
+          double sld_shell,
+          double sld_solvent,
           double theta,
           double phi)
@@ -96,6 +96,6 @@
     const double cyl_y = sin(theta);
 
-    const double sldcs = core_sld - shell_sld;
-    const double sldss = shell_sld- solvent_sld;
+    const double sldcs = sld_core - sld_shell;
+    const double sldss = sld_shell- sld_solvent;
 
     // Compute the angle btw vector q and the
@@ -124,7 +124,7 @@
           double equat_shell,
           double polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld)
+          double sld_core,
+          double sld_shell,
+          double sld_solvent)
 {
     double intensity = core_shell_ellipsoid_kernel(q,
@@ -133,7 +133,7 @@
            equat_shell,
            polar_shell,
-           core_sld,
-           shell_sld,
-           solvent_sld);
+           sld_core,
+           sld_shell,
+           sld_solvent);
 
     return intensity;
@@ -146,7 +146,7 @@
           double equat_shell,
           double polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
+          double sld_core,
+          double sld_shell,
+          double sld_solvent,
           double theta,
           double phi)
@@ -159,7 +159,7 @@
                        equat_shell,
                        polar_shell,
-                       core_sld,
-                       shell_sld,
-                       solvent_sld,
+                       sld_core,
+                       sld_shell,
+                       sld_solvent,
                        theta,
                        phi);
Index: sasmodels/models/core_shell_ellipsoid.py
===================================================================
--- sasmodels/models/core_shell_ellipsoid.py	(revision 5111921153c056b59ff03411e4b7b028d238ae26)
+++ sasmodels/models/core_shell_ellipsoid.py	(revision 29172aa52f4fdcefe534e792b14ec43221514fe1)
@@ -1,5 +1,5 @@
 r"""
 This model provides the form factor, $P(q)$, for a core shell ellipsoid (below)
-where the form factor is normalized by the volume of the cylinder.
+where the form factor is normalized by the volume of the outer [CHECK].
 
 .. math::
@@ -7,5 +7,5 @@
     P(q) = scale * \left<f^2\right>/V + background
 
-where the volume $V = (4/3)\pi(r_{maj}r_{min}^2)$ and the averaging $< >$ is
+where the volume $V = (4/3)\pi(R_{major\_outer}R_{minor\_outer}^2)$ and the averaging $< >$ is
 applied over all orientations for 1D.
 
@@ -22,9 +22,9 @@
 
     P(q) = \frac{scale}{V}\int_0^1
-        \left|F(q,r_{min},r_{max},\alpha)\right|^2d\alpha + background
+        \left|F(q,r_{minor\_core},r_{major\_core},\alpha) + F(q,r_{major\_outer},r_{major\_outer},\alpha)\right|^2d\alpha + background
 
-    \left|F(q,r_{min},r_{max},\alpha)\right|=V\Delta \rho \cdot (3j_1(u)/u)
+    \left|F(q,r_{minor},r_{major},\alpha)\right|=(4\pi/3)r_{major}r_{minor}^2 \Delta \rho \cdot (3j_1(u)/u)
 
-    u = q\left[ r_{maj}^2\alpha ^2 + r_{min}^2(1-\alpha ^2)\right]^{1/2}
+    u = q\left[ r_{major}^2\alpha ^2 + r_{minor}^2(1-\alpha ^2)\right]^{1/2}
 
 where
@@ -43,4 +43,9 @@
 polar core radius, *equat_shell* = $r_{min}$ (or equatorial outer radius),
 and *polar_shell* = $r_{maj}$ (or polar outer radius).
+
+Note:It is the users' responsibility to ensure that shell radii are larger than 
+the core radii, especially if both are polydisperse, in which case the
+core_shell_ellipsoid_xt model may be much better.
+
 
 .. note::
@@ -77,15 +82,16 @@
     single particle scattering amplitude.
     [Parameters]:
-    equat_core = equatorial radius of core,
-    polar_core = polar radius of core,
-    equat_shell = equatorial radius of shell,
-    polar_shell = polar radius (revolution axis) of shell,
-    core_sld = SLD_core
-    shell_sld = SLD_shell
-    solvent_sld = SLD_solvent
+    equat_core = equatorial radius of core, Rminor_core,
+    polar_core = polar radius of core, Rmajor_core,
+    equat_shell = equatorial radius of shell, Rminor_outer,
+    polar_shell = polar radius of shell, Rmajor_outer,
+    sld_core = scattering length density of core,
+    sld_shell = scattering length density of shell,
+    sld_solvent = scattering length density of solvent,
     background = Incoherent bkg
     scale =scale
     Note:It is the users' responsibility to ensure
-    that shell radii are larger than core radii.
+    that shell radii are larger than core radii,
+    especially if both are polydisperse.
     oblate: polar radius < equatorial radius
     prolate :  polar radius > equatorial radius
@@ -97,11 +103,11 @@
 #             ["name", "units", default, [lower, upper], "type", "description"],
 parameters = [
-    ["equat_core",  "Ang",      200,   [0, inf],    "volume",      "Equatorial radius of core"],
-    ["polar_core",  "Ang",       10,   [0, inf],    "volume",      "Polar radius of core"],
-    ["equat_shell", "Ang",      250,   [0, inf],    "volume",      "Equatorial radius of shell"],
-    ["polar_shell", "Ang",       30,   [0, inf],    "volume",      "Polar radius of shell"],
-    ["core_sld",    "1e-6/Ang^2", 2,   [-inf, inf], "",            "Core scattering length density"],
-    ["shell_sld",   "1e-6/Ang^2", 1,   [-inf, inf], "",            "Shell scattering length density"],
-    ["solvent_sld", "1e-6/Ang^2", 6.3, [-inf, inf], "",            "Solvent scattering length density"],
+    ["equat_core",  "Ang",      200,   [0, inf],    "volume",      "Equatorial radius of core, Rminor_core "],
+    ["polar_core",  "Ang",       10,   [0, inf],    "volume",      "Polar radius of core, Rmajor_core"],
+    ["equat_shell", "Ang",      250,   [0, inf],    "volume",      "Equatorial radius of shell, Rminor_outer "],
+    ["polar_shell", "Ang",       30,   [0, inf],    "volume",      "Polar radius of shell, Rmajor_outer"],
+    ["sld_core",    "1e-6/Ang^2", 2,   [-inf, inf], "sld",            "Core scattering length density"],
+    ["sld_shell",   "1e-6/Ang^2", 1,   [-inf, inf], "sld",            "Shell scattering length density"],
+    ["sld_solvent", "1e-6/Ang^2", 6.3, [-inf, inf], "sld",            "Solvent scattering length density"],
     ["theta",       "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation wrt incoming beam"],
     ["phi",         "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation in the plane of the detector"],
@@ -116,7 +122,7 @@
             equat_shell=250.0,
             polar_shell=30.0,
-            core_sld=2.0,
-            shell_sld=1.0,
-            solvent_sld=6.3,
+            sld_core=2.0,
+            sld_shell=1.0,
+            sld_solvent=6.3,
             theta=0,
             phi=0)
@@ -124,7 +130,8 @@
 oldname = 'CoreShellEllipsoidModel'
 
-oldpars = dict(core_sld='sld_core',
-               shell_sld='sld_shell',
-               solvent_sld='sld_solvent',
+oldpars = dict(equat_core='equat_core',polar_core='polar_core',equat_shell='equat_shell',polar_shell='polar_shell',
+               sld_core='sld_core',
+               sld_shell='sld_shell',
+               sld_solvent='sld_solvent',
                theta='axis_theta',
                phi='axis_phi')
@@ -141,7 +148,7 @@
       'equat_shell': 250.0,
       'polar_shell': 30.0,
-      'core_sld': 2.0,
-      'shell_sld': 1.0,
-      'solvent_sld': 6.3,
+      'sld_core': 2.0,
+      'sld_shell': 1.0,
+      'sld_solvent': 6.3,
       'background': 0.001,
       'scale': 1.0,
@@ -155,7 +162,7 @@
       'equat_shell': 54.0,
       'polar_shell': 3.0,
-      'core_sld': 20.0,
-      'shell_sld': 10.0,
-      'solvent_sld': 6.0,
+      'sld_core': 20.0,
+      'sld_shell': 10.0,
+      'sld_solvent': 6.0,
       'background': 0.0,
       'scale': 1.0,
@@ -168,7 +175,7 @@
       'equat_shell': 54.0,
       'polar_shell': 3.0,
-      'core_sld': 20.0,
-      'shell_sld': 10.0,
-      'solvent_sld': 6.0,
+      'sld_core': 20.0,
+      'sld_shell': 10.0,
+      'sld_solvent': 6.0,
       'background': 0.01,
       'scale': 0.01,
Index: sasmodels/models/core_shell_ellipsoid_xt.py
===================================================================
--- sasmodels/models/core_shell_ellipsoid_xt.py	(revision 5111921153c056b59ff03411e4b7b028d238ae26)
+++ sasmodels/models/core_shell_ellipsoid_xt.py	(revision 29172aa52f4fdcefe534e792b14ec43221514fe1)
@@ -1,9 +1,9 @@
 r"""
-An alternative version of $P(q)$ for the core-shell ellipsoid
-(see CoreShellEllipsoidModel), having as parameters the core axial ratio X
-and a shell thickness, which are more often what we would like to determine.
+An alternative version of $P(q)$ for the core_shell_ellipsoid
+having as parameters the core axial ratio X and a shell thickness, 
+which are more often what we would like to determine.
 
 This model is also better behaved when polydispersity is applied than the four
-independent radii in CoreShellEllipsoidModel.
+independent radii in core_shell_ellipsoid model.
 
 Definition
@@ -52,5 +52,5 @@
 ----------
 
-R K Heenan, *Private communication*
+R K Heenan, 2015, reparametrised the core_shell_ellipsoid model
 
 """
@@ -69,10 +69,10 @@
         [Parameters]:
         equat_core = equatorial radius of core,
-        x_core = polar radius of core,
+        x_core = ratio of core polar/equatorial radii,
         t_shell = equatorial radius of outer surface,
-        x_polar_shell = polar radius (revolution axis) of outer surface
-        core_sld = SLD_core
-        shell_sld = SLD_shell
-        solvent_sld = SLD_solvent
+        x_polar_shell = ratio of polar shell thickness to equatorial shell thickness,
+        sld_core = SLD_core
+        sld_shell = SLD_shell
+        sld_solvent = SLD_solvent
         background = Incoherent bkg
         scale =scale
@@ -92,7 +92,7 @@
     ["t_shell",       "Ang",       30,   [0, inf],    "volume",      "Equatorial radius of shell"],
     ["x_polar_shell", "",           1,   [0, inf],    "volume",      "Polar radius of shell"],
-    ["core_sld",      "1e-6/Ang^2", 2,   [-inf, inf], "",            "Core scattering length density"],
-    ["shell_sld",     "1e-6/Ang^2", 1,   [-inf, inf], "",            "Shell scattering length density"],
-    ["solvent_sld",   "1e-6/Ang^2", 6.3, [-inf, inf], "",            "Solvent scattering length density"],
+    ["sld_core",      "1e-6/Ang^2", 2,   [-inf, inf], "",            "Core scattering length density"],
+    ["sld_shell",     "1e-6/Ang^2", 1,   [-inf, inf], "",            "Shell scattering length density"],
+    ["sld_solvent",   "1e-6/Ang^2", 6.3, [-inf, inf], "",            "Solvent scattering length density"],
     ["theta",         "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation wrt incoming beam"],
     ["phi",           "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation in the plane of the detector"],
@@ -108,7 +108,7 @@
             t_shell=30.0,
             x_polar_shell=1.0,
-            core_sld=2.0,
-            shell_sld=1.0,
-            solvent_sld=6.3,
+            sld_core=2.0,
+            sld_shell=1.0,
+            sld_solvent=6.3,
             theta=0,
             phi=0)
@@ -119,7 +119,7 @@
                x_core='X_core',
                t_shell='T_shell',
-               core_sld='sld_core',
-               shell_sld='sld_shell',
-               solvent_sld='sld_solvent',
+               sld_core='sld_core',
+               sld_shell='sld_shell',
+               sld_solvent='sld_solvent',
                theta='axis_theta',
                phi='axis_phi')
@@ -136,7 +136,7 @@
       't_shell': 50.0,
       'x_polar_shell': 0.2,
-      'core_sld': 2.0,
-      'shell_sld': 1.0,
-      'solvent_sld': 6.3,
+      'sld_core': 2.0,
+      'sld_shell': 1.0,
+      'sld_solvent': 6.3,
       'background': 0.001,
       'scale': 1.0,
@@ -150,7 +150,7 @@
       't_shell': 54.0,
       'x_polar_shell': 3.0,
-      'core_sld': 20.0,
-      'shell_sld': 10.0,
-      'solvent_sld': 6.0,
+      'sld_core': 20.0,
+      'sld_shell': 10.0,
+      'sld_solvent': 6.0,
       'background': 0.0,
       'scale': 1.0,
@@ -163,7 +163,7 @@
       't_shell': 54.0,
       'x_polar_shell': 3.0,
-      'core_sld': 20.0,
-      'shell_sld': 10.0,
-      'solvent_sld': 6.0,
+      'sld_core': 20.0,
+      'sld_shell': 10.0,
+      'sld_solvent': 6.0,
       'background': 0.01,
       'scale': 0.01,
Index: sasmodels/models/hardsphere.py
===================================================================
--- sasmodels/models/hardsphere.py	(revision 3bcd03d7e1d8d4938d30f17bf360b83417139376)
+++ sasmodels/models/hardsphere.py	(revision e98c1e008cb289d6be876cc1dfb5c76efa378204)
@@ -58,4 +58,5 @@
 category = "structure-factor"
 structure_factor = True
+single = False
 
 #             ["name", "units", default, [lower, upper], "type","description"],
Index: sasmodels/models/lamellar.py
===================================================================
--- sasmodels/models/lamellar.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ sasmodels/models/lamellar.py	(revision 7c391dd12853b2e9403ba35b776d89aeb6838937)
@@ -5,9 +5,9 @@
 ----------
 
-The scattering intensity $I(q)$ is
+The scattering intensity $I(q)$ for dilute, randomly oriented, "infinitely large" sheets or lamellae is
 
 .. math::
 
-    I(q) = \frac{2\pi P(q)}{\delta q^2}
+    I(q) = scale*\frac{2\pi P(q)}{q^2\delta }
 
 
@@ -16,8 +16,12 @@
 .. math::
 
-    P(q) = \frac{2\Delta\rho^2}{q^2}(1-cos(q\delta))
+   P(q) = \frac{2\Delta\rho^2}{q^2}(1-cos(q\delta)) = \frac{4\Delta\rho^2}{q^2}sin^2(\frac{q\delta}{2})
 
+where $\delta$ is the total layer thickness and $\Delta\rho$ is the scattering length density difference.
 
-where $\delta$ is the bilayer thickness.
+This is the limiting form for a spherical shell of infinitely large radius. Note that the division by $\delta$
+means that $scale$ in sasview is the volume fraction of sheet, $\phi = S\delta$ where $S$ is the area of 
+sheet per unit volume. $S$ is half the Porod surface area per unit volume of a thicker layer (as that would 
+include both faces of the sheet).
 
 The 2D scattering intensity is calculated in the same way as 1D, where
@@ -56,14 +60,11 @@
 
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["sld", "1e-6/Ang^2", 1, [-inf, inf], "",
-               "Layer scattering length density" ],
-              ["solvent_sld", "1e-6/Ang^2", 6, [-inf, inf], "",
-               "Solvent scattering length density" ],
-              ["thickness", "Ang", 50, [0, inf], "volume","Bilayer thickness" ],
+parameters = [ ["thickness", "Ang", 50, [0, inf], "volume","total layer thickness" ],
+               ["sld", "1e-6/Ang^2", 1, [-inf, inf], "","Layer scattering length density" ],
+               ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "","Solvent scattering length density" ],
              ]
 
-
 # No volume normalization despite having a volume parameter
-# This should perhaps be volume normalized?
+# This should perhaps be volume normalized? - it is!
 form_volume = """
     return 1.0;
@@ -71,5 +72,5 @@
 
 Iq = """
-    const double sub = sld - solvent_sld;
+    const double sub = sld - sld_solvent;
     const double qsq = q*q;
     // Original expression
@@ -89,11 +90,11 @@
 
 demo = dict(scale=1, background=0,
-            sld=6, solvent_sld=1,
+            sld=6, sld_solvent=1,
             thickness=40,
             thickness_pd=0.2, thickness_pd_n=40)
 oldname = 'LamellarModel'
-oldpars = dict(sld='sld_bi', solvent_sld='sld_sol', thickness='bi_thick')
+oldpars = dict(sld='sld_bi', sld_solvent='sld_sol', thickness='bi_thick')
 tests = [
-        [ {'scale': 1.0, 'background' : 0.0, 'thickness' : 50.0, 'sld' : 1.0,'solvent_sld' : 6.3, 'thickness_pd' : 0.0, 
+        [ {'scale': 1.0, 'background' : 0.0, 'thickness' : 50.0, 'sld' : 1.0,'sld_solvent' : 6.3, 'thickness_pd' : 0.0, 
            }, [0.001], [882289.54309]]
         ]
Index: smodels/models/lamellarCaille.py
===================================================================
--- sasmodels/models/lamellarCaille.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ 	(revision )
@@ -1,129 +1,0 @@
-r"""
-This model provides the scattering intensity, $I(q) = P(q) S(q)$, for a
-lamellar phase where a random distribution in solution are assumed.
-Here a Caille $S(q)$ is used for the lamellar stacks.
-
-Definition
-----------
-
-The scattering intensity $I(q)$ is
-
-.. math::
-
-    I(q) = 2\pi \frac{P(q)S(q)}{\delta q^2}
-
-The form factor is
-
-.. math::
-
-    P(q) = \frac{2\Delta\rho^2}{q^2}\left(1-\cos q\delta \right)
-
-and the structure factor is
-
-.. math::
-
-    S(q) = 1 + 2 \sum_1^{N-1}\left(1-\frac{n}{N}\right)
-           \cos(qdn)\exp\left(-\frac{2q^2d^2\alpha(n)}{2}\right)
-
-where
-
-.. math::
-    :nowrap:
-
-    \begin{align*}
-    \alpha(n) &= \frac{\eta_{cp}}{4\pi^2} \left(\ln(\pi n)+\gamma_E\right)
-              && \\
-    \gamma_E  &= 0.5772156649
-              && \text{Euler's constant} \\
-    \eta_{cp} &= \frac{q_o^2k_B T}{8\pi\sqrt{K\overline{B}}}
-              && \text{Caille constant}
-    \end{align*}
-
-Here $d$ = (repeat) spacing, $\delta$ = bilayer thickness,
-the contrast $\Delta\rho$ = SLD(headgroup) - SLD(solvent),
-$K$ = smectic bending elasticity, $B$ = compression modulus, and
-$N$ = number of lamellar plates (*n_plates*).
-
-NB: **When the Caille parameter is greater than approximately 0.8 to 1.0, the
-assumptions of the model are incorrect.** And due to a complication of the
-model function, users are responsible for making sure that all the assumptions
-are handled accurately (see the original reference below for more details).
-
-Non-integer numbers of stacks are calculated as a linear combination of
-results for the next lower and higher values.
-
-The 2D scattering intensity is calculated in the same way as 1D, where the
-$q$ vector is defined as
-
-.. math::
-
-    q = \sqrt{q_x^2 + q_y^2}
-
-
-References
-----------
-
-F Nallet, R Laversanne, and D Roux, J. Phys. II France, 3, (1993) 487-502
-
-also in J. Phys. Chem. B, 105, (2001) 11081-11088
-"""
-from numpy import inf
-
-name = "lamellarPS"
-title = "Random lamellar sheet with Caille structure factor"
-description = """\
-    [Random lamellar phase with Caille  structure factor]
-    randomly oriented stacks of infinite sheets
-    with Caille S(Q), having polydisperse spacing.
-    sld = sheet scattering length density
-    sld_solvent = solvent scattering length density
-    background = incoherent background
-    scale = scale factor
-"""
-category = "shape:lamellae"
-
-single = False
-# pylint: disable=bad-whitespace, line-too-long
-#             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [
-    ["thickness",        "Ang",      30.0,  [0, inf],   "volume", "sheet thickness"],
-    ["Nlayers",          "",          20,   [0, inf],   "",       "Number of layers"],
-    ["spacing",          "Ang",      400.,  [0.0,inf],  "volume", "d-spacing of Caille S(Q)"],
-    ["Caille_parameter", "1/Ang^2",    0.1, [0.0,0.8],  "",       "Caille parameter"],
-    ["sld",              "1e-6/Ang^2", 6.3, [-inf,inf], "",       "layer scattering length density"],
-    ["solvent_sld",      "1e-6/Ang^2", 1.0, [-inf,inf], "",       "Solvent scattering length density"],
-    ]
-# pylint: enable=bad-whitespace, line-too-long
-
-source = ["lamellarCaille_kernel.c"]
-
-# No volume normalization despite having a volume parameter
-# This should perhaps be volume normalized?
-form_volume = """
-    return 1.0;
-    """
-
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
-# ER defaults to 0.0
-# VR defaults to 1.0
-
-demo = dict(scale=1, background=0,
-            thickness=67., Nlayers=3.75, spacing=200.,
-            Caille_parameter=0.268, sld=1.0, solvent_sld=6.34,
-            thickness_pd=0.1, thickness_pd_n=100,
-            spacing_pd=0.05, spacing_pd_n=40)
-
-oldname = 'LamellarPSModel'
-oldpars = dict(thickness='delta', Nlayers='N_plates',
-               Caille_parameter='caille',
-               sld='sld_bi', solvent_sld='sld_sol')
-#
-tests = [
-    [{'scale': 1.0, 'background': 0.0, 'thickness': 30., 'Nlayers': 20.0,
-      'spacing': 400., 'Caille_parameter': 0.1, 'sld': 6.3,
-      'solvent_sld': 1.0, 'thickness_pd': 0.0, 'spacing_pd': 0.0},
-     [0.001], [28895.13397]]
-    ]
Index: smodels/models/lamellarCailleHG.py
===================================================================
--- sasmodels/models/lamellarCailleHG.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ 	(revision )
@@ -1,148 +1,0 @@
-# Note: model title and parameter table are inserted automatically
-r"""
-This model provides the scattering intensity, $I(q) = P(q)S(q)$, for a lamellar
-phase where a random distribution in solution are assumed. Here a Caille $S(q)$
-is used for the lamellar stacks.
-
-The scattering intensity $I(q)$ is
-
-.. math::
-
-    I(q) = 2 \pi \frac{P(q)S(q)}{\delta q^2}
-
-
-The form factor $P(q)$ is
-
-.. math::
-
-        P(q) = \frac{4}{q^2}\big\{
-        \Delta\rho_H \left[\sin[q(\delta_H + \delta_T)] - \sin(q\delta_T)\right]
-            + \Delta\rho_T\sin(q\delta_T)\big\}^2
-
-and the structure factor $S(q)$ is
-
-.. math::
-
-    S(q) = 1 + 2 \sum_1^{N-1}\left(1-\frac{n}{N}\right)
-        \cos(qdn)\exp\left(-\frac{2q^2d^2\alpha(n)}{2}\right)
-
-where
-
-.. math::
-    :nowrap:
-
-    \begin{align*}
-    \alpha(n) &= \frac{\eta_{cp}}{4\pi^2} \left(\ln(\pi n)+\gamma_E\right)
-              &&  \\
-    \gamma_E  &= 0.5772156649
-              && \text{Euler's constant} \\
-    \eta_{cp} &= \frac{q_o^2k_B T}{8\pi\sqrt{K\overline{B}}}
-              && \text{Caille constant}
-    \end{align*}
-
-
-$\delta_T$ is the tail length (or *tail_length*), $\delta_H$ is the head
-thickness (or *head_length*), $\Delta\rho_H$ is SLD(headgroup) - SLD(solvent),
-and $\Delta\rho_T$ is SLD(tail) - SLD(headgroup). Here $d$ is (repeat) spacing,
-$K$ is smectic bending elasticity, $B$ is compression modulus, and $N$ is the
-number of lamellar plates (*Nlayers*).
-
-NB: **When the Caille parameter is greater than approximately 0.8 to 1.0, the
-assumptions of the model are incorrect.**  And due to a complication of the
-model function, users are responsible for making sure that all the assumptions
-are handled accurately (see the original reference below for more details).
-
-Non-integer numbers of stacks are calculated as a linear combination of
-results for the next lower and higher values.
-
-The 2D scattering intensity is calculated in the same way as 1D, where
-the $q$ vector is defined as
-
-.. math::
-
-    q = \sqrt{q_x^2 + q_y^2}
-
-
-References
-----------
-
-F Nallet, R Laversanne, and D Roux, J. Phys. II France, 3, (1993) 487-502
-
-also in J. Phys. Chem. B, 105, (2001) 11081-11088
-"""
-from numpy import inf
-
-name = "lamellarCailleHG"
-title = "Random lamellar sheet with Caille structure factor"
-description = """\
-    [Random lamellar phase with Caille  structure factor]
-        randomly oriented stacks of infinite sheets
-        with Caille S(Q), having polydisperse spacing.
-        layer thickness =(H+T+T+H) = 2(Head+Tail)
-        sld = Tail scattering length density
-        sld_head = Head scattering length density
-        sld_solvent = solvent scattering length density
-        background = incoherent background
-        scale = scale factor
-"""
-category = "shape:lamellae"
-
-single = False
-parameters = [
-    #   [ "name", "units", default, [lower, upper], "type",
-    #     "description" ],
-    ["tail_length", "Ang", 10, [0, inf], "volume",
-     "Tail thickness"],
-    ["head_length", "Ang", 2, [0, inf], "volume",
-     "head thickness"],
-    ["Nlayers", "", 30, [0, inf], "",
-     "Number of layers"],
-    ["spacing", "Ang", 40., [0.0, inf], "volume",
-     "d-spacing of Caille S(Q)"],
-    ["Caille_parameter", "", 0.001, [0.0, 0.8], "",
-     "Caille parameter"],
-    ["sld", "1e-6/Ang^2", 0.4, [-inf, inf], "",
-     "Tail scattering length density"],
-    ["head_sld", "1e-6/Ang^2", 2.0, [-inf, inf], "",
-     "Head scattering length density"],
-    ["solvent_sld", "1e-6/Ang^2", 6, [-inf, inf], "",
-     "Solvent scattering length density"],
-    ]
-
-source = ["lamellarCailleHG_kernel.c"]
-
-# No volume normalization despite having a volume parameter
-# This should perhaps be volume normalized?
-form_volume = """
-    return 1.0;
-    """
-
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
-# ER defaults to 0.0
-# VR defaults to 1.0
-
-demo = dict(
-    scale=1, background=0,
-    Nlayers=20, spacing=200., Caille_parameter=0.05,
-    tail_length=15, head_length=10,
-    #sld=-1, head_sld=4.0, solvent_sld=6.0,
-    sld=-1, head_sld=4.1, solvent_sld=6.0,
-    tail_length_pd=0.1, tail_length_pd_n=20,
-    head_length_pd=0.05, head_length_pd_n=30,
-    spacing_pd=0.2, spacing_pd_n=40,
-    )
-
-oldname = 'LamellarPSHGModel'
-
-oldpars = dict(
-    tail_length='deltaT', head_length='deltaH', Nlayers='n_plates',
-    Caille_parameter='caille', sld='sld_tail', head_sld='sld_head',
-    solvent_sld='sld_solvent')
-#
-tests = [[{'scale': 1.0, 'background': 0.0, 'tail_length': 10.0, 'head_length': 2.0,
-           'Nlayers': 30.0, 'spacing': 40., 'Caille_parameter': 0.001, 'sld': 0.4,
-           'head_sld': 2.0, 'solvent_sld': 6.0, 'tail_length_pd': 0.0,
-           'head_length_pd': 0.0, 'spacing_pd': 0.0}, [0.001], [6838238.571488]]]
Index: smodels/models/lamellarCailleHG_kernel.c
===================================================================
--- sasmodels/models/lamellarCailleHG_kernel.c	(revision 12c810f78dbf238faaad5e049615baec671280ae)
+++ 	(revision )
@@ -1,75 +1,0 @@
-/*	LamellarCailleHG kernel - allows for name changes of passed parameters ...
-    Maths identical to LamellarCaille apart from the line for P(Q)
-*/
-
-double Iq(double qval,
-      double tail_length,
-      double head_length,
-      double Nlayers, 
-      double dd,
-      double Cp,
-      double tail_sld,
-      double head_sld,
-      double solvent_sld);
-
-double Iq(double qval,
-      double tail_length,
-      double head_length,
-      double Nlayers, 
-      double dd,
-      double Cp,
-      double tail_sld,
-      double head_sld,
-      double solvent_sld)
-{
-  double NN;   //local variables of coefficient wave
-  double inten,Pq,Sq,alpha,temp,t2;
-  //double dQ, dQDefault, t1, t3;
-  int ii,NNint;
-  // from wikipedia 0.577215664901532860606512090082402431042159335
-  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015
-  //dQDefault = 0.0;    //[=] 1/A, q-resolution, default value
-  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015
-
-  NN = trunc(Nlayers);    //be sure that NN is an integer
-  
-  Pq = (head_sld-solvent_sld)*(sin(qval*(head_length+tail_length))-sin(qval*tail_length)) +
-              (tail_sld-solvent_sld)*sin(qval*tail_length);
-  Pq *= Pq;
-  Pq *= 4.0/(qval*qval);
-
-  NNint = (int)NN;    //cast to an integer for the loop
-  ii=0;
-  Sq = 0.0;
-  for(ii=1;ii<=(NNint-1);ii+=1) {
-
-    //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015
-
-    temp = 0.0;
-    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler);
-    //t1 = 2.0*dQ*dQ*dd*dd*alpha;
-    t2 = 2.0*qval*qval*dd*dd*alpha;
-    //t3 = dQ*dQ*dd*dd*ii*ii;
-
-    temp = 1.0-ii/NN;
-    //temp *= cos(dd*qval*ii/(1.0+t1));
-    temp *= cos(dd*qval*ii);
-    //if (temp < 0) printf("q=%g: ii=%d, cos(dd*q*ii)=cos(%g) < 0\n",qval,ii,dd*qval*ii);
-    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
-    temp *= exp(-t2/2.0 );
-    //temp /= sqrt(1.0+t1);
-
-    Sq += temp;
-  }
-
-  Sq *= 2.0;
-  Sq += 1.0;
-
-  //if (Sq < 0) printf("q=%g: S(q) =%g\n", qval, Sq);
-
-  inten = 2.0*M_PI*Pq*Sq/(dd*qval*qval);
-
-  inten *= 1.0e-04;   // 1/A to 1/cm
-  return(inten);
-}
-
Index: smodels/models/lamellarCaille_kernel.c
===================================================================
--- sasmodels/models/lamellarCaille_kernel.c	(revision 2e35a40faa66788cc99598aba7a610360b33e747)
+++ 	(revision )
@@ -1,69 +1,0 @@
-/*	LamellarCaille kernel - allows for name changes of passed parameters ...
-
-*/
-
-double Iq(double qval,
-      double del,
-      double Nlayers, 
-      double dd,
-	  double Cp, 
-      double sld,
-      double solvent_sld);
-
-double Iq(double qval,
-      double del,
-      double Nlayers, 
-      double dd,
-	  double Cp, 
-      double sld,
-      double solvent_sld)
-{
-  double contr,NN;   //local variables of coefficient wave
-  double inten,Pq,Sq,alpha,temp,t2;
-  //double dQ, dQDefault, t1, t3;
-  int ii,NNint;
-  // from wikipedia 0.577215664901532860606512090082402431042159335
-  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015
-  //dQDefault = 0.0;    //[=] 1/A, q-resolution, default value
-  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015
-
-  NN = trunc(Nlayers);    //be sure that NN is an integer
-  
-  contr = sld - solvent_sld;
-
-  Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del));
-
-  NNint = (int)NN;    //cast to an integer for the loop
-  ii=0;
-  Sq = 0.0;
-  // the vital "=" in ii<=  added March 2015
-  for(ii=1;ii<=(NNint-1);ii+=1) {
-
-    //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015
-
-    temp = 0.0;
-    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler);
-    //t1 = 2.0*dQ*dQ*dd*dd*alpha;
-    t2 = 2.0*qval*qval*dd*dd*alpha;
-    //t3 = dQ*dQ*dd*dd*ii*ii;
-
-    temp = 1.0-ii/NN;
-    //temp *= cos(dd*qval*ii/(1.0+t1));
-    temp *= cos(dd*qval*ii);
-    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
-    temp *= exp(-t2/2.0 );
-    //temp /= sqrt(1.0+t1);
-
-    Sq += temp;
-  }
-
-  Sq *= 2.0;
-  Sq += 1.0;
-
-  inten = 2.0*M_PI*Pq*Sq/(dd*qval*qval);
-
-  inten *= 1.0e-04;   // 1/A to 1/cm
-
-  return(inten);
-}
-
Index: smodels/models/lamellarFFHG.py
===================================================================
--- sasmodels/models/lamellarFFHG.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ 	(revision )
@@ -1,117 +1,0 @@
-# Note: model title and parameter table are inserted automatically
-r"""
-This model provides the scattering intensity, $I(q)$, for a lyotropic lamellar
-phase where a random distribution in solution are assumed. The SLD of the head
-region is taken to be different from the SLD of the tail region.
-
-Definition
-----------
-
-The scattering intensity $I(q)$ is
-
-.. math::
-
-   I(q) = 2\pi\frac{\text{scale}}{2(\delta_H + \delta_T)}  P(q) \frac{1}{q^2}
-
-The form factor $P(q)$ is
-
-.. math::
-
-    P(q) = \frac{4}{q^2}
-        \left\lbrace
-            \Delta \rho_H
-            \left[\sin[q(\delta_H + \delta_T)\ - \sin(q\delta_T)\right]
-            + \Delta\rho_T\sin(q\delta_T)
-        \right\rbrace^2
-
-where $\delta_T$ is *tail_length*, $\delta_H$ is *head_length*,
-$\Delta\rho_H$ is the head contrast (*head_sld* $-$ *solvent_sld*),
-and $\Delta\rho_T$ is tail contrast (*sld* $-$ *solvent_sld*).
-
-The 2D scattering intensity is calculated in the same way as 1D, where
-the $q$ vector is defined as
-
-.. math::
-
-    q = \sqrt{q_x^2 + q_y^2}
-
-
-References
-----------
-
-F Nallet, R Laversanne, and D Roux, J. Phys. II France, 3, (1993) 487-502
-
-also in J. Phys. Chem. B, 105, (2001) 11081-11088
-
-*2014/04/17 - Description reviewed by S King and P Butler.*
-"""
-
-from numpy import inf
-
-name = "lamellar_FFHG"
-title = "Random lamellar phase with Head Groups"
-description = """\
-    [Random lamellar phase with Head Groups]
-        I(q)= 2*pi*P(q)/(2(H+T)*q^(2)), where
-        P(q)= see manual
-        layer thickness =(H+T+T+H) = 2(Head+Tail)
-        sld = Tail scattering length density
-        head_sld = Head scattering length density
-        solvent_sld = solvent scattering length density
-        background = incoherent background
-        scale = scale factor
-"""
-category = "shape:lamellae"
-
-# pylint: disable=bad-whitespace, line-too-long
-#             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["tail_length", "Ang",       15,   [0, inf],  "volume",  "Tail thickness"],
-              ["head_length", "Ang",       10,   [0, inf],  "volume",  "head thickness"],
-              ["sld",         "1e-6/Ang^2", 0.4, [-inf,inf], "",       "Tail scattering length density"],
-              ["head_sld",    "1e-6/Ang^2", 3.0, [-inf,inf], "",       "Head scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2", 6,   [-inf,inf], "",       "Solvent scattering length density"]]
-# pylint: enable=bad-whitespace, line-too-long
-
-# No volume normalization despite having a volume parameter
-# This should perhaps be volume normalized?
-form_volume = """
-    return 1.0;
-    """
-
-Iq = """
-    const double qsq = q*q;
-    const double drh = head_sld - solvent_sld;
-    const double drt = sld - solvent_sld;    //correction 13FEB06 by L.Porcar
-    const double qT = q*tail_length;
-    double Pq, inten;
-    Pq = drh*(sin(q*(head_length+tail_length))-sin(qT)) + drt*sin(qT);
-    Pq *= Pq;
-    Pq *= 4.0/(qsq);
-
-    inten = 2.0e-4*M_PI*Pq/qsq;
-
-    // normalize by the bilayer thickness
-    inten /= 2.0*(head_length+tail_length);
-
-    return inten;
-    """
-
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
-# ER defaults to 0.0
-# VR defaults to 1.0
-
-demo = dict(scale=1, background=0,
-            tail_length=15, head_length=10,
-            sld=0.4, head_sld=3.0, solvent_sld=6.0,
-            tail_length_pd=0.2, tail_length_pd_n=40,
-            head_length_pd=0.01, head_length_pd_n=40)
-
-oldname = 'LamellarFFHGModel'
-oldpars = dict(head_length='h_thickness', tail_length='t_length',
-               sld='sld_tail', head_sld='sld_head', solvent_sld='sld_solvent')
-#
-tests = [[{'scale': 1.0, 'background': 0.0, 'tail_length': 15.0, 'head_length': 10.0,
-           'sld': 0.4, 'head_sld': 3.0, 'solvent_sld': 6.0}, [0.001], [653143.9209]]]
Index: smodels/models/lamellarPC.py
===================================================================
--- sasmodels/models/lamellarPC.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ 	(revision )
@@ -1,153 +1,0 @@
-# Note: model title and parameter table are inserted automatically
-r"""
-This model calculates the scattering from a stack of repeating lamellar
-structures. The stacks of lamellae (infinite in lateral dimension) are
-treated as a paracrystal to account for the repeating spacing. The repeat
-distance is further characterized by a Gaussian polydispersity. **This model
-can be used for large multilamellar vesicles.**
-
-Definition
-----------
-
-In the equations below,
-
-- *scale* is used instead of the mass per area of the bilayer $\Gamma_m$
-  (this corresponds to the volume fraction of the material in the bilayer,
-  *not* the total excluded volume of the paracrystal),
-
-- *sld* $-$ *solvent_sld* is the contrast $\Delta \rho$,
-
-- *thickness* is the layer thickness $t$,
-
-- *Nlayers* is the number of layers $N$,
-
-- *spacing* is the average distance between adjacent layers
-  $\langle D \rangle$, and
-
-- *spacing_polydisp* is the relative standard deviation of the Gaussian
-  layer distance distribution $\sigma_D / \langle D \rangle$.
-
-
-The scattering intensity $I(q)$ is calculated as
-
-.. math::
-
-    I(q) = 2\pi\Delta\rho^2\Gamma_m\frac{P_\text{bil}(q)}{q^2} Z_N(q)
-
-The form factor of the bilayer is approximated as the cross section of an
-infinite, planar bilayer of thickness $t$
-
-.. math::
-
-    P_\text{bil}(q) = \left(\frac{\sin(qt/2)}{qt/2}\right)^2
-
-$Z_N(q)$ describes the interference effects for aggregates
-consisting of more than one bilayer. The equations used are (3-5)
-from the Bergstrom reference:
-
-.. math::
-
-
-    Z_N(q) = \frac{1 - w^2}{1 + w^2 - 2w \cos(q \langle D \rangle)}
-        + x_N S_N + (1 - x_N) S_{N+1}
-
-where
-
-.. math::
-
-    S_N(q) = \frac{a_N}{N}[1 + w^2 - 2 w \cos(q \langle D \rangle)]^2
-
-and
-
-.. math::
-
-    a_N &= 4w^2 - 2(w^3 + w) \cos(q \langle D \rangle) \\
-        &\quad - 4w^{N+2}\cos(Nq \langle D \rangle)
-        + 2 w^{N+3}\cos[(N-1)q \langle D \rangle]
-        + 2w^{N+1}\cos[(N+1)q \langle D \rangle]
-
-for the layer spacing distribution $w = \exp(-\sigma_D^2 q^2/2)$.
-
-Non-integer numbers of stacks are calculated as a linear combination of
-the lower and higher values
-
-.. math::
-
-    N_L = x_N N + (1 - x_N)(N+1)
-
-The 2D scattering intensity is the same as 1D, regardless of the orientation
-of the $q$ vector which is defined as
-
-.. math::
-
-    q = \sqrt{q_x^2 + q_y^2}
-
-
-Reference
----------
-
-M Bergstrom, J S Pedersen, P Schurtenberger, S U Egelhaaf,
-*J. Phys. Chem. B*, 103 (1999) 9888-9897
-
-"""
-
-from numpy import inf
-
-name = "lamellarPC"
-title = "Random lamellar sheet with paracrystal structure factor"
-description = """\
-    [Random lamellar phase with paracrystal structure factor]
-        randomly oriented stacks of infinite sheets
-        with paracrytal S(Q), having polydisperse spacing.
-        sld = sheet scattering length density
-        sld_solvent = solvent scattering length density
-        background = incoherent background
-        scale = scale factor
-"""
-category = "shape:lamellae"
-
-single = False
-
-#             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["thickness", "Ang", 33.0, [0, inf], "volume",
-               "sheet thickness"],
-              ["Nlayers", "", 20, [0, inf], "",
-               "Number of layers"],
-              ["spacing", "Ang", 250., [0.0, inf], "",
-               "d-spacing of paracrystal stack"],
-              ["spacing_polydisp", "Ang", 0.0, [0.0, inf], "",
-               "d-spacing polydispersity"],
-              ["sld", "1e-6/Ang^2", 1.0, [-inf, inf], "",
-               "layer scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2", 6.34, [-inf, inf], "",
-               "Solvent scattering length density"],
-             ]
-
-
-source = ["lamellarPC_kernel.c"]
-
-form_volume = """
-    return 1.0;
-    """
-
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
-
-# ER defaults to 0.0
-# VR defaults to 1.0
-
-demo = dict(scale=1, background=0,
-            thickness=33, Nlayers=20, spacing=250, spacing_polydisp=0.2,
-            sld=1.0, solvent_sld=6.34,
-            thickness_pd=0.2, thickness_pd_n=40)
-
-oldname = 'LamellarPCrystalModel'
-oldpars = dict(spacing_polydisp='pd_spacing', sld='sld_layer',
-               solvent_sld='sld_solvent')
-#
-tests = [
-        [ {'scale': 1.0, 'background' : 0.0, 'thickness' : 33.,'Nlayers' : 20.0, 'spacing' : 250., 'spacing_polydisp' : 0.2,
-            'sld' : 1.0, 'solvent_sld' : 6.34,
-            'thickness_pd' : 0.0, 'thickness_pd_n' : 40 }, [0.001, 0.215268], [21829.3, 0.00487686]]
-        ]
Index: smodels/models/lamellarPC_kernel.c
===================================================================
--- sasmodels/models/lamellarPC_kernel.c	(revision 319ab1478261080f41492caa7ea3b5e7f58409c8)
+++ 	(revision )
@@ -1,86 +1,0 @@
-/*	Lamellar_ParaCrystal - Pedersen's model
-
-*/
-double Iq(double qval,
-      double th,
-      double Nlayers, 
-	  double davg, 
-	  double pd,
-      double sld,
-      double solvent_sld);
-double paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an);
-double paraCryst_an(double ww, double qval, double davg, long Nlayers);
-
-double Iq(double qval,
-      double th,
-      double Nlayers, 
-	  double davg, 
-	  double pd,
-      double sld,
-      double solvent_sld)
-{
-    
-	double inten,contr,xn;
-	double xi,ww,Pbil,Znq,Snq,an;
-	long n1,n2;
-	
-	contr = sld - solvent_sld;
-	//get the fractional part of Nlayers, to determine the "mixing" of N's
-	
-	n1 = (long)trunc(Nlayers);		//rounds towards zero
-	n2 = n1 + 1;
-	xn = (double)n2 - Nlayers;			//fractional contribution of n1
-	
-	ww = exp(-qval*qval*pd*pd*davg*davg/2.0);
-
-	//calculate the n1 contribution
-	an = paraCryst_an(ww,qval,davg,n1);
-	Snq = paraCryst_sn(ww,qval,davg,n1,an);
-	
-	Znq = xn*Snq;
-	
-	//calculate the n2 contribution
-	an = paraCryst_an(ww,qval,davg,n2);
-	Snq = paraCryst_sn(ww,qval,davg,n2,an);
-
-	Znq += (1.0-xn)*Snq;
-	
-	//and the independent contribution
-	Znq += (1.0-ww*ww)/(1.0+ww*ww-2.0*ww*cos(qval*davg));
-	
-	//the limit when Nlayers approaches infinity
-//	Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg))
-	
-	xi = th/2.0;		//use 1/2 the bilayer thickness
-	Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi));
-	
-	inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval);
-	inten *= 1.0e-04;
-//printf("q=%.7e wwm1=%g ww=%.5e an=% 12.5e Snq=% 12.5e Znq=% 12.5e Pbil=% 12.5e\n",qval,wwm1,ww,an,Snq,Znq,Pbil);
-	return(inten);
-}
-
-// functions for the lamellar paracrystal model
-double
-paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an) {
-	
-	double Snq;
-
-	Snq = an/( (double)Nlayers*pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) );
-	
-	return(Snq);
-}
-
-double
-paraCryst_an(double ww, double qval, double davg, long Nlayers) {
-	
-	double an;
-	
-	an = 4.0*ww*ww - 2.0*(ww*ww*ww+ww)*cos(qval*davg);
-	an -= 4.0*pow(ww,(Nlayers+2))*cos((double)Nlayers*qval*davg);
-	an += 2.0*pow(ww,(Nlayers+3))*cos((double)(Nlayers-1)*qval*davg);
-	an += 2.0*pow(ww,(Nlayers+1))*cos((double)(Nlayers+1)*qval*davg);
-	
-	return(an);
-}
-
Index: sasmodels/models/lamellar_hg.py
===================================================================
--- sasmodels/models/lamellar_hg.py	(revision 52a3db3cb286d26bfdebb3ea92ad54cd92a07c7f)
+++ sasmodels/models/lamellar_hg.py	(revision 52a3db3cb286d26bfdebb3ea92ad54cd92a07c7f)
@@ -0,0 +1,122 @@
+# Note: model title and parameter table are inserted automatically
+r"""
+This model provides the scattering intensity, $I(q)$, for a lyotropic lamellar
+phase where a random distribution in solution are assumed. The SLD of the head
+region is taken to be different from the SLD of the tail region.
+
+Definition
+----------
+
+The scattering intensity $I(q)$ is
+
+.. math::
+
+   I(q) = 2\pi\frac{\text{scale}}{2(\delta_H + \delta_T)}  P(q) \frac{1}{q^2}
+
+The form factor $P(q)$ is
+
+.. math::
+
+    P(q) = \frac{4}{q^2}
+        \left\lbrace
+            \Delta \rho_H
+            \left[\sin[q(\delta_H + \delta_T)\ - \sin(q\delta_T)\right]
+            + \Delta\rho_T\sin(q\delta_T)
+        \right\rbrace^2
+
+where $\delta_T$ is *tail_length*, $\delta_H$ is *head_length*,
+$\Delta\rho_H$ is the head contrast (*sld_head* $-$ *sld_solvent*),
+and $\Delta\rho_T$ is tail contrast (*sld* $-$ *sld_solvent*).
+
+The total thickness of the lamellar sheet is $\delta_H$ + $\delta_T$ + $\delta_T$ + $\delta_H$.
+Note that in a non aqueous solvent the chemical "head" group may be the 
+"Tail region" and vice-versa.
+
+The 2D scattering intensity is calculated in the same way as 1D, where
+the $q$ vector is defined as
+
+.. math::
+
+    q = \sqrt{q_x^2 + q_y^2}
+
+
+References
+----------
+
+F Nallet, R Laversanne, and D Roux, J. Phys. II France, 3, (1993) 487-502
+
+also in J. Phys. Chem. B, 105, (2001) 11081-11088
+
+*2014/04/17 - Description reviewed by S King and P Butler.*
+"""
+
+from numpy import inf
+
+name = "lamellar_hg"
+title = "Random lamellar phase with Head and Tail Groups"
+description = """\
+    [Random lamellar phase with Head and Tail Groups]
+        I(q)= 2*pi*P(q)/(2(H+T)*q^(2)), where
+        P(q)= see manual
+        layer thickness =(H+T+T+H) = 2(Head+Tail)
+        sld = Tail scattering length density
+        sld_head = Head scattering length density
+        sld_solvent = solvent scattering length density
+        background = incoherent background
+        scale = scale factor
+"""
+category = "shape:lamellae"
+
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["tail_length", "Ang",       15,   [0, inf],  "volume",  "Tail thickness ( total = H+T+T+H)"],
+              ["head_length", "Ang",       10,   [0, inf],  "volume",  "Head thickness"],
+              ["sld",         "1e-6/Ang^2", 0.4, [-inf,inf], "",       "Tail scattering length density"],
+              ["sld_head",    "1e-6/Ang^2", 3.0, [-inf,inf], "",       "Head scattering length density"],
+              ["sld_solvent", "1e-6/Ang^2", 6,   [-inf,inf], "",       "Solvent scattering length density"]]
+# pylint: enable=bad-whitespace, line-too-long
+
+# No volume normalization despite having a volume parameter
+# This should perhaps be volume normalized?
+form_volume = """
+    return 1.0;
+    """
+
+Iq = """
+    const double qsq = q*q;
+    const double drh = sld_head - sld_solvent;
+    const double drt = sld - sld_solvent;    //correction 13FEB06 by L.Porcar
+    const double qT = q*tail_length;
+    double Pq, inten;
+    Pq = drh*(sin(q*(head_length+tail_length))-sin(qT)) + drt*sin(qT);
+    Pq *= Pq;
+    Pq *= 4.0/(qsq);
+
+    inten = 2.0e-4*M_PI*Pq/qsq;
+
+    // normalize by the bilayer thickness
+    inten /= 2.0*(head_length+tail_length);
+
+    return inten;
+    """
+
+Iqxy = """
+    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
+    """
+
+# ER defaults to 0.0
+# VR defaults to 1.0
+
+demo = dict(scale=1, background=0,
+            tail_length=15, head_length=10,
+            sld=0.4, sld_head=3.0, sld_solvent=6.0,
+            tail_length_pd=0.2, tail_length_pd_n=40,
+            head_length_pd=0.01, head_length_pd_n=40)
+
+oldname = 'LamellarFFHGModel'
+oldpars = dict(head_length='h_thickness', tail_length='t_length',
+               sld='sld_tail', sld_head='sld_head', sld_solvent='sld_solvent')
+#
+tests = [[{'scale': 1.0, 'background': 0.0, 'tail_length': 15.0, 'head_length': 10.0,
+           'sld': 0.4, 'sld_head': 3.0, 'sld_solvent': 6.0}, [0.001], [653143.9209]]]
+# ADDED by: RKH  ON: 18Mar2016  converted from sasview previously, now renaming everything & sorting the docs
Index: sasmodels/models/lamellar_hg_stack_caille.py
===================================================================
--- sasmodels/models/lamellar_hg_stack_caille.py	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
+++ sasmodels/models/lamellar_hg_stack_caille.py	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
@@ -0,0 +1,151 @@
+# Note: model title and parameter table are inserted automatically
+r"""
+This model provides the scattering intensity, $I(q) = P(q)S(q)$, for a lamellar
+phase where a random distribution in solution are assumed. Here a Caille $S(q)$
+is used for the lamellar stacks.
+
+The scattering intensity $I(q)$ is
+
+.. math::
+
+    I(q) = 2 \pi \frac{P(q)S(q)}{q^2\delta }
+
+
+The form factor $P(q)$ is
+
+.. math::
+
+        P(q) = \frac{4}{q^2}\big\{
+        \Delta\rho_H \left[\sin[q(\delta_H + \delta_T)] - \sin(q\delta_T)\right]
+            + \Delta\rho_T\sin(q\delta_T)\big\}^2
+
+and the structure factor $S(q)$ is
+
+.. math::
+
+    S(q) = 1 + 2 \sum_1^{N-1}\left(1-\frac{n}{N}\right)
+        \cos(qdn)\exp\left(-\frac{2q^2d^2\alpha(n)}{2}\right)
+
+where
+
+.. math::
+    :nowrap:
+
+    \begin{align*}
+    \alpha(n) &= \frac{\eta_{cp}}{4\pi^2} \left(\ln(\pi n)+\gamma_E\right)
+              &&  \\
+    \gamma_E  &= 0.5772156649
+              && \text{Euler's constant} \\
+    \eta_{cp} &= \frac{q_o^2k_B T}{8\pi\sqrt{K\overline{B}}}
+              && \text{Caille constant}
+    \end{align*}
+
+
+$\delta_T$ is the tail length (or *tail_length*), $\delta_H$ is the head
+thickness (or *head_length*), $\Delta\rho_H$ is SLD(headgroup) - SLD(solvent),
+and $\Delta\rho_T$ is SLD(tail) - SLD(headgroup). Here $d$ is (repeat) spacing,
+$K$ is smectic bending elasticity, $B$ is compression modulus, and $N$ is the
+number of lamellar plates (*Nlayers*).
+
+NB: **When the Caille parameter is greater than approximately 0.8 to 1.0, the
+assumptions of the model are incorrect.**  And due to a complication of the
+model function, users are responsible for making sure that all the assumptions
+are handled accurately (see the original reference below for more details).
+
+Non-integer numbers of stacks are calculated as a linear combination of
+results for the next lower and higher values.
+
+Be aware that the computations may be very slow.
+
+The 2D scattering intensity is calculated in the same way as 1D, where
+the $q$ vector is defined as
+
+.. math::
+
+    q = \sqrt{q_x^2 + q_y^2}
+
+
+References
+----------
+
+F Nallet, R Laversanne, and D Roux, J. Phys. II France, 3, (1993) 487-502
+
+also in J. Phys. Chem. B, 105, (2001) 11081-11088
+"""
+from numpy import inf
+
+name = "lamellar_hg_stack_caille"
+title = "Random lamellar head/tail/tail/head sheet with Caille structure factor"
+description = """\
+    [Random lamellar phase with Caille  structure factor]
+        randomly oriented stacks of infinite sheets
+        with Caille S(Q), having polydisperse spacing.
+        layer thickness =(H+T+T+H) = 2(Head+Tail)
+        sld = Tail scattering length density
+        sld_head = Head scattering length density
+        sld_solvent = solvent scattering length density
+        background = incoherent background
+        scale = scale factor
+"""
+category = "shape:lamellae"
+
+single = False
+parameters = [
+    #   [ "name", "units", default, [lower, upper], "type",
+    #     "description" ],
+    ["tail_length", "Ang", 10, [0, inf], "volume",
+     "Tail thickness"],
+    ["head_length", "Ang", 2, [0, inf], "volume",
+     "head thickness"],
+    ["Nlayers", "", 30, [0, inf], "",
+     "Number of layers"],
+    ["spacing", "Ang", 40., [0.0, inf], "volume",
+     "d-spacing of Caille S(Q)"],
+    ["Caille_parameter", "", 0.001, [0.0, 0.8], "",
+     "Caille parameter"],
+    ["sld", "1e-6/Ang^2", 0.4, [-inf, inf], "",
+     "Tail scattering length density"],
+    ["sld_head", "1e-6/Ang^2", 2.0, [-inf, inf], "",
+     "Head scattering length density"],
+    ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "",
+     "Solvent scattering length density"],
+    ]
+
+source = ["lamellar_hg_stack_caille_kernel.c"]
+
+# No volume normalization despite having a volume parameter
+# This should perhaps be volume normalized?
+form_volume = """
+    return 1.0;
+    """
+
+Iqxy = """
+    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
+    """
+
+# ER defaults to 0.0
+# VR defaults to 1.0
+
+demo = dict(
+    scale=1, background=0,
+    Nlayers=20, spacing=200., Caille_parameter=0.05,
+    tail_length=15, head_length=10,
+    #sld=-1, sld_head=4.0, sld_solvent=6.0,
+    sld=-1, sld_head=4.1, sld_solvent=6.0,
+    tail_length_pd=0.1, tail_length_pd_n=20,
+    head_length_pd=0.05, head_length_pd_n=30,
+    spacing_pd=0.2, spacing_pd_n=40,
+    )
+
+oldname = 'LamellarPSHGModel'
+
+oldpars = dict(
+    tail_length='deltaT', head_length='deltaH', Nlayers='n_plates',
+    Caille_parameter='caille', sld='sld_tail', sld_head='sld_head',
+    sld_solvent='sld_solvent')
+#
+tests = [[{'scale': 1.0, 'background': 0.0, 'tail_length': 10.0, 'head_length': 2.0,
+           'Nlayers': 30.0, 'spacing': 40., 'Caille_parameter': 0.001, 'sld': 0.4,
+           'sld_head': 2.0, 'sld_solvent': 6.0, 'tail_length_pd': 0.0,
+           'head_length_pd': 0.0, 'spacing_pd': 0.0}, [0.001], [6838238.571488]]]
+# ADDED by: RKH  ON: 18Mar2016  converted from sasview previously, now renaming everything & sorting the docs
Index: sasmodels/models/lamellar_hg_stack_caille_kernel.c
===================================================================
--- sasmodels/models/lamellar_hg_stack_caille_kernel.c	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
+++ sasmodels/models/lamellar_hg_stack_caille_kernel.c	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
@@ -0,0 +1,75 @@
+/*	LamellarCailleHG kernel - allows for name changes of passed parameters ...
+    Maths identical to LamellarCaille apart from the line for P(Q)
+*/
+
+double Iq(double qval,
+      double tail_length,
+      double head_length,
+      double Nlayers, 
+      double dd,
+      double Cp,
+      double tail_sld,
+      double head_sld,
+      double solvent_sld);
+
+double Iq(double qval,
+      double tail_length,
+      double head_length,
+      double Nlayers, 
+      double dd,
+      double Cp,
+      double tail_sld,
+      double head_sld,
+      double solvent_sld)
+{
+  double NN;   //local variables of coefficient wave
+  double inten,Pq,Sq,alpha,temp,t2;
+  //double dQ, dQDefault, t1, t3;
+  int ii,NNint;
+  // from wikipedia 0.577215664901532860606512090082402431042159335
+  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015
+  //dQDefault = 0.0;    //[=] 1/A, q-resolution, default value
+  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015
+
+  NN = trunc(Nlayers);    //be sure that NN is an integer
+  
+  Pq = (head_sld-solvent_sld)*(sin(qval*(head_length+tail_length))-sin(qval*tail_length)) +
+              (tail_sld-solvent_sld)*sin(qval*tail_length);
+  Pq *= Pq;
+  Pq *= 4.0/(qval*qval);
+
+  NNint = (int)NN;    //cast to an integer for the loop
+  ii=0;
+  Sq = 0.0;
+  for(ii=1;ii<=(NNint-1);ii+=1) {
+
+    //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015
+
+    temp = 0.0;
+    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler);
+    //t1 = 2.0*dQ*dQ*dd*dd*alpha;
+    t2 = 2.0*qval*qval*dd*dd*alpha;
+    //t3 = dQ*dQ*dd*dd*ii*ii;
+
+    temp = 1.0-ii/NN;
+    //temp *= cos(dd*qval*ii/(1.0+t1));
+    temp *= cos(dd*qval*ii);
+    //if (temp < 0) printf("q=%g: ii=%d, cos(dd*q*ii)=cos(%g) < 0\n",qval,ii,dd*qval*ii);
+    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
+    temp *= exp(-t2/2.0 );
+    //temp /= sqrt(1.0+t1);
+
+    Sq += temp;
+  }
+
+  Sq *= 2.0;
+  Sq += 1.0;
+
+  //if (Sq < 0) printf("q=%g: S(q) =%g\n", qval, Sq);
+
+  inten = 2.0*M_PI*Pq*Sq/(dd*qval*qval);
+
+  inten *= 1.0e-04;   // 1/A to 1/cm
+  return(inten);
+}
+
Index: sasmodels/models/lamellar_stack_caille.py
===================================================================
--- sasmodels/models/lamellar_stack_caille.py	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
+++ sasmodels/models/lamellar_stack_caille.py	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
@@ -0,0 +1,130 @@
+r"""
+This model provides the scattering intensity, $I(q) = P(q) S(q)$, for a
+lamellar phase where a random distribution in solution are assumed.
+Here a Caille $S(q)$ is used for the lamellar stacks.
+
+Definition
+----------
+
+The scattering intensity $I(q)$ is
+
+.. math::
+
+    I(q) = 2\pi \frac{P(q)S(q)}{q^2\delta }
+
+The form factor is
+
+.. math::
+
+    P(q) = \frac{2\Delta\rho^2}{q^2}\left(1-\cos q\delta \right)
+
+and the structure factor is
+
+.. math::
+
+    S(q) = 1 + 2 \sum_1^{N-1}\left(1-\frac{n}{N}\right)
+           \cos(qdn)\exp\left(-\frac{2q^2d^2\alpha(n)}{2}\right)
+
+where
+
+.. math::
+    :nowrap:
+
+    \begin{align*}
+    \alpha(n) &= \frac{\eta_{cp}}{4\pi^2} \left(\ln(\pi n)+\gamma_E\right)
+              && \\
+    \gamma_E  &= 0.5772156649
+              && \text{Euler's constant} \\
+    \eta_{cp} &= \frac{q_o^2k_B T}{8\pi\sqrt{K\overline{B}}}
+              && \text{Caille constant}
+    \end{align*}
+
+Here $d$ = (repeat) spacing, $\delta$ = bilayer thickness,
+the contrast $\Delta\rho$ = SLD(headgroup) - SLD(solvent),
+$K$ = smectic bending elasticity, $B$ = compression modulus, and
+$N$ = number of lamellar plates (*n_plates*).
+
+NB: **When the Caille parameter is greater than approximately 0.8 to 1.0, the
+assumptions of the model are incorrect.** And due to a complication of the
+model function, users are responsible for making sure that all the assumptions
+are handled accurately (see the original reference below for more details).
+
+Non-integer numbers of stacks are calculated as a linear combination of
+results for the next lower and higher values.
+
+The 2D scattering intensity is calculated in the same way as 1D, where the
+$q$ vector is defined as
+
+.. math::
+
+    q = \sqrt{q_x^2 + q_y^2}
+
+
+References
+----------
+
+F Nallet, R Laversanne, and D Roux, J. Phys. II France, 3, (1993) 487-502
+
+also in J. Phys. Chem. B, 105, (2001) 11081-11088
+"""
+from numpy import inf
+
+name = "lamellar_stack_caille"
+title = "Random lamellar sheet with Caille structure factor"
+description = """\
+    [Random lamellar phase with Caille  structure factor]
+    randomly oriented stacks of infinite sheets
+    with Caille S(Q), having polydisperse spacing.
+    sld = sheet scattering length density
+    sld_solvent = solvent scattering length density
+    background = incoherent background
+    scale = scale factor
+"""
+category = "shape:lamellae"
+
+single = False
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [
+    ["thickness",        "Ang",      30.0,  [0, inf],   "volume", "sheet thickness"],
+    ["Nlayers",          "",          20,   [0, inf],   "",       "Number of layers"],
+    ["spacing",          "Ang",      400.,  [0.0,inf],  "volume", "d-spacing of Caille S(Q)"],
+    ["Caille_parameter", "1/Ang^2",    0.1, [0.0,0.8],  "",       "Caille parameter"],
+    ["sld",              "1e-6/Ang^2", 6.3, [-inf,inf], "",       "layer scattering length density"],
+    ["sld_solvent",      "1e-6/Ang^2", 1.0, [-inf,inf], "",       "Solvent scattering length density"],
+    ]
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lamellar_stack_caille_kernel.c"]
+
+# No volume normalization despite having a volume parameter
+# This should perhaps be volume normalized?
+form_volume = """
+    return 1.0;
+    """
+
+Iqxy = """
+    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
+    """
+
+# ER defaults to 0.0
+# VR defaults to 1.0
+
+demo = dict(scale=1, background=0,
+            thickness=67., Nlayers=3.75, spacing=200.,
+            Caille_parameter=0.268, sld=1.0, sld_solvent=6.34,
+            thickness_pd=0.1, thickness_pd_n=100,
+            spacing_pd=0.05, spacing_pd_n=40)
+
+oldname = 'LamellarPSModel'
+oldpars = dict(thickness='delta', Nlayers='N_plates',
+               Caille_parameter='caille',
+               sld='sld_bi', sld_solvent='sld_sol')
+#
+tests = [
+    [{'scale': 1.0, 'background': 0.0, 'thickness': 30., 'Nlayers': 20.0,
+      'spacing': 400., 'Caille_parameter': 0.1, 'sld': 6.3,
+      'sld_solvent': 1.0, 'thickness_pd': 0.0, 'spacing_pd': 0.0},
+     [0.001], [28895.13397]]
+    ]
+# ADDED by: RKH  ON: 18Mar2016  converted from sasview previously, now renaming everything & sorting the docs
Index: sasmodels/models/lamellar_stack_caille_kernel.c
===================================================================
--- sasmodels/models/lamellar_stack_caille_kernel.c	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
+++ sasmodels/models/lamellar_stack_caille_kernel.c	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
@@ -0,0 +1,69 @@
+/*	LamellarCaille kernel - allows for name changes of passed parameters ...
+
+*/
+
+double Iq(double qval,
+      double del,
+      double Nlayers, 
+      double dd,
+	  double Cp, 
+      double sld,
+      double solvent_sld);
+
+double Iq(double qval,
+      double del,
+      double Nlayers, 
+      double dd,
+	  double Cp, 
+      double sld,
+      double solvent_sld)
+{
+  double contr,NN;   //local variables of coefficient wave
+  double inten,Pq,Sq,alpha,temp,t2;
+  //double dQ, dQDefault, t1, t3;
+  int ii,NNint;
+  // from wikipedia 0.577215664901532860606512090082402431042159335
+  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015
+  //dQDefault = 0.0;    //[=] 1/A, q-resolution, default value
+  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015
+
+  NN = trunc(Nlayers);    //be sure that NN is an integer
+  
+  contr = sld - solvent_sld;
+
+  Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del));
+
+  NNint = (int)NN;    //cast to an integer for the loop
+  ii=0;
+  Sq = 0.0;
+  // the vital "=" in ii<=  added March 2015
+  for(ii=1;ii<=(NNint-1);ii+=1) {
+
+    //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015
+
+    temp = 0.0;
+    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler);
+    //t1 = 2.0*dQ*dQ*dd*dd*alpha;
+    t2 = 2.0*qval*qval*dd*dd*alpha;
+    //t3 = dQ*dQ*dd*dd*ii*ii;
+
+    temp = 1.0-ii/NN;
+    //temp *= cos(dd*qval*ii/(1.0+t1));
+    temp *= cos(dd*qval*ii);
+    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
+    temp *= exp(-t2/2.0 );
+    //temp /= sqrt(1.0+t1);
+
+    Sq += temp;
+  }
+
+  Sq *= 2.0;
+  Sq += 1.0;
+
+  inten = 2.0*M_PI*Pq*Sq/(dd*qval*qval);
+
+  inten *= 1.0e-04;   // 1/A to 1/cm
+
+  return(inten);
+}
+
Index: sasmodels/models/lamellar_stack_paracrystal.py
===================================================================
--- sasmodels/models/lamellar_stack_paracrystal.py	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
+++ sasmodels/models/lamellar_stack_paracrystal.py	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
@@ -0,0 +1,155 @@
+# Note: model title and parameter table are inserted automatically
+r"""
+This model calculates the scattering from a stack of repeating lamellar
+structures. The stacks of lamellae (infinite in lateral dimension) are
+treated as a paracrystal to account for the repeating spacing. The repeat
+distance is further characterized by a Gaussian polydispersity. **This model
+can be used for large multilamellar vesicles.**
+
+Definition
+----------
+
+In the equations below,
+
+- *scale* is used instead of the mass per area of the bilayer $\Gamma_m$
+  (this corresponds to the volume fraction of the material in the bilayer,
+  *not* the total excluded volume of the paracrystal),
+
+- *sld* $-$ *sld_solvent* is the contrast $\Delta \rho$,
+
+- *thickness* is the layer thickness $t$,
+
+- *Nlayers* is the number of layers $N$,
+
+- *spacing* is the average distance between adjacent layers
+  $\langle D \rangle$, and
+
+- *spacing_polydisp* is the relative standard deviation of the Gaussian
+  layer distance distribution $\sigma_D / \langle D \rangle$.
+
+
+The scattering intensity $I(q)$ is calculated as
+
+.. math::
+
+    I(q) = 2\pi\Delta\rho^2\Gamma_m\frac{P_\text{bil}(q)}{q^2} Z_N(q)
+
+The form factor of the bilayer is approximated as the cross section of an
+infinite, planar bilayer of thickness $t$ (compare the equations for the
+lamellar model).
+
+.. math::
+
+    P_\text{bil}(q) = \left(\frac{\sin(qt/2)}{qt/2}\right)^2
+
+$Z_N(q)$ describes the interference effects for aggregates
+consisting of more than one bilayer. The equations used are (3-5)
+from the Bergstrom reference:
+
+.. math::
+
+
+    Z_N(q) = \frac{1 - w^2}{1 + w^2 - 2w \cos(q \langle D \rangle)}
+        + x_N S_N + (1 - x_N) S_{N+1}
+
+where
+
+.. math::
+
+    S_N(q) = \frac{a_N}{N}[1 + w^2 - 2 w \cos(q \langle D \rangle)]^2
+
+and
+
+.. math::
+
+    a_N &= 4w^2 - 2(w^3 + w) \cos(q \langle D \rangle) \\
+        &\quad - 4w^{N+2}\cos(Nq \langle D \rangle)
+        + 2 w^{N+3}\cos[(N-1)q \langle D \rangle]
+        + 2w^{N+1}\cos[(N+1)q \langle D \rangle]
+
+for the layer spacing distribution $w = \exp(-\sigma_D^2 q^2/2)$.
+
+Non-integer numbers of stacks are calculated as a linear combination of
+the lower and higher values
+
+.. math::
+
+    N_L = x_N N + (1 - x_N)(N+1)
+
+The 2D scattering intensity is the same as 1D, regardless of the orientation
+of the $q$ vector which is defined as
+
+.. math::
+
+    q = \sqrt{q_x^2 + q_y^2}
+
+
+Reference
+---------
+
+M Bergstrom, J S Pedersen, P Schurtenberger, S U Egelhaaf,
+*J. Phys. Chem. B*, 103 (1999) 9888-9897
+
+"""
+
+from numpy import inf
+
+name = "lamellar_stack_paracrystal"
+title = "Random lamellar sheet with paracrystal structure factor"
+description = """\
+    [Random lamellar phase with paracrystal structure factor]
+        randomly oriented stacks of infinite sheets
+        with paracrytal S(Q), having polydisperse spacing.
+        sld = sheet scattering length density
+        sld_solvent = solvent scattering length density
+        background = incoherent background
+        scale = scale factor
+"""
+category = "shape:lamellae"
+
+single = False
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["thickness", "Ang", 33.0, [0, inf], "volume",
+               "sheet thickness"],
+              ["Nlayers", "", 20, [0, inf], "",
+               "Number of layers"],
+              ["spacing", "Ang", 250., [0.0, inf], "",
+               "d-spacing of paracrystal stack"],
+              ["spacing_polydisp", "Ang", 0.0, [0.0, inf], "",
+               "d-spacing polydispersity"],
+              ["sld", "1e-6/Ang^2", 1.0, [-inf, inf], "",
+               "layer scattering length density"],
+              ["sld_solvent", "1e-6/Ang^2", 6.34, [-inf, inf], "",
+               "Solvent scattering length density"],
+             ]
+
+
+source = ["lamellar_stack_paracrystal_kernel.c"]
+
+form_volume = """
+    return 1.0;
+    """
+
+Iqxy = """
+    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
+    """
+
+# ER defaults to 0.0
+# VR defaults to 1.0
+
+demo = dict(scale=1, background=0,
+            thickness=33, Nlayers=20, spacing=250, spacing_polydisp=0.2,
+            sld=1.0, sld_solvent=6.34,
+            thickness_pd=0.2, thickness_pd_n=40)
+
+oldname = 'LamellarPCrystalModel'
+oldpars = dict(spacing_polydisp='pd_spacing', sld='sld_layer',
+               sld_solvent='sld_solvent')
+#
+tests = [
+        [ {'scale': 1.0, 'background' : 0.0, 'thickness' : 33.,'Nlayers' : 20.0, 'spacing' : 250., 'spacing_polydisp' : 0.2,
+            'sld' : 1.0, 'sld_solvent' : 6.34,
+            'thickness_pd' : 0.0, 'thickness_pd_n' : 40 }, [0.001, 0.215268], [21829.3, 0.00487686]]
+        ]
+# ADDED by: RKH  ON: 18Mar2016  converted from sasview previously, now renaming everything & sorting the docs
Index: sasmodels/models/lamellar_stack_paracrystal_kernel.c
===================================================================
--- sasmodels/models/lamellar_stack_paracrystal_kernel.c	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
+++ sasmodels/models/lamellar_stack_paracrystal_kernel.c	(revision 6ab4ed8677894e77dfd90dc3fdd5bdac2df2eff8)
@@ -0,0 +1,86 @@
+/*	Lamellar_ParaCrystal - Pedersen's model
+
+*/
+double Iq(double qval,
+      double th,
+      double Nlayers, 
+	  double davg, 
+	  double pd,
+      double sld,
+      double solvent_sld);
+double paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an);
+double paraCryst_an(double ww, double qval, double davg, long Nlayers);
+
+double Iq(double qval,
+      double th,
+      double Nlayers, 
+	  double davg, 
+	  double pd,
+      double sld,
+      double solvent_sld)
+{
+    
+	double inten,contr,xn;
+	double xi,ww,Pbil,Znq,Snq,an;
+	long n1,n2;
+	
+	contr = sld - solvent_sld;
+	//get the fractional part of Nlayers, to determine the "mixing" of N's
+	
+	n1 = (long)trunc(Nlayers);		//rounds towards zero
+	n2 = n1 + 1;
+	xn = (double)n2 - Nlayers;			//fractional contribution of n1
+	
+	ww = exp(-qval*qval*pd*pd*davg*davg/2.0);
+
+	//calculate the n1 contribution
+	an = paraCryst_an(ww,qval,davg,n1);
+	Snq = paraCryst_sn(ww,qval,davg,n1,an);
+	
+	Znq = xn*Snq;
+	
+	//calculate the n2 contribution
+	an = paraCryst_an(ww,qval,davg,n2);
+	Snq = paraCryst_sn(ww,qval,davg,n2,an);
+
+	Znq += (1.0-xn)*Snq;
+	
+	//and the independent contribution
+	Znq += (1.0-ww*ww)/(1.0+ww*ww-2.0*ww*cos(qval*davg));
+	
+	//the limit when Nlayers approaches infinity
+//	Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg))
+	
+	xi = th/2.0;		//use 1/2 the bilayer thickness
+	Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi));
+	
+	inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval);
+	inten *= 1.0e-04;
+//printf("q=%.7e wwm1=%g ww=%.5e an=% 12.5e Snq=% 12.5e Znq=% 12.5e Pbil=% 12.5e\n",qval,wwm1,ww,an,Snq,Znq,Pbil);
+	return(inten);
+}
+
+// functions for the lamellar paracrystal model
+double
+paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an) {
+	
+	double Snq;
+
+	Snq = an/( (double)Nlayers*pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) );
+	
+	return(Snq);
+}
+
+double
+paraCryst_an(double ww, double qval, double davg, long Nlayers) {
+	
+	double an;
+	
+	an = 4.0*ww*ww - 2.0*(ww*ww*ww+ww)*cos(qval*davg);
+	an -= 4.0*pow(ww,(Nlayers+2))*cos((double)Nlayers*qval*davg);
+	an += 2.0*pow(ww,(Nlayers+3))*cos((double)(Nlayers-1)*qval*davg);
+	an += 2.0*pow(ww,(Nlayers+1))*cos((double)(Nlayers+1)*qval*davg);
+	
+	return(an);
+}
+
Index: sasmodels/product.py
===================================================================
--- sasmodels/product.py	(revision 3bcd03d7e1d8d4938d30f17bf360b83417139376)
+++ sasmodels/product.py	(revision d5ba8419b0a797faa3fe41fde9334c5a19282992)
@@ -124,13 +124,15 @@
         # so borrow values from end of p_fixed.  This makes volfraction the
         # first S parameter.
-        start += num_p_fixed - 2
-        par_map['s_fixed'] = np.arange(start, start+num_s_fixed)
+        start += num_p_fixed
+        par_map['s_fixed'] = np.hstack(([start,start],
+                                        np.arange(start, start+num_s_fixed-2)))
         par_map['volfraction'] = num_p_fixed
-        start += num_s_fixed
+        start += num_s_fixed-2
         # vol pars offset from the start of pd pars
         par_map['vol_pars'] = [start+k for k in vol_par_idx]
         par_map['p_pd'] = np.arange(start, start+num_p_pd)
-        start += num_p_pd
-        par_map['s_pd'] = np.arange(start-1, start+num_s_pd)  # should be empty...
+        start += num_p_pd-1
+        par_map['s_pd'] = np.hstack((start,
+                                     np.arange(start, start+num_s_pd-1)))
 
         self.fixed_pars = model_info['partype']['fixed-' + dim]
Index: sasmodels/resolution.py
===================================================================
--- sasmodels/resolution.py	(revision 17bbadde9de73711781f748816b06e6956dc8cfa)
+++ sasmodels/resolution.py	(revision a146eaab015a1387f12891d14d24aa44e466c81c)
@@ -197,5 +197,6 @@
 
 
-    **Algorithm**
+    Definition
+    ----------
 
     We are using the mid-point integration rule to assign weights to each
@@ -441,5 +442,5 @@
     .. math::
 
-        n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n)
+         n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n)
             / (\log q_n - log q_1)
     """
Index: sasmodels/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision fcd7bbda85a603b85101ba7ec6495729497e6b3f)
+++ sasmodels/sasview_model.py	(revision 2622b3fe39474acb34f0bbd16dccdce0561d8a26)
@@ -17,4 +17,5 @@
 from copy import deepcopy
 import warnings
+import collections
 
 import numpy as np
@@ -57,7 +58,8 @@
         ## TODO: reorganize parameter handling
         self.details = dict()
-        self.params = dict()
+        self.params = collections.OrderedDict()
         self.dispersion = dict()
         partype = model.info['partype']
+
         for p in model.info['parameters']:
             self.params[p.name] = p.default
Index: sasmodels/sesans.py
===================================================================
--- sasmodels/sesans.py	(revision 190fc2b167d3d89ca9c3b7d3128bdf8bde9f67d0)
+++ sasmodels/sesans.py	(revision a06430cd1ac1662ddeef978e2f7b97f0579f0945)
@@ -13,5 +13,6 @@
 from numpy import pi, exp
 from scipy.special import jv as besselj
-
+#import direct_model.DataMixin as model
+        
 def make_q(q_max, Rmax):
     r"""
@@ -21,5 +22,121 @@
     q_min = dq = 0.1 * 2*pi / Rmax
     return np.arange(q_min, q_max, dq)
+    
+def make_all_q(data):
+    if not data.needs_all_q:
+        return []
+    elif needs_Iqxy(data):
+        # compute qx, qy
+        Qx, Qy = np.meshgrid(qx, qy)
+        return [Qx, Qy]
+    else:
+        # else only need q
+        return [q]
 
+def transform(data, q_calc, Iq_calc, qmono, Iq_mono):
+    nqmono = len(qmono)
+    if nqmono == 0:
+        result = call_hankel(data, q_calc, Iq_calc)
+    elif nqmono == 1:
+        q = qmono[0]
+        result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono)
+    else:
+        Qx, Qy = [qmono[0], qmono[1]]
+        Qx = np.reshape(Qx, nqx, nqy)
+        Qy = np.reshape(Qy, nqx, nqy)
+        Iq_mono = np.reshape(Iq_mono, nqx, nqy)
+        qx = Qx[0, :]
+        qy = Qy[:, 0]
+        result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono)
+
+    return result
+
+def call_hankel(data, q_calc, Iq_calc):
+    return hankel(data.x, data.lam * 1e-9,
+                  data.sample.thickness / 10,
+                  q_calc, Iq_calc)
+  
+def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono):
+    return hankel(data.x, data.lam * 1e-9,
+                  data.sample.thickness / 10,
+                  q_calc, Iq_calc)
+                  
+def Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono):
+    return hankel(data.x, data.y, data.lam * 1e-9,
+                  data.sample.thickness / 10,
+                  q_calc, Iq_calc)
+                        
+def TotalScatter(model, parameters):  #Work in progress!!
+#    Calls a model with existing model parameters already in place, then integrate the product of q and I(q) from 0 to (4*pi/lambda)
+    allq = np.linspace(0,4*pi/wavelength,1000)
+    allIq = 1
+    integral = allq*allIq
+    
+
+
+def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still
+#==============================================================================
+#     2D Cosine Transform if "wavelength" is a vector
+#==============================================================================
+#allq is the q-space needed to create the total scattering cross-section
+
+    Gprime = np.zeros_like(wavelength, 'd')
+    s = np.zeros_like(wavelength, 'd')
+    sd = np.zeros_like(wavelength, 'd')
+    Gprime = np.zeros_like(wavelength, 'd')
+    f = np.zeros_like(wavelength, 'd')
+    for i, wavelength_i in enumerate(wavelength):
+        z = magfield*wavelength_i
+        allq=np.linspace() #for calculating the Q-range of the  scattering power integral
+        allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow
+        alldq = (allq[1]-allq[0])*1e10
+        sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq)
+        s[i]=1-exp(-sigma)
+        for j, Iqy_j, qy_j in enumerate(qy):
+            for k, Iqz_k, qz_k in enumerate(qz):
+                Iq = np.sqrt(Iqy_j^2+Iqz_k^2)
+                q = np.sqrt(qy_j^2 + qz_k^2)
+                Gintegral = Iq*cos(z*Qz_k)
+                Gprime[i] += Gintegral
+#                sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i]
+#                s[i] += 1-exp(Totalscatter(modelname)*thickness)
+#                For now, work with standard 2-phase scatter
+
+
+                sd[i] += Iq
+        f[i] = 1-s[i]+sd[i]
+        P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i]
+
+
+
+
+def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname):
+#==============================================================================
+#     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS
+#==============================================================================
+#acceptq is the q-space needed to create limited acceptance effect
+    SElength= wavelength*magfield
+    G = np.zeros_like(SElength, 'd')
+    threshold=2*pi*theta/wavelength
+    for i, SElength_i in enumerate(SElength):
+        allq=np.linspace() #for calculating the Q-range of the  scattering power integral
+        allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow
+        alldq = (allq[1]-allq[0])*1e10
+        sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq)
+        s[i]=1-exp(-sigma)
+
+        dq = (q[1]-q[0])*1e10
+        a = (x<threshold)
+        acceptq = a*q
+        acceptIq = a*Iq
+
+        G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq)
+
+#        G[i]=np.sum(integral)
+
+    G *= dq*1e10*2*pi
+
+    P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0]))
+    
 def hankel(SElength, wavelength, thickness, q, Iq):
     r"""
@@ -44,4 +161,7 @@
     """
     G = np.zeros_like(SElength, 'd')
+#==============================================================================
+#     Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS
+#==============================================================================
     for i, SElength_i in enumerate(SElength):
         integral = besselj(0, q*SElength_i)*Iq*q
Index: setup.py
===================================================================
--- setup.py	(revision 3eb331258188803a81fbe7efcc02852d2c6e3e39)
+++ setup.py	(revision f903f0a1f32b4a47b938942c6a78b4aa4d01f065)
@@ -3,5 +3,6 @@
 packages = find_packages(exclude=['contrib', 'docs', 'tests*'])
 package_data = {
-    'sasmodels.models': ['*.c'],
+    'sasmodels.models': ['*.c','lib/*.c'],
+    'sasmodels': ['*.c'],
 }
 required = []
