source: sasmodels/sasmodels/models/stacked_disks.c @ b34fc77

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since b34fc77 was b34fc77, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

minor code tidying to make q transformation regular

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