source: sasmodels/sasmodels/models/stacked_disks.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 d507c3a, checked in by piotr, 8 years ago

Fixed deg→rad conversion

  • Property mode set to 100644
File size: 6.7 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
145    //convert angle degree to radian
146    theta = theta * M_PI/180.0;
147    phi = phi * M_PI/180.0;
148
149    SINCOS(theta, st, ct);
150    SINCOS(phi, sp, cp);
151
152    // silence compiler warnings about unused variable
153    (void) sp;
154
155    // parallelepiped orientation
156    const double cyl_x = ct * cp;
157    const double cyl_y = st;
158
159    // Compute the angle btw vector q and the
160    // axis of the parallelepiped
161    const double cos_val = cyl_x*q_x + cyl_y*q_y;
162
163    // Note: cos(alpha) = 0 and 1 will get an
164    // undefined value from Stackdisc_kern
165    double alpha = acos( cos_val );
166
167    // Call the IGOR library function to get the kernel
168    double d = 2 * layer_thick + core_thick;
169    double halfheight = core_thick/2.0;
170    double answer = _kernel(q,
171                     radius,
172                     core_sld,
173                     layer_sld,
174                     solvent_sld,
175                     halfheight,
176                     layer_thick,
177                     alpha,
178                     sigma_d,
179                     d,
180                     n_stacking);
181
182    answer /= sin(alpha);
183    //convert to [cm-1]
184    answer *= 1.0e-4;
185
186    return answer;
187}
188
189double form_volume(double core_thick,
190                   double layer_thick,
191                   double radius,
192                   double n_stacking){
193    double d = 2 * layer_thick + core_thick;
194    return acos(-1.0) * radius * radius * d * n_stacking;
195}
196
197double Iq(double q,
198          double core_thick,
199          double layer_thick,
200          double radius,
201          double n_stacking,
202          double sigma_d,
203          double core_sld,
204          double layer_sld,
205          double solvent_sld)
206{
207    return stacked_disks_kernel(q,
208                    core_thick,
209                    layer_thick,
210                    radius,
211                    n_stacking,
212                    sigma_d,
213                    core_sld,
214                    layer_sld,
215                    solvent_sld);
216}
217
218// Iqxy is never called since no orientation or magnetic parameters.
219double Iqxy(double qx, double qy,
220            double core_thick,
221            double layer_thick,
222            double radius,
223            double n_stacking,
224            double sigma_d,
225            double core_sld,
226            double layer_sld,
227            double solvent_sld,
228            double theta,
229            double phi)
230{
231    double q = sqrt(qx*qx + qy*qy);
232    return stacked_disks_kernel_2d(q, qx/q, qy/q,
233                    core_thick,
234                    layer_thick,
235                    radius,
236                    n_stacking,
237                    sigma_d,
238                    core_sld,
239                    layer_sld,
240                    solvent_sld,
241                    theta,
242                    phi);
243}
244
Note: See TracBrowser for help on using the repository browser.