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

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

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