Changeset eec5fd6 in sasmodels for sasmodels/models

Timestamp:

Mar 18, 2016 4:58:02 AM (8 years ago)

Author:

richardh

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

ead4bd5

Parents:

4f4e0d5 (diff), 4554131 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'master' of ​https://github.com/SasView/sasmodels

Location:

sasmodels/models

Files:

• bessel.py (modified) (2 diffs)
• lib/j0_cephes.c (modified) (3 diffs)
• lib/j1_cephes.c (modified) (1 diff)
• lib/polevl.c (modified) (1 diff)
• poly_gauss_coil.py (modified) (3 diffs)
• adsorbed_layer.py (modified) (6 diffs)

sasmodels/models/bessel.py

@@ -67 +67 @@
 #Bessel
 parameters = [
+    ["ignored", "", 0.0, [-inf, inf], "", "no parameterless functions"],
              ]
 
-source = ["lib/polevl.c", "lib/j1d.c"]
+source = ["lib/polevl.c", "lib/j1_cephes.c"]
 
 # No volume normalization despite having a volume parameter
@@ -77 +78 @@
 
 Iq = """
-    return j1(q);
+    return 2.0*j1(q)/q;
     """
 

sasmodels/models/lib/j0_cephes.c

@@ -44 +44 @@
  */
 
-/*                                                      y0.c
- *
- *      Bessel function of the second kind, order zero
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, y0();
- *
- * y = y0( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns Bessel function of the second kind, of order
- * zero, of the argument.
- *
- * The domain is divided into the intervals [0, 5] and
- * (5, infinity). In the first interval a rational approximation
- * R(x) is employed to compute
- *   y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.
- * Thus a call to j0() is required.
- *
- * In the second interval, the Hankel asymptotic expansion
- * is employed with two rational functions of degree 6/6
- * and 7/7.
- *
- *
- *
- * ACCURACY:
- *
- *  Absolute error, when y0(x) < 1; else relative error:
- *
- * arithmetic   domain     # trials      peak         rms
- *    DEC       0, 30        9400       7.0e-17     7.9e-18
- *    IEEE      0, 30       30000       1.3e-15     1.6e-16
- *
- */
-
 
 /*
@@ -95 +54 @@
 
 double j0( double );
-
 double j0(double x) {
 
@@ -291 +249 @@
 
     q = 1.0/x;
-    w = sqrtf(q);
+    w = sqrt(q);
 
     p = w * polevl( q, MO, 7);
     w = q*q;
     xn = q * polevl( w, PH, 7) - PIO4F;
-    p = p * cosf(xn + xx);
+    p = p * cos(xn + xx);
     return(p);
 #endif

sasmodels/models/lib/j1_cephes.c

@@ -32 +32 @@
  *    IEEE      0, 30       30000       2.6e-16     1.1e-16
  *
- *
- */
-/*                                                      y1.c
- *
- *      Bessel function of second kind of order one
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, y1();
- *
- * y = y1( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns Bessel function of the second kind of order one
- * of the argument.
- *
- * The domain is divided into the intervals [0, 8] and
- * (8, infinity). In the first interval a 25 term Chebyshev
- * expansion is used, and a call to j1() is required.
- * In the second, the asymptotic trigonometric representation
- * is employed using two rational functions of degree 5/5.
- *
- *
- *
- * ACCURACY:
- *
- *                      Absolute error:
- * arithmetic   domain      # trials      peak         rms
- *    DEC       0, 30       10000       8.6e-17     1.3e-17
- *    IEEE      0, 30       30000       1.0e-15     1.3e-16
- *
- * (error criterion relative when |y1| > 1).
  *
  */

sasmodels/models/lib/polevl.c

@@ -50 +50 @@
 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 */
+
 double polevl( double x, double coef[8], int N );
 double p1evl( double x, double coef[8], int N );

sasmodels/models/poly_gauss_coil.py

