source: sasmodels/sasmodels/models/onion.c @ fdb1487

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since fdb1487 was fdb1487, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

add non-working onion model, named _onion so it is not in the model list

  • Property mode set to 100644
File size: 2.6 KB
Line 
1
2static double
3f_exp(double q, double r, double sld_in, double sld_out,
4    double thickness, double A)
5{
6  const double vol = 4.0/3.0 * M_PI * r * r * r;
7  const double qr = q * r;
8  const double alpha = A * r/thickness;
9  const double bes = sph_j1c(qr);
10  const double B = (sld_out - sld_in)/expm1(A);
11  const double C = sld_in - B;
12  double fun;
13  if (qr == 0.0) {
14    fun = 1.0;
15  } else if (fabs(A) > 0.0) {
16    const double qrsq = qr*qr;
17    const double alphasq = alpha*alpha;
18    const double sumsq = alphasq + qrsq;
19    double sinqr, cosqr;
20    SINCOS(qr, sinqr, cosqr);
21    fun = -3.0(
22            ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq
23                - (alpha*sinqr/qr - cosqr)
24        ) / sumsq;
25  } else {
26    fun = bes;
27  }
28  return vol * (B*fun + C*bes);
29}
30
31static double
32f_linear(double q, double r, double sld, double slope)
33{
34  const double vol = 4.0/3.0 * M_PI * r * r * r;
35  const double qr = q * r;
36  const double bes = sph_j1c(qr);
37  double fun = 0.0;
38  if (qr > 0.0) {
39    const double qrsq = qr*qr;
40    double sinqr, cosqr;
41    SINCOS(qr, sinqr, cosqr);
42    // Jae-He's code seems to simplify to this
43    //     fun = 3.0 * slope * r * (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq);
44    // Rederiving the math, we get the following instead:
45    fun = 3.0 * slope * r * (2.0*cosqr + qr*sinqr)/(qrsq*qrsq);
46  }
47  return vol * (sld*bes + fun);
48}
49
50static double
51f_constant(double q, double r, double sld)
52{
53  const double bes = sph_j1c(q * r);
54  const double vol = 4.0/3.0 * M_PI * r * r * r;
55  return sld * vol * bes;
56}
57
58static double
59form_volume(double core_radius, double n, double thickness[])
60{
61  int i;
62  double r = core_radius;
63  for (i=0; i < n; i++) {
64    r += thickness[i];
65  }
66  return 4.0/3.0 * M_PI * r * r * r;
67}
68
69static double
70Iq(double q, double core_sld, double core_radius, double solvent_sld,
71    double n, double in_sld[], double out_sld[], double thickness[],
72    double A[])
73{
74  int i;
75  r = core_radius;
76  f = f_constant(q, r, core_sld);
77  for (i=0; i<n; i++){
78    const double r0 = r;
79    r += thickness[i];
80    if (r == r0) {
81      // no thickness, so nothing to add
82    } else if (fabs(A[i]) < 1e-16 || sld_out[i] == sld_in[i]) {
83      f -= f_constant(q, r0, sld_in[i]);
84      f += f_constant(q, r, sld_in[i]);
85    } else if (fabs(A[i]) < 1e-4) {
86      const double slope = (sld_out[i] - sld_in[i])/thickness[i];
87      f -= f_linear(q, r0, sld_in[i], slope);
88      f += f_linear(q, r, sld_out[i], slope);
89    } else {
90      f -= f_exp(q, r0, sld_in[i], sld_out[i], thickness[i], A[i]);
91      f += f_exp(q, r, sld_in[i], sld_out[i], thickness[i], A[i]);
92    }
93  }
94  f -= f_constant(q, r, solvent_sld);
95  f2 = f * f * 1.0e-4;
96
97  return f2;
98}
Note: See TracBrowser for help on using the repository browser.