source: sasmodels/sasmodels/models/spherical_sld.c @ 4a82d4d

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

Spherical SLD working expect low q region

  • Property mode set to 100644
File size: 12.6 KB
Line 
1double form_volume(double radius);
2
3double Iq(double q,
4    int n_shells, double thick_inter_0, int func_inter_0, double sld_core_0, double sld_solv,
5    double sld_flat_1, double sld_flat_2, double sld_flat_3, double sld_flat_4, double sld_flat_5,
6    double sld_flat_6, double sld_flat_7, double sld_flat_8, double sld_flat_9, double sld_flat_10,
7    double thick_inter_1, double thick_inter_2, double thick_inter_3, double thick_inter_4, double thick_inter_5,
8    double thick_inter_6, double thick_inter_7, double thick_inter_8, double thick_inter_9, double thick_inter_10,
9    double thick_flat_1, double thick_flat_2, double thick_flat_3, double thick_flat_4, double thick_flat_5,
10    double thick_flat_6, double thick_flat_7, double thick_flat_8, double thick_flat_9, double thick_flat_10,
11    int func_inter_1, int func_inter_2, int func_inter_3, int func_inter_4, int func_inter_5,
12    int func_inter_6, int func_inter_7, int func_inter_8, int func_inter_9, int func_inter_10,
13    double nu_inter_1, double nu_inter_2,double nu_inter_3, double nu_inter_4, double nu_inter_5,
14    double nu_inter_6, double nu_inter_7, double nu_inter_8, double nu_inter_9, double nu_inter_10,
15    int npts_inter, double nu_inter_0, double rad_core_0);
16
17double Iqxy(double qx, double qy,
18    int n_shells, double thick_inter_0, int func_inter_0, double sld_core_0, double sld_solv,
19    double sld_flat_1, double sld_flat_2, double sld_flat_3, double sld_flat_4, double sld_flat_5,
20    double sld_flat_6, double sld_flat_7, double sld_flat_8, double sld_flat_9, double sld_flat_10,
21    double thick_inter_1, double thick_inter_2, double thick_inter_3, double thick_inter_4, double thick_inter_5,
22    double thick_inter_6, double thick_inter_7, double thick_inter_8, double thick_inter_9, double thick_inter_10,
23    double thick_flat_1, double thick_flat_2, double thick_flat_3, double thick_flat_4, double thick_flat_5,
24    double thick_flat_6, double thick_flat_7, double thick_flat_8, double thick_flat_9, double thick_flat_10,
25    int func_inter_1, int func_inter_2, int func_inter_3, int func_inter_4, int func_inter_5,
26    int func_inter_6, int func_inter_7, int func_inter_8, int func_inter_9, int func_inter_10,
27    double nu_inter_1, double nu_inter_2,double nu_inter_3, double nu_inter_4, double nu_inter_5,
28    double nu_inter_6, double nu_inter_7, double nu_inter_8, double nu_inter_9, double nu_inter_10,
29    int npts_inter, double nu_inter_0, double rad_core_0);
30
31//TODO: Check what is for volume for this model
32
33double form_volume(double radius)
34{
35    //This is normalized for each shell. It is spehere but for each shell
36    double vol = 4.0*M_PI/3.0*pow(radius,3);
37    return vol;
38}
39
40
41static double sphere_sld_kernel(double dp[], double q) {
42  int n = dp[0];
43  int i,j,k;
44
45  double scale = dp[1];
46  double thick_inter_core = dp[2];
47  double sld_core = dp[4];
48  double sld_solv = dp[5];
49  double background = dp[6];
50  double npts = dp[57]; //number of sub_layers in each interface
51  double nsl=npts;//21.0; //nsl = Num_sub_layer:  MUST ODD number in double //no other number works now
52  int n_s;
53
54  double sld_i,sld_f,dz,bes,fun,f,vol,qr,r,contr,f2;
55  double sign,slope=0.0;
56  double pi;
57
58  //int* fun_type;
59  //double* sld;
60  //double* thick_inter;
61  //double* thick;
62  //double* fun_coef;
63
64  double total_thick=0.0;
65
66  //fun_type = (int*)malloc((n+2)*sizeof(int));
67  //sld = (double*)malloc((n+2)*sizeof(double));
68  //thick_inter = (double*)malloc((n+2)*sizeof(double));
69  //thick = (double*)malloc((n+2)*sizeof(double));
70  //fun_coef = (double*)malloc((n+2)*sizeof(double));
71
72  //TODO: Solution to avoid mallocs but probablyu can be done better
73  int fun_type[12];
74  double sld[12];
75  double thick_inter[12];
76  double thick[12];
77  double fun_coef[12];
78
79  fun_type[0] = dp[3];
80  fun_coef[0] = fabs(dp[58]);
81  for (i =1; i<=n; i++){
82    sld[i] = dp[i+6];
83    thick_inter[i]= dp[i+16];
84    thick[i] = dp[i+26];
85    fun_type[i] = dp[i+36];
86    fun_coef[i] = fabs(dp[i+46]);
87    total_thick += thick[i] + thick_inter[i];
88  }
89  sld[0] = sld_core;
90  sld[n+1] = sld_solv;
91  thick[0] = dp[59];
92  thick[n+1] = total_thick/5.0;
93  thick_inter[0] = thick_inter_core;
94  thick_inter[n+1] = 0.0;
95  fun_coef[n+1] = 0.0;
96
97  pi = 4.0*atan(1.0);
98  f = 0.0;
99  r = 0.0;
100  vol = 0.0;
101  //vol_pre = 0.0;
102  //vol_sub = 0.0;
103  sld_f = sld_core;
104
105  //floor_nsl = floor(nsl/2.0);
106
107  dz = 0.0;
108  // iteration for # of shells + core + solvent
109  for (i=0;i<=n+1; i++){
110    //iteration for N sub-layers
111    //if (fabs(thick[i]) <= 1e-24){
112    //   continue;
113    //}
114    // iteration for flat and interface
115    for (j=0;j<2;j++){
116      // iteration for sub_shells in the interface
117      // starts from #1 sub-layer
118      for (n_s=1;n_s<=nsl; n_s++){
119        // for solvent, it doesn't have an interface
120        if (i==n+1 && j==1)
121          break;
122        // for flat layers
123        if (j==0){
124          dz = thick[i];
125          sld_i = sld[i];
126          slope = 0.0;
127        }
128        // for interfacial sub_shells
129        else{
130          dz = thick_inter[i]/nsl;
131          // find sld_i at the outer boundary of sub-layer #n_s
132          sld_i = intersldfunc(fun_type[i],nsl, n_s, fun_coef[i], sld[i], sld[i+1]);
133          // calculate slope
134          slope= (sld_i -sld_f)/dz;
135        }
136        contr = sld_f-slope*r;
137        // iteration for the left and right boundary of the shells(or sub_shells)
138        for (k=0; k<2; k++){
139          // At r=0, the contribution to I is zero, so skip it.
140          if ( i == 0 && j == 0 && k == 0){
141            continue;
142          }
143          // On the top of slovent there is no interface; skip it.
144          if (i == n+1 && k == 1){
145            continue;
146          }
147          // At the right side (outer) boundary
148          if ( k == 1){
149            sign = 1.0;
150            r += dz;
151          }
152          // At the left side (inner) boundary
153          else{
154            sign = -1.0;
155          }
156          qr = q * r;
157          fun = 0.0;
158          if(qr == 0.0){
159            // sigular point
160            bes = sign * 1.0;
161          }
162          else{
163            // for flat sub-layer
164            bes = sign *  3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
165            // with linear slope
166            if (fabs(slope) > 0.0 ){
167              fun = sign * 3.0 * r * (2.0*qr*sin(qr)-((qr*qr)-2.0)*cos(qr))/(qr * qr * qr * qr);
168            }
169          }
170          // update total volume
171          vol = 4.0 * pi / 3.0 * r * r * r;
172          // we won't do the following volume correction for now.
173          // substrate empty area of volume
174          //if (k == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && fun_type[i]==0){
175          //  vol_sub += (vol_pre - vol);
176          //}
177          f += vol * (bes * contr + fun * slope);
178        }
179        // remember this sld as sld_f
180        sld_f =sld_i;
181        // no sub-layer iteration (n_s loop) for the flat layer
182        if (j==0)
183          break;
184      }
185    }
186  }
187  //vol += vol_sub;
188  f2 = f * f / vol * 1.0e8;
189  f2 *= scale;
190  f2 += background;
191  //free(fun_type);
192  //free(sld);
193  //free(thick_inter);
194  //free(thick);
195  //free(fun_coef);
196
197  return (f2);
198}
199
200
201/**
202 * Function to evaluate 1D SphereSLD function
203 * @param q: q-value
204 * @return: function value
205 */
206double Iq(double q,
207    int n_shells, double thick_inter_0, int func_inter_0, double sld_core_0, double sld_solv,
208    double sld_flat_1, double sld_flat_2, double sld_flat_3, double sld_flat_4, double sld_flat_5,
209    double sld_flat_6, double sld_flat_7, double sld_flat_8, double sld_flat_9, double sld_flat_10,
210    double thick_inter_1, double thick_inter_2, double thick_inter_3, double thick_inter_4, double thick_inter_5,
211    double thick_inter_6, double thick_inter_7, double thick_inter_8, double thick_inter_9, double thick_inter_10,
212    double thick_flat_1, double thick_flat_2, double thick_flat_3, double thick_flat_4, double thick_flat_5,
213    double thick_flat_6, double thick_flat_7, double thick_flat_8, double thick_flat_9, double thick_flat_10,
214    int func_inter_1, int func_inter_2, int func_inter_3, int func_inter_4, int func_inter_5,
215    int func_inter_6, int func_inter_7, int func_inter_8, int func_inter_9, int func_inter_10,
216    double nu_inter_1, double nu_inter_2,double nu_inter_3, double nu_inter_4, double nu_inter_5,
217    double nu_inter_6, double nu_inter_7, double nu_inter_8, double nu_inter_9, double nu_inter_10,
218    int npts_inter, double nu_inter_0, double rad_core_0) {
219
220    //printf("Number of points %d\n",npts_inter);
221    double intensity;
222    //TODO: Remove this container at later stage. It is only kept to minimize stupid errors now
223    double dp[60];
224    dp[0] = n_shells;
225    //This is scale will also have to be removed at some stage
226    dp[1] = 1.0;
227    dp[2] = thick_inter_0;
228    dp[3] = func_inter_0;
229    dp[4] = sld_core_0;
230    dp[5] = sld_solv;
231    dp[6] = 0.0;
232
233    dp[7] = sld_flat_1;
234    dp[8] = sld_flat_2;
235    dp[9] = sld_flat_3;
236    dp[10] = sld_flat_4;
237    dp[11] = sld_flat_5;
238    dp[12] = sld_flat_6;
239    dp[13] = sld_flat_7;
240    dp[14] = sld_flat_8;
241    dp[15] = sld_flat_9;
242    dp[16] = sld_flat_10;
243
244    dp[17] = thick_inter_1;
245    dp[18] = thick_inter_2;
246    dp[19] = thick_inter_3;
247    dp[20] = thick_inter_4;
248    dp[21] = thick_inter_5;
249    dp[22] = thick_inter_6;
250    dp[23] = thick_inter_7;
251    dp[24] = thick_inter_8;
252    dp[25] = thick_inter_9;
253    dp[26] = thick_inter_10;
254
255    dp[27] = thick_flat_1;
256    dp[28] = thick_flat_2;
257    dp[29] = thick_flat_3;
258    dp[30] = thick_flat_4;
259    dp[31] = thick_flat_5;
260    dp[32] = thick_flat_6;
261    dp[33] = thick_flat_7;
262    dp[34] = thick_flat_8;
263    dp[35] = thick_flat_9;
264    dp[36] = thick_flat_10;
265
266    dp[37] = func_inter_1;
267    dp[38] = func_inter_2;
268    dp[39] = func_inter_3;
269    dp[40] = func_inter_4;
270    dp[41] = func_inter_5;
271    dp[42] = func_inter_6;
272    dp[43] = func_inter_7;
273    dp[44] = func_inter_8;
274    dp[45] = func_inter_9;
275    dp[46] = func_inter_10;
276
277    dp[47] = nu_inter_1;
278    dp[48] = nu_inter_2;
279    dp[49] = nu_inter_3;
280    dp[50] = nu_inter_4;
281    dp[51] = nu_inter_5;
282    dp[52] = nu_inter_6;
283    dp[53] = nu_inter_7;
284    dp[54] = nu_inter_8;
285    dp[55] = nu_inter_9;
286    dp[56] = nu_inter_10;
287
288    dp[57] = npts_inter;
289    dp[58] = nu_inter_0;
290    dp[59] = rad_core_0;
291
292    intensity = sphere_sld_kernel(dp,q);
293    //printf("%10d\n",intensity);
294    return intensity;
295}
296
297/**
298 * Function to evaluate 2D SphereSLD function
299 * @param q_x: value of Q along x
300 * @param q_y: value of Q along y
301 * @return: function value
302 */
303double Iqxy(double qx, double qy,
304    int n_shells, double thick_inter_0, int func_inter_0, double sld_core_0, double sld_solv,
305    double sld_flat_1, double sld_flat_2, double sld_flat_3, double sld_flat_4, double sld_flat_5,
306    double sld_flat_6, double sld_flat_7, double sld_flat_8, double sld_flat_9, double sld_flat_10,
307    double thick_inter_1, double thick_inter_2, double thick_inter_3, double thick_inter_4, double thick_inter_5,
308    double thick_inter_6, double thick_inter_7, double thick_inter_8, double thick_inter_9, double thick_inter_10,
309    double thick_flat_1, double thick_flat_2, double thick_flat_3, double thick_flat_4, double thick_flat_5,
310    double thick_flat_6, double thick_flat_7, double thick_flat_8, double thick_flat_9, double thick_flat_10,
311    int func_inter_1, int func_inter_2, int func_inter_3, int func_inter_4, int func_inter_5,
312    int func_inter_6, int func_inter_7, int func_inter_8, int func_inter_9, int func_inter_10,
313    double nu_inter_1, double nu_inter_2,double nu_inter_3, double nu_inter_4, double nu_inter_5,
314    double nu_inter_6, double nu_inter_7, double nu_inter_8, double nu_inter_9, double nu_inter_10,
315    int npts_inter, double nu_inter_0, double rad_core_0) {
316
317    double q = sqrt(qx*qx + qy*qy);
318    return Iq(q, n_shells, thick_inter_0, func_inter_0, sld_core_0, sld_solv,
319    sld_flat_1, sld_flat_2, sld_flat_3, sld_flat_4, sld_flat_5,
320    sld_flat_6, sld_flat_7, sld_flat_8, sld_flat_9, sld_flat_10,
321    thick_inter_1, thick_inter_2, thick_inter_3, thick_inter_4, thick_inter_5,
322    thick_inter_6, thick_inter_7, thick_inter_8, thick_inter_9, thick_inter_10,
323    thick_flat_1, thick_flat_2, thick_flat_3, thick_flat_4, thick_flat_5,
324    thick_flat_6, thick_flat_7, thick_flat_8, thick_flat_9, thick_flat_10,
325    func_inter_1, func_inter_2, func_inter_3, func_inter_4, func_inter_5,
326    func_inter_6, func_inter_7, func_inter_8, func_inter_9, func_inter_10,
327    nu_inter_1, nu_inter_2, nu_inter_3, nu_inter_4, nu_inter_5,
328    nu_inter_6, nu_inter_7, nu_inter_8, nu_inter_9, nu_inter_10,
329    npts_inter, nu_inter_0, rad_core_0);
330
331}
332
Note: See TracBrowser for help on using the repository browser.