source: sasmodels/sasmodels/models/pringle.c @ 46ed760

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

pringles → pringle name change
  • Property mode set to 100644
File size: 4.8 KB
Line 
1double form_volume(double radius,
2          double thickness);
3
4double Iq(double q,
5          double radius,
6          double thickness,
7          double alpha,
8          double beta,
9          double sld_pringle,
10          double sld_solvent);
11
12double Iqxy(double qx, double qy,
13          double radius,
14          double thickness,
15          double alpha,
16          double beta,
17          double sld_pringle,
18          double sld_solvent);
19
20static
21double pringleC(double radius,
22                double alpha,
23                double beta,
24                double q,
25                double phi,
26                double n) {
27
28    double va, vb;
29    double bessargs, cosarg, bessargcb;
30    double r, retval, yyy;
31
32
33    va = 0;
34    vb = radius;
35
36    // evaluate at Gauss points
37    // remember to index from 0,size-1
38
39    double summ = 0.0;          // initialize integral
40    int ii = 0;
41    do {
42        // Using 76 Gauss points
43        r = (Gauss76Z[ii] * (vb - va) + vb + va) / 2.0;
44
45        bessargs = q*r*sin(phi);
46        cosarg = q*r*r*alpha*cos(phi);
47        bessargcb = q*r*r*beta*cos(phi);
48
49        yyy = Gauss76Wt[ii]*r*cos(cosarg)
50                *sas_JN(n, bessargcb)
51                *sas_JN(2*n, bessargs);
52        summ += yyy;
53
54        ii += 1;
55    } while (ii < N_POINTS_76);                 // end of loop over quadrature points
56    //
57    // calculate value of integral to return
58
59    retval = (vb - va) / 2.0 * summ;
60    retval = retval / pow(r, 2.0);
61
62    return retval;
63}
64
65static
66double pringleS(double radius,
67                double alpha,
68                double beta,
69                double q,
70                double phi,
71                double n) {
72
73    double va, vb, summ;
74    double bessargs, sinarg, bessargcb;
75    double r, retval, yyy;
76    // set up the integration
77    // end points and weights
78
79    va = 0;
80    vb = radius;
81
82    // evaluate at Gauss points
83    // remember to index from 0,size-1
84
85    summ = 0.0;         // initialize integral
86    int ii = 0;
87    do {
88        // Using 76 Gauss points
89        r = (Gauss76Z[ii] * (vb - va) + vb + va) / 2.0;
90
91        bessargs = q*r*sin(phi);
92        sinarg = q*r*r*alpha*cos(phi);
93        bessargcb = q*r*r*beta*cos(phi);
94
95        yyy = Gauss76Wt[ii]*r*sin(sinarg)
96                    *sas_JN(n, bessargcb)
97                    *sas_JN(2*n, bessargs);
98        summ += yyy;
99
100        ii += 1;
101    } while (ii < N_POINTS_76);
102
103    // end of loop over quadrature points
104    //
105    // calculate value of integral to return
106
107    retval = (vb-va)/2.0*summ;
108    retval = retval/pow(r, 2.0);
109
110    return retval;
111}
112
113static
114double _kernel(double thickness,
115               double radius,
116               double alpha,
117               double beta,
118               double q,
119               double phi) {
120
121    const double sincarg = q * thickness * cos(phi) / 2.0;
122    const double sincterm = pow(sin(sincarg) / sincarg, 2.0);
123
124    //calculate sum term from n = -3 to 3
125    double sumterm = 0.0;
126    for (int nn = -3; nn <= 3; nn++) {
127        double powc = pringleC(radius, alpha, beta, q, phi, nn);
128        double pows = pringleS(radius, alpha, beta, q, phi, nn);
129        sumterm += pow(powc, 2.0) + pow(pows, 2.0);
130    }
131    double retval = 4.0 * sin(phi) * sumterm * sincterm;
132
133    return retval;
134
135}
136
137static double pringles_kernel(double q,
138          double radius,
139          double thickness,
140          double alpha,
141          double beta,
142          double sld_pringle,
143          double sld_solvent)
144{
145
146    //upper and lower integration limits
147    const double lolim = 0.0;
148    const double uplim = M_PI / 2.0;
149
150    double summ = 0.0;                  //initialize integral
151
152    double delrho = sld_pringle - sld_solvent; //make contrast term
153
154    for (int i = 0; i < N_POINTS_76; i++) {
155        double phi = (Gauss76Z[i] * (uplim - lolim) + uplim + lolim) / 2.0;
156        summ += Gauss76Wt[i] * _kernel(thickness, radius, alpha, beta, q, phi);
157    }
158
159    double answer = (uplim - lolim) / 2.0 * summ;
160    answer *= delrho*delrho;
161
162    return answer;
163}
164
165double form_volume(double radius,
166        double thickness){
167
168        return 1.0;
169}
170
171double Iq(double q,
172          double radius,
173          double thickness,
174          double alpha,
175          double beta,
176          double sld_pringle,
177          double sld_solvent)
178{
179    const double form = pringles_kernel(q,
180                  radius,
181                  thickness,
182                  alpha,
183                  beta,
184                  sld_pringle,
185                  sld_solvent);
186
187    return 1.0e-4*form*M_PI*radius*radius*thickness;
188}
189
190double Iqxy(double qx, double qy,
191            double radius,
192            double thickness,
193            double alpha,
194            double beta,
195            double sld_pringle,
196            double sld_solvent)
197{
198    double q = sqrt(qx*qx + qy*qy);
199    return Iq(q,
200            radius,
201            thickness,
202            alpha,
203            beta,
204            sld_pringle,
205            sld_solvent);
206}
207
Note: See TracBrowser for help on using the repository browser.