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: explore/J1.py
===================================================================
--- explore/J1.py	(revision aceac39d3abe2b92f7aabc01582ca8c47a704cef)
+++ explore/J1.py	(revision aceac39d3abe2b92f7aabc01582ca8c47a704cef)
@@ -0,0 +1,107 @@
+r"""
+Show numerical precision of $2 J_1(x)/x$.
+"""
+import sys; sys.path.insert(0, '..')
+
+import numpy as np
+try:
+    from mpmath import mp
+except:
+    # CRUFT: mpmath split out into its own package
+    from sympy.mpmath import mp
+#import matplotlib; matplotlib.use('TkAgg')
+import pylab
+
+
+SHOW_DIFF = True # True if show diff rather than function value
+#SHOW_DIFF = False # True if show diff rather than function value
+LINEAR_X = False  # True if q is linearly spaced instead of log spaced
+#LINEAR_X = True # True if q is linearly spaced instead of log spaced
+FUNCTION = "J1(x)"
+
+def mp_fn(vec, bits=500):
+    """
+    Direct calculation using sympy multiprecision library.
+    """
+    with mp.workprec(bits):
+        return [_mp_fn(mp.mpf(x)) for x in vec]
+
+def _mp_fn(x):
+    """
+    Actual function that gets evaluated.  The caller just vectorizes.
+    """
+    #return mp.mpf(2)*mp.j1(x)/x
+    return mp.j1(x)
+
+def np_fn(x, dtype):
+    """
+    Direct calculation using scipy.
+    """
+    from scipy.special import j1 as J1
+    x = np.asarray(x, dtype)
+    #return np.asarray(2, dtype)*J1(x)/x
+    return J1(x)
+
+def sasmodels_fn(x, dtype, platform='ocl'):
+    """
+    Calculation using pade approximant.
+    """
+    from sasmodels import core, data, direct_model
+    model = core.load_model('bessel', dtype=dtype)
+    calculator = direct_model.DirectModel(data.empty_data1D(x), model)
+    #ret = calculator(background=0)
+    #print ret
+    return calculator(background=0)
+
+def plotdiff(x, target, actual, label):
+    """
+    Plot the computed value.
+
+    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
+    """
+    if SHOW_DIFF:
+        err = abs((target-actual)/target)
+        #err = np.clip(err, 0, 1)
+        pylab.loglog(x, err, '-', label=label)
+    else:
+        limits = np.min(target), np.max(target)
+        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
+
+def compare(x, precision, target):
+    r"""
+    Compare the different computation methods using the given precision.
+    """
+    #plotdiff(x, target, mp_fn(x, 11), 'mp 11 bits')
+    plotdiff(x, target, np_fn(x, precision), 'numpy '+precision)
+    plotdiff(x, target, sasmodels_fn(x, precision, 0), 'sasmodels '+precision)
+    pylab.xlabel("qr (1/Ang)")
+    if SHOW_DIFF:
+        pylab.ylabel("relative error")
+    else:
+        pylab.ylabel(FUNCTION)
+        pylab.semilogx(x, target,  '-', label="true value")
+    if LINEAR_X:
+        pylab.xscale('linear')
+
+def main():
+    r"""
+    Compare accuracy of different methods for computing $3 j_1(x)/x$.
+    :return:
+    """
+    if LINEAR_X:
+        qr = np.linspace(1,1000,2000)
+    else:
+        qr = np.logspace(-3,5,400)
+    target = np.asarray(mp_fn(qr), 'double')
+    pylab.subplot(121)
+    compare(qr, 'single', target)
+    pylab.legend(loc='best')
+    pylab.subplot(122)
+    compare(qr, 'double', target)
+    pylab.legend(loc='best')
+    pylab.suptitle(FUNCTION)
+
+if __name__ == "__main__":
+    #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))
+    main()
+    pylab.show()
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision 02e70ffbd2f37e9a91d5d027972b30069bd05c6b)
+++ sasmodels/core.py	(revision d5ba8419b0a797faa3fe41fde9334c5a19282992)
@@ -170,8 +170,7 @@
     """
     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
 
Index: sasmodels/models/bessel.py
===================================================================
--- sasmodels/models/bessel.py	(revision cbd37a7e9b2cd2085f3e1d42c1ab4a48c4db56a7)
+++ sasmodels/models/bessel.py	(revision a5af4e1ca0acba8ff274ff823f33faa5a25b0704)
@@ -78,5 +78,5 @@
 
 Iq = """
