Index: .gitignore
===================================================================
--- .gitignore	(revision 66ebdd668828275e87ada92bedafaa22888a51db)
+++ .gitignore	(revision 4a82d4dd41f712a03b7f365ffe63e370cf57b792)
@@ -1,4 +1,5 @@
 /build/
 /dist/
+/logs/
 /*.csv
 *.pyc
@@ -16,2 +17,4 @@
 /.project
 /.pydevproject
+/.idea
+/sasmodels.egg-info/
Index: doc/_extensions/dollarmath.py
===================================================================
--- doc/_extensions/dollarmath.py	(revision 19dcb93366741a275baa8ad855ce8a310dca78f8)
+++ doc/_extensions/dollarmath.py	(revision 103ea4527aa6bd411ee52e60f2d18b512b56825b)
@@ -12,5 +12,5 @@
 import re
 
-_dollar = re.compile(r"(?:^|(?<=\s))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[.,;\\]))")
+_dollar = re.compile(r"(?:^|(?<=\s|[(]))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[.,;)\\]))")
 _notdollar = re.compile(r"\\[$]")
 
@@ -51,8 +51,15 @@
     assert replace_dollar(u"dollar\$ escape")==u"dollar$ escape"
     assert replace_dollar(u"dollar \$escape\$ too")==u"dollar $escape$ too"
+    assert replace_dollar(u"spaces $in the$ math")==u"spaces :math:`in the` math"
     assert replace_dollar(u"emb\ $ed$\ ed")==u"emb\ :math:`ed`\ ed"
     assert replace_dollar(u"$first$a")==u"$first$a"
     assert replace_dollar(u"a$last$")==u"a$last$"
+    assert replace_dollar(u"$37")==u"$37"
+    assert replace_dollar(u"($37)")==u"($37)"
+    assert replace_dollar(u"$37 - $43")==u"$37 - $43"
+    assert replace_dollar(u"($37, $38)")==u"($37, $38)"
     assert replace_dollar(u"a $mid$dle a")==u"a $mid$dle a"
+    assert replace_dollar(u"a ($in parens$) a")==u"a (:math:`in parens`) a"
+    assert replace_dollar(u"a (again $in parens$) a")==u"a (again :math:`in parens`) a"
 
 if __name__ == "__main__":
Index: doc/_templates/layout.html
===================================================================
--- doc/_templates/layout.html	(revision a98a2cf6961dcaaeea10f5ddb1156173e5df296d)
+++ doc/_templates/layout.html	(revision a98a2cf6961dcaaeea10f5ddb1156173e5df296d)
@@ -0,0 +1,10 @@
+{% extends "!layout.html" %}
+
+{%- block extrahead %}
+ <script type="text/javascript">
+       MathJax.Hub.Config({
+           displayAlign: "left",
+           displayIndent: "2.5em"
+       });
+  </script>
+{% endblock %}
Index: doc/rst_prolog
===================================================================
--- doc/rst_prolog	(revision 19dcb93366741a275baa8ad855ce8a310dca78f8)
+++ doc/rst_prolog	(revision 0a4628d166f9d3d4ed9811ec8882eea0e9675a8e)
@@ -60,2 +60,3 @@
 .. |g/cm^3| replace:: g\ |cdot|\ cm\ :sup:`-3`
 .. |fm^2| replace:: fm\ :sup:`2`
+.. |Ang*cm^-1| replace:: |Ang|\ |cdot|\ cm\ :sup:`-1`
Index: multi_compare.sh
===================================================================
--- multi_compare.sh	(revision 5753e4ebc8ad092cc172a8a013acd92a6a19e86a)
+++ multi_compare.sh	(revision 0a4628d166f9d3d4ed9811ec8882eea0e9675a8e)
@@ -1,3 +1,3 @@
-#!/bin/sh
+#!/bin/bash
 
 sasview=( ../sasview/build/lib.* )
Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision 3964f92c0d8ef50968fdc442ddace0815434fff6)
+++ sasmodels/convert.py	(revision 34d6cabda48bc12af0423cded9c76c0e4bf2572b)
@@ -10,4 +10,5 @@
     'broad_peak',
     'two_lorentzian',
+    "two_power_law",
     'gel_fit',
     'gauss_lorentz_gel',
@@ -151,4 +152,7 @@
         pars['string_thickness_pd_n'] = 0
         pars['number_of_pearls_pd_n'] = 0
+    elif name == 'line':
+        pars['scale'] = 1
+        pars['background'] = 0
     elif name == 'rpa':
         pars['case_num'] = int(pars['case_num'])
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision eafc9fa8b200a3b754e7dafb960c8891361eacef)
+++ sasmodels/generate.py	(revision fa8011eb0e6b527e15b594945876ee39b6a46dbd)
@@ -191,5 +191,5 @@
 The function :func:`make` loads the metadata from the module and returns
 the kernel source.  The function :func:`doc` extracts the doc string
-and adds the parameter table to the top.  The function :func:`sources`
+and adds the parameter table to the top.  The function :func:`model_sources`
 returns a list of files required by the model.
 """
@@ -206,5 +206,5 @@
 import numpy as np
 
-__all__ = ["make", "doc", "sources", "convert_type"]
+#__all__ = ["make", "doc", "model_sources", "convert_type"]
 
 C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c')
@@ -233,4 +233,5 @@
     "degrees": "degree",
     "1/cm": "|cm^-1|",
+    "Ang/cm": "|Ang*cm^-1|",
     "": "None",
     }
Index: sasmodels/kernel_template.c
===================================================================
--- sasmodels/kernel_template.c	(revision 840b85918586f724a23f78252f90ed33ac27fdf7)
+++ sasmodels/kernel_template.c	(revision 960cd80e792d962dac1b7c7326f8409590d8ca81)
@@ -12,40 +12,40 @@
 #ifndef USE_OPENCL
 #  ifdef __cplusplus
-     #include <cstdio>
-     #include <cmath>
-     using namespace std;
-     #if defined(_MSC_VER)
+      #include <cstdio>
+      #include <cmath>
+      using namespace std;
+      #if defined(_MSC_VER)
          #include <limits>
          #include <float.h>
          #define kernel extern "C" __declspec( dllexport )
          inline double trunc(double x) { return x>=0?floor(x):-floor(-x); }
-	     inline double fmin(double x, double y) { return x>y ? y : x; }
-	     inline double fmax(double x, double y) { return x<y ? y : x; }
-	     inline double isnan(double x) { return _isnan(x); }
-	     #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN
-	     static double cephes_expm1(double x) {
-	         // Adapted from the cephes math library.
-	         // Copyright 1984 - 1992 by Stephen L. Moshier
-	         if (x != x || x == 0.0) {
-	             return x; // NaN and +/- 0
-	         } else if (x < -0.5 || x > 0.5) {
-	             return exp(x) - 1.0;
-	         } else {
-	             const double xsq = x*x;
-	             const double p = (((
-	                  +1.2617719307481059087798E-4)*xsq
-                      +3.0299440770744196129956E-2)*xsq
-                      +9.9999999999999999991025E-1);
-                 const double q = ((((
-                      +3.0019850513866445504159E-6)*xsq
-                      +2.5244834034968410419224E-3)*xsq
-                      +2.2726554820815502876593E-1)*xsq
-                      +2.0000000000000000000897E0);
-                 double r = x * p;
-	             r =  r / (q - r);
-	             return r+r;
-	         }
-	     }
-	     #define expm1 cephes_expm1
+         inline double fmin(double x, double y) { return x>y ? y : x; }
+         inline double fmax(double x, double y) { return x<y ? y : x; }
+         inline double isnan(double x) { return _isnan(x); }
+         #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN
+         static double cephes_expm1(double x) {
+            // Adapted from the cephes math library.
+            // Copyright 1984 - 1992 by Stephen L. Moshier
+            if (x != x || x == 0.0) {
+               return x; // NaN and +/- 0
+            } else if (x < -0.5 || x > 0.5) {
+               return exp(x) - 1.0;
+            } else {
+               const double xsq = x*x;
+               const double p = (((
+                  +1.2617719307481059087798E-4)*xsq
+                  +3.0299440770744196129956E-2)*xsq
+                  +9.9999999999999999991025E-1);
+               const double q = ((((
+                  +3.0019850513866445504159E-6)*xsq
+                  +2.5244834034968410419224E-3)*xsq
+                  +2.2726554820815502876593E-1)*xsq
+                  +2.0000000000000000000897E0);
+               double r = x * p;
+               r =  r / (q - r);
+               return r+r;
+             }
+         }
+         #define expm1 cephes_expm1
      #else
          #define kernel extern "C"
@@ -260,6 +260,6 @@
         const double vol_weight = VOLUME_WEIGHT_PRODUCT;
         vol += vol_weight*form_volume(VOLUME_PARAMETERS);
-      #endif
         norm_vol += vol_weight;
+      #endif
       }
       //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,scattering); }
Index: sasmodels/models/be_polyelectrolyte.py
===================================================================
--- sasmodels/models/be_polyelectrolyte.py	(revision 168052c5791506d7dfad00173f5636fc5701ad11)
+++ sasmodels/models/be_polyelectrolyte.py	(revision 0e86967265f7d2cf12cf1cd8abb371e7cc8e8f46)
@@ -127,6 +127,6 @@
     :return:     2D-Intensity
     """
-    iq = Iq(sqrt(qx**2 + qy**2), *args)
-    return iq
+    intensity = Iq(sqrt(qx**2 + qy**2), *args)
+    return intensity
 
 Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values
Index: sasmodels/models/core_shell_bicelle.py
===================================================================
--- sasmodels/models/core_shell_bicelle.py	(revision 8007311db20a23588c62551c8a69f089cf8835e1)
+++ sasmodels/models/core_shell_bicelle.py	(revision fa8011eb0e6b527e15b594945876ee39b6a46dbd)
@@ -7,5 +7,5 @@
 The form factor is normalized by the particle volume.
 
-.. _core-shell-cylinder-geometry:
+.. _core-shell-bicelle-geometry:
 
 .. figure:: img/core_shell_bicelle_geometry.png
Index: sasmodels/models/core_shell_cylinder.py
===================================================================
--- sasmodels/models/core_shell_cylinder.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/core_shell_cylinder.py	(revision f0aa7f88c13775739f654942a5ec343bf210d1f5)
@@ -155,4 +155,7 @@
 
 def ER(radius, thickness, length):
+    """
+        Returns the effective radius used in the S*P calculation
+    """
     radius = radius + thickness
     length = length + 2 * thickness
@@ -161,4 +164,7 @@
 
 def VR(radius, thickness, length):
+    """
+        Returns volume ratio
+    """
     whole = pi * (radius + thickness) ** 2 * (length + 2 * thickness)
     core = pi * radius ** 2 * length
Index: sasmodels/models/core_shell_sphere.c
===================================================================
--- sasmodels/models/core_shell_sphere.c	(revision edc9f8d2c5449514c0206c9ed35db99d865b35b6)
+++ sasmodels/models/core_shell_sphere.c	(revision edc9f8d2c5449514c0206c9ed35db99d865b35b6)
@@ -0,0 +1,31 @@
+double form_volume(double radius, double thickness);
+double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld);
+double Iqxy(double qx, double qy, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld);
+
+
+double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld) {
+    // Core first, then add in shell
+    const double core_qr = q * radius;
+    const double core_contrast = core_sld - shell_sld;
+    const double core_bes = sph_j1c(core_qr);
+    const double core_volume = 4.0 * M_PI / 3.0 * radius * radius * radius;
+    double f = core_volume * core_bes * core_contrast;
+
+    // Now the shell
+    const double shell_qr = q * (radius + thickness);
+    const double shell_contrast = shell_sld - solvent_sld;
+    const double shell_bes = sph_j1c(shell_qr);
+    const double shell_volume = 4.0 * M_PI / 3.0 * pow((radius + thickness), 3);
+    f += shell_volume * shell_bes * shell_contrast;
+    return f * f * 1.0e-4;
+}
+
+double Iqxy(double qx, double qy, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld) {
+    const double q = sqrt(qx*qx+qy*qy);
+    return Iq(q, radius, thickness, core_sld, shell_sld, solvent_sld);
+}
+
+double form_volume(double radius, double thickness)
+{
+    return 4.0 * M_PI / 3.0 * pow((radius + thickness), 3);
+}
Index: sasmodels/models/core_shell_sphere.py
===================================================================
--- sasmodels/models/core_shell_sphere.py	(revision fa8011eb0e6b527e15b594945876ee39b6a46dbd)
+++ sasmodels/models/core_shell_sphere.py	(revision fa8011eb0e6b527e15b594945876ee39b6a46dbd)
@@ -0,0 +1,111 @@
+r"""
+This model provides the form factor, $P(q)$, for a spherical particle with a core-shell structure.
+The form factor is normalized by the particle volume.
+
+Definition
+----------
+
+The 1D scattering intensity is calculated in the following way (Guinier, 1955)
+
+.. math::
+
+    P(q) = \frac{\text{scale}}{V} F^2(q) + \text{background}
+
+where
+
+.. math::
+
+    F^2(q)=\frac{3}{V_s}\left[V_c(\rho_c-\rho_s)\frac{\sin(qr_c)-qr_c\cos(qr_c)}{(qr_c)^3}+
+    V_s(\rho_s-\rho_{solv})\frac{\sin(qr_s)-qr_s\cos(qr_s)}{(qr_s)^3}\right]
+
+where $V_s$ is the volume of the outer shell, $V_c$ is
+the volume of the core, $r_s$ is the radius of the shell, $r_c$ is the radius of the
+core, $\rho_c$ is the scattering length density of the core, $\rho_s$ is the scattering length
+density of the shell, $\rho_{solv}$ is the scattering length density of the solvent.
+
+The 2D scattering intensity is the same as $P(q)$ above, regardless of the
+orientation of the $q$ vector.
+
+NB: The outer most radius (ie, = radius + thickness) is used as the
+effective radius for $S(Q)$ when $P(Q) \cdot S(Q)$ is applied.
+
+Reference
+---------
+
+A Guinier and G Fournet, *Small-Angle Scattering of X-Rays*, John Wiley and Sons, New York, (1955)
+
+Validation
+----------
+
+Validation of our code was done by comparing the output of the 1D model to the output of
+the software provided by NIST (Kline, 2006). Figure 1 shows a comparison of the output of
+our model and the output of the NIST software.
+
+.. figure:: img/core_shell_sphere_1d.jpg
+
+    Comparison of the SasView scattering intensity for a core-shell sphere with
+    the output of the NIST SANS analysis software. The parameters were set to:
+    *scale* = 1.0, *radius* = 60 , *contrast* = 1e-6 |Ang^-2|, and
+    *background* = 0.001 |cm^-1|.
+"""
+
+from numpy import pi, inf
+
+name = "core_shell_sphere"
+title = "Form factor for a monodisperse spherical particle with particle with a core-shell structure."
+description = """
+    F^2(q) = 3/V_s [V_c (core_sld-shell_sld) (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3 
+                   + V_s (shell_sld-solvent_sld) (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3]
+
+            V_s: Volume of the sphere shell
+            V_c: Volume of the sphere core
+            r_s: Shell radius = radius + thickness
+"""
+category = "shape:sphere"
+
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["radius",      "Ang",        60.0, [0, inf],    "volume", "Sphere core radius"],
+              ["thickness",   "Ang",        10.0, [0, inf],    "volume", "Sphere shell thickness"],
+              ["core_sld",    "1e-6/Ang^2", 1.0,  [-inf, inf], "",       "Sphere core scattering length density"],
+              ["shell_sld",   "1e-6/Ang^2", 2.0,  [-inf, inf], "",       "Sphere shell scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2", 3.0,  [-inf, inf],  "",      "Solvent scattering length density"]]
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/sph_j1c.c", "core_shell_sphere.c"]
+
+demo = dict(scale=1, background=0, radius=60, thickness=10,
+            core_sld=1.0, shell_sld=2.0, solvent_sld=0.0)
+
+oldname = 'CoreShellModel'
+oldpars = {}
+
+def ER(radius, thickness):
+    """
+        Equivalent radius
+        @param radius: core radius
+        @param thickness: shell thickness
+    """
+    return radius + thickness
+
+def VR(radius, thickness):
+    """
+        Volume ratio
+        @param radius: core radius
+        @param thickness: shell thickness
+    """
+    whole = 4.0 * pi / 3.0 * pow((radius + thickness), 3)
+    core = 4.0 * pi / 3.0 * radius * radius * radius
+    return whole, whole - core
+
+tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0],
+         [{'radius': 20.0, 'thickness': 10.0}, 'VR', 0.703703704],
+
+         # The SasView test result was 0.00169, with a background of 0.001
+         [{'radius': 60.0,
+           'thickness': 10.0,
+           'core_sld': 1.0,
+           'shell_sld':2.0,
+           'solvent_sld':3.0,
+           'background':0.0
+          }, 0.4, 0.000698838]]
Index: sasmodels/models/correlation_length.py
===================================================================
--- sasmodels/models/correlation_length.py	(revision 5054e805a7c53aa18f659a0c1613f527f05caab3)
+++ sasmodels/models/correlation_length.py	(revision 3e6c5c1168c56a700f8cd12b4025aa6a9764fddf)
@@ -78,15 +78,10 @@
 # names and the target sasview model name.
 oldname = 'CorrLengthModel'
-# pylint: disable=bad-continuation
-oldpars = dict(
-               lorentz_scale='scale_l', porod_scale='scale_p',
+
+oldpars = dict(lorentz_scale='scale_l', porod_scale='scale_p',
                cor_length='length_l', exponent_p='exponent_p',
-               exponent_l='exponent_l'
-               )
+               exponent_l='exponent_l')
 
-tests = [
-         [{}, 0.001, 1009.98],
+tests = [[{}, 0.001, 1009.98],
          [{}, 0.150141, 0.174645],
-         [{}, 0.442528, 0.0203957]
-         ]
-# pylint: enable=bad-continuation
+         [{}, 0.442528, 0.0203957]]
Index: sasmodels/models/cylinder.py
===================================================================
--- sasmodels/models/cylinder.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/cylinder.py	(revision 0c3a22666528577f65591c031dbd899b43980331)
@@ -144,4 +144,7 @@
 
 def ER(radius, length):
+    """
+        Return equivalent radius (ER)
+    """
     ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius))
     return 0.5 * (ddd) ** (1. / 3.)
Index: sasmodels/models/dab.py
===================================================================
--- sasmodels/models/dab.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/dab.py	(revision 94bd8097c6283d5783d3cc4d94cabcc1a41bba78)
@@ -1,4 +1,3 @@
 r"""
-
 Calculates the scattering from a randomly distributed, two-phase system based on
 the Debye-Anderson-Brumberger (DAB) model for such systems. The two-phase system
Index: sasmodels/models/ellipsoid.py
===================================================================
--- sasmodels/models/ellipsoid.py	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
+++ sasmodels/models/ellipsoid.py	(revision 431caaea11d3494f5958adb99c85a4cd92f1290c)
@@ -111,6 +111,7 @@
 ----------
 
-L A Feigin and D I Svergun. *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum,
-New York, 1987.
+L A Feigin and D I Svergun.
+*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*,
+Plenum Press, New York, 1987.
 """
 
Index: sasmodels/models/elliptical_cylinder.c
===================================================================
--- sasmodels/models/elliptical_cylinder.c	(revision a8b3cdb60df4f2e22c9b3c5c3110bb51b7dba2c6)
+++ sasmodels/models/elliptical_cylinder.c	(revision a8b3cdb60df4f2e22c9b3c5c3110bb51b7dba2c6)
@@ -0,0 +1,170 @@
+double form_volume(double r_minor, double r_ratio, double length);
+double Iq(double q, double r_minor, double r_ratio, double length,
+          double sld, double solvent_sld);
+double Iqxy(double qx, double qy, double r_minor, double r_ratio, double length,
+            double sld, double solvent_sld, double theta, double phi, double psi);
+
+
+double _elliptical_cylinder_kernel(double q, double r_minor, double r_ratio, double theta);
+
+double _elliptical_cylinder_kernel(double q, double r_minor, double r_ratio, double theta)
+{
+    // This is the function LAMBDA1^2 in Feigin's notation
+    // q is the q-value for the calculation (1/A)
+    // r_minor is the transformed radius"a" in Feigin's notation
+    // r_ratio is the ratio (major radius)/(minor radius) of the Ellipsoid [=] ---
+    // theta is the dummy variable of the integration
+
+    double retval,arg;
+
+    arg = q*r_minor*sqrt((1.0+r_ratio*r_ratio)/2+(1.0-r_ratio*r_ratio)*cos(theta)/2);
+    if (arg == 0.0){
+        retval = 1.0;
+    }else{
+        retval = 2.0*NR_BessJ1(arg)/arg;
+    }
+    return retval*retval ;
+}
+
+
+double form_volume(double r_minor, double r_ratio, double length)
+{
+    return M_PI * r_minor * r_minor * r_ratio * length;
+}
+
+double Iq(double q, double r_minor, double r_ratio, double length,
+          double sld, double solvent_sld) {
+
+    const int nordi=76; //order of integration
+    const int nordj=20;
+    double va,vb;       //upper and lower integration limits
+    double summ,zi,yyy,answer;         //running tally of integration
+    double summj,vaj,vbj,zij,arg,si;            //for the inner integration
+
+    // orientational average limits
+    va = 0.0;
+    vb = 1.0;
+    // inner integral limits
+    vaj=0.0;
+    vbj=M_PI;
+
+    //initialize integral
+    summ = 0.0;
+
+    const double delrho = sld - solvent_sld;
+
+    for(int i=0;i<nordi;i++) {
+        //setup inner integral over the ellipsoidal cross-section
+        summj=0;
+        zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;     //the "x" dummy
+        arg = r_minor*sqrt(1.0-zi*zi);
+        for(int j=0;j<nordj;j++) {
+            //20 gauss points for the inner integral
+            zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0;        //the "y" dummy
+            yyy = Gauss20Wt[j] * _elliptical_cylinder_kernel(q, arg, r_ratio, zij);
+            summj += yyy;
+        }
+        //now calculate the value of the inner integral
+        answer = (vbj-vaj)/2.0*summj;
+        //divide integral by Pi
+        answer /= M_PI;
+
+        //now calculate outer integral
+        arg = q*length*zi/2.0;
+        if (arg == 0.0){
+            si = 1.0;
+        }else{
+            si = sin(arg) * sin(arg) / arg / arg;
+        }
+        yyy = Gauss76Wt[i] * answer * si;
+        summ += yyy;
+    }
+
+    answer = (vb-va)/2.0*summ;
+    // Multiply by contrast^2
+    answer *= delrho*delrho;
+
+    const double vol = form_volume(r_minor, r_ratio, length);
+    return answer*vol*vol*1e-4;
+}
+
+
+double Iqxy(double qx, double qy, double r_minor, double r_ratio, double length,
+            double sld, double solvent_sld, double theta, double phi, double psi) {
+    const double _theta = theta * M_PI / 180.0;
+    const double _phi = phi * M_PI / 180.0;
+    const double _psi = psi * M_PI / 180.0;
+    const double q = sqrt(qx*qx+qy*qy);
+    const double q_x = qx/q;
+    const double q_y = qy/q;
+
+    //Cylinder orientation
+    double cyl_x = cos(_theta) * cos(_phi);
+    double cyl_y = sin(_theta);
+
+    //cyl_z = -cos(_theta) * sin(_phi);
+
+    // q vector
+    //q_z = 0;
+
+    // Note: cos(alpha) = 0 and 1 will get an
+    // undefined value from CylKernel
+    //alpha = acos( cos_val );
+
+    //ellipse orientation:
+    // the elliptical corss section was transformed and projected
+    // into the detector plane already through sin(alpha)and furthermore psi remains as same
+    // on the detector plane.
+    // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt
+    // the wave vector q.
+
+    //x- y- component on the detector plane.
+    const double ella_x =  -cos(_phi)*sin(_psi) * sin(_theta)+sin(_phi)*cos(_psi);
+    const double ella_y =  sin(_psi)*cos(_theta);
+    const double ellb_x =  -sin(_theta)*cos(_psi)*cos(_phi)-sin(_psi)*sin(_phi);
+    const double ellb_y =  cos(_theta)*cos(_psi);
+
+    // Compute the angle btw vector q and the
+    // axis of the cylinder
+    double cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
+
+    // calculate the axis of the ellipse wrt q-coord.
+    double cos_nu = ella_x*q_x + ella_y*q_y;
+    double cos_mu = ellb_x*q_x + ellb_y*q_y;
+
+    // The following test should always pass
+    if (fabs(cos_val)>1.0) {
+      //printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
+      cos_val = 1.0;
+    }
+    if (fabs(cos_nu)>1.0) {
+      //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n");
+      cos_nu = 1.0;
+    }
+    if (fabs(cos_mu)>1.0) {
+      //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n");
+      cos_mu = 1.0;
+    }
+
+    const double r_major = r_ratio * r_minor;
+    const double qr = q*sqrt( r_major*r_major*cos_nu*cos_nu + r_minor*r_minor*cos_mu*cos_mu );
+    const double qL = q*length*cos_val/2.0;
+
+    double Be;
+    if (qr==0){
+      Be = 0.5;
+    }else{
+      Be = NR_BessJ1(qr)/qr;
+    }
+
+    double Si;
+    if (qL==0){
+      Si = 1.0;
+    }else{
+      Si = sin(qL)/qL;
+    }
+
+    const double k = 2.0 * Be * Si;
+    const double vol = form_volume(r_minor, r_ratio, length);
+    return (sld - solvent_sld) * (sld - solvent_sld) * k * k *vol*vol*1.0e-4;
+}
Index: sasmodels/models/elliptical_cylinder.py
===================================================================
--- sasmodels/models/elliptical_cylinder.py	(revision fa8011eb0e6b527e15b594945876ee39b6a46dbd)
+++ sasmodels/models/elliptical_cylinder.py	(revision fa8011eb0e6b527e15b594945876ee39b6a46dbd)
@@ -0,0 +1,146 @@
+# pylint: disable=line-too-long
+r"""
+This function calculates the scattering from an elliptical cylinder.
+
+Definition for 2D (orientated system)
+-------------------------------------
+
+The angles |theta| and |phi| define the orientation of the axis of the cylinder. The angle |bigpsi| is defined as the
+orientation of the major axis of the ellipse with respect to the vector *Q*\ . A gaussian polydispersity can be added
+to any of the orientation angles, and also for the minor radius and the ratio of the ellipse radii.
+
+.. figure:: img/elliptical_cylinder_geometry.gif
+
+    *a* = *r_minor* and |nu|\ :sub:`n` = $r_ratio$ (i.e., $r_major / r_minor$).
+
+The function calculated is
+
+.. math::
+
+    I(\mathbf{q})=\frac{1}{V_{cyl}}\int{d\psi}\int{d\phi}\int{p(\theta,\phi,\psi)F^2(\mathbf{q},\alpha,\psi)\sin(\theta)d\theta}
+
+with the functions
+
+.. math::
+
+    F(\mathbf{q},\alpha,\psi)=2\frac{J_1(a)\sin(b)}{ab}
+    \\
+    a = \mathbf{q}\sin(\alpha)\left[ r^2_{major}\sin^2(\psi)+r^2_{minor}\cos(\psi) \right]^{1/2}
+    \\
+    b=\mathbf{q}\frac{L}{2}\cos(\alpha)
+
+and the angle |bigpsi| is defined as the orientation of the major axis of the ellipse with respect to the vector $\vec q$ .
+The angle $\alpha$ is the angle between the axis of the cylinder and $\vec q$.
+
+
+Definition for 1D (no preferred orientation)
+--------------------------------------------
+
+The form factor is averaged over all possible orientation before normalized by the particle volume
+
+.. math::
+    P(q) = scale  <F^2> / V
+
+To provide easy access to the orientation of the elliptical cylinder, we define the axis of the cylinder using two
+angles |theta|, |phi| and |bigpsi|. As for the case of the cylinder, the angles |theta| and |phi| are defined on
+Figure 2 of CylinderModel. The angle |bigpsi| is the rotational angle around its own long_c axis against the *q* plane.
+For example, |bigpsi| = 0 when the *r_minor* axis is parallel to the *x*\ -axis of the detector.
+
+All angle parameters are valid and given only for 2D calculation; ie, an oriented system.
+
+.. figure:: img/elliptical_cylinder_geometry_2d.jpg
+
+    Definition of angles for 2D
+
+.. figure:: img/core_shell_bicelle_fig2.jpg
+
+    Examples of the angles for oriented elliptical cylinders against the detector plane.
+
+NB: The 2nd virial coefficient of the cylinder is calculated based on the averaged radius (= sqrt(*r_minor*\ :sup:`2` \* *r_ratio*))
+and length values, and used as the effective radius for *S(Q)* when *P(Q)* \* *S(Q)* is applied.
+
+
+.. figure:: img/elliptical_cylinder_comparison_1d.jpg
+
+    1D plot using the default values (w/1000 data point).
+
+Validation
+----------
+
+Validation of our code was done by comparing the output of the 1D calculation to the angular average of the output of
+the 2D calculation over all possible angles. The figure below shows the comparison where the solid dot refers to
+averaged 2D values while the line represents the result of the 1D calculation (for the 2D averaging, values of 76, 180,
+and 76 degrees are taken for the angles of |theta|, |phi|, and |bigpsi| respectively).
+
+.. figure:: img/elliptical_cylinder_validation_1d.gif
+
+    Comparison between 1D and averaged 2D.
+
+In the 2D average, more binning in the angle |phi| is necessary to get the proper result. The following figure shows
+the results of the averaging by varying the number of angular bins.
+
+.. figure:: img/elliptical_cylinder_averaging.gif
+
+    The intensities averaged from 2D over different numbers of bins and angles.
+
+Reference
+---------
+
+L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum,
+New York, (1987)
+"""
+
+import math
+from numpy import pi, inf
+
+name = "elliptical_cylinder"
+title = "Form factor for an elliptical cylinder."
+description = """
+    Form factor for an elliptical cylinder.
+    See L A Feigin and D I Svergun, Structure Analysis by Small-Angle X-Ray and Neutron Scattering, Plenum, New York, (1987).
+"""
+category = "shape:cylinder"
+
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["r_minor",     "Ang",        20.0,  [0, inf],    "volume",      "Ellipse minor radius"],
+              ["r_ratio",     "",           1.5,   [1, inf],    "volume",      "Ratio of major radius over minor radius"],
+              ["length",      "Ang",        400.0, [1, inf],    "volume",      "Length of the cylinder"],
+              ["sld",         "1e-6/Ang^2", 4.0,   [-inf, inf], "",            "Cylinder scattering length density"],
+              ["solvent_sld", "1e-6/Ang^2", 1.0,   [-inf, inf], "",            "Solvent scattering length density"],
+              ["theta",       "degrees",    90.0,  [-360, 360], "orientation", "In plane angle"],
+              ["phi",         "degrees",    0,     [-360, 360], "orientation", "Out of plane angle"],
+              ["psi",         "degrees",    0,     [-360, 360], "orientation", "Major axis angle relative to Q"]]
+
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/nr_bess_j1.c", "lib/gauss76.c", "lib/gauss20.c", "elliptical_cylinder.c"]
+
+demo = dict(scale=1, background=0, r_minor=100, r_ratio=1.5, length=400.0,
+            sld=4.0, solvent_sld=1.0, theta=10.0, phi=20, psi=30, theta_pd=10, phi_pd=2, psi_pd=3)
+
+oldname = 'EllipticalCylinderModel'
+oldpars = dict(theta='cyl_theta', phi='cyl_phi', psi='cyl_psi', sld='sldCyl', solvent_sld='sldSolv')
+
+def ER(r_minor, r_ratio, length):
+    """
+        Equivalent radius
+        @param r_minor: Ellipse minor radius
+        @param r_ratio: Ratio of major radius over minor radius
+        @param length: Length of the cylinder
+    """
+    radius = math.sqrt(r_minor * r_minor * r_ratio)
+    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius))
+    return 0.5 * (ddd) ** (1. / 3.)
+
+tests = [[{'r_minor': 20.0, 'r_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024],
+         [{'r_minor': 20.0, 'r_ratio': 1.2, 'length':300.0}, 'VR', 1],
+
+         # The SasView test result was 0.00169, with a background of 0.001
+         [{'r_minor': 20.0,
+           'r_ratio': 1.5,
+           'sld': 4.0,
+           'length':400.0,
+           'solvent_sld':1.0,
+           'background':0.0
+          }, 0.001, 675.504402]]
Index: sasmodels/models/fcc.py
===================================================================
--- sasmodels/models/fcc.py	(revision 9aac25dcf63daf48040d6a470ddda4ccb13abbd6)
+++ sasmodels/models/fcc.py	(revision c0ccea8f892e537b9ac3a573d9883875986f4b15)
@@ -114,4 +114,5 @@
 single = False
 
+# pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type","description"],
 parameters = [["dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"],
@@ -124,4 +125,5 @@
               ["psi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"]
              ]
+# pylint: enable=bad-whitespace, line-too-long
 
 source = ["lib/J1.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc.c"]
Index: sasmodels/models/guinier.py
===================================================================
--- sasmodels/models/guinier.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/guinier.py	(revision 723cebed4612ad396b30ca46c6f1e5e801142cba)
@@ -36,7 +36,5 @@
 
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [
-              ["rg", "Ang", 60.0, [0, inf], "", "Radius of Gyration"],
-              ]
+parameters = [["rg", "Ang", 60.0, [0, inf], "", "Radius of Gyration"]]
 
 Iq = """
@@ -51,5 +49,5 @@
 
 # parameters for demo
-demo = dict(scale=1.0,rg=60.0)
+demo = dict(scale=1.0, rg=60.0)
 
 # For testing against the old sasview models, include the converted parameter
@@ -59,5 +57,3 @@
 
 # parameters for unit tests
-tests = [
-         [{'rg' : 31.5}, 0.005, 0.991756]
-         ]
+tests = [[{'rg' : 31.5}, 0.005, 0.991756]]
Index: sasmodels/models/guinier_porod.py
===================================================================
--- sasmodels/models/guinier_porod.py	(revision fa8011eb0e6b527e15b594945876ee39b6a46dbd)
+++ sasmodels/models/guinier_porod.py	(revision fa8011eb0e6b527e15b594945876ee39b6a46dbd)
@@ -0,0 +1,130 @@
+# pylint: disable=line-too-long
+r"""
+Calculates the scattering for a generalized Guinier/power law object.
+This is an empirical model that can be used to determine the size
+and dimensionality of scattering objects, including asymmetric objects
+such as rods or platelets, and shapes intermediate between spheres
+and rods or between rods and platelets.
+
+Definition
+----------
+
+The following functional form is used
+
+.. math::
+    I(q) = \frac{G}{Q^s} \ \exp{\left[ \frac{-Q^2R_g^2}{3-s} \right]} \textrm{ for } Q \leq Q_1
+    \\
+    I(q) = D / Q^m \textrm{ for } Q \geq Q_1
+
+This is based on the generalized Guinier law for such elongated objects
+(see the Glatter reference below). For 3D globular objects (such as spheres), $s = 0$
+and one recovers the standard Guinier formula.
+For 2D symmetry (such as for rods) $s = 1$,
+and for 1D symmetry (such as for lamellae or platelets) $s = 2$.
+A dimensionality parameter ($3-s$) is thus defined, and is 3 for spherical objects,
+2 for rods, and 1 for plates.
+
+Enforcing the continuity of the Guinier and Porod functions and their derivatives yields
+
+.. math::
+    Q_1 = \frac{1}{R_g} \sqrt{(m-s)(3-s)/2}
+
+and
+
+.. math::
+    D = G \ \exp{ \left[ \frac{-Q_1^2 R_g^2}{3-s} \right]} \ Q_1^{m-s}
+      = \frac{G}{R_g^{m-s}} \ \exp{\left[  -\frac{m-s}{2} \right]} \left( \frac{(m-s)(3-s)}{2} \right)^{\frac{m-s}{2}}
+
+
+Note that the radius-of-gyration for a sphere of radius R is given by $R_g = R \sqrt(3/5)$.
+
+The cross-sectional radius-of-gyration for a randomly oriented cylinder
+of radius R is given by $R_g = R / \sqrt(2)$.
+
+The cross-sectional radius-of-gyration of a randomly oriented lamella
+of thickness $T$ is given by $R_g = T / \sqrt(12)$.
+
+For 2D data: 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}
+
+.. figure:: img/guinier_porod_model.jpg
+
+    Guinier-Porod model for $R_g=100$ |Ang|, $s=1$, $m=3$, and $background=0.1$.
+
+
+Reference
+---------
+
+A Guinier, G Fournet, Small-Angle Scattering of X-Rays, John Wiley and Sons, New York, (1955)
+
+O Glatter, O Kratky, Small-Angle X-Ray Scattering, Academic Press (1982)
+Check out Chapter 4 on Data Treatment, pages 155-156.
+"""
+
+from numpy import inf, sqrt, power, exp
+
+name = "guinier_porod"
+title = "Guinier-Porod function"
+description = """\
+         I(q) = scale/q^s* exp ( - R_g^2 q^2 / (3-s) ) for q<= ql
+         = scale/q^m*exp((-ql^2*Rg^2)/(3-s))*ql^(m-s) for q>=ql
+                        where ql = sqrt((m-s)(3-s)/2)/Rg.
+                        List of parameters:
+                        scale = Guinier Scale
+                        s = Dimension Variable
+                        Rg = Radius of Gyration [A] 
+                        m = Porod Exponent
+                        background  = Background [1/cm]"""
+
+category = "shape-independent"
+
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["rg", "Ang", 60.0, [0, inf], "", "Radius of gyration"],
+              ["s",  "",    1.0,  [0, inf], "", "Dimension variable"],
+              ["m",  "",    3.0,  [0, inf], "", "Porod exponent"]]
+# pylint: enable=bad-whitespace, line-too-long
+
+# pylint: disable=C0103
+def Iq(q, rg, s, m):
+    """
+    @param q: Input q-value
+    """
+    n = 3.0 - s
+
+    # Take care of the singular points
+    if rg <= 0.0:
+        return 0.0
+    if (n-3.0+m) <= 0.0:
+        return 0.0
+
+    # Do the calculation and return the function value
+    q1 = sqrt((n-3.0+m)*n/2.0)/rg
+    if q < q1:
+        iq = (1.0/power(q, (3.0-n)))*exp((-q*q*rg*rg)/n)
+    else:
+        iq = (1.0/power(q, m))*exp(-(n-3.0+m)/2.0)*power(((n-3.0+m)*n/2.0),
+                                                         ((n-3.0+m)/2.0))/power(rg, (n-3.0+m))
+    return iq
+
+Iq.vectorized = False  # Iq accepts an array of q values
+
+def Iqxy(qx, qy, *args):
+    """
+    @param qx:   Input q_x-value
+    @param qy:   Input q_y-value
+    @param args: Remaining arguments
+    """
+    return Iq(sqrt(qx ** 2 + qy ** 2), *args)
+
+Iqxy.vectorized = False # Iqxy accepts an array of qx, qy values
+
+demo = dict(scale=1.5, background=0.5, rg=60, s=1.0, m=3.0)
+
+oldname = "GuinierPorodModel"
+oldpars = dict(scale='scale', background='background', s='dim', m='m', rg='rg')
+
+tests = [[{'scale': 1.5, 'background':0.5}, 0.04, 5.290096890253155]]
Index: sasmodels/models/hollow_cylinder.py
===================================================================
--- sasmodels/models/hollow_cylinder.py	(revision ec2ca99b9e3d90e4b91919c523a94d9855ca9c3c)
+++ sasmodels/models/hollow_cylinder.py	(revision e0fd91396c863eb30f45111903316f0bde64c127)
@@ -93,4 +93,5 @@
 source = ["lib/J1.c", "lib/gauss76.c", "hollow_cylinder.c"]
 
+# pylint: disable=W0613
 def ER(radius, core_radius, length):
     """
Index: sasmodels/models/lamellarCailleHG.py
===================================================================
--- sasmodels/models/lamellarCailleHG.py	(revision 13ed84c3c71257d2df9c2ed4b4e98c0bfa4d7d69)
+++ sasmodels/models/lamellarCailleHG.py	(revision 652a78a9a04424b9e84f32b660bf775a67d37463)
@@ -146,7 +146,6 @@
     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]]
-        ]
+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: sasmodels/models/lamellarFFHG.py
===================================================================
--- sasmodels/models/lamellarFFHG.py	(revision 7f47777b1828197366e17d650a92f794a3433fcc)
+++ sasmodels/models/lamellarFFHG.py	(revision 3eb6b901aa1e301f7d79b984a69d3d37509b2eb0)
@@ -68,16 +68,12 @@
 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"],
-             ]
+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
@@ -113,8 +109,8 @@
 
 demo = dict(scale=1, background=0,
-            tail_length=15,head_length=10,
+            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)
+            tail_length_pd=0.2, tail_length_pd_n=40,
+            head_length_pd=0.01, head_length_pd_n=40)
 
 oldname = 'LamellarFFHGModel'
@@ -122,8 +118,4 @@
                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]]
-        ]
-
-
+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: sasmodels/models/lib/gauss76.c
===================================================================
--- sasmodels/models/lib/gauss76.c	(revision 994d77f9b713691e41106d1f9af9661dfb2dc419)
+++ sasmodels/models/lib/gauss76.c	(revision 66d119fbbe121774089e17e3086b4809c9cf2a80)
@@ -7,7 +7,8 @@
  *
  */
