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

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 36cbbe74 was 36cbbe74, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

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