source: sasview/sansmodels/src/sans/models/c_extensions/onion.c @ ef70686

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since ef70686 was 3be94e8, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Re #5 fixing problem with malloc
  • Property mode set to 100644
File size: 8.1 KB
Line 
1/**
2 * Scattering model for a onion
3 */
4
5#include <math.h>
6#include "onion.h"
7#include <stdio.h>
8#include <stdlib.h>
9// some details can be found in sld_cal.c
10double so_kernel(double dp[], double q) {
11        int n = dp[0];
12        double scale = dp[1];
13        double rad_core0 = dp[2];
14        double sld_core0 = dp[3];
15        double sld_solv = dp[4];
16        double background = dp[5];
17        int i,j;
18  double bes,fun,alpha,f,vol,vol_pre,vol_sub,qr,r,contr,f2;
19  double sign;
20  double pi;
21  double r0 = 0.0;
22
23  double *sld_out;
24  double *slope;
25  double *sld_in;
26  double *thick;
27  double *A;
28  int *fun_type;
29
30  sld_out = (double*)malloc((n+2)*sizeof(double));
31  slope = (double*)malloc((n+2)*sizeof(double));
32  sld_in = (double*)malloc((n+2)*sizeof(double));
33  thick = (double*)malloc((n+2)*sizeof(double));
34  A = (double*)malloc((n+2)*sizeof(double));
35  fun_type = (int*)malloc((n+2)*sizeof(int));
36
37        for (i =1; i<=n; i++){
38                sld_out[i] = dp[i+5];
39                sld_in[i] = dp[i+15];
40                A[i] = dp[i+25];
41                thick[i] = dp[i+35];
42                fun_type[i] = dp[i+45];
43        }
44        sld_out[0] = sld_core0;
45        sld_out[n+1] = sld_solv;
46        sld_in[0] = sld_core0;
47        sld_in[n+1] = sld_solv;
48        thick[0] = rad_core0;
49        thick[n+1] = 1e+10;
50        A[0] = 0.0;
51        A[n+1] = 0.0;
52        fun_type[0] = 0;
53        fun_type[n+1] = 0;
54
55
56    pi = 4.0*atan(1.0);
57    f = 0.0;
58    r = 0.0;
59    vol = 0.0;
60    vol_pre = 0.0;
61    vol_sub = 0.0;
62
63    for (i =0; i<= n+1; i++){
64        if (thick[i] == 0.0){
65                continue;
66        }
67        if (fun_type[i]== 0 ){
68                slope[i] = 0.0;
69                A[i] = 0.0;
70        }
71        vol_pre = vol;
72        switch(fun_type[i]){
73            case 2 :
74                                        r0 = r;
75                                        if (A[i] == 0.0){
76                                                slope[i] = 0.0;
77                                        }
78                                        else{
79                                                slope[i] = (sld_out[i]-sld_in[i])/(exp(A[i])-1.0);
80                                        }
81                    for (j=0; j<2; j++){
82                        if ( i == 0 && j == 0){
83                            continue;
84                            }
85                        if (i == n+1 && j == 1){
86                            continue;
87                            }
88                        if ( j == 1){
89                            sign = 1.0;
90                            r += thick[i];
91                            }
92                        else{
93                            sign = -1.0;
94                                }
95                        qr = q * r;
96                        alpha = A[i] * r/thick[i];
97                        fun = 0.0;
98                        if(qr == 0.0){
99                            fun = sign * 1.0;
100                            bes = sign * 1.0;
101                            }
102                        else{
103                            if (fabs(A[i]) > 0.0 ){
104                                fun = 3.0 * ((alpha*alpha - qr * qr) * sin(qr) - 2.0 * alpha * qr * cos(qr))/ ((alpha * alpha + qr * qr) * (alpha * alpha + qr * qr) * qr);
105                                fun = fun - 3.0 * (alpha * sin(qr) - qr * cos(qr)) / ((alpha * alpha + qr * qr) * qr);
106                                fun = - sign *fun;
107                                bes = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
108                                }
109                            else {
110                                fun = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
111                                bes = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
112                                }
113                            }
114                        contr = slope[i]*exp(A[i]*(r-r0)/thick[i]);
115
116                        vol = 4.0 * pi / 3.0 * r * r * r;
117                        //if (j == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && A[i]==0.0){
118                        //      vol_sub += (vol_pre - vol);
119                        //}
120                        f += vol * (contr * (fun) + (sld_in[i]-slope[i]) * bes);
121                        }
122                        break;
123            default :
124                                                if (fun_type[i]==0){
125                                                        slope[i] = 0.0;
126                                                }
127                                                else{
128                                                        slope[i]= (sld_out[i] -sld_in[i])/thick[i];
129                                                }
130                        contr = sld_in[i]-slope[i]*r;
131                        for (j=0; j<2; j++){
132                            if ( i == 0 && j == 0){
133                                continue;
134                                }
135                            if (i == n+1 && j == 1){
136                                continue;
137                                }
138                            if ( j == 1){
139                                sign = 1.0;
140                                r += thick[i];
141                                }
142                            else{
143                                sign = -1.0;
144                                                                }
145
146                            qr = q * r;
147                            fun = 0.0;
148                            if(qr == 0.0){
149                                bes = sign * 1.0;
150                                }
151                            else{
152                                bes = sign *  3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
153                                if (fabs(slope[i]) > 0.0 ){
154                                    fun = sign * 3.0 * r * (2.0*qr*sin(qr)-((qr*qr)-2.0)*cos(qr))/(qr * qr * qr * qr);
155                                    }
156                                }
157                            vol = 4.0 * pi / 3.0 * r * r * r;
158                            //if (j == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && fun_type[i]==0){
159                            //  vol_sub += (vol_pre - vol);
160                            //}
161                            f += vol * (bes * contr + fun * slope[i]);
162                            }
163                            break;
164                                }
165
166        }
167    //vol += vol_sub;
168    f2 = f * f / vol * 1.0e8;
169        f2 *= scale;
170        f2 += background;
171
172  free(sld_out);
173  free(slope);
174  free(sld_in);
175  free(thick);
176  free(A);
177  free(fun_type);
178
179    return (f2);
180}
181/**
182 * Function to evaluate 1D scattering function
183 * @param pars: parameters of the sphere
184 * @param q: q-value
185 * @return: function value
186 */
187double onion_analytical_1D(OnionParameters *pars, double q) {
188        double dp[56];
189
190        dp[0] = pars->n_shells;
191        dp[1] = pars->scale;
192        dp[2] = pars->rad_core0;
193        dp[3] = pars->sld_core0;
194        dp[4] = pars->sld_solv;
195        dp[5] = pars->background;
196
197        dp[6] = pars->sld_out_shell1;
198        dp[7] = pars->sld_out_shell2;
199        dp[8] = pars->sld_out_shell3;
200        dp[9] = pars->sld_out_shell4;
201        dp[10] = pars->sld_out_shell5;
202        dp[11] = pars->sld_out_shell6;
203        dp[12] = pars->sld_out_shell7;
204        dp[13] = pars->sld_out_shell8;
205        dp[14] = pars->sld_out_shell9;
206        dp[15] = pars->sld_out_shell10;
207
208        dp[16] = pars->sld_in_shell1;
209        dp[17] = pars->sld_in_shell2;
210        dp[18] = pars->sld_in_shell3;
211        dp[19] = pars->sld_in_shell4;
212        dp[20] = pars->sld_in_shell5;
213        dp[21] = pars->sld_in_shell6;
214        dp[22] = pars->sld_in_shell7;
215        dp[23] = pars->sld_in_shell8;
216        dp[24] = pars->sld_in_shell9;
217        dp[25] = pars->sld_in_shell10;
218
219        dp[26] = pars->A_shell1;
220        dp[27] = pars->A_shell2;
221        dp[28] = pars->A_shell3;
222        dp[29] = pars->A_shell4;
223        dp[30] = pars->A_shell5;
224        dp[31] = pars->A_shell6;
225        dp[32] = pars->A_shell7;
226        dp[33] = pars->A_shell8;
227        dp[34] = pars->A_shell9;
228        dp[35] = pars->A_shell10;
229
230        dp[36] = pars->thick_shell1;
231        dp[37] = pars->thick_shell2;
232        dp[38] = pars->thick_shell3;
233        dp[39] = pars->thick_shell4;
234        dp[40] = pars->thick_shell5;
235        dp[41] = pars->thick_shell6;
236        dp[42] = pars->thick_shell7;
237        dp[43] = pars->thick_shell8;
238        dp[44] = pars->thick_shell9;
239        dp[45] = pars->thick_shell10;
240
241        dp[46] = pars->func_shell1;
242        dp[47] = pars->func_shell2;
243        dp[48] = pars->func_shell3;
244        dp[49] = pars->func_shell4;
245        dp[50] = pars->func_shell5;
246        dp[51] = pars->func_shell6;
247        dp[52] = pars->func_shell7;
248        dp[53] = pars->func_shell8;
249        dp[54] = pars->func_shell9;
250        dp[55] = pars->func_shell10;
251
252        return so_kernel(dp, q);
253}
254
255/**
256 * Function to evaluate 2D scattering function
257 * @param pars: parameters of the sphere
258 * @param q: q-value
259 * @return: function value
260 */
261double onion_analytical_2D(OnionParameters *pars, double q, double phi) {
262        return onion_analytical_1D(pars,q);
263}
264
265double onion_analytical_2DXY(OnionParameters *pars, double qx, double qy){
266        return onion_analytical_1D(pars,sqrt(qx*qx+qy*qy));
267}
Note: See TracBrowser for help on using the repository browser.