source: sasmodels/sasmodels/models/pearl_necklace.c @ 2c74c11

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

implicit Iqxy; fix divide by 0 for q=0
  • Property mode set to 100644
File size: 3.9 KB
Line 
1double _pearl_necklace_kernel(double q, double radius, double edge_separation,
2        double thick_string, double num_pearls, double sld_pearl,
3        double sld_string, double sld_solv);
4double form_volume(double radius, double edge_separation,
5        double string_thickness, double number_of_pearls);
6
7double Iq(double q, double radius, double edge_separation,
8        double string_thickness, double number_of_pearls, double sld, 
9        double string_sld, double solvent_sld);
10
11// From Igor library
12double _pearl_necklace_kernel(double q, double radius, double edge_separation, double thick_string,
13        double num_pearls, double sld_pearl, double sld_string, double sld_solv)
14{
15        //relative slds
16        double contrast_pearl = sld_pearl - sld_solv;
17        double contrast_string = sld_string - sld_solv;
18       
19        // number of string segments
20        num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls
21        double num_strings = num_pearls - 1.0;
22       
23        //Pi
24        double pi = 4.0*atan(1.0);
25       
26        // center to center distance between the neighboring pearls
27        double A_s = edge_separation + 2.0 * radius;
28       
29        // Repeated Calculations
30        double sincasq = sinc(q*A_s);
31        double oneminussinc = 1 - sincasq;
32        double q_r = q * radius;
33        double q_edge = q * edge_separation;
34       
35        // each volume
36        double string_vol = edge_separation * pi * thick_string * thick_string / 4.0;
37        double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius;
38
39        //total volume
40        double tot_vol;
41        //each masses
42        double m_r= contrast_string * string_vol;
43        double m_s= contrast_pearl * pearl_vol;
44        double psi, gamma, beta;
45        //form factors
46        double sss, srr, srs; //cross
47        double srr_1, srr_2, srr_3;
48        double form_factor;
49        tot_vol = num_strings * string_vol;
50        tot_vol += num_pearls * pearl_vol;
51
52        //sine functions of a pearl
53        psi = sin(q_r);
54        psi -= q_r * cos(q_r);
55        psi *= 3.0;
56        psi /= q_r * q_r * q_r;
57
58        // Note take only 20 terms in Si series: 10 terms may be enough though.
59        gamma = Si(q_edge);
60        gamma /= (q_edge);
61        beta = Si(q * (A_s - radius));
62        beta -= Si(q_r);
63        beta /= q_edge;
64
65        // form factor for num_pearls
66        sss = 1.0 - pow(sincasq, num_pearls);
67        sss /= oneminussinc * oneminussinc;
68        sss *= -sincasq;
69        sss -= num_pearls / 2.0;
70        sss += num_pearls / oneminussinc;
71        sss *= 2.0 * m_s * psi * m_s * psi;
72
73        // form factor for num_strings (like thin rods)
74        srr_1 = -sinc(q_edge/2.0) * sinc(q_edge/2.0);
75
76        srr_1 += 2.0 * gamma;
77        srr_1 *= num_strings;
78        srr_2 = 2.0/oneminussinc;
79        srr_2 *= num_strings;
80        srr_2 *= beta * beta;
81        srr_3 = 1.0 - pow(sincasq, num_strings);
82        srr_3 /= oneminussinc * oneminussinc;
83        srr_3 *= beta * beta;
84        srr_3 *= -2.0;
85
86        // total srr
87        srr = srr_1 + srr_2 + srr_3;
88        srr *= m_r * m_r;
89
90        // form factor for correlations
91        srs = 1.0;
92        srs -= pow(sincasq, num_strings);
93        srs /= oneminussinc * oneminussinc;
94        srs *= -sincasq;
95        srs += num_strings/oneminussinc;
96        srs *= 4.0;
97        srs *= (m_r * m_s * beta * psi);
98
99        form_factor = sss + srr + srs;
100        form_factor /= (tot_vol * 1.0e4); // norm by volume and A^-1 to cm^-1
101
102        return (form_factor);
103}
104
105double form_volume(double radius, double edge_separation,
106        double string_thickness, double number_of_pearls)
107{
108        double total_vol;
109
110        double pi = 4.0*atan(1.0);
111        double number_of_strings = number_of_pearls - 1.0;
112       
113        double string_vol = edge_separation * pi * string_thickness * string_thickness / 4.0;
114        double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius;
115
116        total_vol = number_of_strings * string_vol;
117        total_vol += number_of_pearls * pearl_vol;
118
119        return(total_vol);
120}
121
122double Iq(double q, double radius, double edge_separation,
123        double string_thickness, double number_of_pearls, double sld, 
124        double string_sld, double solvent_sld)
125{
126        double value, tot_vol;
127       
128        if (string_thickness >= radius || number_of_pearls <= 0) {
129                return NAN;
130        }
131       
132        value = _pearl_necklace_kernel(q, radius, edge_separation, string_thickness,
133                number_of_pearls, sld, string_sld, solvent_sld);
134        tot_vol = form_volume(radius, edge_separation, string_thickness, number_of_pearls);
135
136        return value*tot_vol;
137}
Note: See TracBrowser for help on using the repository browser.