Changeset d119f34 in sasmodels for sasmodels/models/onion.c

Timestamp:

Aug 2, 2016 10:13:29 AM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

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:

0dd4199, 3a45c2c

Parents:

f95556f

Message:

onion: now matches sasview

File:

• sasmodels/models/onion.c (modified) (2 diffs)

sasmodels/models/onion.c

@@ -2 +2 @@
 static double
 f_exp(double q, double r, double sld_in, double sld_out,
-    double thickness, double A)
-{
-  const double vol = M_4PI_3 * cube(r);
-  const double qr = q * r;
-  const double alpha = A * r/thickness;
-  const double bes = sph_j1c(qr);
-  const double B = (sld_out - sld_in)/expm1(A);
-  const double C = sld_in - B;
-  double fun;
-  if (qr == 0.0) {
-    fun = 1.0;
-  } else if (fabs(A) > 0.0) {
-    const double qrsq = qr*qr;
-    const double alphasq = alpha*alpha;
-    const double sumsq = alphasq + qrsq;
-    double sinqr, cosqr;
-    SINCOS(qr, sinqr, cosqr);
-    fun = -3.0*(
-            ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq
-                - (alpha*sinqr/qr - cosqr)
-        ) / sumsq;
-  } else {
-    fun = bes;
-  }
-  return vol * (B*fun + C*bes);
-}
-
-static double
-f_linear(double q, double r, double sld, double slope)
+    double thickness, double A, double side)
 {
   const double vol = M_4PI_3 * cube(r);
   const double qr = q * r;
   const double bes = sph_j1c(qr);
-  double fun = 0.0;
-  if (qr > 0.0) {
-    const double qrsq = qr*qr;
+  const double alpha = A * r/thickness;
+  double result;
+  if (qr == 0.0) {
+    result = 1.0;
+  } else if (fabs(A) > 0.0) {
+    const double qrsq = qr * qr;
+    const double alphasq = alpha * alpha;
+    const double sumsq = alphasq + qrsq;
     double sinqr, cosqr;
     SINCOS(qr, sinqr, cosqr);
-    // Jae-He's code seems to simplify to this
-    //     fun = 3.0 * slope * r * (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq);
-    // Rederiving the math, we get the following instead:
-    fun = 3.0 * slope * r * (2.0*cosqr + qr*sinqr)/(qrsq*qrsq);
+    const double t1 = (alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr;
+    const double t2 = alpha*sinqr/qr - cosqr;
+    const double fun = -3.0*(t1/sumsq - t2)/sumsq;
+    const double slope = (sld_out - sld_in)/expm1(A);
+    const double contrast = slope*exp(A*side);
+    result = contrast*fun + (sld_in-slope)*bes;
+  } else {
+    result = sld_in*bes;
   }
-  return vol * (sld*bes + fun);
-}
-
-static double
-f_constant(double q, double r, double sld)
-{
-  const double bes = sph_j1c(q * r);
-  const double vol = M_4PI_3 * cube(r);
-  return sld * vol * bes;
+  return vol * result;
 }
 
@@ -69 +42 @@
 static double
 Iq(double q, double sld_core, double core_radius, double sld_solvent,
-    double n, double sld_in[], double sld_out[], double thickness[],
+    double n_shells, double sld_in[], double sld_out[], double thickness[],
     double A[])
 {
-  int i;
-  double r = core_radius;
-  double f = f_constant(q, r, sld_core);
-  for (i=0; i<n; i++){
-    const double r0 = r;
-    r += thickness[i];
-    if (r == r0) {
-      // no thickness, so nothing to add
-    } else if (fabs(A[i]) < 1.0e-16 || sld_out[i] == sld_in[i]) {
-      f -= f_constant(q, r0, sld_in[i]);
-      f += f_constant(q, r, sld_in[i]);
-    } else if (fabs(A[i]) < 1.0e-4) {
-      const double slope = (sld_out[i] - sld_in[i])/thickness[i];
-      f -= f_linear(q, r0, sld_in[i], slope);
-      f += f_linear(q, r, sld_out[i], slope);
-    } else {
-      f -= f_exp(q, r0, sld_in[i], sld_out[i], thickness[i], A[i]);
-      f += f_exp(q, r, sld_in[i], sld_out[i], thickness[i], A[i]);
-    }
+  int n = (int)(n_shells+0.5);
+  double r_out = core_radius;
+  double f = f_exp(q, r_out, sld_core, 0.0, 0.0, 0.0, 0.0);
+  for (int i=0; i < n; i++){
+    const double r_in = r_out;
+    r_out += thickness[i];
+    f -= f_exp(q, r_in, sld_in[i], sld_out[i], thickness[i], A[i], 0.0);
+    f += f_exp(q, r_out, sld_in[i], sld_out[i], thickness[i], A[i], 1.0);
   }
-  f -= f_constant(q, r, sld_solvent);
+  f -= f_exp(q, r_out, sld_solvent, 0.0, 0.0, 0.0, 0.0);
   const double f2 = f * f * 1.0e-4;