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 a8b3cdb60df4f2e22c9b3c5c3110bb51b7dba2c6)
+++ sasmodels/models/elliptical_cylinder.py	(revision a8b3cdb60df4f2e22c9b3c5c3110bb51b7dba2c6)
@@ -0,0 +1,160 @@
+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.
+
+.. image:: img/elliptical_cylinder_geometry.gif
+
+    *Figure.* *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 *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
+
+The returned value is scaled to units of |cm^-1|.
+
+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.
+
+.. image:: img/elliptical_cylinder_geometry_2d.jpg
+
+    *Figure. Definition of angles for 2D*
+
+.. image:: img/core_shell_bicelle_fig2.jpg
+
+    *Figure. 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.
+
+
+.. image:: img/elliptical_cylinder_comparison_1d.jpg
+
+    *Figure. 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).
+
+.. image:: img/elliptical_cylinder_validation_1d.gif
+
+    *Figure. 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.
+
+.. image:: img/elliptical_cylinder_averaging.gif
+
+    *Figure. 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 = """
+    F^2(q) = 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)
+    
+In 1D we have
+
+.. math::
+    P(q) = scale  <F^2> / V
+"""
+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
+    """
+    r = math.sqrt(r_minor * r_minor * r_ratio )
+    ddd = 0.75 * r * (2 * r * length + (length + r) * (length + pi * r))
+    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/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);
+}