@@ -50 +50 @@
 """
 
-from numpy import inf, sqrt, power
+from numpy import inf, sqrt, exp, power
 
 name =  "poly_gauss_coil"
@@ -69 +69 @@
 def Iq(q, i_zero, radius_gyration, polydispersity):
     # pylint: disable = missing-docstring
-    # need to trap the case of the polydispersity being 1 (ie, monodispersity)
     u = polydispersity - 1.0
-    if polydispersity == 1:
-       minusoneonu = -1.0 / u
-    else:
-       minusoneonu = -1.0 / u
     z = ((q * radius_gyration) * (q * radius_gyration)) / (1.0 + 2.0 * u)
     if (q == 0).any():
-       inten = i_zero
+        inten = i_zero
     else:
-       inten = i_zero * 2.0 * (power((1.0 + u * z),minusoneonu) + z - 1.0 ) / ((1.0 + u) * (z * z))
+    # need to trap the case of the polydispersity being 1 (ie, monodispersity!)
+        if polydispersity == 1:
+            inten = i_zero * 2.0 * (exp(-z) + z - 1.0 ) / (z * z)
+        else:
+            minusoneonu = -1.0 / u
+            inten = i_zero * 2.0 * (power((1.0 + u * z),minusoneonu) + z - 1.0 ) / ((1.0 + u) * (z * z))
     return inten
-Iq.vectorized =  True  # Iq accepts an array of q values
+#Iq.vectorized =  True  # Iq accepts an array of q values
 
 def Iqxy(qx, qy, *args):
@@ -100 +100 @@
                background = 'background')
 
+# these unit test values taken from SasView 3.1.2
 tests =  [
-    [{'scale': 70.0, 'radius_gyration': 75.0, 'polydispersity': 2.0, 'background': 0.0},
+    [{'scale': 1.0, 'i_zero': 70.0, 'radius_gyration': 75.0, 'polydispersity': 2.0, 'background': 0.0},
      [0.0106939, 0.469418], [57.6405, 0.169016]],
     ]

sasmodels/models/adsorbed_layer.py

@@ -6 +6 @@
 
 r"""
-This model describes the scattering from a layer of surfactant or polymer adsorbed on spherical particles under the conditions that (i) the particles (cores) are contrast-matched to the dispersion medium, (ii) *S(Q)* ~ 1 (ie, the particle volume fraction is dilute), (iii) the particle radius is >> layer thickness (ie, the interface is locally flat), and (iv) scattering from excess unadsorbed adsorbate in the bulk medium is absent or has been corrected for.
+This model describes the scattering from a layer of surfactant or polymer adsorbed on large, smooth, notionally spherical particles under the conditions that (i) the particles (cores) are contrast-matched to the dispersion medium, (ii) *S(Q)* ~ 1 (ie, the particle volume fraction is dilute), (iii) the particle radius is >> layer thickness (ie, the interface is locally flat), and (iv) scattering from excess unadsorbed adsorbate in the bulk medium is absent or has been corrected for.
 
 Unlike many other core-shell models, this model does not assume any form for the density distribution of the adsorbed species normal to the interface (cf, a core-shell model normally assumes the density distribution to be a homogeneous step-function). For comparison, if the thickness of a (traditional core-shell like) step function distribution is *t*, the second moment about the mean of the density distribution (ie, the distance of the centre-of-mass of the distribution from the interface), |sigma| =  sqrt((*t* :sup:`2` )/12).
