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

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since d507c3a was f12357f, checked in by krzywon, 8 years ago

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