+#define N_POINTS_76 76
 
 // Gaussians
-constant double Gauss76Wt[76]={
+constant double Gauss76Wt[N_POINTS_76]={
 	.00126779163408536,		//0
 	.00294910295364247,
@@ -88,5 +89,5 @@
 };
 
-constant double Gauss76Z[76]={
+constant double Gauss76Z[N_POINTS_76]={
 	-.999505948362153,		//0
 	-.997397786355355,
Index: sasmodels/models/lib/nr_bess_j1.c
===================================================================
--- sasmodels/models/lib/nr_bess_j1.c	(revision a8b3cdb60df4f2e22c9b3c5c3110bb51b7dba2c6)
+++ sasmodels/models/lib/nr_bess_j1.c	(revision a8b3cdb60df4f2e22c9b3c5c3110bb51b7dba2c6)
@@ -0,0 +1,28 @@
+double NR_BessJ1(double x);
+double NR_BessJ1(double x)
+{
+    double ax,z;
+    double xx,y,ans,ans1,ans2;
+
+    if ((ax=fabs(x)) < 8.0) {
+        y=x*x;
+        ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
+                                                  +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
+        ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
+                                               +y*(99447.43394+y*(376.9991397+y*1.0))));
+        ans=ans1/ans2;
+    } else {
+        z=8.0/ax;
+        y=z*z;
+        xx=ax-2.356194491;
+        ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
+                                   +y*(0.2457520174e-5+y*(-0.240337019e-6))));
+        ans2=0.04687499995+y*(-0.2002690873e-3
+                              +y*(0.8449199096e-5+y*(-0.88228987e-6
+                                                     +y*0.105787412e-6)));
+        ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
+        if (x < 0.0) ans = -ans;
+    }
+
+    return(ans);
+}
Index: sasmodels/models/line.py
===================================================================
--- sasmodels/models/line.py	(revision fa8011eb0e6b527e15b594945876ee39b6a46dbd)
+++ sasmodels/models/line.py	(revision fa8011eb0e6b527e15b594945876ee39b6a46dbd)
@@ -0,0 +1,107 @@
+r"""
+This model calculates intensity using simple linear function
+
+Definition
+----------
+
+The scattering intensity $I(q)$ is calculated as
+
+.. math::
+
+    I(q) = A + B \cdot q
+
+.. note::
+    For 2D plots intensity has different definition than other shape independent models
+
+.. math::
+    I(q) = I(qx) \cdot I(qy)
+
+References
+----------
+
+None.
+
+"""
+from numpy import inf
+
+name = "line"
+title = "Line model"
+description = """\
+      I(q) = A + B*q
+
+      List of default parameters:
+      A = intercept
+      B = slope
+      """
+category = "shape-independent"
+
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type", "description"],
+parameters = [["intercept", "1/cm",   1.0, [-inf, inf], "", "intercept in linear model"],
+              ["slope",     "Ang/cm", 1.0, [-inf, inf], "", "slope in linear model"],
+             ]
+# pylint: enable=bad-whitespace, line-too-long
+
+def Iq(q, intercept, slope):
+    """
+    :param q:           Input q-value
+    :param intercept:   Intrecept in linear model
+    :param slope:       Slope in linear model
+    :return:            Calculated Intensity
+    """
+    inten = intercept + slope*q
+    # TODO: In SasView code additional formula for list has been specified.
+    # if inten(x) = intercept + slope*x:
+    # then if q is a list, Iq=inten(x[0]*math.cos(x[1]))*inten(x[0]*math.sin(x[1]))
+    return inten
+
+Iq.vectorized = True # Iq accepts an array of q values
+
+def Iqxy(qx, qy, *args):
+    """
+    :param qx:   Input q_x-value
+    :param qy:   Input q_y-value
+    :param args: Remaining arguments
+    :return:     2D-Intensity
+    """
+    # TODO: SasView documention lists 2D intensity as Iq(qx)*Iq(qy) but code says:
+    # return self._line(x[1])
+    return Iq(qx, *args)*Iq(qy, *args)
+
+Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values
+
+demo = dict(scale=1.0, background=0, intercept=1.0, slope=1.0)
+
+oldname = "LineModel"
+oldpars = dict(intercept='A', slope='B', background=None, scale=None)
+
+tests = [
+
+    [{'intercept':   1.0,
+      'slope': 1.0,
+     }, 1.0, 2.0],
+
+    [{'intercept':   1.0,
+      'slope': 1.0,
+     }, 0.0, 1.0],
+
+    [{'intercept':   1.0,
+      'slope': 1.0,
+     }, 0.4, 1.4],
+
+    [{'intercept':   1.0,
+      'slope': 1.0,
+     }, 1.3, 2.3],
+
+    [{'intercept':   1.0,
+      'slope': 1.0,
+     }, 0.5, 1.5],
+
+    [{'intercept':   1.0,
+      'slope': 1.0,
+     }, [0.4, 0.5], [1.4, 1.5]],
+
+    [{'intercept':   1.0,
+      'slope': 1.0,
+     }, (1.3, 1.57), 5.911],
+]
Index: sasmodels/models/linear_pearls.py
===================================================================
--- sasmodels/models/linear_pearls.py	(revision 577912b439fe119a0663629c075b6b67323269ca)
+++ sasmodels/models/linear_pearls.py	(revision 32c743d8fe636333bd6fa7aa99f46c623b72890c)
Index: sasmodels/models/lorentz.py
===================================================================
--- sasmodels/models/lorentz.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/lorentz.py	(revision eb5901bce0a70fc35f2edf395a5b02199feab0cb)
@@ -56,5 +56,5 @@
 
 # parameters for demo
