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

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

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