Index: sasmodels/models/flexible_cylinder_ex.c
===================================================================
--- sasmodels/models/flexible_cylinder_ex.c	(revision 504abeee8bc5d67a14d7598227542961d3d748ed)
+++ sasmodels/models/flexible_cylinder_ex.c	(revision 504abeee8bc5d67a14d7598227542961d3d748ed)
@@ -0,0 +1,103 @@
+double form_volume(double length, double kuhn_length, double radius);
+double Iq(double q, double length, double kuhn_length, double radius,
+          double axis_ratio, double sld, double solvent_sld);
+double Iqxy(double qx, double qy, double length, double kuhn_length,
+            double radius, double axis_ratio, double sld, double solvent_sld);
+double flexible_cylinder_ex_kernel(double q, double length, double kuhn_length,
+                                double radius, double axis_ratio, double sld,
+                                double solvent_sld);
+double elliptical_crosssection(double q, double a, double b);
+
+double form_volume(double length, double kuhn_length, double radius)
+{
+    return 1.0;
+}
+
+double
+elliptical_crosssection(double q, double a, double b)
+{
+    double uplim,lolim,Pi,summ,arg,zi,yyy,answer;
+    int i,nord=76;
+
+    Pi = 4.0*atan(1.0);
+    lolim=0.0;
+    uplim=Pi/2.0;
+    summ=0.0;
+
+    for(i=0;i<nord;i++) {
+		zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
+		arg = q*sqrt(a*a*sin(zi)*sin(zi)+b*b*cos(zi)*cos(zi));
+		yyy = pow((2.0 * J1(arg) / arg),2);
+		yyy *= Gauss76Wt[i];
+		summ += yyy;
+    }
+    answer = (uplim-lolim)/2.0*summ;
+    answer *= 2.0/Pi;
+    return(answer);
+
+}
+
+double flexible_cylinder_ex_kernel(double q,
+          double length,
+          double kuhn_length,
+          double radius,
+          double axis_ratio,
+          double sld,
+          double solvent_sld)
+{
+
+	double Pi,flex,crossSect, cont;
+
+	Pi = 4.0*atan(1.0);
+	cont = sld - solvent_sld;
+	crossSect = elliptical_crosssection(q,radius,(radius*axis_ratio));
+
+	flex = Sk_WR(q,length,kuhn_length);
+	flex *= crossSect;
+	flex *= Pi*radius*radius*axis_ratio*axis_ratio*length;
+	flex *= cont*cont;
+	flex *= 1.0e-4;
+
+	return flex;
+}
+
+double Iq(double q,
+          double length,
+          double kuhn_length,
+          double radius,
+          double axis_ratio,
+          double sld,
+          double solvent_sld)
+{
+
+	double result = flexible_cylinder_ex_kernel(q,
+	                length,
+	                kuhn_length,
+	                radius,
+	                axis_ratio,
+	                sld,
+	                solvent_sld);
+
+	return result;
+}
+
+double Iqxy(double qx, double qy,
+            double length,
+            double kuhn_length,
+            double radius,
+            double axis_ratio,
+            double sld,
+            double solvent_sld)
+{
+	double q;
+	q = sqrt(qx*qx+qy*qy);
+	double result = flexible_cylinder_ex_kernel(q,
+	                length,
+	                kuhn_length,
+	                radius,
+	                axis_ratio,
+	                sld,
+	                solvent_sld);
+
+	return result;
+}
Index: sasmodels/models/flexible_cylinder_ex.py
===================================================================
--- sasmodels/models/flexible_cylinder_ex.py	(revision 504abeee8bc5d67a14d7598227542961d3d748ed)
+++ sasmodels/models/flexible_cylinder_ex.py	(revision 504abeee8bc5d67a14d7598227542961d3d748ed)
@@ -0,0 +1,166 @@
+r"""
+This model calculates the form factor for a flexible cylinder with an
+elliptical cross section and a uniform scattering length density.
+The non-negligible diameter of the cylinder is included by accounting
+for excluded volume interactions within the walk of a single cylinder.
+The form factor is normalized by the particle volume such that
+
+.. math::
+
+    P(q) = \text{scale} \left<F^2\right>/V + \text{background}
+
+where the averaging $\left<\ldots\right>$ is over all possible orientations
+of the flexible cylinder.
+
+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}
+
+
+Definitions
+-----------
+
+The function calculated is from the reference given below. From that paper,
+"Method 3 With Excluded Volume" is used.
+The model is a parameterization of simulations of a discrete representation of
+the worm-like chain model of Kratky and Porod applied in the pseudo-continuous
+limit. See equations (13, 26-27) in the original reference for the details.
+
+.. note::
+
+    There are several typos in the original reference that have been corrected
+    by WRC. Details of the corrections are in the reference below. Most notably
+
+    - Equation (13): the term $(1 - w(QR))$ should swap position with $w(QR)$
+
+    - Equations (23) and (24) are incorrect; WRC has entered these into
+      Mathematica and solved analytically. The results were then converted to
+      code.
+
+    - Equation (27) should be $q0 = max(a3/sqrt(RgSquare),3)$ instead of
+      $max(a3*b/sqrt(RgSquare),3)$
+
+    - The scattering function is negative for a range of parameter values and
+      q-values that are experimentally accessible. A correction function has been
+      added to give the proper behavior.
+
+.. figure:: img/flexible_cylinder_ex_geometry.jpg
+
+
+The chain of contour length, $L$, (the total length) can be described as a chain
+of some number of locally stiff segments of length $l_p$, the persistence length
+(the length along the cylinder over which the flexible cylinder can be considered
+a rigid rod).
+The Kuhn length $(b = 2*l_p)$ is also used to describe the stiffness of a chain.
+
+The cross section of the cylinder is elliptical, with minor radius $a$ .
+The major radius is larger, so of course, **the axis ratio (parameter 5) must be
+greater than one.** Simple constraints should be applied during curve fitting to
+maintain this inequality.
+
+The returned value is in units of $cm^-1$, on absolute scale.
+
+In the parameters, the $sldCyl$ and $sldSolv$ represent the SLD of the
+chain/cylinder and solvent respectively. The *scale*, and the contrast are both
+multiplicative factors in the model and are perfectly correlated. One or both of
+these parameters must be held fixed during model fitting.
+
+**No inter-cylinder interference effects are included in this calculation.**
+
+.. figure:: img/flexible_cylinder_ex_1d.jpg
+
+    1D plot using the default values (w/1000 data point).
+
+References
+----------
+
+J S Pedersen and P Schurtenberger. *Scattering functions of semiflexible polymers
+with and without excluded volume effects.* Macromolecules, 29 (1996) 7602-7612
+
+Correction of the formula can be found in
+
+W R Chen, P D Butler and L J Magid, *Incorporating Intermicellar Interactions in
+the Fitting of SANS Data from Cationic Wormlike Micelles.* Langmuir,
+22(15) 2006 6539-6548
+"""
+from numpy import inf
+
+name = "flexible_cylinder_ex"
+title = "Flexible cylinder wth an elliptical cross section and a uniform " \
+        "scattering length density."
+description = """Note : scale and contrast=sldCyl-sldSolv are both multiplicative
+        factors in the
+        model and are perfectly correlated. One or
+        both of these parameters must be held fixed
+        during model fitting.
+        """
+
+category = "shape:cylinder"
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type", "description"],
+parameters = [
+    ["length",      "Ang",       1000.0, [0, inf],    "volume", "Length of the flexible cylinder"],
+    ["kuhn_length", "Ang",        100.0, [0, inf],    "volume", "Kuhn length of the flexible cylinder"],
+    ["radius",      "Ang",         20.0, [0, inf],    "volume", "Radius of the flexible cylinder"],
+    ["axis_ratio",  "",             1.5, [0, inf],    "",       "Axis_ratio (major_radius/radius"],
+    ["sld",         "1e-6/Ang^2",   1.0, [-inf, inf], "",       "Cylinder scattering length density"],
+    ["solvent_sld", "1e-6/Ang^2",   6.3, [-inf, inf], "",       "Solvent scattering length density"],
+    ]
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/J1.c", "lib/gauss76.c", "lib/wrc_cyl.c", "flexible_cylinder_ex.c"]
+
+demo = dict(scale=1.0, background=0.0001,
+            length=1000.0,
+            kuhn_length=100.0,
+            radius=20.0,
+            axis_ratio=1.5,
+            sld=1.0,
+            solvent_sld=6.3)
+
+oldname = 'FlexCylEllipXModel'
+oldpars = dict(sld='sldCyl', solvent_sld='sldSolv')
+
+
+tests = [
+    # Accuracy tests based on content in test/utest_other_models.py
+    # Currently fails in OCL
+    # [{'length':     1000.0,
+    #  'kuhn_length': 100.0,
+    #  'radius':       20.0,
+    #  'axis_ratio':    1.5,
+    #  'sld':           1.0,
+    #  'solvent_sld':   6.3,
+    #  'background':    0.0001,
+    # }, 0.001, 3509.2187],
+
+    # Additional tests with larger range of parameters
+    [{'length':     1000.0,
+      'kuhn_length': 100.0,
+      'radius':       20.0,
+      'axis_ratio':    1.5,
+      'sld':           1.0,
+      'solvent_sld':   6.3,
+      'background':    0.0001,
+     }, 1.0, 0.00223819],
+    [{'length':        10.0,
+      'kuhn_length': 800.0,
+      'radius':        2.0,
+      'axis_ratio':    0.5,
+      'sld':           6.0,
+      'solvent_sld':  12.3,
+      'background':    0.001,
+     }, 0.1, 0.390281],
+    [{'length':        100.0,
+      'kuhn_length': 800.0,
+      'radius':       50.0,
+      'axis_ratio':    4.5,
+      'sld':           0.1,
+      'solvent_sld':   5.1,
+      'background':    0.0,
+     }, 1.0, 0.0016338264790]
+    ]
+
Index: sasmodels/models/lib/wrc_cyl.c
===================================================================
--- sasmodels/models/lib/wrc_cyl.c	(revision f94d8a27ee285627fab8769c64a054dbc75e774f)
+++ sasmodels/models/lib/wrc_cyl.c	(revision 504abeee8bc5d67a14d7598227542961d3d748ed)
@@ -297,9 +297,25 @@
 }
 
-
+static double
+Sdebye_kernel(double arg)
+{
+    // ORIGINAL
+    double result = 2.0*(exp(-arg) + arg -1.0)/(pow((arg),2));
+
+    // CONVERSION 1 from http://herbie.uwplse.org/
+    //
+    // exhibits discontinuity - needs more investigation
+    //double a1 = 1.0/6.0;
+    //double a2 = 1.0/72.0;
+    //double a3 = 1.0/24.0;
+    //double result = pow((1.0 - a1*arg - (a2+a3)*arg*arg), 2);
+
+    return result;
+}
 static double
 Sdebye(double q, double L, double b)
 {
-    return 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2));
+    double arg = u_WR(q,L,b);
+    return Sdebye_kernel(arg);
 }
 
@@ -308,5 +324,7 @@
 Sdebye1(double q, double L, double b)
 {
-    return 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) );
+    double arg = u1(q,L,b);
+    return Sdebye_kernel(arg);
+
 }
 
