Index: sasmodels/models/core_shell_ellipsoid.c
===================================================================
--- sasmodels/models/core_shell_ellipsoid.c	(revision 81dd61963c319ce8b1c38d14a6f3056a517fccdd)
+++ sasmodels/models/core_shell_ellipsoid.c	(revision 81dd61963c319ce8b1c38d14a6f3056a517fccdd)
@@ -0,0 +1,175 @@
+double form_volume(double equat_core,
+                   double polar_core,
+                   double equat_shell,
+                   double polar_shell);
+double Iq(double q,
+          double equat_core,
+          double polar_core,
+          double equat_shell,
+          double polar_shell,
+          double core_sld,
+          double shell_sld,
+          double solvent_sld);
+
+
+double Iqxy(double qx, double qy,
+          double equat_core,
+          double polar_core,
+          double equat_shell,
+          double polar_shell,
+          double core_sld,
+          double shell_sld,
+          double solvent_sld,
+          double theta,
+          double phi);
+
+
+double form_volume(double equat_core,
+                   double polar_core,
+                   double equat_shell,
+                   double polar_shell)
+{
+    double vol = 4.0*M_PI/3.0*equat_shell*equat_shell*polar_shell;
+    return vol;
+}
+
+static double
+core_shell_ellipsoid_kernel(double q,
+          double equat_core,
+          double polar_core,
+          double equat_shell,
+          double polar_shell,
+          double core_sld,
+          double shell_sld,
+          double solvent_sld)
+{
+	double delpc,delps;
+	double uplim,lolim;		//upper and lower integration limits
+	double summ,zi,yyy,answer; //running tally of integration
+
+	lolim = 0.0;
+	uplim = 1.0;
+
+	summ = 0.0;	 //initialize intergral
+
+	delpc = core_sld - shell_sld; //core - shell
+	delps = shell_sld - solvent_sld; //shell - solvent
+
+	for(int i=0;i<76;i++) {
+		zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
+		yyy = Gauss76Wt[i] * gfn4(zi,
+		                          equat_core,
+		                          polar_core,
+		                          equat_shell,
+		                          polar_shell,
+		                          delpc,
+		                          delps,
+		                          q);
+		summ += yyy;
+	}
+
+	answer = (uplim-lolim)/2.0*summ;
+
+	//convert to [cm-1]
+	answer *= 1.0e-4;
+
+	return answer;
+}
+
+static double
+core_shell_ellipsoid_kernel_2d(double q, double q_x, double q_y,
+          double equat_core,
+          double polar_core,
+          double equat_shell,
+          double polar_shell,
+          double core_sld,
+          double shell_sld,
+          double solvent_sld,
+          double theta,
+          double phi)
+{
+    double cyl_x, cyl_y;
+    double cos_val;
+    double answer;
+    double sldcs,sldss;
+
+    //convert angle degree to radian
+    theta = theta * M_PI/180.0;
+    phi = phi * M_PI/180.0;
+
+
+    // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis.
+    cyl_x = cos(theta) * cos(phi);
+    cyl_y = sin(theta);
+
+    sldcs = core_sld - shell_sld;
+    sldss = shell_sld- solvent_sld;
+
+    // Compute the angle btw vector q and the
+    // axis of the cylinder
+    cos_val = cyl_x*q_x + cyl_y*q_y;
+
+    // Call the IGOR library function to get the kernel: MUST use gfn4 not gf2 because of the def of params.
+    answer = gfn4(cos_val,
+                  equat_core,
+                  polar_core,
+                  equat_shell,
+                  polar_shell,
+                  sldcs,
+                  sldss,
+                  q);
+
+    //convert to [cm-1]
+    answer *= 1.0e-4;
+
+    return answer;
+}
+
+double Iq(double q,
+          double equat_core,
+          double polar_core,
+          double equat_shell,
+          double polar_shell,
+          double core_sld,
+          double shell_sld,
+          double solvent_sld)
+{
+    double intensity = core_shell_ellipsoid_kernel(q,
+           equat_core,
+           polar_core,
+           equat_shell,
+           polar_shell,
+           core_sld,
+           shell_sld,
+           solvent_sld);
+
+    return intensity;
+}
+
+
+double Iqxy(double qx, double qy,
+          double equat_core,
+          double polar_core,
+          double equat_shell,
+          double polar_shell,
+          double core_sld,
+          double shell_sld,
+          double solvent_sld,
+          double theta,
+          double phi)
+{
+    double q;
+    q = sqrt(qx*qx+qy*qy);
+    double intensity = core_shell_ellipsoid_kernel_2d(q, qx/q, qy/q,
+                       equat_core,
+                       polar_core,
+                       equat_shell,
+                       polar_shell,
+                       core_sld,
+                       shell_sld,
+                       solvent_sld,
+                       theta,
+                       phi);
+
+    return intensity;
+}
Index: sasmodels/models/core_shell_ellipsoid.py
===================================================================
--- sasmodels/models/core_shell_ellipsoid.py	(revision 81dd61963c319ce8b1c38d14a6f3056a517fccdd)
+++ sasmodels/models/core_shell_ellipsoid.py	(revision 81dd61963c319ce8b1c38d14a6f3056a517fccdd)
@@ -0,0 +1,178 @@
+r"""
+This model provides the form factor, $P(q)$, for a core shell ellipsoid (below)
+where the form factor is normalized by the volume of the cylinder.
+
+.. math::
+
+    P(q) = scale * \left<f^2\right>/V + background
+
+where the volume $V = (4/3)\pi(r_{maj}r_{min}^2)$ and the averaging $< >$ is
+applied over all orientations for 1D.
+
+.. figure:: img/core_shell_ellipsoid_fig1.gif
+
+    The returned value is in units of $cm^{-1}$, on absolute scale.
+
+Definition
+----------
+
+The form factor calculated is
+
+.. math::
+
+    P(q) = \frac{scale}{V}\int_0^1
+        \left|F(q,r_{min},r_{max},\alpha)\right|^2d\alpha + background
+
+    \left|F(q,r_{min},r_{max},\alpha)\right|=V\Delta \rho \cdot (3j_1(u)/u)
+
+    u = q\left[ r_{maj}^2\alpha ^2 + r_{min}^2(1-\alpha ^2)\right]^{1/2}
+
+where
+
+.. math::
+
+    j_1(u)=(\sin x - x \cos x)/x^2
+
+To provide easy access to the orientation of the core-shell ellipsoid,
+we define the axis of the solid ellipsoid using two angles $\theta$ and $\phi$.
+These angles are defined on Figure 2 of the CylinderModel.
+The contrast is defined as SLD(core) - SLD(shell) and SLD(shell) - SLD(solvent).
+
+In the parameters, *equat_core* = equatorial core radius, *polar_core* =
+polar core radius, *equat_shell* = $r_{min}$ (or equatorial outer radius),
+and *polar_shell* = $r_{maj}$ (or polar outer radius).
+
+.. note::
+    The 2nd virial coefficient of the solid ellipsoid is calculated based on
+    the *radius_a* (= *polar_shell)* and *radius_b (= equat_shell)* values,
+    and used as the effective radius for *S(Q)* when $P(Q) * S(Q)$ is applied.
+
+.. figure:: img/core_shell_ellipsoid_1d.jpg
+
+    1D plot using the default values (w/200 data point).
+
+.. figure:: img/core_shell_ellipsoid_fig2.jpg
+
+    The angles for oriented core_shell_ellipsoid.
+
+Our model uses the form factor calculations implemented in a c-library provided
+by the NIST Center for Neutron Research (Kline, 2006).
+
+References
+----------
+
+M Kotlarchyk, S H Chen, *J. Chem. Phys.*, 79 (1983) 2461
+
+S J Berr, *Phys. Chem.*, 91 (1987) 4760
+
+"""
+
+from numpy import inf, sin, cos, pi
+
+name = "core_shell_ellipsoid"
+title = "Form factor for an spheroid ellipsoid particle with a core shell structure."
+description = """
+    [SpheroidCoreShellModel] Calculates the form factor for an spheroid
+    ellipsoid particle with a core_shell structure.
+    The form factor is averaged over all possible
+    orientations of the ellipsoid such that P(q)
+    = scale*<f^2>/Vol + bkg, where f is the
+    single particle scattering amplitude.
+    [Parameters]:
+    equat_core = equatorial radius of core,
+    polar_core = polar radius of core,
+    equat_shell = equatorial radius of shell,
+    polar_shell = polar radius (revolution axis) of shell,
+    core_sld = SLD_core
+    shell_sld = SLD_shell
+    solvent_sld = SLD_solvent
+    background = Incoherent bkg
+    scale =scale
+    Note:It is the users' responsibility to ensure
+    that shell radii are larger than core radii.
+    oblate: polar radius < equatorial radius
+    prolate :  polar radius > equatorial radius
+    """
+category = "shape:ellipsoid"
+
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type", "description"],
+parameters = [
+    ["equat_core",  "Ang",      200,   [0, inf],    "volume",      "Equatorial radius of core"],
+    ["polar_core",  "Ang",       10,   [0, inf],    "volume",      "Polar radius of core"],
+    ["equat_shell", "Ang",      250,   [0, inf],    "volume",      "Equatorial radius of shell"],
+    ["polar_shell", "Ang",       30,   [0, inf],    "volume",      "Polar radius of shell"],
+    ["core_sld",    "1e-6/Ang^2", 2,   [-inf, inf], "",            "Core scattering length density"],
+    ["shell_sld",   "1e-6/Ang^2", 1,   [-inf, inf], "",            "Shell scattering length density"],
+    ["solvent_sld", "1e-6/Ang^2", 6.3, [-inf, inf], "",            "Solvent scattering length density"],
+    ["theta",       "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation wrt incoming beam"],
+    ["phi",         "degrees",    0,   [-inf, inf], "orientation", "Oblate orientation in the plane of the detector"],
+    ]
+# pylint: enable=bad-whitespace, line-too-long
+
+source = ["lib/gfn.c", "lib/gauss76.c", "core_shell_ellipsoid.c"]
+
+demo = dict(scale=1, background=0.001,
+            equat_core=200.0,
+            polar_core=10.0,
+            equat_shell=250.0,
+            polar_shell=30.0,
+            core_sld=2.0,
+            shell_sld=1.0,
+            solvent_sld=6.3,
+            theta=0,
+            phi=0)
+
+oldname = 'CoreShellEllipsoidModel'
+
+oldpars = dict(core_sld='sld_core',
+               shell_sld='sld_shell',
+               solvent_sld='sld_solvent',
+               theta='axis_theta',
+               phi='axis_phi')
+
+q = 0.1
+phi = pi/6
+qx = q*cos(phi)
+qy = q*sin(phi)
+
+tests = [
+    # Accuracy tests based on content in test/utest_other_models.py
+    [{'equat_core': 200.0,
+      'polar_core': 20.0,
+      'equat_shell': 250.0,
+      'polar_shell': 30.0,
+      'core_sld': 2.0,
+      'shell_sld': 1.0,
+      'solvent_sld': 6.3,
+      'background': 0.001,
+      'scale': 1.0,
+     }, 1.0, 0.00189402],
+
+    # Additional tests with larger range of parameters
+    [{'background': 0.01}, 0.1, 8.86741],
+
+    [{'equat_core': 20.0,
+      'polar_core': 200.0,
+      'equat_shell': 54.0,
+      'polar_shell': 3.0,
+      'core_sld': 20.0,
+      'shell_sld': 10.0,
+      'solvent_sld': 6.0,
+      'background': 0.0,
+      'scale': 1.0,
+     }, 0.01, 26150.4],
+
+    [{'background': 0.001}, (0.4, 0.5), 0.00170471],
+
+    [{'equat_core': 20.0,
+      'polar_core': 200.0,
+      'equat_shell': 54.0,
+      'polar_shell': 3.0,
+      'core_sld': 20.0,
+      'shell_sld': 10.0,
+      'solvent_sld': 6.0,
+      'background': 0.01,
+      'scale': 0.01,
+     }, (qx, qy), 0.105764],
+    ]
Index: sasmodels/models/lib/gfn.c
===================================================================
--- sasmodels/models/lib/gfn.c	(revision 81dd61963c319ce8b1c38d14a6f3056a517fccdd)
+++ sasmodels/models/lib/gfn.c	(revision 81dd61963c319ce8b1c38d14a6f3056a517fccdd)
@@ -0,0 +1,40 @@
+//
+//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN
+//                       BY (53) & (58-59) IN CHEN AND
+//                       KOTLARCHYK REFERENCE
+//
+//       <OBLATE ELLIPSOID>
+// function gfn4 for oblate ellipsoids
+double
+gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
+{
+	// local variables
+	double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi;
+
+	Pi = 4.0*atan(1.0);
+	pi43=4.0/3.0*Pi;
+  	aa = crmaj;
+ 	bb = crmin;
+ 	u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx));
+ 	ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx));
+   	uq = sqrt(u2)*qq;
+ 	ut= sqrt(ut2)*qq;
+	vc = pi43*aa*aa*bb;
+   	vt = pi43*trmaj*trmaj*trmin;
+   	if (uq == 0.0){
+   		siq = 1.0/3.0;
+   	}else{
+   		siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
+   	}
+   	if (ut == 0.0){
+   		sit = 1.0/3.0;
+   	}else{
+   		sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
+   	}
+   	gfnc = 3.0*siq*vc*delpc;
+  	gfnt = 3.0*sit*vt*delps;
+  	tgfn = gfnc+gfnt;
+  	gfn4 = tgfn*tgfn;
+
+  	return (gfn4);
+}
