source: sasmodels/sasmodels/models/pearl_necklace.c @ 9b5fd42

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 9b5fd42 was 9b5fd42, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

remove compiler warning from pearl necklace

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