source: sasmodels/sasmodels/models/onion.c @ 1e7b0db0

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 1e7b0db0 was 9762341, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

onion: update variables to match parameter names

  • Property mode set to 100644
File size: 1.7 KB
Line 
1
2static double
3f_exp(double q, double r, double sld_in, double sld_out,
4    double thickness, double A, double side)
5{
6  const double vol = M_4PI_3 * cube(r);
7  const double qr = q * r;
8  const double bes = sph_j1c(qr);
9  const double alpha = A * r/thickness;
10  double result;
11  if (qr == 0.0) {
12    result = 1.0;
13  } else if (fabs(A) > 0.0) {
14    const double qrsq = qr * qr;
15    const double alphasq = alpha * alpha;
16    const double sumsq = alphasq + qrsq;
17    double sinqr, cosqr;
18    SINCOS(qr, sinqr, cosqr);
19    const double t1 = (alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr;
20    const double t2 = alpha*sinqr/qr - cosqr;
21    const double fun = -3.0*(t1/sumsq - t2)/sumsq;
22    const double slope = (sld_out - sld_in)/expm1(A);
23    const double contrast = slope*exp(A*side);
24    result = contrast*fun + (sld_in-slope)*bes;
25  } else {
26    result = sld_in*bes;
27  }
28  return vol * result;
29}
30
31static double
32form_volume(double radius_core, double n, double thickness[])
33{
34  int i;
35  double r = radius_core;
36  for (i=0; i < n; i++) {
37    r += thickness[i];
38  }
39  return M_4PI_3*cube(r);
40}
41
42static double
43Iq(double q, double sld_core, double radius_core, double sld_solvent,
44    double n_shells, double sld_in[], double sld_out[], double thickness[],
45    double A[])
46{
47  int n = (int)(n_shells+0.5);
48  double r_out = radius_core;
49  double f = f_exp(q, r_out, sld_core, 0.0, 0.0, 0.0, 0.0);
50  for (int i=0; i < n; i++){
51    const double r_in = r_out;
52    r_out += thickness[i];
53    f -= f_exp(q, r_in, sld_in[i], sld_out[i], thickness[i], A[i], 0.0);
54    f += f_exp(q, r_out, sld_in[i], sld_out[i], thickness[i], A[i], 1.0);
55  }
56  f -= f_exp(q, r_out, sld_solvent, 0.0, 0.0, 0.0, 0.0);
57  const double f2 = f * f * 1.0e-4;
58
59  return f2;
60}
Note: See TracBrowser for help on using the repository browser.