Index: .gitignore
===================================================================
--- .gitignore	(revision df4dc86e6a7ca1d5dd8961a1f3a5f0f54d3f8e37)
+++ .gitignore	(revision 946cdc8e3a248b1d6d56a371c9aef64f499bced8)
@@ -1,4 +1,5 @@
 /build/
 /dist/
+/*.csv
 *.pyc
 *.cl
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision 9c117a293e54d65d07bc4902eeab78dc58a6a8f8)
+++ sasmodels/bumps_model.py	(revision 9c117a293e54d65d07bc4902eeab78dc58a6a8f8)
@@ -17,5 +17,5 @@
 try:
     from .kernelcl import load_model as _loader
-except ImportError,exc:
+except RuntimeError,exc:
     import warnings
     warnings.warn(str(exc))
Index: sasmodels/direct_model.py
===================================================================
--- sasmodels/direct_model.py	(revision ae7b97b7619696aff6e6f1d4549b539e6ecf8a1d)
+++ sasmodels/direct_model.py	(revision ae7b97b7619696aff6e6f1d4549b539e6ecf8a1d)
@@ -0,0 +1,96 @@
+import warnings
+
+import numpy as np
+
+from . import models
+from . import weights
+
+try:
+    from .kernelcl import load_model
+except ImportError,exc:
+    warnings.warn(str(exc))
+    warnings.warn("using ctypes instead")
+    from .kerneldll import load_model
+
+def load_model_definition(model_name):
+    __import__('sasmodels.models.'+model_name)
+    model_definition = getattr(models, model_name, None)
+    return model_definition
+
+# load_model is imported above.  It looks like the following
+#def load_model(model_definition, dtype='single):
+#    if kerneldll:
+#        if source is newer than compiled: compile
+#        load dll
+#        return kernel
+#    elif kernelcl:
+#        compile source on context
+#        return kernel
+
+
+def make_kernel(model, q_vectors):
+    """
+    Return a computation kernel from the model definition and the q input.
+    """
+    input = model.make_input(q_vectors)
+    return model(input)
+
+def get_weights(kernel, pars, name):
+    """
+    Generate the distribution for parameter *name* given the parameter values
+    in *pars*.
+
+    Searches for "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma"
+    """
+    relative = name in kernel.info['partype']['pd-rel']
+    limits = kernel.info['limits']
+    disperser = pars.get(name+'_pd_type', 'gaussian')
+    value = pars.get(name)
+    npts = pars.get(name+'_pd_n', 0)
+    width = pars.get(name+'_pd', 0.0)
+    nsigma = pars.get(name+'_pd_nsigma', 3.0)
+    v,w = weights.get_weights(
+        disperser, npts, width, nsigma,
+        value, limits[name], relative)
+    return v,w/np.sum(w)
+
+def call_kernel(kernel, pars):
+    fixed_pars = [pars.get(name, kernel.info['defaults'][name])
+                  for name in kernel.fixed_pars]
+    pd_pars = [get_weights(kernel, pars, name) for name in kernel.pd_pars]
+    return kernel(fixed_pars, pd_pars)
+
+class DirectModel:
+    def __init__(self, name, q_vectors, dtype='single'):
+        self.model_definition = load_model_definition(name)
+        self.model = load_model(self.model_definition, dtype=dtype)
+        q_vectors = [np.ascontiguousarray(q,dtype=dtype) for q in q_vectors]
+        self.kernel = make_kernel(self.model, q_vectors)
+    def __call__(self, pars):
+        return call_kernel(self.kernel, pars)
+
+def demo():
+    import sys
+    if len(sys.argv) < 3:
+        print "usage: python -m sasmodels.direct_model modelname (q|qx,qy) par=val ..."
+        sys.exit(1)
+    model_name = sys.argv[1]
+    values = [float(v) for v in sys.argv[2].split(',')]
+    if len(values) == 1:
+        q = values[0]
+        q_vectors = [[q]]
+    elif len(values) == 2:
+        qx,qy = values
+        q_vectors = [[qx],[qy]]
+    else:
+        print "use q or qx,qy"
+        sys.exit(1)
+    model = DirectModel(model_name, q_vectors)
+    pars = dict((k,float(v))
+                for pair in sys.argv[3:]
+                for k,v in [pair.split('=')])
+    Iq = model(pars)
+    print Iq[0]
+
+if __name__ == "__main__":
+    demo()
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision 3271e207340eb2fd5f01b1d8ba963ab4fe9ede8c)
+++ sasmodels/generate.py	(revision ae7b97b7619696aff6e6f1d4549b539e6ecf8a1d)
@@ -731,4 +731,5 @@
     info['limits'] = dict((p[0],p[3]) for p in info['parameters'])
     info['partype'] = categorize_parameters(info['parameters'])
+    info['defaults'] = dict((p[0],p[2]) for p in info['parameters'])
 
     source = make_model(info)
@@ -757,10 +758,15 @@
     print "time:",toc()
 
-def demo():
-    from .models import cylinder
-    source, info = make(cylinder)
-    #print doc(cylinder)
-    print source
+def main():
+    if len(sys.argv) <= 1:
+        print "usage: python -m sasmodels.generate modelname"
+    else:
+        name = sys.argv[1]
+        import sasmodels.models
+        __import__('sasmodels.models.'+name)
+        model = getattr(sasmodels.models, name)
+        source, info = make(model); print source
+        #print doc(model)
 
 if __name__ == "__main__":
-    demo()
+    main()
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision 87fce00c768607daaeafbbd67ca19eb27459bc2e)
+++ sasmodels/kernelcl.py	(revision da32ec3c8f278d6e8b397d2dc433724bdd4b0d2f)
@@ -30,7 +30,14 @@
 try:
     import pyopencl as cl
-except RuntimeError,exc:
+except ImportError,exc:
     warnings.warn(str(exc))
     raise RuntimeError("OpenCL not available")
+
+try:
+    context = cl.create_some_context(interactive=False)
+    del context
+except cl.RuntimeError, exc:
+    warnings.warn(str(exc))
+    raise RuntimeError("OpenCl not available")
 
 from pyopencl import mem_flags as mf
Index: sasmodels/kerneldll.py
===================================================================
--- sasmodels/kerneldll.py	(revision 87fce00c768607daaeafbbd67ca19eb27459bc2e)
+++ sasmodels/kerneldll.py	(revision afb2c7807a7e3f42111b09c19299e1e75edcf0a5)
@@ -24,5 +24,6 @@
     #COMPILE = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG /Tp%(source)s /link /DLL /INCREMENTAL:NO /MANIFEST /OUT:%(output)s"
     #COMPILE = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG /Tp%(source)s /openmp /link /DLL /INCREMENTAL:NO /MANIFEST /OUT:%(output)s"
-    COMPILE = "gcc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %(source)s -o %(output)s -lm"
+    #COMPILE = "gcc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %(source)s -o %(output)s -lm"
+    COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
 else:
     COMPILE = "cc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %(source)s -o %(output)s -lm"
Index: sasmodels/models/bcc.py
===================================================================
--- sasmodels/models/bcc.py	(revision c95dc9084eac92b9289ad0551b8d6b7bcdfa4240)
+++ sasmodels/models/bcc.py	(revision e166cb94a55915b1239c64a45fdab6799a81625f)
@@ -80,5 +80,5 @@
 from numpy import pi, inf
 
-name = "BCCparacrystal"
+name = "bcc_paracrystal"
 title = "Body-centred cubic lattic with paracrystalline distortion"
 description = """
Index: sasmodels/models/broad_peak.py
===================================================================
--- sasmodels/models/broad_peak.py	(revision 636adb6f79af5036379849690e91a2cfdf1b5c02)
+++ sasmodels/models/broad_peak.py	(revision 636adb6f79af5036379849690e91a2cfdf1b5c02)
@@ -0,0 +1,115 @@
+r"""
+This model calculates an empirical functional form for SAS data characterized
+by a broad scattering peak. Many SAS spectra are characterized by a broad peak
+even though they are from amorphous soft materials. For example, soft systems
+that show a SAS peak include copolymers, polyelectrolytes, multiphase systems,
+layered structures, etc.
+
+The d-spacing corresponding to the broad peak is a characteristic distance 
+between the scattering inhomogeneities (such as in lamellar, cylindrical, or 
+spherical morphologies, or for bicontinuous structures).
+
+The returned value is scaled to units of |cm^-1|, absolute scale.
+
+Definition
+----------
+
+The scattering intensity *I(q)* is calculated as
+
+.. image:: img/image174.jpg
+
+Here the peak position is related to the d-spacing as *Q0* = 2|pi| / *d0*.
+
+For 2D data: The 2D scattering intensity is calculated in the same way as 1D,
+where the *q* vector is defined as
+
+.. image:: img/image040.gif
+
+=====================  =========  =============
+Parameter name          Units     Default value
+=====================  =========  =============
+lorentz_scale(=C)       None      10
+porod_scale  (=A)       None      1e-05
+lorentz_length (= |xi| )  |Ang|     50
+peak_pos  (=Q0)         |Ang^-1|  0.1
+porod_exp  (=n)         None      2
+lorentz_exp (=m)        None      3
+Background (=B)        |cm^-1|   0.1
+==================  ========  =============
+
+.. image:: img/image175.jpg
+
+*Figure. 1D plot using the default values (w/200 data point).*
+
+REFERENCE
+
+None.
+
+*2013/09/09 - Description reviewed by King, S and Parker, P.*
+
+"""
+
+import numpy as np
+from numpy import pi, inf, sin, cos, sqrt, exp, log
+
+name = "broad_peak"
+title = "Broad Lorentzian type peak on top of a power law decay"
+description = """\
+      I(q) = scale_p/pow(q,exponent)+scale_l/
+      (1.0 + pow((fabs(q-q_peak)*length_l),exponent_l) )+ background
+
+      List of default parameters:
+      porod_scale = Porod term scaling
+      porod_exp = Porod exponent
+      lorentz_scale = Lorentzian term scaling
+      lorentz_length = Lorentzian screening length [A]
+      peak_pos = peak location [1/A]
+      lorentz_exp = Lorentzian exponent
+      background = Incoherent background"""
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type",
+#     "description" ],
+
+    [ "porod_scale", "", 1.0e-05, [-inf,inf], "",
+      "Power law scale factor" ],
+    [ "porod_exp", "", 3.0, [-inf,inf], "",
+      "Exponent of power law" ],
+    [ "lorentz_scale", "", 10.0, [-inf,inf], "",
+      "Scale factor for broad Lorentzian peak" ],
+    [ "lorentz_length", "Ang",  50.0, [-inf, inf], "",
+      "Lorentzian screening length" ],
+    [ "peak_pos", "1/Ang",  0.1, [-inf, inf], "",
+      "Peak postion in q" ],
+    [ "lorentz_exp", "",  2.0, [-inf, inf], "",
+      "exponent of Lorentz function" ],
+    ]
+
+
+def form_volume():
+    return 1
+
+def Iq(q, porod_scale, porod_exp, lorentz_scale, lorentz_length, peak_pos, lorentz_exp):
+    inten = porod_scale/pow(q,porod_exp) + lorentz_scale/(1.0 \
+        + pow((math.fabs(q-peak_pos)*lorentz_length),lorentz_exp))
+    return inten  
+
+# FOR VECTORIZED VERSION, UNCOMMENT THE NEXT LINE
+Iq.vectorized = True
+
+def Iqxy(qx, qy, porod_scale, porod_exp, lorentz_scale, lorentz_length, peak_pos, lorentz_exp):
+    return Iq(sqrt(qx**2 + qy**2), porod_scale, porod_exp, lorentz_scale, lorentz_length, peak_pos, lorentz_exp)
+
+# FOR VECTORIZED VERSION, UNCOMMENT THE NEXT LINE
+Iqxy.vectorized = True
+
+
+demo = dict(
+    scale=1, background=0,
+    porod_scale=1.0e-05, porod_exp=3,
+    lorentz_scale=10,lorentz_length=50, peak_pos=0.1, lorentz_exp=2,
+    )
+oldname = "BroadPeakModel"
+oldpars = dict(porod_scale='scale_p', porod_exp='exponent_p', 
+        lorentz_scale='scale_l', lorentz_length='length_l', peak_pos='q_peak', 
+        lorentz_exp='exponent_l')
Index: sasmodels/models/fcc.c
===================================================================
--- sasmodels/models/fcc.c	(revision 3271e207340eb2fd5f01b1d8ba963ab4fe9ede8c)
+++ sasmodels/models/fcc.c	(revision e7b3d7b015b4e57d299b8d962e6f11825c0bad69)
@@ -5,10 +5,10 @@
     double theta, double phi, double psi);
 
-double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
-double _BCCeval(double Theta, double Phi, double temp1, double temp3);
+double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
+double _FCCeval(double Theta, double Phi, double temp1, double temp3);
 double _sphereform(double q, double radius, double sld, double solvent_sld);
 
 
-double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
+double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
 
 	const double Da = d_factor*dnn;
@@ -16,9 +16,9 @@
 	const double temp3 = q*dnn;
 
-	double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
+	double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
 	return(retVal);
 }
 
-double _BCCeval(double Theta, double Phi, double temp1, double temp3) {
+double _FCCeval(double Theta, double Phi, double temp1, double temp3) {
 
 	double temp6,temp7,temp8,temp9,temp10;
@@ -26,9 +26,9 @@
 
 	temp6 = sin(Theta);
-	temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta);
-	temp8 = -1.0*sin(Theta)*cos(Phi)-sin(Theta)*sin(Phi)+cos(Theta);
-	temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)-cos(Theta);
+	temp7 = sin(Theta)*sin(Phi)+cos(Theta);
+	temp8 = -1.0*sin(Theta)*cos(Phi)+cos(Theta);
+	temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi);
 	temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
-	result = pow(1.0-(temp10*temp10),3)*temp6/((1.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10)));
+	result = pow((1.0-(temp10*temp10)),3)*temp6/((1.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10)));
 
 	return (result);
@@ -54,6 +54,6 @@
 
 	//Volume fraction calculated from lattice symmetry and sphere radius
-	const double s1 = dnn/sqrt(0.75);
-	const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
+	const double s1 = dnn*sqrt(2.0);
+	const double latticescale = 4.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
 
     const double va = 0.0;
@@ -71,5 +71,5 @@
 			//20 gauss points for the inner integral
 			double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;		//the inner dummy is theta
-			double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi);
+			double yyy = Gauss150Wt[j] * _FCC_Integrand(q,dnn,d_factor,ztheta,zphi);
 			summj += yyy;
 		}
@@ -140,13 +140,13 @@
     // The following test should always pass
     if (fabs(cos_val_b3)>1.0) {
-      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
       cos_val_b3 = 1.0;
     }
     if (fabs(cos_val_b2)>1.0) {
-      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
       cos_val_b2 = 1.0;
     }
     if (fabs(cos_val_b1)>1.0) {
-      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
       cos_val_b1 = 1.0;
     }
Index: sasmodels/models/fcc.py
===================================================================
--- sasmodels/models/fcc.py	(revision 3271e207340eb2fd5f01b1d8ba963ab4fe9ede8c)
+++ sasmodels/models/fcc.py	(revision e7b3d7b015b4e57d299b8d962e6f11825c0bad69)
@@ -1,7 +1,7 @@
-#bcc paracrystal model
+#fcc paracrystal model
 #note model title and parameter table are automatically inserted
 #note - calculation requires double precision
-"""
-Calculates the scattering from a **body-centered cubic lattice** with paracrystalline distortion. Thermal vibrations
+r"""
+Calculates the scattering from a **face-centered cubic lattice** with paracrystalline distortion. Thermal vibrations
 are considered to be negligible, and the size of the paracrystal is infinitely large. Paracrystalline distortion is
 assumed to be isotropic and characterized by a Gaussian distribution.
@@ -14,14 +14,14 @@
 The scattering intensity *I(q)* is calculated as
 
-.. image:: img/image167.jpg
+.. image:: img/image158.jpg
 
 where *scale* is the volume fraction of spheres, *Vp* is the volume of the primary particle, *V(lattice)* is a volume
 correction for the crystal structure, *P(q)* is the form factor of the sphere (normalized), and *Z(q)* is the
-paracrystalline structure factor for a body-centered cubic structure.
+paracrystalline structure factor for a face-centered cubic structure.
 
-Equation (1) of the 1990 reference is used to calculate *Z(q)*, using equations (29)-(31) from the 1987 paper for
+Equation (1) of the 1990 reference is used to calculate *Z(q)*, using equations (23)-(25) from the 1987 paper for
 *Z1*\ , *Z2*\ , and *Z3*\ .
 
-The lattice correction (the occupied volume of the lattice) for a body-centered cubic structure of particles of radius
+The lattice correction (the occupied volume of the lattice) for a face-centered cubic structure of particles of radius
 *R* and nearest neighbor separation *D* is
 
@@ -34,7 +34,7 @@
 where *g* is a fractional distortion based on the nearest neighbor distance.
 
-The body-centered cubic lattice is
+The face-centered cubic lattice is
 
-.. image:: img/image168.jpg
+.. image:: img/image161.jpg
 
 For a crystal, diffraction peaks appear at reduced q-values given by
@@ -42,8 +42,8 @@
 .. image:: img/image162.jpg
 
-where for a body-centered cubic lattice, only reflections where (\ *h* + *k* + *l*\ ) = even are allowed and
-reflections where (\ *h* + *k* + *l*\ ) = odd are forbidden. Thus the peak positions correspond to (just the first 5)
+where for a face-centered cubic lattice *h*\ , *k*\ , *l* all odd or all even are allowed and reflections where
+*h*\ , *k*\ , *l* are mixed odd/even are forbidden. Thus the peak positions correspond to (just the first 5)
 
-.. image:: img/image169.jpg
+.. image:: img/image163.jpg
 
 **NB: The calculation of** *Z(q)* **is a double numerical integral that must be carried out with a high density of**
@@ -52,8 +52,8 @@
 makes a triple integral. Very, very slow. Go get lunch!
 
-This example dataset is produced using 200 data points, *qmin* = 0.001 |Ang^-1|, *qmax* = 0.1 |Ang^-1| and the above
+This example dataset is produced using 200 data points, *qmin* = 0.01 |Ang^-1|, *qmax* = 0.1 |Ang^-1| and the above
 default values.
 
-.. image:: img/image170.jpg
+.. image:: img/image164.jpg
 
 *Figure. 1D plot in the linear scale using the default values (w/200 data point).*
@@ -65,5 +65,5 @@
 .. image:: img/image165.gif
 
-.. image:: img/image171.jpg
+.. image:: img/image166.jpg
 
 *Figure. 2D plot using the default values (w/200X200 pixels).*
@@ -80,8 +80,8 @@
 from numpy import pi, inf
 
-name = "BCCparacrystal"
-title = "Body-centred cubic lattic with paracrystalline distortion"
+name = "fcc_paracrystal"
+title = "Face-centred cubic lattic with paracrystalline distortion"
 description = """
-    Calculates the scattering from a **body-centered cubic lattice** with paracrystalline distortion. Thermal vibrations
+    Calculates the scattering from a **face-centered cubic lattice** with paracrystalline distortion. Thermal vibrations
     are considered to be negligible, and the size of the paracrystal is infinitely large. Paracrystalline distortion is
     assumed to be isotropic and characterized by a Gaussian distribution.
@@ -100,5 +100,5 @@
     ]
 
-source = [ "lib/J1.c", "lib/gauss150.c", "bcc.c" ]
+source = [ "lib/J1.c", "lib/gauss150.c", "fcc.c" ]
 
 def ER(radius, length):
@@ -119,5 +119,5 @@
 # For testing against the old sasview models, include the converted parameter
 # names and the target sasview model name.
-oldname='BCCrystalModel'
+oldname='FCCrystalModel'
 oldpars=dict(sld='sldSph',
              solvent_sld='sldSolv')
Index: sasmodels/models/hardsphere.py
===================================================================
--- sasmodels/models/hardsphere.py	(revision 301e09699a18f9efe29522bd73400f85dc156d6b)
+++ sasmodels/models/hardsphere.py	(revision 301e09699a18f9efe29522bd73400f85dc156d6b)
@@ -0,0 +1,107 @@
+# Note: model title and parameter table are inserted automatically
+r"""This calculates the interparticle structure factor for monodisperse spherical particles interacting through hard
+sphere (excluded volume) interactions.
+
+The calculation uses the Percus-Yevick closure where the interparticle potential is
+
+.. image:: img/HardSphere_223.PNG
+
+where *r* is the distance from the center of the sphere of a radius *R*.
+
+For a 2D plot, the wave transfer is defined as
+
+.. math::
+
+    Q = \sqrt{Q_x^2 + Q_y^2}
+
+
+==============  ========  =============
+Parameter name  Units     Default value
+==============  ========  =============
+effect_radius   |Ang|     50.0
+volfraction     None      0.2
+==============  ========  =============
+
+.. image:: img/HardSphere_224.jpg
+
+*Figure. 1D plot using the default values (in linear scale).*
+
+REFERENCE
+
+J K Percus, J Yevick, *J. Phys. Rev.*, 110, (1958) 1
+"""
+
+from numpy import pi, inf
+
+name = "hardsphere"
+title = "Hard sphere structure factor, with Percus-Yevick closure"
+description = """\
+	[Hard sphere structure factor, with Percus-Yevick closure]
+        Interparticle S(Q) for random, non-interacting spheres.
+		May be a reasonable approximation for other shapes of
+		particles that freely rotate, and for moderately polydisperse
+        systems. Though strictly the maths needs to be modified - 
+		which sasview does not do yet.
+		effect_radius is the hard sphere radius
+		volfraction is the volume fraction occupied by the spheres.
+"""
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type",
+#     "description" ],
+    [ "effect_radius", "Ang",  50.0, [0, inf], "volume",
+      "effective radius of hard sphere" ],
+    [ "volfraction", "",  0.2, [0, 0.74], "",
+      "volume fraction of hard spheres" ],
+     ]
+
+# No volume normalization despite having a volume parameter
+# This should perhaps be volume normalized?
+form_volume = """
+    return 1.0;
+    """
+
+Iq = """
+	double denom,dnum,alpha,beta,gamm,a,asq,ath,afor,rca,rsa;
+	double calp,cbeta,cgam,prefac,c,vstruc;
+	double struc;
+	
+	//  compute constants
+	denom = pow((1.0-volfraction),4);
+	dnum = pow((1.0 + 2.0*volfraction),2);
+	alpha = dnum/denom;
+	beta = -6.0*volfraction*pow((1.0 + volfraction/2.0),2)/denom;
+	gamm = 0.50*volfraction*dnum/denom;
+	//
+	//  calculate the structure factor
+	//     
+	a = 2.0*q*effect_radius;
+	asq = a*a;
+	ath = asq*a;
+	afor = ath*a;
+	SINCOS(a,rsa,rca);
+	//rca = cos(a);
+	//rsa = sin(a);
+	calp = alpha*(rsa/asq - rca/a);
+	cbeta = beta*(2.0*rsa/asq - (asq - 2.0)*rca/ath - 2.0/ath);
+	cgam = gamm*(-rca/a + (4.0/a)*((3.0*asq - 6.0)*rca/afor + (asq - 6.0)*rsa/ath + 6.0/afor));
+	prefac = -24.0*volfraction/a;
+	c = prefac*(calp + cbeta + cgam);
+	vstruc = 1.0/(1.0-c);
+	struc = vstruc;
+	
+	return(struc);
+   """
+
+Iqxy = """
+    // never called since no orientation or magnetic parameters.
+    return -1.0;
+    """
+
+# ER defaults to 0.0
+# VR defaults to 1.0
+
+demo = dict(effect_radius = 200,volfraction = 0.2,effect_radius_pd = 0.1,effect_radius_pd_n = 40)
+oldname = 'HardsphereStructure'
+oldpars = dict()
+
Index: sasmodels/models/lamellarCaille.py
===================================================================
--- sasmodels/models/lamellarCaille.py	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
+++ sasmodels/models/lamellarCaille.py	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
@@ -0,0 +1,120 @@
+# 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
+
+.. image:: img/lamellarCaille_139.PNG
+
+The form factor is
+
+.. image:: img/lamellarCaille_134.PNG
+
+and the structure factor is
+
+.. image:: img/lamellarCaille_.PNG
+
+where
+
+.. image:: img/lamellarCaille_.PNG
+
+Here *d* = (repeat) spacing, |delta| = bilayer thickness, the contrast |drho| = 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}
+
+The returned value is in units of |cm^-1|, on absolute scale.
+
+==============  ========  =============
+Parameter name  Units     Default value
+==============  ========  =============
+background      |cm^-1|   0.0
+contrast        |Ang^-2|  5e-06
+scale           None      1
+delta           |Ang|     30
+n_plates        None      20
+spacing         |Ang|     400
+caille          |Ang^-2|  0.1
+==============  ========  =============
+
+.. image:: img/lamellarPS_142.jpg
+
+*Figure. 1D plot using the default values (w/6000 data point).*
+
+Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research
+(Kline, 2006).
+
+REFERENCE
+
+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 pi, 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
+"""
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type",
+#     "description" ],
+    [ "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", "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" ],
+    ]
+
+source = [ "lamellarCaille_kernel.c"]
+
+# No volume normalization despite having a volume parameter
+# This should perhaps be volume normalized?
+form_volume = """
+    return 1.0;
+    """
+
+Iqxy = """
+    // never called since no orientation or magnetic parameters.
+    return -1.0;
+    """
+
+# 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')
Index: sasmodels/models/lamellarCailleHG.py
===================================================================
--- sasmodels/models/lamellarCailleHG.py	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
+++ sasmodels/models/lamellarCailleHG.py	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
@@ -0,0 +1,134 @@
+# 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
+
+.. image:: img/lamellarCailleHG_139.PNG
+
+The form factor is
+
+.. image:: img/lamellarCailleHG_143.PNG
+
+and the structure factor is
+
+.. image:: img/lamellarCailleHG_140.PNG
+
+where
+
+.. image:: img/lamellarCailleHG_141.PNG
+
+where |delta|\ T = tail length (or *tail_length*), |delta|\ H = head thickness (or *h_thickness*),
+|drho|\ H = SLD(headgroup) - SLD(solvent), and |drho|\ T = SLD(tail) - SLD(headgroup).
+Here *d* = (repeat) spacing, *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}
+
+The returned value is in units of |cm^-1|, on absolute scale.
+
+==============  ========  =============
+Parameter name  Units     Default value
+==============  ========  =============
+background      |cm^-1|   0.001
+sld_head        |Ang^-2|  2e-06
+scale           None      1
+sld_solvent     |Ang^-2|  6e-06
+deltaH          |Ang|     2
+deltaT          |Ang|     10
+sld_tail        |Ang^-2|  0
+n_plates        None      30
+spacing         |Ang|     40
+caille          |Ang^-2|  0.001
+==============  ========  =============
+
+.. image:: img/lamellarCailleHG_142.jpg
+
+*Figure. 1D plot using the default values (w/6000 data point).*
+
+Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research
+(Kline, 2006).
+
+REFERENCE
+
+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 pi, 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
+"""
+
+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 = """
+    // never called since no orientation or magnetic parameters.
+    return -1.0;
+    """
+
+# 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,
+		tail_length_pd= 0.1, tail_length_pd_n=20,
+		head_length_pd= 0.05, head_length_pd_n=30,
+		spacing_pd= 0, 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')
Index: sasmodels/models/lamellarCailleHG_kernel.c
===================================================================
--- sasmodels/models/lamellarCailleHG_kernel.c	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
+++ sasmodels/models/lamellarCailleHG_kernel.c	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
@@ -0,0 +1,71 @@
+/*	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);
+
+  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);
+    //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/lamellarCaille_kernel.c
===================================================================
--- sasmodels/models/lamellarCaille_kernel.c	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
+++ sasmodels/models/lamellarCaille_kernel.c	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
@@ -0,0 +1,68 @@
+/*	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;
+  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/lamellarFFHG.py
===================================================================
--- sasmodels/models/lamellarFFHG.py	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
+++ sasmodels/models/lamellarFFHG.py	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
@@ -0,0 +1,128 @@
+# 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.
+
+*2.1.31.1. Definition*
+
+The scattering intensity *I(q)* is
+
+.. math::
+
+    I(Q) = 2\pi{P(Q) \over (2(|delta|\ H +|delta|\ T) Q^2)
+
+The form factor is
+
+.. image:: img/lamellarFFHG_.jpg
+
+where |delta|\ T = tail length (or *tail_length*), |delta|\ H = head thickness (or *h_thickness*),
+|drho|\ H = SLD(headgroup) - SLD(solvent), and |drho|\ T = SLD(tail) - SLD(solvent).
+
+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}
+
+The returned value is in units of |cm^-1|, on absolute scale. In the parameters, *sld_tail* = SLD of the tail group,
+and *sld_head* = SLD of the head group.
+
+==============  ========  =============
+Parameter name  Units     Default value
+==============  ========  =============
+background      |cm^-1|   0.0
+head_sld        |Ang^-2|  3e-06
+scale            None      1
+solvent_sld     |Ang^-2|  6e-06
+head_length     |Ang|     10
+tail_length     |Ang|     15
+sld (tail)      |Ang^-2|  0
+==============  ========  =============
+
+.. image:: img/lamellarFFHG_138.jpg
+
+*Figure. 1D plot using the default values (w/1000 data point).*
+
+Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research
+(Kline, 2006).
+
+REFERENCE
+
+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 pi, 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
+		sld_head = Head scattering length density
+		sld_solvent = solvent scattering length density
+		background = incoherent background
+		scale = scale factor
+"""
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type",
+#     "description" ],
+    [ "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" ],
+    ]
+
+# 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;	
+	
+	return inten /= 2.0*(head_length+tail_length);	//normalize by the bilayer thickness	
+		
+    """
+	
+Iqxy = """
+    // never called since no orientation or magnetic parameters.
+    return -1.0;
+    """
+
+# 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', sld='sld_tail', head_sld='sld_head', solvent_sld='sld_solvent')
+
Index: sasmodels/models/lamellarPC.py
===================================================================
--- sasmodels/models/lamellarPC.py	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
+++ sasmodels/models/lamellarPC.py	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
@@ -0,0 +1,117 @@
+# 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.**
+
+*2.1.33.1. Definition*
+
+The scattering intensity *I(q)* is calculated as
+
+.. image:: img/image145.jpg
+
+The form factor of the bilayer is approximated as the cross section of an infinite, planar bilayer of thickness *t*
+
+.. image:: img/image146.jpg
+
+Here, the scale factor is used instead of the mass per area of the bilayer (*G*). The scale factor is the volume
+fraction of the material in the bilayer, *not* the total excluded volume of the paracrystal. *Z*\ :sub:`N`\ *(q)*
+describes the interference effects for aggregates consisting of more than one bilayer. The equations used are (3-5)
+from the Bergstrom reference below.
+
+Non-integer numbers of stacks are calculated as a linear combination of the lower and higher values
+
+.. image:: img/image147.jpg
+
+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}
+
+The parameters of the model are *Nlayers* = no. of layers, and *pd_spacing* = polydispersity of spacing.
+
+==============  ========  =============
+Parameter name  Units     Default value
+==============  ========  =============
+background      |cm^-1|   0
+scale           None      1
+Nlayers         None      20
+pd_spacing      None      0.2
+sld_layer       |Ang^-2|  1e-6
+sld_solvent     |Ang^-2|  6.34e-6
+spacing         |Ang|     250
+thickness       |Ang|     33
+==============  ========  =============
+
+.. image:: img/image148.jpg
+
+*Figure. 1D plot using the default values above (w/20000 data point).*
+
+Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research
+(Kline, 2006).
+
+REFERENCE
+
+M Bergstrom, J S Pedersen, P Schurtenberger, S U Egelhaaf, *J. Phys. Chem. B*, 103 (1999) 9888-9897
+
+"""
+
+from numpy import pi, 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
+"""
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type",
+#     "description" ],
+    [ "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 of paracrystal stack" ],
+    [ "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"]
+
+# No volume normalization despite having a volume parameter
+# This should perhaps be volume normalized?
+form_volume = """
+    return 1.0;
+    """
+
+Iqxy = """
+    // never called since no orientation or magnetic parameters.
+    return -1.0;
+    """
+
+# 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')
+
+
Index: sasmodels/models/lamellarPC_kernel.c
===================================================================
--- sasmodels/models/lamellarPC_kernel.c	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
+++ sasmodels/models/lamellarPC_kernel.c	(revision dc02af0f5a8fa1ee637e4aeeee5e5fe5c977ce74)
@@ -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 = 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;
+	
+	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);
+}
+
