source: sasmodels/sasmodels/models/stacked_disks.c @ 797a8e3

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 797a8e3 was 3ac4e1b, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

stacked_disks: code tidying
  • Property mode set to 100644
File size: 5.6 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_J1c(besarg1);
56    //const double be2 = sas_J1c(besarg2);
57    const double be2 = be1;
58    const double si1 = sinc(sinarg1);
59    const double si2 = sinc(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    double debye_arg = -0.5*square(qd_cos_alpha*sigma_dnn);
74    double sq=0.0;
75    for (int kk=1; kk<n_stacking; kk++) {
76        sq += (n_stacking-kk) * cos(qd_cos_alpha*kk) * exp(debye_arg*kk);
77    }
78    // end of loop for S(q)
79    sq = 1.0 + 2.0*sq/n_stacking;
80
81    return pq * sq;
82}
83
84
85static
86double stacked_disks_kernel(double q,
87                            double thick_core,
88                            double thick_layer,
89                            double radius,
90                            double n_stacking,
91                            double sigma_dnn,
92                            double core_sld,
93                            double layer_sld,
94                            double solvent_sld)
95{
96/*    StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks
97like clay platelets that are not exfoliated
98*/
99    double summ = 0.0;    //initialize integral
100
101    double d = 2.0*thick_layer+thick_core;
102    double halfheight = 0.5*thick_core;
103
104    for(int i=0; i<N_POINTS_76; i++) {
105        double zi = (Gauss76Z[i] + 1.0)*M_PI_4;
106        double sin_alpha, cos_alpha; // slots to hold sincos function output
107        SINCOS(zi, sin_alpha, cos_alpha);
108        double yyy = _kernel(q,
109                           radius,
110                           core_sld,
111                           layer_sld,
112                           solvent_sld,
113                           halfheight,
114                           thick_layer,
115                           sin_alpha,
116                           cos_alpha,
117                           sigma_dnn,
118                           d,
119                           n_stacking);
120        summ += Gauss76Wt[i] * yyy * sin_alpha;
121    }
122
123    double answer = M_PI_4*summ;
124
125    //Convert to [cm-1]
126    return 1.0e-4*answer;
127}
128
129double form_volume(double thick_core,
130                   double thick_layer,
131                   double radius,
132                   double n_stacking){
133    double d = 2.0 * thick_layer + thick_core;
134    return M_PI * radius * radius * d * n_stacking;
135}
136
137double Iq(double q,
138          double thick_core,
139          double thick_layer,
140          double radius,
141          double n_stacking,
142          double sigma_dnn,
143          double core_sld,
144          double layer_sld,
145          double solvent_sld)
146{
147    return stacked_disks_kernel(q,
148                    thick_core,
149                    thick_layer,
150                    radius,
151                    n_stacking,
152                    sigma_dnn,
153                    core_sld,
154                    layer_sld,
155                    solvent_sld);
156}
157
158
159double
160Iqxy(double qx, double qy,
161     double thick_core,
162     double thick_layer,
163     double radius,
164     double n_stacking,
165     double sigma_dnn,
166     double core_sld,
167     double layer_sld,
168     double solvent_sld,
169     double theta,
170     double phi)
171{
172    double q, sin_alpha, cos_alpha;
173    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
174
175    double d = 2.0 * thick_layer + thick_core;
176    double halfheight = 0.5*thick_core;
177    double answer = _kernel(q,
178                     radius,
179                     core_sld,
180                     layer_sld,
181                     solvent_sld,
182                     halfheight,
183                     thick_layer,
184                     sin_alpha,
185                     cos_alpha,
186                     sigma_dnn,
187                     d,
188                     n_stacking);
189
190    //convert to [cm-1]
191    answer *= 1.0e-4;
192
193    return answer;
194}
195
Note: See TracBrowser for help on using the repository browser.