Index: sasmodels/models/binary_hard_sphere.c
===================================================================
--- sasmodels/models/binary_hard_sphere.c	(revision 115d0f015d008c59c827203a96ee8894b7ec49b4)
+++ sasmodels/models/binary_hard_sphere.c	(revision 115d0f015d008c59c827203a96ee8894b7ec49b4)
@@ -0,0 +1,193 @@
+double form_volume(void);
+
+double Iq(double q, 
+    double lg_radius, double sm_radius,
+    double lg_vol_frac, double sm_vol_frac,
+    double lg_sld, double sm_sld, double solvent_sld
+    );
+    
+double Iqxy(double qx, double qy,
+    double lg_radius, double sm_radius,
+    double lg_vol_frac, double sm_vol_frac,
+    double lg_sld, double sm_sld, double solvent_sld
+    );
+
+int calculate_psfs(double qval,
+    double r2, double nf2,
+    double aa, double phi,
+    double *s11, double *s22, double *s12
+    );
+
+double form_volume(void)
+{
+    return 1.0;
+}
+
+double Iq(double q,
+    double lg_radius, double sm_radius,
+    double lg_vol_frac, double sm_vol_frac,
+    double lg_sld, double sm_sld, double solvent_sld)
+    
+{
+    double r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten;       //my local names
+    double psf11,psf12,psf22;
+    double phi1,phi2,phr,a3;
+    double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2;
+    int err;
+    
+    r2 = lg_radius;
+    r1 = sm_radius;
+    phi2 = lg_vol_frac;
+    phi1 = sm_vol_frac;
+    rho2 = lg_sld;
+    rho1 = sm_sld;
+    rhos = solvent_sld;
+    
+    
+    phi = phi1 + phi2;
+    aa = r1/r2;
+    //calculate the number fraction of larger spheres (eqn 2 in reference)
+    a3=aa*aa*aa;
+    phr=phi2/phi;
+    nf2 = phr*a3/(1.0-phr+phr*a3);
+    // calculate the PSF's here
+    err = calculate_psfs(q,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
+    
+    // /* do form factor calculations  */
+    
+    v1 = 4.0*M_PI/3.0*r1*r1*r1;
+    v2 = 4.0*M_PI/3.0*r2*r2*r2;
+    
+    n1 = phi1/v1;
+    n2 = phi2/v2;
+    
+    qr1 = r1*q;
+    qr2 = r2*q;
+
+    //if (qr1 == 0){
+        //sc1 = 1.0/3.0;
+    //}else{
+        //sc1 = (sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1;
+    //}
+    //if (qr2 == 0){
+        //sc2 = 1.0/3.0;
+    //}else{
+        //sc2 = (sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2;
+    //}
+    sc1 = sph_j1c(qr1);
+    sc2 = sph_j1c(qr2);
+    b1 = r1*r1*r1*(rho1-rhos)*4.0/3.0*M_PI*sc1;
+    b2 = r2*r2*r2*(rho2-rhos)*4.0/3.0*M_PI*sc2;
+    inten = n1*b1*b1*psf11;
+    inten += sqrt(n1*n2)*2.0*b1*b2*psf12;
+    inten += n2*b2*b2*psf22;
+    ///* convert I(1/A) to (1/cm)  */
+    inten *= 1.0e8;
+    ///*convert rho^2 in 10^-6A to A*/
+    inten *= 1.0e-12;    
+    return(inten);
+}
+
+
+double Iqxy(double qx, double qy,
+    double lg_radius, double sm_radius,
+    double lg_vol_frac, double sm_vol_frac,
+    double lg_sld, double sm_sld, double solvent_sld)
+    
+{
+    double q = sqrt(qx*qx + qy*qy);
+    return Iq(q,
+        lg_radius, sm_radius,
+        lg_vol_frac, sm_vol_frac,
+        lg_sld, sm_sld, solvent_sld);
+}
+
+int calculate_psfs(double qval,
+    double r2, double nf2,
+    double aa, double phi,
+    double *s11, double *s22, double *s12)
+
+{
+    //  variable qval,r2,nf2,aa,phi,&s11,&s22,&s12
+    
+    //   calculate constant terms
+    double s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4;
+    double a1,a2i,a2,b1,b2,b12,gm1,gm12;
+    double err=0.0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3;
+    double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13;
+    double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1;
+    
+    s2 = 2.0*r2;
+    s1 = aa*s2;
+    v = phi;
+    a3 = aa*aa*aa;
+    v1=((1.-nf2)*a3/(nf2+(1.-nf2)*a3))*v;
+    v2=(nf2/(nf2+(1.-nf2)*a3))*v;
+    g11=((1.+.5*v)+1.5*v2*(aa-1.))/(1.-v)/(1.-v);
+    g22=((1.+.5*v)+1.5*v1*(1./aa-1.))/(1.-v)/(1.-v);
+    g12=((1.+.5*v)+1.5*(1.-aa)*(v1-v2)/(1.+aa))/(1.-v)/(1.-v);
+    wmv = 1/(1.-v);
+    wmv3 = wmv*wmv*wmv;
+    wmv4 = wmv*wmv3;
+    a1=3.*wmv4*((v1+a3*v2)*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2))) + ((v1+a3*v2)*(1.+2.*v)+(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)-3.*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*wmv3;
+    a2i=((v1+a3*v2)*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*3*wmv4 + ((v1+a3*v2)*(1.+2.*v)+a3*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*aa-3.*v1*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*wmv3;
+    a2=a2i/a3;
+    b1=-6.*(v1*g11*g11+.25*v2*(1.+aa)*(1.+aa)*aa*g12*g12);
+    b2=-6.*(v2*g22*g22+.25*v1/a3*(1.+aa)*(1.+aa)*g12*g12);
+    b12=-3.*aa*(1.+aa)*(v1*g11/aa/aa+v2*g22)*g12;
+    gm1=(v1*a1+a3*v2*a2)*.5;
+    gm12=2.*gm1*(1.-aa)/aa;
+    //c  
+    //c   calculate the direct correlation functions and print results
+    //c
+    //  do 20 j=1,npts
+    
+    yy=qval*s2;
+    //c   calculate direct correlation functions
+    //c   ----c11
+    ay=aa*yy;
+    ay2 = ay*ay;
+    ay3 = ay*ay*ay;
+    t1=a1*(sin(ay)-ay*cos(ay));
+    t2=b1*(2.*ay*sin(ay)-(ay2-2.)*cos(ay)-2.)/ay;
+    t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3;
+    f11=24.*v1*(t1+t2+t3)/ay3;
+    
+    //c ------c22
+    y2=yy*yy;
+    y3=yy*y2;
+    tt1=a2*(sin(yy)-yy*cos(yy));
+    tt2=b2*(2.*yy*sin(yy)-(y2-2.)*cos(yy)-2.)/yy;
+    tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3;
+    f22=24.*v2*(tt1+tt2+tt3)/y3;
+    
+    //c   -----c12
+    yl=.5*yy*(1.-aa);
+    yl3=yl*yl*yl;
+    wma3 = (1.-aa)*(1.-aa)*(1.-aa);
+    y1=aa*yy;
+    y13 = y1*y1*y1;
+    ttt1=3.*wma3*v*sqrt(nf2)*sqrt(1.-nf2)*a1*(sin(yl)-yl*cos(yl))/((nf2+(1.-nf2)*a3)*yl3);
+    t21=b12*(2.*y1*cos(y1)+(y1*y1-2.)*sin(y1));
+    t22=gm12*((3.*y1*y1-6.)*cos(y1)+(y1*y1*y1-6.*y1)*sin(y1)+6.)/y1;
+    t23=gm1*((4.*y13-24.*y1)*cos(y1)+(y13*y1-12.*y1*y1+24.)*sin(y1))/(y1*y1);
+    t31=b12*(2.*y1*sin(y1)-(y1*y1-2.)*cos(y1)-2.);
+    t32=gm12*((3.*y1*y1-6.)*sin(y1)-(y1*y1*y1-6.*y1)*cos(y1))/y1;
+    t33=gm1*((4.*y13-24.*y1)*sin(y1)-(y13*y1-12.*y1*y1+24.)*cos(y1)+24.)/(y1*y1);
+    t41=cos(yl)*((sin(y1)-y1*cos(y1))/(y1*y1) + (1.-aa)/(2.*aa)*(1.-cos(y1))/y1);
+    t42=sin(yl)*((cos(y1)+y1*sin(y1)-1.)/(y1*y1) + (1.-aa)/(2.*aa)*sin(y1)/y1);
+    ttt2=sin(yl)*(t21+t22+t23)/(y13*y1);
+    ttt3=cos(yl)*(t31+t32+t33)/(y13*y1);
+    ttt4=a1*(t41+t42)/y1;
+    f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3);
+    
+    c11=f11;
+    c22=f22;
+    c12=f12;
+    *s11=1./(1.+c11-(c12)*c12/(1.+c22));
+    *s22=1./(1.+c22-(c12)*c12/(1.+c11)); 
+    *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));   
+    
+    return(err);
+}
+
Index: sasmodels/models/binary_hard_sphere.py
===================================================================
--- sasmodels/models/binary_hard_sphere.py	(revision 115d0f015d008c59c827203a96ee8894b7ec49b4)
+++ sasmodels/models/binary_hard_sphere.py	(revision 115d0f015d008c59c827203a96ee8894b7ec49b4)
@@ -0,0 +1,128 @@
+r"""
+
+The binary hard sphere model provides the scattering intensity, for binary
+mixture of hard spheres including hard sphere interaction between those
+particles, using rhw Percus-Yevick closure. The calculation is an exact
+multi-component solution that properly accounts for the 3 partial structure
+factors.
+
+Definition
+----------
+
+.. math::
+
+    \begin{eqnarray}
+    I(q) = (1-x)f_1^2(q) S_{11}(q) + 2[x(1-x)]^{1/2} f_1(q)f_2(q)S_{12}(q) + \\
+    x\,f_2^2(q)S_{22}(q) \\
+    \end{eqnarray}
+
+where $S_{ij}$ are the partial structure factors and $f_i$ are the scattering
+amplitudes of the particles. The subscript 1 is for the smaller particle and 2
+is for the larger. The number fraction of the larger particle,
+($x = n2/(n1+n2)$, where $n$ = the number density) is internally calculated
+based on the diameter ratio and the volume fractions.
+
+.. math::
+
+    \begin{eqnarray}
+    x &=& \frac{(\phi_2 / \phi)\alpha^3}{(1-(\phi_2/\phi) + (\phi_2/\phi) \\
+    \alpha^3)} \phi &=& \phi_1 + \phi_2 = \text{total volume fraction} \\
+    \alpha &=& R_1/R_2 = \text{size ratio}
+    \end{eqnarray}
+
+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}
+
+
+**NOTE 1:** The volume fractions and the scattering contrasts are loosely
+correlated, so holding as many parameters fixed to known values during fitting
+will improve the robustness of the fit.
+
+**NOTE 2:** Since the calculation uses the Percus-Yevick closure, all of the
+limitations of that closure relation apply here. Specifically, one should be
+wary of results for (total) volume fractions greater than approximately 40%.
+Depending on the size ratios or number fractions, the limit on total volume
+fraction may be lower.
+
+**NOTE 3:** The heavy arithmatic operations also mean that at present the
+function is poorly behaved at very low qr.  In some cases very large qr may
+also be poorly behaved.  These should however be outside any useful region of
+qr.
+
+The code for this model is based originally on a c-library implementation by the
+NIST Center for Neutron Research (Kline, 2006).
+
+See the references for details.
+
+References
+----------
+
+N W Ashcroft and D C Langreth, *Physical Review*, 156 (1967) 685-692
+[Errata found in *Phys. Rev.* 166 (1968) 934]
+
+S R Kline, *J Appl. Cryst.*, 39 (2006) 895
+
+**Author:** N/A **on:**
+
+**Modified by:** Paul Butler **on:** March 18, 2016
+
+**Reviewed by:** Paul Butler **on:** March 18, 2016
+"""
+
+import numpy as np
+from numpy import pi, inf
+
+name = "binary_hard_sphere"
+title = "binary mixture of hard spheres with hard sphere interactions."
+description = """Describes the scattering from a mixture of two distinct
+monodisperse, hard sphere particles.
+        [Parameters];
+        radius_lg: large radius of binary hard sphere,
+        radius_sm: small radius of binary hard sphere,
+        volfraction_lg: volume fraction of large spheres,
+        volfraction_sm: volume fraction of small spheres,
+        sld_lg: large sphere  scattering length density,
+        sld_sm: small sphere scattering length density,
+        sld_solvent: solvent scattering length density.
+"""
+category = "shape:sphere"
+
+#             ["name", "units", default, [lower, upper], "type", "description"],
+parameters = [["radius_lg", "Ang", 100, [0, inf], "",
+               "radius of large particle"],
+              ["radius_sm", "Ang", 25, [0, inf], "",
+               "radius of small particle"],
+              ["volfraction_lg", "", 0.1, [0, 1], "",
+               "volume fraction of large particle"],
+              ["volfraction_sm", "", 0.2, [0, 1], "",
+               "volume fraction of small particle"],
+              ["sld_lg", "1e-6/Ang^2", 3.5, [-inf, inf], "",
+               "scattering length density of large particle"],
+              ["sld_sm", "1e-6/Ang^2", 0.5, [-inf, inf], "",
+               "scattering length density of small particle"],
+              ["sld_solvent", "1e-6/Ang^2", 6.36, [-inf, inf], "",
+               "Solvent scattering length density"],
+             ]
+
+source = ["lib/sph_j1c.c", "binary_hard_sphere.c"]
+
+# parameters for demo and documentation
+demo = dict(scale=1, background=0,
+            sld_lg=3.5, sld_sm=0.5, sld_solvent=6.36,
+            radius_lg=100, radius_sm=20,
+            volfraction_lg=0.1, volfraction_sm=0.2)
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname = 'BinaryHSModel'
+oldpars = dict(volfraction_lg='vol_frac_ls', volfraction_sm='vol_frac_ss',
+               sld_lg='ls_sld', sld_sm='ss_sld', sld_solvent='solvent_sld',
+               radius_lg='l_radius', radius_sm='s_radius')
+
+# NOTE: test results taken from values returned by SasView 3.1.2
+tests = [[{'scale':1.0}, 0.001, 25.8927262013]]
+
Index: sasmodels/models/hardsphere.py
===================================================================
--- sasmodels/models/hardsphere.py	(revision 3bcd03d7e1d8d4938d30f17bf360b83417139376)
+++ sasmodels/models/hardsphere.py	(revision e98c1e008cb289d6be876cc1dfb5c76efa378204)
@@ -58,4 +58,5 @@
 category = "structure-factor"
 structure_factor = True
+single = False
 
 #             ["name", "units", default, [lower, upper], "type","description"],
