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);
+}
+
