source: sasmodels/sasmodels/models/stacked_disks.c @ 0d6e865

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0d6e865 was 0d6e865, checked in by dirk, 8 years ago

Rewriting selected models in spherical coordinates. This fixes ticket #492 for these models.
  • Property mode set to 100644
File size: 5.7 KB
RevLine 
[a807206]1double form_volume(double thick_core,
2                   double thick_layer,
[66d119f]3                   double radius,
4                   double n_stacking);
5
6double Iq(double q,
[a807206]7          double thick_core,
8          double thick_layer,
[66d119f]9          double radius,
10          double n_stacking,
[a807206]11          double sigma_dnn,
[66d119f]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,
[a807206]23               double thick_layer,
[66d119f]24               double zi,
[a807206]25               double sigma_dnn,
[66d119f]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);
[a807206]40        const double sinarg2 = qq*(halfheight+thick_layer)*cos(zi);
[66d119f]41
[2c74c11]42        const double be1 = sas_J1c(besarg1);
[43b7eea]43        const double be2 = sas_J1c(besarg2);
[66d119f]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;
[a807206]50        const double totald=2.0*(thick_layer+halfheight);
[66d119f]51
[43b7eea]52        const double t1 = area*(2.0*halfheight)*dr1*(si1)*(be1);
53        const double t2 = area*dr2*(totald*si2-2.0*halfheight*si1)*(be2);
[66d119f]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) {
[a807206]61                double dexpt=qq*cos(zi)*qq*cos(zi)*d*d*sigma_dnn*sigma_dnn*kk/2.0;
[66d119f]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,
[a807206]76                            double thick_core,
77                            double thick_layer,
[66d119f]78                            double radius,
79                            double n_stacking,
[a807206]80                            double sigma_dnn,
[66d119f]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
[a807206]90        double d=2.0*thick_layer+thick_core;
91        double halfheight = thick_core/2.0;
[66d119f]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,
[a807206]102                                   thick_layer,
[66d119f]103                                   zi,
[a807206]104                                   sigma_dnn,
[66d119f]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,
[a807206]119                            double thick_core,
120                            double thick_layer,
[66d119f]121                            double radius,
122                            double n_stacking,
[a807206]123                            double sigma_dnn,
[66d119f]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;
[d507c3a]132
133    //convert angle degree to radian
134    theta = theta * M_PI/180.0;
135    phi = phi * M_PI/180.0;
136
[66d119f]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
[0d6e865]144    const double cyl_x = st * cp;
145    const double cyl_y = st * sp;
[66d119f]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
[a807206]156    double d = 2 * thick_layer + thick_core;
157    double halfheight = thick_core/2.0;
[66d119f]158    double answer = _kernel(q,
159                     radius,
160                     core_sld,
161                     layer_sld,
162                     solvent_sld,
163                     halfheight,
[a807206]164                     thick_layer,
[66d119f]165                     alpha,
[a807206]166                     sigma_dnn,
[66d119f]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
[a807206]177double form_volume(double thick_core,
178                   double thick_layer,
[66d119f]179                   double radius,
180                   double n_stacking){
[a807206]181    double d = 2 * thick_layer + thick_core;
[66d119f]182    return acos(-1.0) * radius * radius * d * n_stacking;
183}
184
185double Iq(double q,
[a807206]186          double thick_core,
187          double thick_layer,
[66d119f]188          double radius,
189          double n_stacking,
[a807206]190          double sigma_dnn,
[66d119f]191          double core_sld,
192          double layer_sld,
193          double solvent_sld)
194{
195    return stacked_disks_kernel(q,
[a807206]196                    thick_core,
197                    thick_layer,
[66d119f]198                    radius,
199                    n_stacking,
[a807206]200                    sigma_dnn,
[66d119f]201                    core_sld,
202                    layer_sld,
203                    solvent_sld);
204}
Note: See TracBrowser for help on using the repository browser.