more edits, add unit test core_multi_shell, lib funcs for ellipsoid

sasmodels/models/core_multi_shell.py

@@ -50,7 +50,7 @@
 from __future__ import division
 
 import numpy as np
-from numpy import inf
+from numpy import pi, inf
 
 name = "core_multi_shell"
 title = "This model provides the scattering from a spherical core with 1 to 10 \
@@ -146,3 +146,11 @@ def profile(sld_core, radius, sld_solvent, n, sld, thickness):
     rho.append(sld_solvent)
 
     return np.asarray(z), np.asarray(rho)
+
+# RKH 22Sept2023, copy in tests from core_shell_sphere, first tests the volume calc, second I(Q=0.4)
+# TODO needs more tests with further layers and S(Q)
+tests = [
+    [{'radius': 20.0, 'n':1, 'thickness': [10.0] }, 0.1, None, None, 30.0, 4.*pi/3*30**3, 1.0],
+    [{'radius': 60.0, 'n':1, 'thickness': [10.0], 'sld_core': 1.0, 'sld': [2.0],
+      'sld_solvent': 3.0, 'background': 0.001}, 0.4, .001698838]
+]

sasmodels/models/core_multi_shell_cylinder.py

@@ -2,7 +2,8 @@
 Definition
 ----------
 
