Index: sasmodels/models/flexible_cylinder.c
===================================================================
--- sasmodels/models/flexible_cylinder.c	(revision f94d8a27ee285627fab8769c64a054dbc75e774f)
+++ sasmodels/models/flexible_cylinder.c	(revision f94d8a27ee285627fab8769c64a054dbc75e774f)
@@ -0,0 +1,64 @@
+double form_volume(double length, double kuhn_length, double radius);
+double Iq(double q, double length, double kuhn_length, double radius,
+          double sld, double solvent_sld);
+double Iqxy(double qx, double qy, double length, double kuhn_length,
+            double radius, double sld, double solvent_sld);
+double flexible_cylinder_kernel(double q, double length, double kuhn_length,
+                                double radius, double sld, double solvent_sld);
+
+
+double form_volume(double length, double kuhn_length, double radius)
+{
+
+      return 0.0;
+}
+
+double flexible_cylinder_kernel(double q,
+          double length,
+          double kuhn_length,
+          double radius,
+          double sld,
+          double solvent_sld)
+{
+
+	double Pi = 4.0*atan(1.0);
+
+	double cont = sld-solvent_sld;
+	double qr = q*radius;
+	double flex = Sk_WR(q,length,kuhn_length);
+    double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr);
+
+	flex *= crossSect;
+	flex *= Pi*radius*radius*length;
+	flex *= cont*cont;
+	flex *= 1.0e-4;
+
+	return flex;
+}
+
+double Iq(double q,
+          double length,
+          double kuhn_length,
+          double radius,
+          double sld,
+          double solvent_sld)
+{
+
+	double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);
+
+	return result;
+}
+
+double Iqxy(double qx, double qy,
+            double length,
+            double kuhn_length,
+            double radius,
+            double sld,
+            double solvent_sld)
+{
+	double q;
+	q = sqrt(qx*qx+qy*qy);
+	double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);
+
+	return result;
+}
Index: sasmodels/models/flexible_cylinder.py
===================================================================
--- sasmodels/models/flexible_cylinder.py	(revision f94d8a27ee285627fab8769c64a054dbc75e774f)
+++ sasmodels/models/flexible_cylinder.py	(revision f94d8a27ee285627fab8769c64a054dbc75e774f)
@@ -0,0 +1,129 @@
+r"""
+This model provides the form factor, $P(q)$, for a flexible cylinder where the form factor
+is normalized by the volume of the cylinder.
+**Inter-cylinder interactions are NOT provided for.**
+
+.. math::
+
+    P(q) = \text{scale} \left<F^2\right>/V + \text{background}
+
+where the averaging $\left<\ldots\right>$ is applied only for the 1D calculation
+
+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
+-----------
+
+.. figure:: img/flexible_cylinder_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 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.
+
+
+.. figure:: img/flexible_cylinder_1d.jpg
+
+    1D plot using the default values (w/1000 data point).
+
+
+Our model uses the form factor calculations implemented in a c-library provided by the NIST Center
+for Neutron Research (Kline, 2006).
+
+
+From the reference:
+
+    'Method 3 With Excluded Volume' is used.
+    The model is a parametrization of simulations of a discrete representation
+    of the worm-like chain model of Kratky and Porod applied in the pseudocontinuous limit.
+    See equations (13,26-27) in the original reference for the details.
+
+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"
+title = "Flexible cylinder where the form factor is normalized by the volume of the cylinder."
+description = """Note : scale and contrast=sld-solvent_sld 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"
+
+#             ["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"],
+              ["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"],
+             ]
+
+source = ["lib/J1.c", "lib/wrc_cyl.c", "flexible_cylinder.c"]
+
+demo = dict(scale=1.0, background=0.0001,
+            length=1000.0,
+            kuhn_length=100.0,
+            radius=20.0,
+            sld=1.0,
+            solvent_sld=6.3)
+
+oldname = 'FlexibleCylinderModel'
+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,
+         #  '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,
+           'sld':           1.0,
+           'solvent_sld':   6.3,
+           'background':    0.0001,
+           }, 1.0, 0.000595345],
+        [{'length':        10.0,
+           'kuhn_length': 800.0,
+           'radius':        2.0,
+           'sld':           6.0,
+           'solvent_sld':  12.3,
+           'background':    0.001,
+           }, 0.1, 1.55228],
+       [{'length':        100.0,
+           'kuhn_length': 800.0,
+           'radius':       50.0,
+           'sld':           0.1,
+           'solvent_sld':   5.1,
+           'background':    0.0,
+           }, 1.0, 0.000938456]
+         ]
+
Index: sasmodels/models/lib/wrc_cyl.c
===================================================================
--- sasmodels/models/lib/wrc_cyl.c	(revision f94d8a27ee285627fab8769c64a054dbc75e774f)
+++ sasmodels/models/lib/wrc_cyl.c	(revision f94d8a27ee285627fab8769c64a054dbc75e774f)
@@ -0,0 +1,415 @@
+/*
+    Functions for WRC implementation of flexible cylinders
+*/
+static double
+AlphaSquare(double x)
+{
+    return pow( (1.0 + (x/3.12)*(x/3.12) +
+         (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );
+}
+
+//
+static double
+Rgsquarezero(double q, double L, double b)
+{
+    return (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) -
+         0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b))));
+}
+
+//
+static double
+Rgsquareshort(double q, double L, double b)
+{
+    return AlphaSquare(L/b) * Rgsquarezero(q,L,b);
+}
+
+//
+static double
+Rgsquare(double q, double L, double b)
+{
+    return AlphaSquare(L/b)*L*b/6.0;
+}
+
+static double
+sech_WR(double x)
+{
+    return(1/cosh(x));
+}
+
+static double
+a1long(double q, double L, double b, double p1, double p2, double q0)
+{
+    double yy,C,C1,C2,C3,C4,C5,miu,Rg2;
+    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15;
+    double E,pi;
+    double b2,b3,b4,q02,q03,q04,q05,Rg22;
+
+    pi = 4.0*atan(1.0);
+    E = 2.718281828459045091;
+
+    if( L/b > 10.0) {
+        C = 3.06/pow((L/b),0.44);
+    } else {
+        C = 1.0;
+    }
+
+    C1 = 1.22;
+    C2 = 0.4288;
+    C3 = -1.651;
+    C4 = 1.523;
+    C5 = 0.1477;
+    miu = 0.585;
+
+    Rg2 = Rgsquare(q,L,b);
+    Rg22 = Rg2*Rg2;
+    b2 = b*b;
+    b3 = b*b*b;
+    b4 = b3*b;
+    q02 = q0*q0;
+    q03 = q0*q0*q0;
+    q04 = q03*q0;
+    q05 = q04*q0;
+
+    t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 +
+        (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))));
+
+    t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
+        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) -
+        tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
+
+    t3 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
+        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
+        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));
+
+    t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
+
+    t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) -
+        b*p2*pow(q0,((-1.0) - p1 - p2))));
+
+    t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) +
+        (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) +
+        (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 +
+        (7.0*b2)/(15.0*q02*Rg2)))*Rg2)/b)));
+
+    t7 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
+        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
+        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) +
+        (sqrt(Rg2)*q0)/b)/C5),2));
+
+    t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
+        (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));
+
+    t9 = (2.0*b4*(((2.0*q0*Rg2)/b -
+        (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) -
+        tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
+
+    t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
+        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) +
+        (sqrt(Rg2)*q0)/b)/C5))))));
+
+    t11 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) -
+        3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) -
+        2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) -
+        1.0/miu)))/miu));
+
+    t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
+
+    t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 +
+        (7.0*b2)/(15.0*q02* Rg2))) +
+        (7.0*b2)/(15.0*q02*Rg2))));
+
+    t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
+        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) +
+        (sqrt(Rg2)*q0)/b)/C5))))));
+
+    t15 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
+        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
+        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));
+
+
+    yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) +
+        1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*
+        (((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -
+        t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) -
+        b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t13/L +
+        t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) +
+        (sqrt(Rg2)*q0)/b)/C5)))))))))));
+
+    return (yy);
+}
+
+static double
+a2long(double q, double L, double b, double p1, double p2, double q0)
+{
+    double yy,C1,C2,C3,C4,C5,miu,C,Rg2;
+    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi;
+    double E,b2,b3,b4,q02,q03,q04,q05,Rg22;
+
+    pi = 4.0*atan(1.0);
+    E = 2.718281828459045091;
+    if( L/b > 10.0) {
+        C = 3.06/pow((L/b),0.44);
+    } else {
+        C = 1.0;
+    }
+
+    C1 = 1.22;
+    C2 = 0.4288;
+    C3 = -1.651;
+    C4 = 1.523;
+    C5 = 0.1477;
+    miu = 0.585;
+
+    Rg2 = Rgsquare(q,L,b);
+    Rg22 = Rg2*Rg2;
+    b2 = b*b;
+    b3 = b*b*b;
+    b4 = b3*b;
+    q02 = q0*q0;
+    q03 = q0*q0*q0;
+    q04 = q03*q0;
+    q05 = q04*q0;
+
+    t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) -
+        b*p2*pow(q0,((-1.0) - p1 - p2)) ));
+
+    t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) +
+        (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) +
+        (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 +
+        (7*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L;
+
+    t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
+        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
+        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*
+        pow(sech_WR(((-C4) +(sqrt(Rg2)*q0)/b)/C5),2.0))/(2.0*C5);
+
+    t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
+        (q02*Rg2)/b2))*pow(sech_WR(((-C4) +
+        (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22);
+
+    t5 = (2.0*b4*(((2.0*q0*Rg2)/b -
+        (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*
+        ((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) +
+        (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
+
+    t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
+        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) +
+        (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22);
+
+    t7 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) -
+        3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),
+        ((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),
+        ((-1.0) - 1.0/miu)))/miu));
+
+    t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
+
+    t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 +
+        (7.0    *b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L;
+
+    t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) +
+        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) +
+        (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
+
+    yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) +
+        t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*
+        (((-((b*pi)/(L*q0))) + t9 + t10 +
+        1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
+        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
+        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*
+        ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))));
+
+    return (yy);
+}
+
+//
+static double
+a1short(double q, double L, double b, double p1short, double p2short, double q0)
+{
+    double yy,Rg2_sh;
+    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3;
+    double pi;
+
+    E = 2.718281828459045091;
+    pi = 4.0*atan(1.0);
+    Rg2_sh = Rgsquareshort(q,L,b);
+    Rg2_sh2 = Rg2_sh*Rg2_sh;
+    b3 = b*b*b;
+    t1 = ((q0*q0*Rg2_sh)/(b*b));
+    Et1 = pow(E,t1);
+    Emt1 =pow(E,-t1);
+    q02 = q0*q0;
+    q0p = pow(q0,(-4.0 + p1short) );
+
+    yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*
+        ((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short +
+        2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh -
+        2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 +
+        Et1*p2short*pi*q02*q0*Rg2_sh2))))));
+
+    return(yy);
+}
+
+static double
+a2short(double q, double L, double b, double p1short, double p2short, double q0)
+{
+    double yy,Rg2_sh;
+    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p;
+    double pi;
+
+    E = 2.718281828459045091;
+    pi = 4.0*atan(1.0);
+    Rg2_sh = Rgsquareshort(q,L,b);
+    Rg2_sh2 = Rg2_sh*Rg2_sh;
+    t1 = ((q0*q0*Rg2_sh)/(b*b));
+    Et1 = pow(E,t1);
+    Emt1 =pow(E,-t1);
+    q02 = q0*q0;
+    q0p = pow(q0,(-4.0 + p2short) );
+
+    //E is the number e
+    yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*
+        ((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short +
+        2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh -
+        2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 +
+        Et1*p1short*pi*q02*q0*Rg2_sh2)))))));
+
+    return (yy);
+}
+
+//WR named this w (too generic)
+static double
+w_WR(double x)
+{
+    return 0.5*(1 + tanh((x - 1.523)/0.1477));
+}
+
+//
+static double
+u1(double q, double L, double b)
+{
+    return Rgsquareshort(q,L,b)*q*q;
+}
+
+static double
+u_WR(double q, double L, double b)
+{
+    return Rgsquare(q,L,b)*q*q;
+}
+
+
+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));
+}
+
+//
+static double
+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) );
+}
+
+//
+static double
+Sexv(double q, double L, double b)
+{
+    double yy,C1,C2,C3,miu,Rg2;
+    C1=1.22;
+    C2=0.4288;
+    C3=-1.651;
+    miu = 0.585;
+
+    Rg2 = Rgsquare(q,L,b);
+
+    yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) +
+        w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) +
+        C2*pow((q*sqrt(Rg2)),(-2.0/miu)) +
+        C3*pow((q*sqrt(Rg2)),(-3.0/miu)));
+
+    return (yy);
+}
+
+
+static double
+Sexvnew(double q, double L, double b)
+{
+    double yy,C1,C2,C3,miu;
+    double del=1.05,C_star2,Rg2;
+
+    C1=1.22;
+    C2=0.4288;
+    C3=-1.651;
+    miu = 0.585;
+
+    //calculating the derivative to decide on the corection (cutoff) term?
+    // I have modified this from WRs original code
+
+    if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) {
+        C_star2 = 0.0;
+    } else {
+        C_star2 = 1.0;
+    }
+
+    Rg2 = Rgsquare(q,L,b);
+
+    yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) +
+        C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) +
+        C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu)));
+
+    return (yy);
+}
+
+double Sk_WR(double q, double L, double b)
+{
+  //
+  double p1,p2,p1short,p2short,q0;
+  double C,ans,q0short,Sexvmodify,pi;
+
+  pi = 4.0*atan(1.0);
+
+  p1 = 4.12;
+  p2 = 4.42;
+  p1short = 5.36;
+  p2short = 5.62;
+  q0 = 3.1;
+
+  q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);
+
+
+  if(L/b > 10.0) {
+    C = 3.06/pow((L/b),0.44);
+  } else {
+    C = 1.0;
+  }
+  //
+
+  //
+  if( L > 4*b ) { // Longer Chains
+    if (q*b <= 3.1) {   //Modified by Yun on Oct. 15,
+
+      Sexvmodify = Sexvnew(q, L, b);
+
+      ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) -
+            (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L);
+
+    } else { //q(i)*b > 3.1
+      ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) +
+            a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L);
+    }
+  } else { //L <= 4*b Shorter Chains
+    if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
+      if (q*b<=0.01) {
+        ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0;
+      } else {
+        ans = Sdebye1(q,L,b);
+      }
+    } else {  //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)
+      ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) +
+            a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) +
+            pi/(q*L);
+    }
+  }
+
+  return(ans);
+}