-demo = dict(scale=1.0,background=0.0,cor_length=50.0)
+demo = dict(scale=1.0, background=0.0, cor_length=50.0)
 
 # For testing against the old sasview models, include the converted parameter
@@ -64,5 +64,3 @@
 
 # parameters for unit tests
-tests = [
-         [{'cor_length' : 250},0.01,0.137931]
-         ]
+tests = [[{'cor_length': 250}, 0.01, 0.137931]]
Index: sasmodels/models/multi_shell.c
===================================================================
--- sasmodels/models/multi_shell.c	(revision d51ea74056d9370690a3247b3f61a63cf942f189)
+++ sasmodels/models/multi_shell.c	(revision d51ea74056d9370690a3247b3f61a63cf942f189)
@@ -0,0 +1,76 @@
+static
+double multi_shell_kernel(double q,
+          double core_radius,
+          double s_thickness,
+          double w_thickness,
+          double core_sld,
+          double shell_sld,
+          double n_pairs)
+{
+    //calculate with a loop, two shells at a time
+    int ii = 0;
+    double fval = 0.0;
+    double voli = 0.0;
+    const double sldi = core_sld-shell_sld;
+
+    do {
+        double ri = core_radius + (double)ii*(s_thickness + w_thickness);
+
+        // layer 1
+        voli = 4.0*M_PI/3.0*ri*ri*ri;
+        fval += voli*sldi*sph_j1c(ri*q);
+
+        ri += s_thickness;
+
+        // layer 2
+        voli = 4.0*M_PI/3.0*ri*ri*ri;
+        fval -= voli*sldi*sph_j1c(ri*q);
+
+        //do 2 layers at a time
+        ii += 1;
+
+    } while(ii <= n_pairs-1);  //change to make 0 < n_pairs < 2 correspond to
+                               //unilamellar vesicles (C. Glinka, 11/24/03)
+
+    fval *= 1.0e-4*fval/voli;
+
+    return(fval);
+}
+
+static
+double Iq(double q,
+          double core_radius,
+          double s_thickness,
+          double w_thickness,
+          double core_sld,
+          double shell_sld,
+          double n_pairs)
+{
+    return multi_shell_kernel(q,
+           core_radius,
+           s_thickness,
+           w_thickness,
+           core_sld,
+           shell_sld,
+           n_pairs);
+}
+
+static
+double Iqxy(double qx, double qy,
+          double core_radius,
+          double s_thickness,
+          double w_thickness,
+          double core_sld,
+          double shell_sld,
+          double n_pairs)
+{
+    double q = sqrt(qx*qx + qy*qy);
+
+    return multi_shell_kernel(q,
+           core_radius,
+           s_thickness,
+           w_thickness,
+           core_sld,
+           shell_sld,
+           n_pairs);
+}
Index: sasmodels/models/multi_shell.py
===================================================================
--- sasmodels/models/multi_shell.py	(revision d51ea74056d9370690a3247b3f61a63cf942f189)
+++ sasmodels/models/multi_shell.py	(revision d51ea74056d9370690a3247b3f61a63cf942f189)
@@ -0,0 +1,117 @@
+r"""
+This model provides the form factor, *P(q)*, for a multi-lamellar vesicle
+with *N* shells where the core is filled with solvent and the shells are
+interleaved with layers of solvent. For *N = 1*, this returns the VesicleModel.
+
+For information about polarised and magnetic scattering, click here_.
+
+.. _here: polar_mag_help.html
+
+Definition
+----------
+
+This model is a trivial extension of the CoreShell function to a larger number
+of shells. See the core_shell_ function for a diagram and documentation.
+
+.. _core_shell: core_shell_sphere.html
+
+Be careful! The SLDs and scale can be highly correlated. Hold as many of these
+parameters fixed as possible.
+
+.. figure:: img/multi_shell_fig1.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}
+
+.. note:
+    The outer most radius
+    $core_radius + n_pairs * s_thickness + (n_pairs - 1) * w_thickness$
+    is used as the effective radius for *S(Q)* when $P(Q) * S(Q)$ is applied.
+
+
+.. figure:: img/multi_shell_1d.jpg
+
+    1D plot using the default values (with 200 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
+---------
+B Cabane, *Small Angle Scattering Methods*,
+in *Surfactant Solutions: New Methods of Investigation*,
+Ch.2, Surfactant Science Series Vol. 22, Ed. R Zana and M Dekker,
+New York, (1987).
+
+"""
+
+from numpy import inf
+
+name = "multi_shell"
+title = "Multi shell model"
+description = """
+    MultiShell (Sphere) Model (or Multilamellar Vesicles): Model parameters;
+    scale : scale factor
+    core_radius : Core radius of the multishell
+    s_thickness: shell thickness
+    w_thickness: water thickness
+    core_sld: core scattering length density
+    shell_sld: shell scattering length density
+    n_pairs:number of pairs of water/shell
+    background: incoherent background
+        """
+category = "shape:sphere"
+
+# pylint: disable=bad-whitespace, line-too-long
+#   ["name", "units", default, [lower, upper], "type","description"],
+parameters = [
+    ["core_radius", "Ang",        60.0, [0.0, inf],  "", "Core radius of the multishell"],
+    ["s_thickness", "Ang",        10.0, [0.0, inf],  "", "Shell thickness"],
+    ["w_thickness", "Ang",        10.0, [0.0, inf],  "", "Water thickness"],
+    ["core_sld",    "1e-6/Ang^2",  6.4, [-inf, inf], "", "Core scattering length density"],
+    ["shell_sld",   "1e-6/Ang^2",  0.4, [-inf, inf], "", "Shell scattering length density"],
+    ["n_pairs",     "",            2.0, [1.0, inf],  "", "Number of pairs of water and shell"],
+    ]
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/sph_j1c.c", "multi_shell.c"]
+
+polydispersity = ["core_radius", "n_pairs"]
+
+demo = dict(scale=1, background=0,
+            core_radius=60.0,
+            s_thickness=10.0,
+            w_thickness=10.0,
+            core_sld=6.4,
+            shell_sld=0.4,
+            n_pairs=2.0)
+
+oldname = 'MultiShellModel'
+oldpars = dict()
+
+tests = [
+    # Accuracy tests based on content in test/utest_other_models.py
+    [{'core_radius': 60.0,
+      's_thickness': 10.0,
+      'w_thickness': 10.0,
+      'core_sld': 6.4,
+      'shell_sld': 0.4,
+      'n_pairs': 2.0,
+      'scale': 1.0,
+      'background': 0.001,
+     }, 0.001, 2442.81],
+
+    [{'core_radius': 60.0,
+      's_thickness': 10.0,
+      'w_thickness': 10.0,
+      'core_sld': 6.4,
+      'shell_sld': 0.4,
+      'n_pairs': 2.0,
+      'scale': 1.0,
+      'background': 0.001,
+     }, (0.001, 0.30903), 1.61873],
+    ]
Index: sasmodels/models/peak_lorentz.py
===================================================================
--- sasmodels/models/peak_lorentz.py	(revision 14ba6f663691e07e3ad9e93eadb16aedf8b2caa6)
+++ sasmodels/models/peak_lorentz.py	(revision 04b0b30a637f374651ae93ca42ae86785daeb735)
@@ -11,5 +11,6 @@
     I(q) = \frac{scale}{\bigl(1+\bigl(\frac{q-q_0}{B}\bigr)^2\bigr)} + background
 