-Core multi-shell cylinders or discs.
+Core multi-shell cylinders or discs, may be used for multi-walled tubes or
+a stck of discs.
 
 There must be a minumum of ONE shell (May set its sld to match solvent, 
 thickness1 and face1 to zero or very small so the shell does not contribute 
@@ -14,8 +15,8 @@
 
 2d scattering is so far only minimally tested.
 
-Use of S(Q) with this I(Q) currently not giving correct I(Q) in sasview, 
-(though passes unit tests), due to a more general sasview v5 bug.
+Use of S(Q) with this model currently may not give the correct I(Q) though 
+passes unit tests), due to a more general sasview v5 bug, so check carefully.
 
 The sld profile plot show profiles along both radius and half length 
 simultaneously! (A simple edit of the py code will change which is 
@@ -28,6 +29,14 @@
 
 $Router = radius\_core + thickness1 + thickness2 + thickness3 + ...$
 
+Each new shell must completely enclose the previous ones. Setting face1, 
+face2 etc. to very small values will give a multi-walled tube. By setting 
+some sld\_shell to match sld\_solvent then a hollow tube or multilammellar 
+tube may be generated. Conversely setting thickness1, thickness2 etc. to 
+very small values will generate a stack of discs, then by setting some 
+sld\_shell to match sld\_solvent a multilammellar disc stack.
+
+
 .. math::
 
 Scattered intensity is calculated by a numerical integration over angle 
@@ -83,10 +92,11 @@
 
 An approximation for the effects of "Gaussian interfacial roughness" $\sigma$
 is included, by multiplying $I(Q)$ by
-$\exp\left \{ -\frac{1}{2}Q^2\sigma^2 \right \}$. This applies, in some way, to
+$\exp\left \{ -\frac{1}{2}Q^2\sigma^2 \right \}$. This applies to
 all interfaces in the model not just the external ones. (Note that for a one
 dimensional system convolution of the scattering length density profile with
 a Gaussian of standard deviation $\sigma$ does exactly this multiplication.)
+To work well $\sigma$ must be much smaller than the particle size. 
 Leave $\sigma$ set to zero for the usual sharp interfaces.
 There is some debate as to whether the factor of 1/2 is needed or not.
 
@@ -120,8 +130,8 @@
 ----------------------------
 
 * **Author:** Richard Heenan **Date:** April 2021
-* **Last Modified by:** Richard Heenan **Date:** April 2021
-* **Last Reviewed by:** Richard Heenan **Date:** April 2021
+* **Last Modified by:** Richard Heenan **Date:** Sept 2023
+* **Last Reviewed by:**   **Date:** 
 """
 
 from __future__ import division

sasmodels/models/core_shell_ellipsoid.c

@@ -4,26 +4,7 @@
 //     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN
 //                       BY (53) & (58-59) IN CHEN AND
 //                       KOTLARCHYK REFERENCE
-//
-//       <OBLATE ELLIPSOID>
-static double
-_cs_ellipsoid_kernel(double qab, double qc,
-    double equat_core, double polar_core,
-    double equat_shell, double polar_shell,
-    double sld_core_shell, double sld_shell_solvent)
-{
-    const double qr_core = sqrt(square(equat_core*qab) + square(polar_core*qc));
-    const double si_core = sas_3j1x_x(qr_core);
-    const double volume_core = M_4PI_3*equat_core*equat_core*polar_core;
-    const double fq_core = si_core*volume_core*sld_core_shell;
-
-    const double qr_shell = sqrt(square(equat_shell*qab) + square(polar_shell*qc));
-    const double si_shell = sas_3j1x_x(qr_shell);
-    const double volume_shell = M_4PI_3*equat_shell*equat_shell*polar_shell;
-    const double fq_shell = si_shell*volume_shell*sld_shell_solvent;
-
-    return fq_core + fq_shell;
-}
+// see functions in /lib/cs_ellipsoid_funcs.c
 
 static double
 form_volume(double radius_equat_core,
@@ -44,28 +25,6 @@ radius_from_volume(double radius_equat_core, double x_core, double thick_shell,
     return cbrt(volume_ellipsoid/M_4PI_3);
 }
 
-static double
-radius_from_curvature(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell)
-{
-    // Trivial cases
-    if (1.0 == x_core && 1.0 == x_polar_shell) return radius_equat_core + thick_shell;
-    if ((radius_equat_core + thick_shell)*(radius_equat_core*x_core + thick_shell*x_polar_shell) == 0.)  return 0.;
-
-    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
-    const double radius_equat_tot = radius_equat_core + thick_shell;
-    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell;
-    const double ratio = (radius_polar_tot < radius_equat_tot
-                          ? radius_polar_tot / radius_equat_tot
-                          : radius_equat_tot / radius_polar_tot);
-    const double e1 = sqrt(1.0 - ratio*ratio);
-    const double b1 = 1.0 + asin(e1) / (e1 * ratio);
-    const double bL = (1.0 + e1) / (1.0 - e1);
-    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL);
-    const double delta = 0.75 * b1 * b2;
-    const double ddd = 2.0 * (delta + 1.0) * radius_polar_tot * radius_equat_tot * radius_equat_tot;
-    return 0.5 * cbrt(ddd);
-}
-
 static double
 radius_effective(int mode, double radius_equat_core, double x_core, double thick_shell, double x_polar_shell)
 {

sasmodels/models/core_shell_ellipsoid.py

@@ -100,7 +100,8 @@
 ----------------------------
 
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
-* **Last Modified by:** Richard Heenan (reparametrised model) **Date:** 2015
+* **Modified by:** Richard Heenan (reparametrised model) **Date:** 2015
+* **Last Modified by:** Richard Heenan (share functions with _tied odel) **Date:** 2023
 * **Last Reviewed by:** Steve King **Date:** March 27, 2019
 """
 
@@ -149,7 +150,7 @@
     ]
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"]
+source = ["lib/sas_3j1x_x.c", "lib/gauss76.c","lib/cs_ellipsoid_funcs.c", "core_shell_ellipsoid.c"]
 have_Fq = True
 radius_effective_modes = [
     "average outer curvature", "equivalent volume sphere",

sasmodels/models/core_shell_ellipsoid_tied.py

@@ -152,7 +152,7 @@
 
 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "lib/cubic_solve_reparam3.c", "lib/cs_ellipsoid_funcs.c",
  "core_shell_ellipsoid_tied.c"]
-# TODO get core_shell_ellipsoid to use the new cs/_ellipsoid_funcs
+
 have_Fq = True
 radius_effective_modes = [
     "average outer curvature", "equivalent volume sphere",