source: sasmodels/sasmodels/models/stacked_disks.c @ 6c3e266

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 6c3e266 was 6c3e266, checked in by wojciech, 6 years ago

Code clean-up - removed commented line

  • Property mode set to 100644
File size: 5.9 KB
Line 
1double form_volume(double thick_core,
2                   double thick_layer,
3                   double radius,
4                   double n_stacking);
5
6double Iq(double q,
7          double thick_core,
8          double thick_layer,
9          double radius,
10          double n_stacking,
11          double sigma_dnn,
12          double core_sld,
13          double layer_sld,
14          double solvent_sld);
15
16double Iqxy(double qx, double qy,
17          double thick_core,
18          double thick_layer,
19          double radius,
20          double n_stacking,
21          double sigma_dnn,
22          double core_sld,
23          double layer_sld,
24          double solvent_sld,
25          double theta,
26          double phi);
27
28static
29double _kernel(double q,
30               double radius,
31               double core_sld,
32               double layer_sld,
33               double solvent_sld,
34               double halfheight,
35               double thick_layer,
36               double sin_alpha,
37               double cos_alpha,
38               double sigma_dnn,
39               double d,
40               double n_stacking)
41
42{
43    // q is the q-value for the calculation (1/A)
44    // radius is the core radius of the cylinder (A)
45    // *_sld are the respective SLD's
46    // halfheight is the *Half* CORE-LENGTH of the cylinder = L (A)
47    // zi is the dummy variable for the integration (x in Feigin's notation)
48
49    const double besarg1 = q*radius*sin_alpha;
50    //const double besarg2 = q*radius*sin_alpha;
51
52    const double sinarg1 = q*halfheight*cos_alpha;
53    const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha;
54
55    const double be1 = sas_2J1x_x(besarg1);
56    //const double be2 = sas_2J1x_x(besarg2);
57    const double be2 = be1;
58    const double si1 = sas_sinx_x(sinarg1);
59    const double si2 = sas_sinx_x(sinarg2);
60
61    const double dr1 = core_sld - solvent_sld;
62    const double dr2 = layer_sld - solvent_sld;
63    const double area = M_PI*radius*radius;
64    const double totald = 2.0*(thick_layer + halfheight);
65
66    const double t1 = area * (2.0*halfheight) * dr1 * si1 * be1;
67    const double t2 = area * dr2 * (totald*si2 - 2.0*halfheight*si1) * be2;
68
69    double pq = square(t1 + t2);
70
71    // loop for the structure factor S(q)
72    double qd_cos_alpha = q*d*cos_alpha;
73    //d*cos_alpha is the projection of d onto q (in other words the component
74    //of d that is parallel to q.
75    double debye_arg = -0.5*square(qd_cos_alpha*sigma_dnn);
76    double sq=0.0;
77    for (int kk=1; kk<n_stacking; kk++) {
78        sq += (n_stacking-kk) * cos(qd_cos_alpha*kk) * exp(debye_arg*kk);
79    }
80    // end of loop for S(q)
81    sq = 1.0 + 2.0*sq/n_stacking;
82
83    return pq * sq * n_stacking;
84    // volume normalization should be per disk not per stack but form_volume
85    // is per stack so correct here for now.  Could change form_volume but
86    // if one ever wants to use P*S we need the ER based on the total volume
87}
88
89
90static
91double stacked_disks_kernel(double q,
92                            double thick_core,
93                            double thick_layer,
94                            double radius,
95                            double n_stacking,
96                            double sigma_dnn,
97                            double core_sld,
98                            double layer_sld,
99                            double solvent_sld)
100{
101/*    StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks
102like clay platelets that are not exfoliated
103*/
104    double summ = 0.0;    //initialize integral
105
106    double d = 2.0*thick_layer+thick_core;
107    double halfheight = 0.5*thick_core;
108
109    for(int i=0; i<N_POINTS_76; i++) {
110        double zi = (Gauss76Z[i] + 1.0)*M_PI_4;
111        double sin_alpha, cos_alpha; // slots to hold sincos function output
112        SINCOS(zi, sin_alpha, cos_alpha);
113        double yyy = _kernel(q,
114                           radius,
115                           core_sld,
116                           layer_sld,
117                           solvent_sld,
118                           halfheight,
119                           thick_layer,
120                           sin_alpha,
121                           cos_alpha,
122                           sigma_dnn,
123                           d,
124                           n_stacking);
125        summ += Gauss76Wt[i] * yyy * sin_alpha;
126    }
127
128    double answer = M_PI_4*summ;
129
130    //Convert to [cm-1]
131    return 1.0e-4*answer;
132}
133
134double form_volume(double thick_core,
135                   double thick_layer,
136                   double radius,
137                   double n_stacking){
138    double d = 2.0 * thick_layer + thick_core;
139    return M_PI * radius * radius * d * n_stacking;
140}
141
142double Iq(double q,
143          double thick_core,
144          double thick_layer,
145          double radius,
146          double n_stacking,
147          double sigma_dnn,
148          double core_sld,
149          double layer_sld,
150          double solvent_sld)
151{
152    return stacked_disks_kernel(q,
153                    thick_core,
154                    thick_layer,
155                    radius,
156                    n_stacking,
157                    sigma_dnn,
158                    core_sld,
159                    layer_sld,
160                    solvent_sld);
161}
162
163
164double
165Iqxy(double qx, double qy,
166     double thick_core,
167     double thick_layer,
168     double radius,
169     double n_stacking,
170     double sigma_dnn,
171     double core_sld,
172     double layer_sld,
173     double solvent_sld,
174     double theta,
175     double phi)
176{
177    double q, sin_alpha, cos_alpha;
178    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
179
180    double d = 2.0 * thick_layer + thick_core;
181    double halfheight = 0.5*thick_core;
182    double answer = _kernel(q,
183                     radius,
184                     core_sld,
185                     layer_sld,
186                     solvent_sld,
187                     halfheight,
188                     thick_layer,
189                     sin_alpha,
190                     cos_alpha,
191                     sigma_dnn,
192                     d,
193                     n_stacking);
194
195    //convert to [cm-1]
196    answer *= 1.0e-4;
197
198    return answer;
199}
200
Note: See TracBrowser for help on using the repository browser.