-with the peak having height of $I_0$ centered at $q_0$ and having a HWHM (half-width half-maximum) of B.
+with the peak having height of $I_0$ centered at $q_0$ and having
+a HWHM (half-width half-maximum) of B.
 
 For 2D data the scattering intensity is calculated in the same way as 1D,
@@ -55,4 +56,7 @@
 
 def Iq(q, peak_pos, peak_hwhm):
+    """
+        Return I(q)
+    """
     inten = (1/(1+((q-peak_pos)/peak_hwhm)**2))
     return inten
@@ -60,4 +64,7 @@
 
 def Iqxy(qx, qy, *args):
+    """
+        Return I(qx, qy)
+    """
     return Iq(sqrt(qx ** 2 + qy ** 2), *args)
 Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
Index: sasmodels/models/porod.py
===================================================================
--- sasmodels/models/porod.py	(revision 1ce6e8231452d582e846afe023a22e987fe25c95)
+++ sasmodels/models/porod.py	(revision 1ce6e8231452d582e846afe023a22e987fe25c95)
@@ -0,0 +1,58 @@
+r"""
+This model fits the Porod function
+
+.. math::
+
+    I(q) = C/q^4
+    \\
+    C = 2\pi (\Delta\rho)^2 S_v
+
+to the data directly without any need for linearisation (cf. Log I(q) vs Log q).
+
+Here $C$ is the scale factor and $S_v$ is the specific surface area (ie, surface area / volume)
+of the sample, and $\Delta\rho$ is the contrast factor.
+
+For 2D data: 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}
+
+"""
+
+from numpy import sqrt, power
+
+name = "porod"
+title = "Porod function"
+description = """\
+          I(q) = scale/q^4 + background
+"""
+
+category = "shape-independent"
+
+parameters = []
+
+def Iq(q):
+    """
+    @param q: Input q-value
+    """
+    return 1.0/power(q, 4)
+
+Iq.vectorized = True  # Iq accepts an array of q values
+
+def Iqxy(qx, qy, *args):
+    """
+    @param qx:   Input q_x-value
+    @param qy:   Input q_y-value
+    @param args: Remaining arguments
+    """
+    return Iq(sqrt(qx ** 2 + qy ** 2), *args)
+
+Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
+
+demo = dict(scale=1.5, background=0.5)
+
+oldname = "PorodModel"
+oldpars = dict(scale='scale', background='background')
+
+tests = [[{'scale': 0.00001, 'background':0.01}, 0.04, 3.916250]]
Index: sasmodels/models/rpa.py
===================================================================
--- sasmodels/models/rpa.py	(revision 82c299fa51ab091ee8ead0a235d4177ccfeecb99)
+++ sasmodels/models/rpa.py	(revision fa8011eb0e6b527e15b594945876ee39b6a46dbd)
@@ -43,5 +43,5 @@
 component.
 
