Round to nearest integer for star polymer arms.

sasmodels/models/star_polymer.c

@@ -2,21 +2,21 @@ double form_volume(void);
 
 double Iq(double q, double rg_sq, double arms);
 
-static double star_polymer_kernel(double q, double rg_sq, double arms)
+static double star_polymer_kernel(double q, double rg_sq, int n_arms)
 {
     // Given the kernel function:
     //    I = 2(v + (e^-v - 1) + (a-1)/2 (e^-v - 1)^2) / av^2
     // At q = 0, u_2 and v are 0 so term1 = term2 = 0, hence 0/0 = nan
-    // We can clean up the values around low q using Taylor expansions for
+    // We can clean up the values at low q using Taylor expansions for
     //    T1 = (e^-v - 1)/v = [-1 +1/2 -1/6 +1/24 -1/120 +1/720 ...]
     //    T2 = [(e^-v - 1)/v]^2 = [+1 -1 +7/12 -1/4 +31/360 -1/40 ...]
     // Rewriting
     //    I = 2/av + 2/av (e^-v - 1)/v + (a-1)/a ((e^v - 1)/v)^2
     // Since the first term of T1 is -1 this cancels the leading 2/av.
     // Muliplying the remaining terms by 2/v yields T1' = 2*T1[1:]
-    //    I = 1/a T1'(v) + (1 - 1/a) T2(v)
+    //    I = 1/a T1'(v) + (1 - 1/a) T2(v) = (T1'(v) - T2(v))/a + T2(v)
     //    T1' = [+1 -1/3 +1/12 -1/60 +1/360 -1/2520]
-    // Substituting v=0 we get lim[q->0] I = 1/a + (1 - 1/a) = 1 as expected.
+    // Substituting v=0 we get lim[q->0] I = (1 - 1)/a + 1 = 1 as expected.
     // After testing it turns out that term2 is numerically stable, so we
     // only use the Taylor series for T1', and protect T2 against v = 0.
 
@@ -32,16 +32,15 @@ static double star_polymer_kernel(double q, double rg_sq, double arms)
 
     // Note: arms should be one or more
     const double u_sq = rg_sq * q * q;
-    const double v = u_sq * arms / (3.0 * arms - 2.0);
+    const double v = (u_sq * n_arms) / (3 * n_arms - 2);
     double term1;
     if (u_sq < _SPLIM) {
         term1 = 1. + v*(-1./3. + v*(1./12. + v*(-1./60. + v*(1./360. + v*(-1./2520)))));
     } else {
         term1 = 2.0 * (v + expm1(-v)) / (v * v);
     }
     double term2 = square(expm1(-v)/v); // we trap q == 0 above, so v > 0 here
-    return term1/arms + term2*(1.0 - 1.0/arms);
-
+    return (term1 - term2)/n_arms + term2;
 }
 
 double form_volume(void)
@@ -51,5 +50,9 @@ double form_volume(void)
 
 double Iq(double q, double rg_sq, double arms)
 {
-    return star_polymer_kernel(q, rg_sq, arms);
+    // TODO: round or truncate?
+    // rounding: stacked_disk pearl_necklace multilayer_vesicle linear_pearls
+    //      lamellar_stack_caille lamellar_hg_stack_caille
+    // truncation: lamellar_stack_paracrystal
+    return star_polymer_kernel(q, rg_sq, (int)(arms + 0.5));
 }

sasmodels/models/star_polymer.py

@@ -78,7 +78,7 @@
 # pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type","description"],
 parameters = [["rg_squared", "Ang^2", 100.0, [0.0, inf], "", "Ensemble radius of gyration SQUARED of the full polymer"],
-              ["arms",    "",      3,   [1.0, 6.0], "", "Number of arms in the model"],
+              ["arms",    "",      3,   [1.0, 6.0], "", "Number of arms in the model (integer)"],
              ]
 # pylint: enable=bad-whitespace, line-too-long
 
@@ -95,13 +95,17 @@ def random():
     return pars
 
 tests = [[{'rg_squared': 2.0,
-           'arms':    3.3,
-          }, 0.5, 0.851646091108],
+           'arms':    3.3, # should equal the value for arms=3 in the next test
+          }, 0.5, 0.8518871750393033],
+
+         [{'rg_squared': 2.0,
+           'arms':    3.0,
+          }, 0.5, 0.8518871750393033],
 
          [{'rg_squared':    1.0,
            'arms':       2.0,
            'background': 1.8,
           }, 1.0, 2.53575888234],
         ]
 # 2016-03-23 RKH edited docs, would this better use rg not rg^2 ? Numerical noise at extremely small q.rg
-# 2025-01-17 PAK Fixed numerical noise.
+# 2025-01-17 PAK Fixed numerical noise. Use nearest integer for number of arms.