source: sasmodels/sasmodels/models/stacked_disks.c @ 66d119f

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

Converted StackedDisks?
  • Property mode set to 100644
File size: 6.6 KB
Line 
1double form_volume(double core_thick,
2                   double layer_thick,
3                   double radius,
4                   double n_stacking);
5
6double Iq(double q,
7          double core_thick,
8          double layer_thick,
9          double radius,
10          double n_stacking,
11          double sigma_d,
12          double core_sld,
13          double layer_sld,
14          double solvent_sld);
15
16double Iqxy(double qx, double qy,
17            double core_thick,
18            double layer_thick,
19            double radius,
20            double n_stacking,
21            double sigma_d,
22            double core_sld,
23            double layer_sld,
24            double solvent_sld,
25            double theta,
26            double phi);
27
28static
29double _kernel(double qq,
30               double radius,
31               double core_sld,
32               double layer_sld,
33               double solvent_sld,
34               double halfheight,
35               double layer_thick,
36               double zi,
37               double sigma_d,
38               double d,
39               double n_stacking)
40
41{
42        // qq is the q-value for the calculation (1/A)
43        // radius is the core radius of the cylinder (A)
44        // *_sld are the respective SLD's
45        // halfheight is the *Half* CORE-LENGTH of the cylinder = L (A)
46        // zi is the dummy variable for the integration (x in Feigin's notation)
47
48        const double besarg1 = qq*radius*sin(zi);
49        const double besarg2 = qq*radius*sin(zi);
50
51        const double sinarg1 = qq*halfheight*cos(zi);
52        const double sinarg2 = qq*(halfheight+layer_thick)*cos(zi);
53
54    const double be1 = J1(besarg1)/besarg1;
55        const double be2 = J1(besarg2)/besarg2;
56        const double si1 = sin(sinarg1)/sinarg1;
57        const double si2 = sin(sinarg2)/sinarg2;
58
59        const double dr1 = (core_sld-solvent_sld);
60        const double dr2 = (layer_sld-solvent_sld);
61        const double area = M_PI*radius*radius;
62        const double totald=2.0*(layer_thick+halfheight);
63
64        const double t1 = 2.0*area*(2.0*halfheight)*dr1*(si1)*(be1);
65        const double t2 = 2.0*area*dr2*(totald*si2-2.0*halfheight*si1)*(be2);
66
67
68        double retval =((t1+t2)*(t1+t2))*sin(zi);
69
70        // loop for the structure facture S(q)
71        double sqq=0.0;
72        for(int kk=1;kk<n_stacking;kk+=1) {
73                double dexpt=qq*cos(zi)*qq*cos(zi)*d*d*sigma_d*sigma_d*kk/2.0;
74                sqq=sqq+(n_stacking-kk)*cos(qq*cos(zi)*d*kk)*exp(-1.*dexpt);
75        }
76
77        // end of loop for S(q)
78        sqq=1.0+2.0*sqq/n_stacking;
79
80        retval *= sqq;
81
82        return(retval);
83}
84
85
86static
87double stacked_disks_kernel(double q,
88                            double core_thick,
89                            double layer_thick,
90                            double radius,
91                            double n_stacking,
92                            double sigma_d,
93                            double core_sld,
94                            double layer_sld,
95                            double solvent_sld)
96{
97/*      StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks
98like clay platelets that are not exfoliated
99*/
100        double summ = 0.0;      //initialize integral
101
102        double d=2.0*layer_thick+core_thick;
103        double halfheight = core_thick/2.0;
104
105        for(int i=0;i<N_POINTS_76;i++) {
106                double zi = (Gauss76Z[i] + 1.0)*M_PI/4.0;
107                double yyy = Gauss76Wt[i] *
108                    _kernel(q,
109                                   radius,
110                                   core_sld,
111                                   layer_sld,
112                                   solvent_sld,
113                                   halfheight,
114                                   layer_thick,
115                                   zi,
116                                   sigma_d,
117                                   d,
118                                   n_stacking);
119                summ += yyy;
120        }
121
122        double answer = M_PI/4.0*summ;
123
124        //Convert to [cm-1]
125        answer *= 1.0e-4;
126
127        return answer;
128}
129
130static double stacked_disks_kernel_2d(double q, double q_x, double q_y,
131                            double core_thick,
132                            double layer_thick,
133                            double radius,
134                            double n_stacking,
135                            double sigma_d,
136                            double core_sld,
137                            double layer_sld,
138                            double solvent_sld,
139                            double theta,
140                            double phi)
141{
142
143    double ct, st, cp, sp;
144    SINCOS(theta, st, ct);
145    SINCOS(phi, sp, cp);
146
147    // silence compiler warnings about unused variable
148    (void) sp;
149
150    // parallelepiped orientation
151    const double cyl_x = ct * cp;
152    const double cyl_y = st;
153
154    // Compute the angle btw vector q and the
155    // axis of the parallelepiped
156    const double cos_val = cyl_x*q_x + cyl_y*q_y;
157
158    // Note: cos(alpha) = 0 and 1 will get an
159    // undefined value from Stackdisc_kern
160    double alpha = acos( cos_val );
161
162    // Call the IGOR library function to get the kernel
163    double d = 2 * layer_thick + core_thick;
164    double halfheight = core_thick/2.0;
165    double answer = _kernel(q,
166                     radius,
167                     core_sld,
168                     layer_sld,
169                     solvent_sld,
170                     halfheight,
171                     layer_thick,
172                     alpha,
173                     sigma_d,
174                     d,
175                     n_stacking);
176
177    answer /= sin(alpha);
178    //convert to [cm-1]
179    answer *= 1.0e-4;
180
181    return answer;
182}
183
184double form_volume(double core_thick,
185                   double layer_thick,
186                   double radius,
187                   double n_stacking){
188    double d = 2 * layer_thick + core_thick;
189    return acos(-1.0) * radius * radius * d * n_stacking;
190}
191
192double Iq(double q,
193          double core_thick,
194          double layer_thick,
195          double radius,
196          double n_stacking,
197          double sigma_d,
198          double core_sld,
199          double layer_sld,
200          double solvent_sld)
201{
202    return stacked_disks_kernel(q,
203                    core_thick,
204                    layer_thick,
205                    radius,
206                    n_stacking,
207                    sigma_d,
208                    core_sld,
209                    layer_sld,
210                    solvent_sld);
211}
212
213// Iqxy is never called since no orientation or magnetic parameters.
214double Iqxy(double qx, double qy,
215            double core_thick,
216            double layer_thick,
217            double radius,
218            double n_stacking,
219            double sigma_d,
220            double core_sld,
221            double layer_sld,
222            double solvent_sld,
223            double theta,
224            double phi)
225{
226    double q = sqrt(qx*qx + qy*qy);
227    return stacked_disks_kernel_2d(q, qx/q, qy/q,
228                    core_thick,
229                    layer_thick,
230                    radius,
231                    n_stacking,
232                    sigma_d,
233                    core_sld,
234                    layer_sld,
235                    solvent_sld,
236                    theta,
237                    phi);
238}
239
Note: See TracBrowser for help on using the repository browser.