-.. figure:: img/image215.jpg
+.. figure:: img/rpa_1d.jpg
 
     1D plot using the default values (w/500 data points).
@@ -89,5 +89,5 @@
 #   ["name", "units", default, [lower, upper], "type","description"],
 parameters = [
-    ["case_num", CASES, 0, [0, 10], "", "Component organization" ],
+    ["case_num", CASES, 0, [0, 10], "", "Component organization"],
 
     ["Na", "", 1000.0, [1, inf], "", "Degree of polymerization"],
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
+++ sasmodels/models/sphere.py	(revision c6915512ab98df810af8905ba552572737ceb9fa)
@@ -100,4 +100,7 @@
 
 def ER(radius):
+    """
+        Return equivalent radius (ER)
+    """
     return radius
 
Index: sasmodels/models/spherepy.py
===================================================================
--- sasmodels/models/spherepy.py	(revision eb69cce00f0781ad8f3eec0eb5965a8576d815b5)
+++ sasmodels/models/spherepy.py	(revision d2950f4360cecbd3a71cc4eb7a857316a173acf4)
@@ -118,5 +118,5 @@
     g[low] = sqrt(1 - dlow2 / 4.) * (1 + dlow2 / 8.) + dlow2 / 2.*(1 - dlow2 / 16.) * log(dlow / (2. + sqrt(4. - dlow2)))
     return g
-sesans.vectorized = True  # sesans accepts and array of z values
+sesans.vectorized = True  # sesans accepts an array of z values
 
 def ER(radius):
Index: sasmodels/models/stacked_disks.c
===================================================================
--- sasmodels/models/stacked_disks.c	(revision d507c3a946bb1b2103ca169efe83263b3955b43d)
+++ sasmodels/models/stacked_disks.c	(revision d507c3a946bb1b2103ca169efe83263b3955b43d)
@@ -0,0 +1,244 @@
+double form_volume(double core_thick,
+                   double layer_thick,
+                   double radius,
+                   double n_stacking);
+
+double Iq(double q,
+          double core_thick,
+          double layer_thick,
+          double radius,
+          double n_stacking,
+          double sigma_d,
+          double core_sld,
+          double layer_sld,
+          double solvent_sld);
+
+double Iqxy(double qx, double qy,
+            double core_thick,
+            double layer_thick,
+            double radius,
+            double n_stacking,
+            double sigma_d,
+            double core_sld,
+            double layer_sld,
+            double solvent_sld,
+            double theta,
+            double phi);
+
+static
+double _kernel(double qq,
+               double radius,
+               double core_sld,
+               double layer_sld,
+               double solvent_sld,
+               double halfheight,
+               double layer_thick,
+               double zi,
+               double sigma_d,
+               double d,
+               double n_stacking)
+
+{
+	// qq is the q-value for the calculation (1/A)
+	// radius is the core radius of the cylinder (A)
+	// *_sld are the respective SLD's
+	// halfheight is the *Half* CORE-LENGTH of the cylinder = L (A)
+	// zi is the dummy variable for the integration (x in Feigin's notation)
+
+	const double besarg1 = qq*radius*sin(zi);
+	const double besarg2 = qq*radius*sin(zi);
+
+	const double sinarg1 = qq*halfheight*cos(zi);
+	const double sinarg2 = qq*(halfheight+layer_thick)*cos(zi);
+
+    const double be1 = J1(besarg1)/besarg1;
+	const double be2 = J1(besarg2)/besarg2;
+	const double si1 = sin(sinarg1)/sinarg1;
+	const double si2 = sin(sinarg2)/sinarg2;
+
+	const double dr1 = (core_sld-solvent_sld);
+	const double dr2 = (layer_sld-solvent_sld);
+	const double area = M_PI*radius*radius;
+	const double totald=2.0*(layer_thick+halfheight);
+
+	const double t1 = 2.0*area*(2.0*halfheight)*dr1*(si1)*(be1);
+	const double t2 = 2.0*area*dr2*(totald*si2-2.0*halfheight*si1)*(be2);
+
+
+	double retval =((t1+t2)*(t1+t2))*sin(zi);
+
+	// loop for the structure facture S(q)
+	double sqq=0.0;
+	for(int kk=1;kk<n_stacking;kk+=1) {
+		double dexpt=qq*cos(zi)*qq*cos(zi)*d*d*sigma_d*sigma_d*kk/2.0;
+		sqq=sqq+(n_stacking-kk)*cos(qq*cos(zi)*d*kk)*exp(-1.*dexpt);
+	}
+
+	// end of loop for S(q)
+	sqq=1.0+2.0*sqq/n_stacking;
+
+	retval *= sqq;
+
+	return(retval);
+}
+
+
+static
+double stacked_disks_kernel(double q,
+                            double core_thick,
+                            double layer_thick,
+                            double radius,
+                            double n_stacking,
+                            double sigma_d,
+                            double core_sld,
+                            double layer_sld,
+                            double solvent_sld)
+{
+/*	StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks
+like clay platelets that are not exfoliated
+*/
+	double summ = 0.0;	//initialize integral
+
+	double d=2.0*layer_thick+core_thick;
+	double halfheight = core_thick/2.0;
+
+	for(int i=0;i<N_POINTS_76;i++) {
+		double zi = (Gauss76Z[i] + 1.0)*M_PI/4.0;
+		double yyy = Gauss76Wt[i] *
+                    _kernel(q,
+		                   radius,
+		                   core_sld,
+		                   layer_sld,
+		                   solvent_sld,
+		                   halfheight,
+		                   layer_thick,
+		                   zi,
+		                   sigma_d,
+		                   d,
+		                   n_stacking);
+		summ += yyy;
+	}
+
+	double answer = M_PI/4.0*summ;
+
+	//Convert to [cm-1]
+	answer *= 1.0e-4;
+
+	return answer;
+}
+
+static double stacked_disks_kernel_2d(double q, double q_x, double q_y,
+                            double core_thick,
+                            double layer_thick,
+                            double radius,
+                            double n_stacking,
+                            double sigma_d,
+                            double core_sld,
+                            double layer_sld,
+                            double solvent_sld,
+                            double theta,
+                            double phi)
+{
+
+    double ct, st, cp, sp;
+
+    //convert angle degree to radian
+    theta = theta * M_PI/180.0;
+    phi = phi * M_PI/180.0;
+
+    SINCOS(theta, st, ct);
+    SINCOS(phi, sp, cp);
+
+    // silence compiler warnings about unused variable
+    (void) sp;
+
+    // parallelepiped orientation
+    const double cyl_x = ct * cp;
+    const double cyl_y = st;
+
+    // Compute the angle btw vector q and the
+    // axis of the parallelepiped
+    const double cos_val = cyl_x*q_x + cyl_y*q_y;
+
+    // Note: cos(alpha) = 0 and 1 will get an
+    // undefined value from Stackdisc_kern
+    double alpha = acos( cos_val );
+
+    // Call the IGOR library function to get the kernel
+    double d = 2 * layer_thick + core_thick;
+    double halfheight = core_thick/2.0;
+    double answer = _kernel(q,
+                     radius,
+                     core_sld,
+                     layer_sld,
+                     solvent_sld,
+                     halfheight,
+                     layer_thick,
+                     alpha,
+                     sigma_d,
+                     d,
+                     n_stacking);
+
+    answer /= sin(alpha);
+    //convert to [cm-1]
+    answer *= 1.0e-4;
+
+    return answer;
+}
+
+double form_volume(double core_thick,
+                   double layer_thick,
+                   double radius,
+                   double n_stacking){
+    double d = 2 * layer_thick + core_thick;
+    return acos(-1.0) * radius * radius * d * n_stacking;
+}
+
+double Iq(double q,
+          double core_thick,
+          double layer_thick,
+          double radius,
+          double n_stacking,
+          double sigma_d,
+          double core_sld,
+          double layer_sld,
+          double solvent_sld)
+{
+    return stacked_disks_kernel(q,
+                    core_thick,
+                    layer_thick,
+                    radius,
+                    n_stacking,
+                    sigma_d,
+                    core_sld,
+                    layer_sld,
+                    solvent_sld);
+}
+
+// Iqxy is never called since no orientation or magnetic parameters.
+double Iqxy(double qx, double qy,
+            double core_thick,
+            double layer_thick,
+            double radius,
+            double n_stacking,
+            double sigma_d,
+            double core_sld,
+            double layer_sld,
+            double solvent_sld,
+            double theta,
+            double phi)
+{
+    double q = sqrt(qx*qx + qy*qy);
+    return stacked_disks_kernel_2d(q, qx/q, qy/q,
+                    core_thick,
+                    layer_thick,
+                    radius,
+                    n_stacking,
+                    sigma_d,
+                    core_sld,
+                    layer_sld,
+                    solvent_sld,
+                    theta,
+                    phi);
+}
+
Index: sasmodels/models/stacked_disks.py
===================================================================
--- sasmodels/models/stacked_disks.py	(revision d507c3a946bb1b2103ca169efe83263b3955b43d)
+++ sasmodels/models/stacked_disks.py	(revision d507c3a946bb1b2103ca169efe83263b3955b43d)
@@ -0,0 +1,220 @@
+r"""
+This model provides the form factor, $P(q)$, for stacked discs (tactoids)
+with a core/layer structure where the form factor is normalized by the volume
+of the cylinder. Assuming the next neighbor distance (d-spacing) in a stack
+of parallel discs obeys a Gaussian distribution, a structure factor $S(q)$
+proposed by Kratky and Porod in 1949 is used in this function.
+
+Note that the resolution smearing calculation uses 76 Gauss quadrature points
+to properly smear the model since the function is HIGHLY oscillatory,
+especially around the q-values that correspond to the repeat distance of
+the layers.
+
+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}
+
+Definition
+----------
+
+.. figure:: img/stacked_disks_fig1.gif
+
+The scattered intensity $I(q)$ is calculated as
+
+.. math::
+
+    I(q) = N\int_{0}^{\pi /2}\left[ \Delta \rho_t
+    \left( V_tf_t(q) - V_cf_c(q)\right) + \Delta \rho_cV_cf_c(q)
+    \right]^2S(q)\sin{\alpha d\alpha} + background
+
+where the contrast
+
+.. math::
+
+    \Delta \rho_i = \rho_i - \rho_{solvent}
+
+and *N* is the number of discs per unit volume,
+$\alpha$ is the angle between the axis of the disc and *q*,
+and $V_t$ and $V_c$ are the total volume and the core volume of
+a single disc, respectively.
+
+.. math::
+
+    \left\langle f_{t}^2(q)\right\rangle_{\alpha} =
+    \int_{0}^{\pi/2}\left[
+    \left(\frac{sin(q(d+h)\cos{\alpha})}{q(d+h)\cos{\alpha}}\right)
+    \left(\frac{2J_1(qR\sin{\alpha})}{qR\sin{\alpha}} \right)
+    \right]^2 \sin{\alpha d\alpha}
+
+    \left\langle f_{c}^2(q)\right\rangle_{\alpha} =
+    \int_{0}^{\pi/2}\left[
+    \left(\frac{sin(qh)\cos{\alpha})}{qh\cos{\alpha}}\right)
+    \left(\frac{2J_1(qR\sin{\alpha})}{qR\sin{\alpha}} \right)
+    \right]^2 \sin{\alpha d\alpha}
+
+where *d* = thickness of the layer (layer_thick),
+*2h* = core thickness (core_thick), and *R* = radius of the disc (radius).
+
+.. math::
+
+    S(q) = 1 + \frac{1}{2}\sum_{k=1}^n(n-k)\cos{(kDq\cos{\alpha})}
+    exp\left[ -k(q\cos{\alpha})^2\sigma_D/2\right]
+
+where *n* = the total number of the disc stacked (n_stacking),
+*D* = the next neighbor center-to-center distance (d-spacing),
+and $\sigma_D$ = the Gaussian standard deviation of the d-spacing (sigma_d).
+
+To provide easy access to the orientation of the stacked disks, we define
+the axis of the cylinder using two angles $\theta$ and $\varphi$.
+These angles are defined on Figure 2 of cylinder_model.
+
+.. note::
+    The 2nd virial coefficient of the cylinder is calculated based on the
+    *radius* and *length* = *n_stacking* * (*core_thick* + 2 * *layer_thick*)
+    values, and used as the effective radius for $S(Q)$ when $P(Q) * S(Q)$
+    is applied.
+
+.. figure:: img/stacked_disks_1d.jpg
+
+    1D plot using the default values (w/1000 data point).
+
+.. figure:: img/stacked_disks_fig2.jpg
+
+    Examples of the angles for oriented stacked disks against
+    the detector plane.
+
+.. figure:: img/stacked_disks_fig3.jpg
+
+    Examples of the angles for oriented pp against the detector plane.
+
+Our model uses the form factor calculations implemented in a c-library provided
+by the NIST Center for Neutron Research (Kline, 2006)
+
+Reference
+---------
+
+A Guinier and G Fournet, *Small-Angle Scattering of X-Rays*, John Wiley and Sons, New York, 1955
+
+O Kratky and G Porod, *J. Colloid Science*, 4, (1949) 35
+
+J S Higgins and H C Benoit, *Polymers and Neutron Scattering*, Clarendon, Oxford, 1994
+
+"""
+
+from numpy import inf
+
+name = "stacked_disks"
+title = ""
+description = """\
+    One layer of disk consists of a core, a top layer, and a bottom layer.
+    radius =  the radius of the disk
+    core_thick = thickness of the core
+    layer_thick = thickness of a layer
+    core_sld = the SLD of the core
+    layer_sld = the SLD of the layers
+    n_stacking = the number of the disks
+    sigma_d =  Gaussian STD of d-spacing
+    solvent_sld = the SLD of the solvent
+    """
+category = "shape:cylinder"
+
+# pylint: disable=bad-whitespace, line-too-long
+#   ["name", "units", default, [lower, upper], "type","description"],
+parameters = [
+    ["core_thick",  "Ang",        10.0, [0, inf],    "volume",      "Thickness of the core disk"],
+    ["layer_thick", "Ang",        10.0, [0, inf],    "volume",      "Thickness of the stacked disk"],
+    ["radius",      "Ang",        15.0, [0, inf],    "volume",      "Radius of the stacked disk"],
+    ["n_stacking",  "",            1.0, [0, inf],    "volume",      "Number of stacking"],
+    ["sigma_d",     "Ang",         0,   [0, inf],    "",            "GSD of disks sigma_d"],
+    ["core_sld",    "1e-6/Ang^2",  4,   [-inf, inf], "",            "Core scattering length density"],
+    ["layer_sld",   "1e-6/Ang^2",  0.0, [-inf, inf], "",            "Layer scattering length density"],
+    ["solvent_sld", "1e-6/Ang^2",  5.0, [-inf, inf], "",            "Solvent scattering length density"],
+    ["theta",       "degrees",     0,   [-inf, inf], "orientation", "Orientation of the stacked disk axis w/respect incoming beam"],
+    ["phi",         "degrees",     0,   [-inf, inf], "orientation", "Orientation of the stacked disk in the plane of the detector"],
+    ]
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/gauss76.c", "lib/J1.c", "stacked_disks.c"]
+
+demo = dict(background=0.001,
+            scale=0.01,
+            core_thick=10.0,
+            layer_thick=10.0,
+            radius=15.0,
+            n_stacking=1,
+            sigma_d=0,
+            core_sld=4,
+            layer_sld=0.0,
+            solvent_sld=5.0,
+            theta=0,
+            phi=0)
+
+
+oldname = 'StackedDisksModel'
+
+oldpars = dict(theta='axis_theta',
+               phi='axis_phi')
+
+tests = [
+    # Accuracy tests based on content in test/utest_extra_models.py
+    [{'core_thick': 10.0,
+      'layer_thick': 15.0,
+      'radius': 3000.0,
+      'n_stacking': 1.0,
+      'sigma_d': 0.0,
+      'core_sld': 4.0,
+      'layer_sld': -0.4,
+      'solvent_sd': 5.0,
+      'theta': 0.0,
+      'phi': 0.0,
+      'scale': 0.01,
+      'background': 0.001,
+     }, 0.001, 5075.12],
+
+    [{'core_thick': 10.0,
+      'layer_thick': 15.0,
+      'radius': 3000.0,
+      'n_stacking': 1.0,
+      'sigma_d': 0.0,
+      'core_sld': 4.0,
+      'layer_sld': -0.4,
+      'solvent_sd': 5.0,
+      'theta': 0.0,
+      'phi': 0.0,
+      'scale': 0.01,
+      'background': 0.001,
+     }, [0.001, 90.0], [5075.12, 0.001]],
+
+    [{'core_thick': 10.0,
+      'layer_thick': 15.0,
+      'radius': 3000.0,
+      'n_stacking': 1.0,
+      'sigma_d': 0.0,
+      'core_sld': 4.0,
+      'layer_sld': -0.4,
+      'solvent_sd': 5.0,
+      'theta': 0.0,
+      'phi': 0.0,
+      'scale': 0.01,
+      'background': 0.001,
+     }, ([0.4, 0.5]), [0.00105074, 0.00121761]],
+
+    [{'core_thick': 10.0,
+      'layer_thick': 15.0,
+      'radius': 3000.0,
+      'n_stacking': 1.0,
+      'sigma_d': 0.0,
+      'core_sld': 4.0,
+      'layer_sld': -0.4,
+      'solvent_sd': 5.0,
+      'theta': 0.0,
+      'phi': 0.0,
+      'scale': 0.01,
+      'background': 0.001,
+     }, ([1.3, 1.57]), [0.0010039, 0.0010038]],
+    ]
+
+
Index: sasmodels/models/triaxial_ellipsoid.py
===================================================================
--- sasmodels/models/triaxial_ellipsoid.py	(revision 9c461c713e9961727442c844be3e6354d6548c2f)
+++ sasmodels/models/triaxial_ellipsoid.py	(revision 469e763f1e5fc119d32baa75e7d38058e275d079)
@@ -82,6 +82,6 @@
 ----------
 
-L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum,
-New York, 1987.
+L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray
+and Neutron Scattering*, Plenum, New York, 1987.
 """
 
@@ -120,4 +120,7 @@
 
 def ER(req_minor, req_major, rpolar):
+    """
+        Returns the effective radius used in the S*P calculation
+    """
     import numpy as np
     from .ellipsoid import ER as ellipsoid_ER
Index: sasmodels/models/two_power_law.py
===================================================================
--- sasmodels/models/two_power_law.py	(revision 7e1d090d6f1579784a2df1a534eda3abdc9c2d5a)
+++ sasmodels/models/two_power_law.py	(revision 7e1d090d6f1579784a2df1a534eda3abdc9c2d5a)
@@ -0,0 +1,163 @@
+r"""
+This model calculates an empirical functional form for SAS data characterized
+by two power laws.
+
+Definition
+----------
+
+The scattering intensity $I(q)$ is calculated as
+
+.. math::
+
+    I(q) = \begin{cases}
+    A q^{-m1} + \text{background} & q <= qc \\
+    C q^{-m2} + \text{background} & q > qc
+    \end{cases}
+
+where $qc$ = the location of the crossover from one slope to the other,
+$A$ = the scaling coefficent that sets the overall intensity of the lower Q power law region,
+$m1$ = power law exponent at low Q,
+$m2$ = power law exponent at high Q
+
+The scaling of the second power law region (coefficent C) is then automatically scaled to
+match the first by following formula
+
+.. math::
+    C = \frac{A}{qc^{-m1} qc^{-m2}}
+
+.. note::
+    Be sure to enter the power law exponents as positive values!
+
+For 2D data the 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}
+
+
+.. figure:: img/two_power_law_1d.jpg
+
+    1D plot using the default values (with 500 data point).
+
+References
+----------
+
+None.
+
+"""
+
+from numpy import power
+from numpy import sqrt
+from numpy import inf
+from numpy import concatenate
+
+name = "two_power_law"
+title = "Two Power Law intensity calculation"
+description = """
+            I(q) = coef_A*pow(qval,-1.0*power1) + background for q<=qc
+            =C*pow(qval,-1.0*power2) + background for q>qc
+            where C=coef_A*pow(qc,-1.0*power1)/pow(qc,-1.0*power2).
+
+            coef_A = scaling coefficent
+            qc = crossover location [1/A]
+            power_1 (=m1) = power law exponent at low Q
+            power_2 (=m2) = power law exponent at high Q
+            background = Incoherent background [1/cm]
+        """
+category = "shape-independent"
+
+# pylint: disable=bad-whitespace, line-too-long
+#            ["name", "units", default, [lower, upper], "type", "description"],
+parameters = [["coefficent_1",  "",         1.0, [-inf, inf], "", "coefficent A in low Q region"],
+              ["crossover",     "1/Ang",    0.04,[0, inf],    "", "crossover location"],
+              ["power_1",       "",         1.0, [0, inf],    "", "power law exponent at low Q"],
+              ["power_2",       "",         4.0, [0, inf],    "", "power law exponent at high Q"],
+             ]
+# pylint: enable=bad-whitespace, line-too-long
+
+
+def Iq(q,
+       coefficent_1=1.0,
+       crossover=0.04,
+       power_1=1.0,
+       power_2=4.0,
+      ):
+
+    """
+    :param q:                   Input q-value (float or [float, float])
+    :param coefficent_1:        Scaling coefficent at low Q
+    :param crossover:           Crossover location
+    :param power_1:             Exponent of power law function at low Q
+    :param power_2:             Exponent of power law function at high Q
+    :return:                    Calculated intensity
+    """
+
+    #Two sub vectors are created to treat crossover values
+    q_lower = q[q <= crossover]
+    q_upper = q[q > crossover]
+    coefficent_2 = coefficent_1*power(crossover, -1.0*power_1)/power(crossover, -1.0*power_2)
+    intensity_lower = coefficent_1*power(q_lower, -1.0*power_1)
+    intensity_upper = coefficent_2*power(q_upper, -1.0*power_2)
+    intensity = concatenate((intensity_lower, intensity_upper), axis=0)
+
+    return intensity
+
+Iq.vectorized = True  # Iq accepts an array of q values
+
+def Iqxy(qx, qy, *args):
+    """
+    :param qx:   Input q_x-value
+    :param qy:   Input q_y-value
+    :param args: Remaining arguments
+    :return:     2D-Intensity
+    """
+
+    return Iq(sqrt(qx**2 + qy**2), *args)
+
+Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values
+
+demo = dict(scale=1, background=0.0,
+            coefficent_1=1.0,
+            crossover=0.04,
+            power_1=1.0,
+            power_2=4.0)
+
+oldname = "TwoPowerLawModel"
+oldpars = dict(coefficent_1='coef_A',
+               crossover='qc',
+               power_1='power1',
+               power_2='power2',
+               background='background')
+
+tests = [
+    # Accuracy tests based on content in test/utest_extra_models.py
+    [{'coefficent_1':     1.0,
+      'crossover':  0.04,
+      'power_1':    1.0,
+      'power_2':    4.0,
+      'background': 0.0,
+     }, 0.001, 1000],
+
+    [{'coefficent_1':     1.0,
+      'crossover':  0.04,
+      'power_1':    1.0,
+      'power_2':    4.0,
+      'background': 0.0,
+     }, 0.150141, 0.125945],
+
+    [{'coeffcent_1':    1.0,
+      'crossover':  0.04,
+      'power_1':    1.0,
+      'power_2':    4.0,
+      'background': 0.0,
+     }, 0.442528, 0.00166884],
+
+    [{'coeffcent_1':    1.0,
+      'crossover':  0.04,
+      'power_1':    1.0,
+      'power_2':    4.0,
+      'background': 0.0,
+     }, (0.442528, 0.00166884), 0.00166884],
+
+]
Index: sasmodels/models/vesicle.py
===================================================================
--- sasmodels/models/vesicle.py	(revision 216fa6d00e42770a4a028c2953a73319ad05e924)
+++ sasmodels/models/vesicle.py	(revision fa8011eb0e6b527e15b594945876ee39b6a46dbd)
@@ -24,5 +24,5 @@
 is a flat background level (due for example to incoherent scattering in the
 case of neutrons), and $j_1$ is the spherical bessel function
-$j_1 = (sin(x) - x cos(x))/ x^2$.
+$j_1 = (\sin(x) - x \cos(x))/ x^2$.
 
 The functional form is identical to a "typical" core-shell structure, except
@@ -35,5 +35,7 @@
 thickness = $R_{\text{tot}} - R_{\text{core}}$.
 
-.. figure: img/vesicle_geometry.jpg
+.. figure:: img/vesicle_geometry.jpg
+
+    Vesicle geometry.
 
 The 2D scattering intensity is the same as *P(q)* above, regardless of the
@@ -48,9 +50,9 @@
 radius for *S(Q)* when *P(Q)* \* *S(Q)* is applied.
 
-.. image:: img/vesicle_1d.jpg
+.. figure:: img/vesicle_1d.jpg
 
-*Figure. 1D plot using the default values given in the table
-(w/200 data point). Polydispersity and instrumental resolution normally
-will smear out most of the rapidly oscillating features.*
+    1D plot using the default values given in the table (w/200 data point).
+    Polydispersity and instrumental resolution normally will smear out most
+    of the rapidly oscillating features.
 
 REFERENCE
@@ -60,5 +62,4 @@
 """
 
-import numpy as np
 from numpy import pi, inf
 
@@ -131,5 +132,5 @@
 # NOTE: test results taken from values returned by SasView 3.1.2
 tests = [[{}, 0.0010005303255, 17139.8268799],
-         [{}, 0.200027832249, 0.130387268704 ],
+         [{}, 0.200027832249, 0.130387268704],
          [{}, 'ER', 130.],
          [{}, 'VR', 0.54483386436],
Index: setup.py
===================================================================
--- setup.py	(revision 040575fa2fdfdbba541ff27f43057fb6c1354010)
+++ setup.py	(revision 3eb331258188803a81fbe7efcc02852d2c6e3e39)
@@ -11,5 +11,5 @@
     version = "1.0.0a",
     description = "sasmodels package",
-    long_description=open('README.rst').read(),
+    long_description=open('README.md').read(),
     author = "SasView Collaboration",
     author_email = "management@sasview.org",
