Index: sasmodels/models/barbell.c
===================================================================
--- sasmodels/models/barbell.c	(revision 58f41fe82f2ac143c9a5830fff81950cc382d40c)
+++ sasmodels/models/barbell.c	(revision 58f41fe82f2ac143c9a5830fff81950cc382d40c)
@@ -0,0 +1,124 @@
+double form_volume(double bell_radius, double radius, double length);
+double Iq(double q, double sld, double solvent_sld, double bell_radius, double radius, double length);
+double Iqxy(double qx, double qy, double sld, double solvent_sld,
+    double bell_radius, double radius, double length, double theta, double phi);
+
+//barbell kernel - same as dumbell
+double _bell_kernel(double q, double h, double bell_radius,
+    double length, double sin_alpha, double cos_alpha);
+double _bell_kernel(double q, double h, double bell_radius,
+    double length, double sin_alpha, double cos_alpha)
+{
+    const double upper = 1.0;
+    const double lower = -1.0*h/bell_radius;
+
+    double total = 0.0;
+    for (int i = 0; i < 76; i++){
+        const double t = 0.5*(Gauss76Z[i]*(upper-lower)+upper+lower);
+	    const double arg1 = q*cos_alpha*(bell_radius*t+h+length*0.5);
+	    const double arg2 = q*bell_radius*sin_alpha*sqrt(1.0-t*t);
+
+        const double be = (arg2 == 0.0 ? 0.5 :J1(arg2)/arg2);
+
+	    double Fq = cos(arg1)*(1.0-t*t)*be;
+
+	    total += Gauss76Wt[i] * Fq;
+    }
+    const double integral = 0.5*(upper-lower)*total;
+    return 4.0*M_PI*bell_radius*bell_radius*bell_radius*integral;
+}
+
+double form_volume(double bell_radius,
+        double radius,
+        double length)
+{
+
+    // bell radius should never be less than radius when this is called
+    const double hdist = sqrt(bell_radius*bell_radius - radius*radius);
+    const double p1 = 2*bell_radius*bell_radius*bell_radius/3.0;
+    const double p2 = bell_radius*bell_radius*hdist;
+    const double p3 = hdist*hdist*hdist/3.0;
+
+    return M_PI*radius*radius*length + 2.0*M_PI*(p1+p2+p3);
+}
+
+double Iq(double q, double sld,
+    double solvent_sld,
+    double bell_radius,
+    double radius,
+    double length)
+{
+    double sn, cn; // slots to hold sincos function output
+
+    if (bell_radius < radius) return -1.0;
+
+    const double lower = 0.0;
+    const double upper = M_PI_2;
+    const double h = sqrt(fabs(bell_radius*bell_radius-radius*radius));
+    double total = 0.0;
+    for (int i = 0; i < 76; i++){
+        const double alpha= 0.5*(Gauss76Z[i]*(upper-lower) + upper + lower);
+        SINCOS(alpha, sn, cn);
+
+        const double bell_Fq = _bell_kernel(q, h, bell_radius, length, sn, cn);
+
+        const double arg1 = q*length*0.5*cn;
+        const double arg2 = q*radius*sn;
+        // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1
+        const double be = (arg2 == 0.0 ? 0.5 :J1(arg2)/arg2);
+        const double si = (arg1 == 0.0 ? 1.0 :sin(arg1)/arg1);
+        const double cyl_Fq = M_PI*radius*radius*length*si*2.0*be;
+
+        const double Aq = cyl_Fq + bell_Fq;
+        total += Gauss76Wt[i] * Aq * Aq * sn;
+    }
+
+    const double form = total*(upper-lower)*0.5;
+
+    //Contrast and volume normalization
+    const double s = (sld - solvent_sld);
+    return form*1.0e-4*s*s;
+}
+
+
+
+double Iqxy(double qx, double qy,
+    double sld,
+    double solvent_sld,
+    double bell_radius,
+    double radius,
+    double length,
+    double theta,
+    double phi)
+{
+     double sn, cn; // slots to hold sincos function output
+
+    // Exclude invalid inputs.
+    if (bell_radius < radius) return -1.0;
+
+    // Compute angle alpha between q and the cylinder axis
+    SINCOS(theta*M_PI_180, sn, cn);
+    // # The following correction factor exists in sasview, but it can't be
+    // # right, so we are leaving it out for now.
+    const double q = sqrt(qx*qx+qy*qy);
+    const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q);
+    const double alpha = acos(cos_val); // rod angle relative to q
+    SINCOS(alpha, sn, cn);
+
+    const double h = sqrt(bell_radius*bell_radius - radius*radius); // negative h
+    const double bell_Fq = _bell_kernel(q, h, bell_radius, length, sn, cn);
+
+    const double besarg = q*radius*sn;
+    const double siarg = q*0.5*length*cn;
+    // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1
+    const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg);
+    const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg);
+    const double cyl_Fq = M_PI*radius*radius*length*2.0*bj*si;
+
+    // Volume weighted average F(q)
+    const double Aq = cyl_Fq + bell_Fq;
+
+    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1
+    const double s = (sld - solvent_sld);
+    return 1.0e-4 * Aq * Aq * s * s; // form_volume(radius, cap_radius, length);
+}
Index: sasmodels/models/barbell.py
===================================================================
--- sasmodels/models/barbell.py	(revision 58f41fe82f2ac143c9a5830fff81950cc382d40c)
+++ sasmodels/models/barbell.py	(revision 58f41fe82f2ac143c9a5830fff81950cc382d40c)
@@ -0,0 +1,159 @@
+#barbell model
+# cylinder model
+# Note: model title and parameter table are inserted automatically
+r"""
+
+Calculates the scattering from a barbell-shaped cylinder (This model simply becomes the DumBellModel when the length of
+the cylinder, *L*, is set to zero). That is, a sphereocylinder with spherical end caps that have a radius larger than
+that of the cylinder and the center of the end cap radius lies outside of the cylinder. All dimensions of the BarBell
+are considered to be monodisperse. See the diagram for the details of the geometry and restrictions on parameter values.
+
+Definition
+----------
+
+The returned value is scaled to units of |cm^-1|\ |sr^-1|, absolute scale.
+
+The barbell geometry is defined as
+
+.. image:: img/image105.jpg
+
+where *r* is the radius of the cylinder. All other parameters are as defined in the diagram.
+
+Since the end cap radius
+*R* >= *r* and by definition for this geometry *h* < 0, *h* is then defined by *r* and *R* as
+
+*h* = -1 \* sqrt(*R*\ :sup:`2` - *r*\ :sup:`2`)
+
+The scattered intensity *I(q)* is calculated as
+
+.. math::
+
+    I(Q) = \frac{(\Delta \rho)^2}{V} \left< A^2(Q)\right>
+
+where the amplitude *A(q)* is given as
+
+.. math::
+
+    A(Q) =&\ \pi r^2L
+        {\sin\left(\tfrac12 QL\cos\theta\right)
+            \over \tfrac12 QL\cos\theta}
+        {2 J_1(Qr\sin\theta) \over Qr\sin\theta} \\
+        &\ + 4 \pi R^3 \int_{-h/R}^1 dt
+        \cos\left[ Q\cos\theta
+            \left(Rt + h + {\tfrac12} L\right)\right]
+        \times (1-t^2)
+        {J_1\left[QR\sin\theta \left(1-t^2\right)^{1/2}\right]
+             \over QR\sin\theta \left(1-t^2\right)^{1/2}}
+
+The < > brackets denote an average of the structure over all orientations. <*A* :sup:`2`\ *(q)*> is then the form
+factor, *P(q)*. The scale factor is equivalent to the volume fraction of cylinders, each of volume, *V*. Contrast is
+the difference of scattering length densities of the cylinder and the surrounding solvent.
+
+The volume of the barbell is
+
+.. math::
+
+    V = \pi r_c^2 L + 2\pi\left(\tfrac23R^3 + R^2h-\tfrac13h^3\right)
+
+
+and its radius-of-gyration is
+
+.. math::
+
+    R_g^2 =&\ \left[ \tfrac{12}{5}R^5
+        + R^4\left(6h+\tfrac32 L\right)
+        + R^2\left(4h^2 + L^2 + 4Lh\right)
+        + R^2\left(3Lh^2 + \tfrac32 L^2h\right) \right. \\
+        &\ \left. + \tfrac25 h^5 - \tfrac12 Lh^4 - \tfrac12 L^2h^3
+        + \tfrac14 L^3r^2 + \tfrac32 Lr^4 \right]
+        \left( 4R^3 6R^2h - 2h^3 + 3r^2L \right)^{-1}
+
+**The requirement that** *R* >= *r* **is not enforced in the model!** It is up to you to restrict this during analysis.
+
+This example dataset is produced by running the Macro PlotBarbell(), using 200 data points, *qmin* = 0.001 |Ang^-1|,
+*qmax* = 0.7 |Ang^-1| and the following default values
+
+==============  ========  =============
+Parameter name  Units     Default value
+==============  ========  =============
+scale           None      1.0
+len_bar         |Ang|     400.0
+rad_bar         |Ang|     20.0
+rad_bell        |Ang|     40.0
+sld_barbell     |Ang^-2|  1.0e-006
+sld_solv        |Ang^-2|  6.3e-006
+background      |cm^-1|   0
+==============  ========  =============
+
+.. image:: img/image110.jpg
+
+*Figure. 1D plot using the default values (w/256 data point).*
+
+For 2D data: The 2D scattering intensity is calculated similar to the 2D cylinder model. For example, for
+|theta| = 45 deg and |phi| = 0 deg with default values for other parameters
+
+.. image:: img/image111.jpg
+
+*Figure. 2D plot (w/(256X265) data points).*
+
+.. image:: img/image061.jpg
+
+*Figure. Examples of the angles for oriented pp against the detector plane.*
+
+.. image:: img/image062.jpg
+
+Figure. Definition of the angles for oriented 2D barbells.
+
+REFERENCE
+
+H Kaya, *J. Appl. Cryst.*, 37 (2004) 37 223-230
+
+H Kaya and N R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda and errata)
+
+"""
+from numpy import pi, inf
+
+name = "barbell"
+title = "Cylinder with spherical end caps"
+description = """
+        Calculates the scattering from a barbell-shaped cylinder. That is a sphereocylinder with spherical end caps that have a radius larger than that of the cylinder and the center of the end cap
+		radius lies outside of the cylinder.
+		Note: As the length of cylinder(bar) -->0,it becomes a dumbbell. And when rad_bar = rad_bell, it is a spherocylinder.
+		It must be that rad_bar <(=) rad_bell.
+"""
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type","description" ],
+    [ "sld", "1e-6/Ang^2", 4, [-inf,inf], "", "Barbell scattering length density" ],
+    [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "","Solvent scattering length density" ],
+    [ "bell_radius", "Ang",  40, [0, inf], "volume","Spherical bell radius" ],
+    [ "radius", "Ang",  20, [0, inf], "volume","Cylindrical bar radius" ],
+    [ "length", "Ang",  400, [0, inf], "volume","Cylinder bar length" ],
+    [ "theta", "degrees", 60, [-inf, inf], "orientation","In plane angle" ],
+    [ "phi", "degrees", 60, [-inf, inf], "orientation","Out of plane angle" ],
+    ]
+
+source = [ "lib/J1.c", "lib/gauss76.c", "barbell.c" ]
+
+def ER(radius, length):
+    ddd = 0.75*radius*(2*radius*length + (length+radius)*(length+pi*radius))
+    return 0.5 * (ddd)**(1./3.)
+
+# parameters for demo
+demo = dict(
+    scale=1, background=0,
+    sld=6, solvent_sld=1,
+    bell_radius = 40, radius=20, length=400,
+    theta=60, phi=60,
+    radius_pd=.2, radius_pd_n=0,
+    length_pd=.2,length_pd_n=0,
+    theta_pd=15, theta_pd_n=0,
+    phi_pd=15, phi_pd_n=0,
+    )
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname='BarBellModel'
+oldpars=dict(sld='sld_barbell',
+             solvent_sld='sld_solv', bell_radius='rad_bell',
+             radius='rad_bar',length='len_bar')