-    return 2.0*j1(q)/q;
+    return J1(q);
     """
 
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/hardsphere.py
===================================================================
--- sasmodels/models/hardsphere.py	(revision 3bcd03d7e1d8d4938d30f17bf360b83417139376)
+++ sasmodels/models/hardsphere.py	(revision 934f9066e07d143504e27d267915514971b7ceb0)
@@ -58,4 +58,5 @@
 category = "structure-factor"
 structure_factor = True
+single = False
 
 #             ["name", "units", default, [lower, upper], "type","description"],
@@ -74,4 +75,13 @@
 Iq = """
       double D,A,B,G,X,X2,X4,S,C,FF,HARDSPH;
+      // these are c compiler instructions, can also put normal code inside the "if else" structure 
+      #if FLOAT_SIZE > 4
+      // double precision    orig had 0.2, don't call the variable cutoff as PAK already has one called that! Must use UPPERCASE name please.
+      //  0.05 better, 0.1 OK
+      #define CUTOFFHS 0.05
+      #else
+      // 0.1 bad, 0.2 OK, 0.3 good, 0.4 better, 0.8 no good
+      #define CUTOFFHS 0.4  
+      #endif
 
       if(fabs(radius_effective) < 1.E-12) {
@@ -96,6 +106,6 @@
       G=0.5*volfraction*A;
 
-      if(X < 0.2) {
-      // RKH Feb 2016, use Taylor series expansion for small X, IT IS VERY PICKY ABOUT THE X CUT OFF VALUE, ought to be lower in double. 
+      if(X < CUTOFFHS) {
+      // RKH Feb 2016, use Taylor series expansion for small X
       // else no obvious way to rearrange the equations to avoid needing a very high number of significant figures. 
       // Series expansion found using Mathematica software. Numerical test in .xls showed terms to X^2 are sufficient 
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/models/lib/j1_cephes.c
===================================================================
--- sasmodels/models/lib/j1_cephes.c	(revision e2af2a953458fc0c3f5a1c0c5754c8455115eac8)
+++ sasmodels/models/lib/j1_cephes.c	(revision fad5dc1e719d8dec9f7d9592debc76d14942f224)
@@ -39,21 +39,6 @@
 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
 */
-double j1(double );
-
-double j1(double x) {
-
-//Cephes double pression function
-#if FLOAT_SIZE>4
-
-    double w, z, p, q, xn;
-
-    const double DR1 = 5.78318596294678452118E0;
-    const double DR2 = 3.04712623436620863991E1;
-    const double Z1 = 1.46819706421238932572E1;
-    const double Z2 = 4.92184563216946036703E1;
-    const double THPIO4 =  2.35619449019234492885;
-    const double SQ2OPI = 0.79788456080286535588;
-
-    double RP[8] = {
+
+constant double RP[8] = {
     -8.99971225705559398224E8,
     4.52228297998194034323E11,
@@ -63,20 +48,18 @@
     0.0,
     0.0,
-    0.0
-    };
-
-    double RQ[8] = {
-    /* 1.00000000000000000000E0,*/
-    6.20836478118054335476E2,
-    2.56987256757748830383E5,
-    8.35146791431949253037E7,
-    2.21511595479792499675E10,
-    4.74914122079991414898E12,
-    7.84369607876235854894E14,
-    8.95222336184627338078E16,
-    5.32278620332680085395E18,
-    };
-
-    double PP[8] = {
+    0.0 };
+
+constant double RQ[8] = {
+        6.20836478118054335476E2,
+        2.56987256757748830383E5,
+        8.35146791431949253037E7,
+        2.21511595479792499675E10,
+        4.74914122079991414898E12,
+        7.84369607876235854894E14,
+        8.95222336184627338078E16,
+        5.32278620332680085395E18
+    };
+
+constant double PP[8] = {
     7.62125616208173112003E-4,
     7.31397056940917570436E-2,
@@ -86,7 +69,8 @@
     5.21451598682361504063E0,
     1.00000000000000000254E0,
-    0.0,
-    };
-    double PQ[8] = {
+    0.0} ;
+
+
+constant double PQ[8] = {
     5.71323128072548699714E-4,
     6.88455908754495404082E-2,
@@ -96,8 +80,7 @@
     5.20982848682361821619E0,
     9.99999999999999997461E-1,
-    0.0,
-    };
-
-    double QP[8] = {
+    0.0 };
+
+constant double QP[8] = {
     5.10862594750176621635E-2,
     4.98213872951233449420E0,
@@ -107,9 +90,7 @@
     5.97489612400613639965E2,
     2.11688757100572135698E2,
-    2.52070205858023719784E1,
-    };
-
-    double QQ[8] = {
-    /* 1.00000000000000000000E0,*/
+    2.52070205858023719784E1 };
+
+constant double QQ[8] = {
     7.42373277035675149943E1,
     1.05644886038262816351E3,
@@ -119,41 +100,7 @@
     2.82619278517639096600E3,
     3.36093607810698293419E2,
-    0.0,
-    };
-
-    w = x;
-    if( x < 0 )
-	    w = -x;
-
-    if( w <= 5.0 )
-	{
-	    z = x * x;
-	    w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 );
-	    w = w * x * (z - Z1) * (z - Z2);
-	    return( w );
-	}
-
-    w = 5.0/x;
-    z = w * w;
-    p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
-    q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
-    xn = x - THPIO4;
-
-    double sn, cn;
-    SINCOS(xn, sn, cn);
-    p = p * cn - w * q * sn;
-
-    return( p * SQ2OPI / sqrt(x) );
-
-
-//Single precission version of cephes
-#else
-    double xx, w, z, p, q, xn;
-
-    const double Z1 = 1.46819706421238932572E1;
-    const double THPIO4F =  2.35619449019234492885;    /* 3*pi/4 */
-
-
-    double JP[8] = {
+    0.0 };
+
+constant double JP[8] = {
         -4.878788132172128E-009,
         6.009061827883699E-007,
@@ -166,5 +113,5 @@
     };
 
-    double MO1[8] = {
+constant double MO1[8] = {
         6.913942741265801E-002,
         -2.284801500053359E-001,
@@ -177,5 +124,5 @@
     };
 
-    double PH1[8] = {
+constant double PH1[8] = {
         -4.497014141919556E+001,
         5.073465654089319E+001,
@@ -188,4 +135,51 @@
     };
 
+double J1(double x) {
+
+//Cephes double pression function
+#if FLOAT_SIZE>4
+
+    double w, z, p, q, xn;
+
+    const double Z1 = 1.46819706421238932572E1;
+    const double Z2 = 4.92184563216946036703E1;
+    const double THPIO4 =  2.35619449019234492885;
+    const double SQ2OPI = 0.79788456080286535588;
+
+    w = x;
+    if( x < 0 )
+	    w = -x;
+
+    if( w <= 5.0 )
+	{
+	    z = x * x;
+	    w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 );
+	    w = w * x * (z - Z1) * (z - Z2);
+	    return( w );
+	}
+
+    w = 5.0/x;
+    z = w * w;
+
+    p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
+    q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
+
+    xn = x - THPIO4;
+
+    double sn, cn;
+    SINCOS(xn, sn, cn);
+    p = p * cn - w * q * sn;
+
+    return( p * SQ2OPI / sqrt(x) );
+
+
+//Single precission version of cephes
+#else
+    double xx, w, z, p, q, xn;
+
+    const double Z1 = 1.46819706421238932572E1;
+    const double THPIO4F =  2.35619449019234492885;    /* 3*pi/4 */
+
+
     xx = x;
     if( xx < 0 )
@@ -211,2 +205,6 @@
 }
 
+
+/*double J1c(double x) {
+    return
+}*/
Index: sasmodels/models/lib/polevl.c
===================================================================
--- sasmodels/models/lib/polevl.c	(revision e2af2a953458fc0c3f5a1c0c5754c8455115eac8)
+++ sasmodels/models/lib/polevl.c	(revision fad5dc1e719d8dec9f7d9592debc76d14942f224)
@@ -51,8 +51,6 @@
 */
 
-double polevl( double x, double coef[8], int N );
-double p1evl( double x, double coef[8], int N );
 
-double polevl( double x, double coef[8], int N ) {
+double polevl( double x, constant double *coef, int N ) {
 
     int i = 0;
@@ -74,5 +72,5 @@
  */
 
-double p1evl( double x, double coef[8], int N ) {
+double p1evl( double x, constant double *coef, int N ) {
     int i=0;
     double ans = x+coef[i];
@@ -84,4 +82,3 @@
 
     return( ans );
-
 }
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/sasview_model.py
===================================================================
--- sasmodels/sasview_model.py	(revision 2622b3fe39474acb34f0bbd16dccdce0561d8a26)
+++ sasmodels/sasview_model.py	(revision 08376e707b88225bb8c73e3058792f08aabb9b07)
@@ -23,9 +23,12 @@
 from . import core
 
-def make_class(model_info, dtype='single', namestyle='name'):
+def standard_models():
+    return [make_class(model_name) for model_name in core.list_models()]
+
+def make_class(model_name, namestyle='name'):
     """
     Load the sasview model defined in *kernel_module*.
 
-    Returns a class that can be used directly as a sasview model.
+    Returns a class that can be used directly as a sasview model.t
 
     Defaults to using the new name for a model.  Setting
@@ -33,9 +36,9 @@
     compatible with SasView.
     """
-    model = core.build_model(model_info, dtype=dtype)
+    model_info = core.load_model_info(model_name)
     def __init__(self, multfactor=1):
-        SasviewModel.__init__(self, model)
-    attrs = dict(__init__=__init__)
-    ConstructedModel = type(model.info[namestyle], (SasviewModel,), attrs)
+        SasviewModel.__init__(self)
+    attrs = dict(__init__=__init__, _model_info=model_info)
+    ConstructedModel = type(model_info[namestyle], (SasviewModel,), attrs)
     return ConstructedModel
 
@@ -44,11 +47,11 @@
     Sasview wrapper for opencl/ctypes model.
     """
-    def __init__(self, model):
-        """Initialization"""
-        self._model = model
-
-        self.name = model.info['name']
-        self.oldname = model.info['oldname']
-        self.description = model.info['description']
+    def __init__(self):
+        self._kernel = None
+        model_info = self._model_info
+
+        self.name = model_info['name']
+        self.oldname = model_info['oldname']
+        self.description = model_info['description']
         self.category = None
         self.multiplicity_info = None
@@ -60,7 +63,7 @@
         self.params = collections.OrderedDict()
         self.dispersion = dict()
-        partype = model.info['partype']
-
-        for p in model.info['parameters']:
+        partype = model_info['partype']
+
+        for p in model_info['parameters']:
             self.params[p.name] = p.default
             self.details[p.name] = [p.units] + p.limits
@@ -83,8 +86,8 @@
 
         ## independent parameter name and unit [string]
-        self.input_name = model.info.get("input_name", "Q")
-        self.input_unit = model.info.get("input_unit", "A^{-1}")
-        self.output_name = model.info.get("output_name", "Intensity")
-        self.output_unit = model.info.get("output_unit", "cm^{-1}")
+        self.input_name = model_info.get("input_name", "Q")
+        self.input_unit = model_info.get("input_unit", "A^{-1}")
+        self.output_name = model_info.get("output_name", "Intensity")
+        self.output_unit = model_info.get("output_unit", "cm^{-1}")
 
         ## _persistency_dict is used by sas.perspectives.fitting.basepage
@@ -95,4 +98,17 @@
         ## New fields introduced for opencl rewrite
         self.cutoff = 1e-5
+
+    def __get_state__(self):
+        state = self.__dict__.copy()
+        model_id = self._model_info['id']
+        state.pop('_kernel')
+        # May need to reload model info on set state since it has pointers
+        # to python implementations of Iq, etc.
+        #state.pop('_model_info')
+        return state
+
+    def __set_state__(self, state):
+        self.__dict__ = state
+        self._kernel = None
 
     def __str__(self):
@@ -187,5 +203,5 @@
         # TODO: fix test so that parameter order doesn't matter
         ret = ['%s.%s' % (d.lower(), p)
-               for d in self._model.info['partype']['pd-2d']
+               for d in self._model_info['partype']['pd-2d']
                for p in ('npts', 'nsigmas', 'width')]
         #print(ret)
@@ -261,5 +277,5 @@
             # Check whether we have a list of ndarrays [qx,qy]
             qx, qy = qdist
-            partype = self._model.info['partype']
+            partype = self._model_info['partype']
             if not partype['orientation'] and not partype['magnetic']:
                 return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2))
@@ -284,6 +300,8 @@
         to the card for each evaluation.
         """
+        if self._kernel is None:
+            self._kernel = core.build_model(self._model_info)
         q_vectors = [np.asarray(q) for q in args]
-        fn = self._model(q_vectors)
+        fn = self._kernel(q_vectors)
         pars = [self.params[v] for v in fn.fixed_pars]
         pd_pars = [self._get_weights(p) for p in fn.pd_pars]
@@ -299,5 +317,5 @@
         :return: the value of the effective radius
         """
-        ER = self._model.info.get('ER', None)
+        ER = self._model_info.get('ER', None)
         if ER is None:
             return 1.0
@@ -314,5 +332,5 @@
         :return: the value of the volf ratio
         """
-        VR = self._model.info.get('VR', None)
+        VR = self._model_info.get('VR', None)
         if VR is None:
             return 1.0
@@ -358,5 +376,5 @@
         parameter set in the vector.
         """
-        pars = self._model.info['partype']['volume']
+        pars = self._model_info['partype']['volume']
         return core.dispersion_mesh([self._get_weights(p) for p in pars])
 
@@ -368,6 +386,6 @@
         from . import weights
 
-        relative = self._model.info['partype']['pd-rel']
-        limits = self._model.info['limits']
+        relative = self._model_info['partype']['pd-rel']
+        limits = self._model_info['limits']
         dis = self.dispersion[par]
         value, weight = weights.get_weights(
@@ -375,2 +393,3 @@
             self.params[par], limits[par], par in relative)
         return value, weight / np.sum(weight)
+
