Index: CylinderModel.py
===================================================================
--- CylinderModel.py	(revision 32c160aab45a38358aa878cf4622c2cec2d87725)
+++ CylinderModel.py	(revision a1542aaea8c3c4b7e87ee6444cd821536e45ba02)
@@ -5,4 +5,4 @@
 """
 from sasmodels.sasview_model import make_class
-from sasmodels.models import cylinder_clone as cylinder
+from sasmodels.models import cylinder
 CylinderModel = make_class(cylinder, dtype='single')
Index: compare.py
===================================================================
--- compare.py	(revision 373d1b6113e17053b78f2b86735a454ab02666cf)
+++ compare.py	(revision 373d1b6113e17053b78f2b86735a454ab02666cf)
@@ -10,5 +10,7 @@
 
 from sasmodels.bumps_model import BumpsModel, plot_data, tic
-from sasmodels import kernelcl, kerneldll
+try: from sasmodels import kernelcl
+except: from sasmodels import kerneldll as kernelcl
+from sasmodels import kerneldll
 from sasmodels.convert import revert_model
 
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision 250fa25ad675ea9a85cddbc91dd32cd7a1255fb0)
+++ sasmodels/kernelcl.py	(revision 87fce00c768607daaeafbbd67ca19eb27459bc2e)
@@ -30,5 +30,5 @@
 try:
     import pyopencl as cl
-except ImportError,exc:
+except RuntimeError,exc:
     warnings.warn(str(exc))
     raise RuntimeError("OpenCL not available")
Index: sasmodels/kerneldll.py
===================================================================
--- sasmodels/kerneldll.py	(revision f786ff3daee63fedbada60cec1eb19964c370615)
+++ sasmodels/kerneldll.py	(revision 87fce00c768607daaeafbbd67ca19eb27459bc2e)
@@ -23,6 +23,6 @@
     # Can't find VCOMP90.DLL (don't know why), so remove openmp support from windows compiler build
     #COMPILE = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG /Tp%(source)s /link /DLL /INCREMENTAL:NO /MANIFEST /OUT:%(output)s"
-    COMPILE = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG /Tp%(source)s /openmp /link /DLL /INCREMENTAL:NO /MANIFEST /OUT:%(output)s"
-    #COMPILE = "gcc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %(source)s -o %(output)s -lm"
+    #COMPILE = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG /Tp%(source)s /openmp /link /DLL /INCREMENTAL:NO /MANIFEST /OUT:%(output)s"
+    COMPILE = "gcc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %(source)s -o %(output)s -lm"
 else:
     COMPILE = "cc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %(source)s -o %(output)s -lm"
Index: sasmodels/models/barbell.c
===================================================================
--- sasmodels/models/barbell.c	(revision 0706431576a8920e01fbe34cc6c63d2af1648d77)
+++ sasmodels/models/barbell.c	(revision 0706431576a8920e01fbe34cc6c63d2af1648d77)
@@ -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);
+
+	    const 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.0*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(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; //form_volume(bell_radius,radius,length);
+}
+
+
+
+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 + sn*qy;
+    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)/sn;
+
+    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 0706431576a8920e01fbe34cc6c63d2af1648d77)
+++ sasmodels/models/barbell.py	(revision 0706431576a8920e01fbe34cc6c63d2af1648d77)
@@ -0,0 +1,158 @@
+#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/orientation.jpg
+
+Figure. Definition of the angles for oriented 2D barbells.
+
+.. image:: img/orientation2.jpg
+
+*Figure. Examples of the angles for oriented pp against the detector plane.*
+
+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):
+    return 1.0
+
+# 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=5,
+    length_pd=.2,length_pd_n=5,
+    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')
Index: sasmodels/models/bcc.c
===================================================================
--- sasmodels/models/bcc.c	(revision c95dc9084eac92b9289ad0551b8d6b7bcdfa4240)
+++ sasmodels/models/bcc.c	(revision c95dc9084eac92b9289ad0551b8d6b7bcdfa4240)
@@ -0,0 +1,175 @@
+double form_volume(double radius);
+double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld);
+double Iqxy(double qx, double qy, double dnn,
+    double d_factor, double radius,double sld, double solvent_sld,
+    double theta, double phi, double psi);
+
+double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
+double _BCCeval(double Theta, double Phi, double temp1, double temp3);
+double _sphereform(double q, double radius, double sld, double solvent_sld);
+
+
+double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
+
+	const double Da = d_factor*dnn;
+	const double temp1 = q*q*Da*Da;
+	const double temp3 = q*dnn;
+
+	double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
+	return(retVal);
+}
+
+double _BCCeval(double Theta, double Phi, double temp1, double temp3) {
+
+	double temp6,temp7,temp8,temp9,temp10;
+	double result;
+
+	temp6 = sin(Theta);
+	temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta);
+	temp8 = -1.0*sin(Theta)*cos(Phi)-sin(Theta)*sin(Phi)+cos(Theta);
+	temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)-cos(Theta);
+	temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
+	result = pow(1.0-(temp10*temp10),3)*temp6/((1.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10)));
+
+	return (result);
+}
+
+double _sphereform(double q, double radius, double sld, double solvent_sld){
+    const double qr = q*radius;
+    double sn, cn;
+    SINCOS(qr, sn, cn);
+    const double bes = (qr == 0.0 ? 1.0 : 3.0*(sn-qr*cn)/(qr*qr*qr));
+    const double fq = bes * (sld - solvent_sld)*form_volume(radius);
+    return 1.0e-4*fq*fq;
+}
+
+double form_volume(double radius){
+    return 1.333333333333333*M_PI*radius*radius*radius;
+}
+
+
+double Iq(double q, double dnn,
+  double d_factor, double radius,
+  double sld, double solvent_sld){
+
+	//Volume fraction calculated from lattice symmetry and sphere radius
+	const double s1 = dnn/sqrt(0.75);
+	const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
+
+    const double va = 0.0;
+    const double vb = 2.0*M_PI;
+    const double vaj = 0.0;
+    const double vbj = M_PI;
+
+    double summ = 0.0;
+    double answer = 0.0;
+	for(int i=0; i<150; i++) {
+		//setup inner integral over the ellipsoidal cross-section
+		double summj=0.0;
+		const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;		//the outer dummy is phi
+		for(int j=0;j<150;j++) {
+			//20 gauss points for the inner integral
+			double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;		//the inner dummy is theta
+			double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi);
+			summj += yyy;
+		}
+		//now calculate the value of the inner integral
+		double answer = (vbj-vaj)/2.0*summj;
+
+		//now calculate outer integral
+		summ = summ+(Gauss150Wt[i] * answer);
+	}		//final scaling is done at the end of the function, after the NT_FP64 case
+
+	answer = (vb-va)/2.0*summ;
+	answer = answer*_sphereform(q,radius,sld,solvent_sld)*latticescale;
+
+    return answer;
+
+
+}
+
+
+double Iqxy(double qx, double qy, double dnn,
+    double d_factor, double radius,double sld, double solvent_sld,
+    double theta, double phi, double psi){
+
+  double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
+  double q_z;
+  double cos_val_b3, cos_val_b2, cos_val_b1;
+  double a1_dot_q, a2_dot_q,a3_dot_q;
+  double answer;
+  double Zq, Fkq, Fkq_2;
+
+  //convert to q and make scaled values
+  double q = sqrt(qx*qx+qy*qy);
+  double q_x = qx/q;
+  double q_y = qy/q;
+
+  //convert angle degree to radian
+  theta = theta * M_PI_180;
+  phi = phi * M_PI_180;
+  psi = psi * M_PI_180;
+
+  const double Da = d_factor*dnn;
+  const double qDa_2 = q*q*Da*Da;
+  const double s1 = dnn/sqrt(0.75);
+
+
+  //the occupied volume of the lattice
+  const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
+  // q vector
+  q_z = 0.0; // for SANS; assuming qz is negligible
+  /// Angles here are respect to detector coordinate
+  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
+    // b3 axis orientation
+    b3_x = cos(theta) * cos(phi);
+    b3_y = sin(theta);
+    //b3_z = -cos(theta) * sin(phi);
+    cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
+
+    //alpha = acos(cos_val_b3);
+    // b1 axis orientation
+    b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+    b1_y = sin(psi)*cos(theta);
+    cos_val_b1 = b1_x*q_x + b1_y*q_y;
+    // b2 axis orientation
+    b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+  	b2_y = cos(theta)*cos(psi);
+    cos_val_b2 = b2_x*q_x + b2_y*q_y;
+
+    // The following test should always pass
+    if (fabs(cos_val_b3)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b3 = 1.0;
+    }
+    if (fabs(cos_val_b2)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b2 = 1.0;
+    }
+    if (fabs(cos_val_b1)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b1 = 1.0;
+    }
+    // Compute the angle btw vector q and the a3 axis
+    a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
+
+    // a1 axis
+    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
+
+    // a2 axis
+    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
+
+
+    // Get Fkq and Fkq_2
+    Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
+    Fkq_2 = Fkq*Fkq;
+    // Call Zq=Z1*Z2*Z3
+    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
+
+  // Use SphereForm directly from libigor
+  answer = _sphereform(q,radius,sld,solvent_sld)*Zq*latticescale;
+
+  return answer;
+ }
Index: sasmodels/models/bcc.py
===================================================================
--- sasmodels/models/bcc.py	(revision c95dc9084eac92b9289ad0551b8d6b7bcdfa4240)
+++ sasmodels/models/bcc.py	(revision c95dc9084eac92b9289ad0551b8d6b7bcdfa4240)
@@ -0,0 +1,123 @@
+#bcc paracrystal model
+#note model title and parameter table are automatically inserted
+#note - calculation requires double precision
+"""
+Calculates the scattering from a **body-centered cubic lattice** with paracrystalline distortion. Thermal vibrations
+are considered to be negligible, and the size of the paracrystal is infinitely large. Paracrystalline distortion is
+assumed to be isotropic and characterized by a Gaussian distribution.
+
+The returned value is scaled to units of |cm^-1|\ |sr^-1|, absolute scale.
+
+Definition
+----------
+
+The scattering intensity *I(q)* is calculated as
+
+.. image:: img/image167.jpg
+
+where *scale* is the volume fraction of spheres, *Vp* is the volume of the primary particle, *V(lattice)* is a volume
+correction for the crystal structure, *P(q)* is the form factor of the sphere (normalized), and *Z(q)* is the
+paracrystalline structure factor for a body-centered cubic structure.
+
+Equation (1) of the 1990 reference is used to calculate *Z(q)*, using equations (29)-(31) from the 1987 paper for
+*Z1*\ , *Z2*\ , and *Z3*\ .
+
+The lattice correction (the occupied volume of the lattice) for a body-centered cubic structure of particles of radius
+*R* and nearest neighbor separation *D* is
+
+.. image:: img/image159.jpg
+
+The distortion factor (one standard deviation) of the paracrystal is included in the calculation of *Z(q)*
+
+.. image:: img/image160.jpg
+
+where *g* is a fractional distortion based on the nearest neighbor distance.
+
+The body-centered cubic lattice is
+
+.. image:: img/image168.jpg
+
+For a crystal, diffraction peaks appear at reduced q-values given by
+
+.. image:: img/image162.jpg
+
+where for a body-centered cubic lattice, only reflections where (\ *h* + *k* + *l*\ ) = even are allowed and
+reflections where (\ *h* + *k* + *l*\ ) = odd are forbidden. Thus the peak positions correspond to (just the first 5)
+
+.. image:: img/image169.jpg
+
+**NB: The calculation of** *Z(q)* **is a double numerical integral that must be carried out with a high density of**
+**points to properly capture the sharp peaks of the paracrystalline scattering.** So be warned that the calculation is
+SLOW. Go get some coffee. Fitting of any experimental data must be resolution smeared for any meaningful fit. This
+makes a triple integral. Very, very slow. Go get lunch!
+
+This example dataset is produced using 200 data points, *qmin* = 0.001 |Ang^-1|, *qmax* = 0.1 |Ang^-1| and the above
+default values.
+
+.. image:: img/image170.jpg
+
+*Figure. 1D plot in the linear scale using the default values (w/200 data point).*
+
+The 2D (Anisotropic model) is based on the reference below where *I(q)* is approximated for 1d scattering. Thus the
+scattering pattern for 2D may not be accurate. Note that we are not responsible for any incorrectness of the 2D model
+computation.
+
+.. image:: img/image165.gif
+
+.. image:: img/image171.jpg
+
+*Figure. 2D plot using the default values (w/200X200 pixels).*
+
+REFERENCE
+
+Hideki Matsuoka et. al. *Physical Review B*, 36 (1987) 1754-1765
+(Original Paper)
+
+Hideki Matsuoka et. al. *Physical Review B*, 41 (1990) 3854 -3856
+(Corrections to FCC and BCC lattice structure calculation)
+"""
+
+from numpy import pi, inf
+
+name = "BCCparacrystal"
+title = "Body-centred cubic lattic with paracrystalline distortion"
+description = """
+    Calculates the scattering from a **body-centered cubic lattice** with paracrystalline distortion. Thermal vibrations
+    are considered to be negligible, and the size of the paracrystal is infinitely large. Paracrystalline distortion is
+    assumed to be isotropic and characterized by a Gaussian distribution.
+    """
+
+parameters = [
+#   [ "name", "units", default, [lower, upper], "type","description" ],
+    [ "dnn", "Ang", 220, [-inf,inf],"","Nearest neighbour distance"],
+    [ "d_factor", "", 0.06,[-inf,inf],"","Paracrystal distortion factor" ],
+    [ "radius", "Ang",  40, [0, inf], "volume","Particle radius" ],
+    [ "sld", "1e-6/Ang^2", 4, [-inf,inf], "", "Particle scattering length density" ],
+    [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "","Solvent scattering length density" ],
+    [ "theta", "degrees", 60, [-inf, inf], "orientation","In plane angle" ],
+    [ "phi", "degrees", 60, [-inf, inf], "orientation","Out of plane angle" ],
+    [ "psi", "degrees", 60, [-inf,inf], "orientation","Out of plane angle"]
+    ]
+
+source = [ "lib/J1.c", "lib/gauss150.c", "bcc.c" ]
+
+def ER(radius, length):
+    return 0
+
+# parameters for demo
+demo = dict(
+    scale=1, background=0,
+    dnn=220, d_factor=0.06, sld=4, solvent_sld=1,
+    radius=40,
+    theta=60, phi=60, psi=60,
+    radius_pd=.2, radius_pd_n=0.2,
+    theta_pd=15, theta_pd_n=0,
+    phi_pd=15, phi_pd_n=0,
+    psi_pd=15, psi_pd_n=0,
+    )
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname='BCCrystalModel'
+oldpars=dict(sld='sldSph',
+             solvent_sld='sldSolv')
Index: smodels/models/cylinder_clone.c
===================================================================
--- sasmodels/models/cylinder_clone.c	(revision 994d77f9b713691e41106d1f9af9661dfb2dc419)
+++ 	(revision )
@@ -1,71 +1,0 @@
-double form_volume(double radius, double length);
-double Iq(double q, double sld, double solvent_sld, double radius, double length);
-double Iqxy(double qx, double qy, double sld, double solvent_sld, double radius, double length, double theta, double phi);
-
-
-// twovd = 2 * volume * delta_rho
-// besarg = q * R * sin(alpha)
-// siarg = q * L/2 * cos(alpha)
-double _cyl(double twovd, double besarg, double siarg, double alpha);
-double _cyl(double twovd, double besarg, double siarg, double alpha)
-{
-    const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg);
-    const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg);
-    return twovd*si*bj;
-}
-
-double form_volume(double radius, double length)
-{
-    return M_PI*radius*radius*length;
-}
-double Iq(double q,
-    double sldCyl,
-    double sldSolv,
-    double radius,
-    double length)
-{
-    const double qr = q*radius;
-    const double qh = q*0.5*length;
-    const double twovd = 2.0*(sldCyl-sldSolv)*form_volume(radius, length);
-    double total = 0.0;
-    // double lower=0, upper=M_PI_2;
-    for (int i=0; i<76 ;i++) {
-        // translate a point in [-1,1] to a point in [lower,upper]
-        //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
-        const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2);
-        double sn, cn;
-        SINCOS(alpha, sn, cn);
-        const double fq = _cyl(twovd, qr*sn, qh*cn, alpha);
-        total += Gauss76Wt[i] * fq * fq * sn;
-    }
-    // translate dx in [-1,1] to dx in [lower,upper]
-    //const double form = (upper-lower)/2.0*total;
-    return 1.0e8 * total * M_PI_4;
-}
-
-double Iqxy(double qx, double qy,
-    double sldCyl,
-    double sldSolv,
-    double radius,
-    double length,
-    double cyl_theta,
-    double cyl_phi)
-{
-    double sn, cn; // slots to hold sincos function output
-
-    // Compute angle alpha between q and the cylinder axis
-    SINCOS(cyl_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 spherical_integration = fabs(cn)*M_PI_2;
-    const double q = sqrt(qx*qx+qy*qy);
-    const double cos_val = cn*cos(cyl_phi*M_PI_180)*(qx/q) + sn*(qy/q);
-    const double alpha = acos(cos_val);
-
-    const double qr = q*radius;
-    const double qh = q*0.5*length;
-    const double twovd = 2.0*(sldCyl-sldSolv)*form_volume(radius, length);
-    SINCOS(alpha, sn, cn);
-    const double fq = _cyl(twovd, qr*sn, qh*cn, alpha);
-    return 1.0e8 * fq * fq * spherical_integration;
-}
Index: smodels/models/cylinder_clone.py
===================================================================
--- sasmodels/models/cylinder_clone.py	(revision 19dcb93366741a275baa8ad855ce8a310dca78f8)
+++ 	(revision )
@@ -1,43 +1,0 @@
-
-from numpy import pi, inf
-
-name = "CylinderModel"
-title = "Cylinder with uniform scattering length density"
-description = """\
-         f(q)= 2*(sldCyl - sldSolv)*V*sin(qLcos(alpha/2))
-                /[qLcos(alpha/2)]*J1(qRsin(alpha/2))/[qRsin(alpha)]
-
-                P(q,alpha)= scale/V*f(q)^(2)+background
-                V: Volume of the cylinder
-                R: Radius of the cylinder
-                L: Length of the cylinder
-                J1: The bessel function
-                alpha: angle betweenthe axis of the
-                cylinder and the q-vector for 1D
-                :the ouput is P(q)=scale/V*integral
-                from pi/2 to zero of...
-                f(q)^(2)*sin(alpha)*dalpha+ bkg
-"""
-
-parameters = [
-    #   [ "name", "units", default, [lower, upper], "type",
-    #     "description" ],
-    [ "sldCyl", "1/Ang^2", 4e-6, [-inf,inf], "",
-      "Cylinder scattering length density" ],
-    [ "sldSolv", "1/Ang^2", 1e-6, [-inf,inf], "",
-      "Solvent scattering length density" ],
-    [ "radius", "Ang",  20, [0, inf], "volume",
-      "Cylinder radius" ],
-    [ "length", "Ang",  400, [0, inf], "volume",
-      "Cylinder length" ],
-    [ "cyl_theta", "degrees", 60, [-inf, inf], "orientation",
-      "In plane angle" ],
-    [ "cyl_phi", "degrees", 60, [-inf, inf], "orientation",
-      "Out of plane angle" ],
-    ]
-source = [ "lib/J1.c", "lib/gauss76.c", "cylinder_clone.c" ]
-
-def ER(radius, length):
-    ddd = 0.75*radius*(2*radius*length + (length+radius)*(length+pi*radius))
-    return 0.5 * (ddd)**(1./3.)
-
