source: sasmodels/sasmodels/models/pearl_necklace.c @ 2126131

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

pearl necklace: code cleanup; use sph_j1c; force num_pearls to integer in form_volume

  • Property mode set to 100644
File size: 3.3 KB
Line 
1double form_volume(double radius, double edge_sep,
2    double thick_string, double num_pearls);
3double Iq(double q, double radius, double edge_sep,
4    double thick_string, double num_pearls, double sld, 
5    double string_sld, double solvent_sld);
6
7#define INVALID(v) (v.thick_string >= v.radius || v.num_pearls <= 0)
8
9// From Igor library
10static double
11_pearl_necklace_kernel(double q, double radius, double edge_sep, double thick_string,
12    double num_pearls, double sld_pearl, double sld_string, double sld_solv)
13{
14    // number of string segments
15    num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls
16    const double num_strings = num_pearls - 1.0;
17
18    //each masses: contrast * volume
19    const double contrast_pearl = sld_pearl - sld_solv;
20    const double contrast_string = sld_string - sld_solv;
21    const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string;
22    const double pearl_vol = M_4PI_3 * radius * radius * radius;
23    const double m_string = contrast_string * string_vol;
24    const double m_pearl = contrast_pearl * pearl_vol;
25
26    // center to center distance between the neighboring pearls
27    const double A_s = edge_sep + 2.0 * radius;
28
29    //sine functions of a pearl
30    // Note: lim_(q->0) Si(q*a)/(q*b) = a/b
31    // So therefore:
32    //    beta = q==0. ? 1.0 : (Si(q*(A_s-radius)) - Si(q*radius))/q_edge;
33    //    gamma = q==0. ? 1.0 : Si(q_edge)/q_edge;
34    // But there is a 1/(1-sinc) term below which blows up so don't bother
35    const double q_edge = q * edge_sep;
36    const double beta = (Si(q*(A_s-radius)) - Si(q*radius)) / q_edge;
37    const double gamma = Si(q_edge) / q_edge;
38    const double psi = sph_j1c(q*radius);
39
40    // Precomputed sinc terms
41    const double si = sinc(q*A_s);
42    const double omsi = 1.0 - si;
43    const double pow_si = pow(si, num_pearls);
44
45    // form factor for num_pearls
46    const double sss = 2.0*square(m_pearl*psi) * (
47        - si * (1.0 - pow_si) / (omsi*omsi)
48        + num_pearls / omsi
49        - 0.5 * num_pearls
50        );
51
52    // form factor for num_strings (like thin rods)
53    const double srr = m_string * m_string * (
54        - 2.0 * (1.0 - pow_si/si)*beta*beta / (omsi*omsi)
55        + 2.0 * num_strings*beta*beta / omsi
56        + num_strings * (2.0*gamma - square(sinc(q_edge/2.0)))
57        );
58
59    // form factor for correlations
60    const double srs = 4.0 * m_string * m_pearl * beta * psi * (
61        - si * (1.0 - pow_si/si) / (omsi*omsi)
62        + num_strings / omsi
63        );
64
65    const double form = sss + srr + srs;
66
67    return 1.0e-4 * form;
68}
69
70double form_volume(double radius, double edge_sep,
71    double thick_string, double num_pearls)
72{
73    num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls
74
75    const double num_strings = num_pearls - 1.0;
76    const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string;
77    const double pearl_vol = M_4PI_3 * radius * radius * radius;
78    const double volume = num_strings*string_vol + num_pearls*pearl_vol;
79
80    return volume;
81}
82
83double Iq(double q, double radius, double edge_sep,
84    double thick_string, double num_pearls, double sld, 
85    double string_sld, double solvent_sld)
86{
87    const double form = _pearl_necklace_kernel(q, radius, edge_sep,
88        thick_string, num_pearls, sld, string_sld, solvent_sld);
89
90    return form;
91}
Note: See TracBrowser for help on using the repository browser.