Changeset 7e923c2 in sasmodels for sasmodels/models

Timestamp:

Aug 7, 2018 10:45:52 PM (6 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

1a7ddc9

Parents:

c036ddb (diff), e9b17b18 (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' into beta_approx

Location:

sasmodels/models

Files:

• guinier.py (modified) (4 diffs)
• core_shell_sphere.py (modified) (1 diff)
• ellipsoid.c (modified) (2 diffs)
• ellipsoid.py (modified) (1 diff)
• lib/sphere_form.c (modified) (1 diff)
• sphere.c (added)
• sphere.py (modified) (1 diff)

sasmodels/models/guinier.py

@@ -7 +7 @@
 .. math::
 
-    I(q) = \text{scale} \cdot \exp{\left[ \frac{-Q^2R_g^2}{3} \right]}
+    I(q) = \text{scale} \cdot \exp{\left[ \frac{-Q^2 R_g^2 }{3} \right]}
             + \text{background}
 
@@ -19 +19 @@
 
 .. math:: q = \sqrt{q_x^2 + q_y^2}
+
+In scattering, the radius of gyration $R_g$ quantifies the objects's
+distribution of SLD (not mass density, as in mechanics) from the objects's
+SLD centre of mass. It is defined by
+
+.. math:: R_g^2 = \frac{\sum_i\rho_i\left(r_i-r_0\right)^2}{\sum_i\rho_i}
+
+where $r_0$ denotes the object's SLD centre of mass and $\rho_i$ is the SLD at
+a point $i$.
+
+Notice that $R_g^2$ may be negative (since SLD can be negative), which happens
+when a form factor $P(Q)$ is increasing with $Q$ rather than decreasing. This
+can occur for core/shell particles, hollow particles, or for composite
+particles with domains of different SLDs in a solvent with an SLD close to the
+average match point. (Alternatively, this might be regarded as there being an
+internal inter-domain "structure factor" within a single particle which gives
+rise to a peak in the scattering).
+
+To specify a negative value of $R_g^2$ in SasView, simply give $R_g$ a negative
+value ($R_g^2$ will be evaluated as $R_g |R_g|$). Note that the physical radius 
+of gyration, of the exterior of the particle, will still be large and positive. 
+It is only the apparent size from the small $Q$ data that will give a small or 
+negative value of $R_g^2$.
 
 References
@@ -42 +65 @@
 
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["rg", "Ang", 60.0, [0, inf], "", "Radius of Gyration"]]
+parameters = [["rg", "Ang", 60.0, [-inf, inf], "", "Radius of Gyration"]]
 
 Iq = """
-    double exponent = rg*rg*q*q/3.0;
+    double exponent = fabs(rg)*rg*q*q/3.0;
     double value = exp(-exponent);
     return value;
@@ -66 +89 @@
 
 # parameters for demo
-demo = dict(scale=1.0, rg=60.0)
+demo = dict(scale=1.0,  background=0.001, rg=60.0 )
 
 # parameters for unit tests

sasmodels/models/core_shell_sphere.py

@@ -58 +58 @@
 title = "Form factor for a monodisperse spherical particle with particle with a core-shell structure."
 description = """
-    F^2(q) = 3/V_s [V_c (sld_core-sld_shell) (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3
-                   + V_s (sld_shell-sld_solvent) (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3]
+    F(q) = [V_c (sld_core-sld_shell) 3 (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3
+            + V_s (sld_shell-sld_solvent) 3 (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3]
 
             V_s: Volume of the sphere shell

sasmodels/models/ellipsoid.c

@@ -5 +5 @@
 }
 
+/* Fq overrides Iq
 static  double
 Iq(double q,
@@ -37 +38 @@
     return 1.0e-4 * s * s * form;
 }
+*/
+
+static void
+Fq(double q,
+    double *F1,
+    double *F2,
+    double sld,
+    double sld_solvent,
+    double radius_polar,
+    double radius_equatorial)
+{
+    // Using ratio v = Rp/Re, we can implement the form given in Guinier (1955)
+    //     i(h) = int_0^pi/2 Phi^2(h a sqrt(cos^2 + v^2 sin^2) cos dT
+    //          = int_0^pi/2 Phi^2(h a sqrt((1-sin^2) + v^2 sin^2) cos dT
+    //          = int_0^pi/2 Phi^2(h a sqrt(1 + sin^2(v^2-1)) cos dT
+    // u-substitution of
+    //     u = sin, du = cos dT
+    //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du
+    const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0;
+    // translate a point in [-1,1] to a point in [0, 1]
+    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2;
+    const double zm = 0.5;
+    const double zb = 0.5;
+    double total_F2 = 0.0;
+    double total_F1 = 0.0;
+    for (int i=0;i<GAUSS_N;i++) {
+        const double u = GAUSS_Z[i]*zm + zb;
+        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one);
+        const double f = sas_3j1x_x(q*r);
+        total_F2 += GAUSS_W[i] * f * f;
+        total_F1 += GAUSS_W[i] * f;
+    }
+    // translate dx in [-1,1] to dx in [lower,upper]
+    const double form_squared_avg = total_F2*zm;
+    const double form_avg = total_F1*zm;
+    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
+    *F1 = 1e-2 * s * form_avg;
+    *F2 = 1e-4 * s * s * form_squared_avg;
+}
 
 static double

sasmodels/models/ellipsoid.py

@@ -161 +161 @@
              ]
 
+
 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"]
+
+have_Fq = True
 
 def ER(radius_polar, radius_equatorial):

sasmodels/models/lib/sphere_form.c

@@ -13 +13 @@
     return 1.0e-4*square(contrast * fq);
 }
+

sasmodels/models/sphere.py

@@ -69 +69 @@
 source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c"]
 
-# No volume normalization despite having a volume parameter
-# This should perhaps be volume normalized?
-form_volume = """
+c_code = """
+static double form_volume(double radius)
+{
     return sphere_volume(radius);
-    """
+}
 
-Iq = """
-    return sphere_form(q, radius, sld, sld_solvent);
-    """
+static void Fq(double q, double *F1,double *F2, double sld, double solvent_sld, double radius)
+{
+    const double fq = sas_3j1x_x(q*radius);
+    const double contrast = (sld - solvent_sld);
+    const double form = 1e-2 * contrast * sphere_volume(radius) * fq;
+    *F1 = form;
+    *F2 = form*form;
+}
+"""
+
+# TODO: figure this out by inspection
+have_Fq = True
 
 def ER(radius):