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

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 75fbd17 was 0164899a, checked in by Jae Cho <jhjcho@…>, 14 years ago

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