@@ -34 +34 @@
 
 description =  """
-    Evaluates the scattering from particles
+    Evaluates the scattering from large particles
     with an adsorbed layer of surfactant or
     polymer, independent of the form of the
@@ -42 +42 @@
 
 #             ["name", "units", default, [lower, upper], "type", "description"],
-parameters =  [["second_moment", "Ang", 23.0, [0.0, inf], "", "Second moment"],
-              ["adsorbed_amount", "mg/m2", 1.9, [0.0, inf], "", "Adsorbed amount"],
-              ["density_poly", "g/cm3", 0.7, [0.0, inf], "", "Polymer density"],
-              ["radius", "Ang", 500.0, [0.0, inf], "", "Particle radius"],
-              ["volfraction", "None", 0.14, [0.0, inf], "", "Particle vol fraction"],
-              ["polymer_sld", "1/Ang^2", 1.5e-06, [-inf, inf], "", "Polymer SLD"],
-              ["solvent_sld", "1/Ang^2", 6.3e-06, [-inf, inf], "", "Solvent SLD"]]
+parameters =  [["second_moment", "Ang", 23.0, [0.0, inf], "", "Second moment of polymer distribution"],
+              ["adsorbed_amount", "mg/m2", 1.9, [0.0, inf], "", "Adsorbed amount of polymer"],
+              ["density_shell", "g/cm3", 0.7, [0.0, inf], "", "Bulk density of polymer in the shell"],
+              ["radius", "Ang", 500.0, [0.0, inf], "", "Core particle radius"],
+              ["volfraction", "None", 0.14, [0.0, inf], "", "Core particle volume fraction"],
+              ["sld_shell", "1e-6/Ang^2", 1.5, [-inf, inf], "sld", "Polymer shell SLD"],
+              ["sld_solvent", "1e-6/Ang^2", 6.3, [-inf, inf], "sld", "Solvent SLD"]]
 
 # NB: Scale and Background are implicit parameters on every model
-def Iq(q, second_moment, adsorbed_amount, density_poly, radius, 
-        volfraction, polymer_sld, solvent_sld):
+def Iq(q, second_moment, adsorbed_amount, density_shell, radius, 
+        volfraction, sld_shell, sld_solvent):
     # pylint: disable = missing-docstring
-    deltarhosqrd =  (polymer_sld - solvent_sld) * (polymer_sld - solvent_sld)
-    numerator =  6.0 * pi * volfraction * (adsorbed_amount * adsorbed_amount)
-    denominator =  (q * q) * (density_poly * density_poly) * radius
-    eterm =  exp(-1.0 * (q * q) * (second_moment * second_moment))
-    #scale by 10^10 for units conversion to cm^-1
-    inten =  1.0e+10 * deltarhosqrd * ((numerator / denominator) * eterm)
+#    deltarhosqrd =  (sld_shell - sld_solvent) * (sld_shell - sld_solvent)
+#    numerator =  6.0 * pi * volfraction * (adsorbed_amount * adsorbed_amount)
+#    denominator =  (q * q) * (density_shell * density_shell) * radius
+#    eterm =  exp(-1.0 * (q * q) * (second_moment * second_moment))
+#    #scale by 10^-2 for units conversion to cm^-1
+#    inten =  1.0e-02 * deltarhosqrd * ((numerator / denominator) * eterm)
+    aa =  (sld_shell - sld_solvent) * adsorbed_amount / q / density_shell 
+    bb = q * second_moment
+    #scale by 10^-2 for units conversion to cm^-1
+    inten =  6.0e-02 * pi * volfraction * aa * aa * exp(-bb * bb) / radius
     return inten
 Iq.vectorized =  True  # Iq accepts an array of q values
@@ -71 +75 @@
             second_moment = 23.0,
             adsorbed_amount = 1.9,
-            density_poly = 0.7,
+            density_shell = 0.7,
             radius = 500.0,
             volfraction = 0.14,
-            polymer_sld = 1.5e-06,
-            solvent_sld = 6.3e-06,
+            sld_shell = 1.5,
+            sld_solvent = 6.3,
             background = 0.0)
 
@@ -82 +86 @@
                second_moment = 'second_moment',
                adsorbed_amount = 'ads_amount',
-               density_poly = 'density_poly',
+               density_shell = 'density_poly',
                radius = 'radius_core',
                volfraction = 'volf_cores',
-               polymer_sld = 'sld_poly',
-               solvent_sld = 'sld_solv',
+               sld_shell = 'sld_poly',
+               sld_solvent = 'sld_solv',
                background = 'background')
 
@@ -92 +96 @@
 tests =  [
     [{'scale': 1.0, 'second_moment': 23.0, 'adsorbed_amount': 1.9, 
-     'density_poly': 0.7, 'radius': 500.0, 'volfraction': 0.14, 
-     'polymer_sld': 1.5e-06, 'solvent_sld': 6.3e-06, 'background': 0.0},
+     'density_shell': 0.7, 'radius': 500.0, 'volfraction': 0.14, 
+     'sld_shell': 1.5, 'sld_solvent': 6.3, 'background': 0.0},
      [0.0106939, 0.469418], [73.741, 9.65391e-53]],
     ]
+# ADDED by: SMK  ON: 16Mar2016  convert from sasview, check vs SANDRA, 18Mar2016 RKH some edits & renaming