source: sasview/sansmodels/src/sans/models/c_extensions/refl_adv.c @ d760c0a

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

Re #4 A few more warnings

  • Property mode set to 100644
File size: 7.0 KB
Line 
1// The original code, of which work was not DANSE funded,
2// was provided by J. Cho.
3/**
4 * NR model Parratt method
5 */
6#include <math.h>
7#include "refl_adv.h"
8#include "libmultifunc/librefl.h"
9#include <stdio.h>
10#include <stdlib.h>
11
12#define lamda 4.62
13
14
15double re_adv_kernel(double dp[], double q) {
16        int n = dp[0];
17        int i,j;
18        double nsl;
19
20        double scale = dp[1];
21        double thick_inter_sub = dp[2];
22        double sld_sub = dp[4];
23        double sld_super = dp[5];
24        double background = dp[6];
25        double npts = dp[69]; //number of sub_layers in each interface
26
27        double total_thick=0.0;
28
29  int n_s;
30  double sld_i,sldim_i,dz,phi,R,ko2;
31  double pi;
32
33  int* fun_type;
34  double* sld;
35  double* sld_im;
36  double* thick_inter;
37  double* thick;
38  double* fun_coef;
39  complex  phi1,alpha,alpha2,kn,fnm,fnp,rn,Xn,nn,nn2,an,nnp1,one,two,n_sub,n_sup,knp1,Xnp1;
40
41  fun_type = (int*)malloc((n+2)*sizeof(int));
42  sld = (double*)malloc((n+2)*sizeof(double));
43  sld_im = (double*)malloc((n+2)*sizeof(double));
44  thick_inter = (double*)malloc((n+2)*sizeof(double));
45  thick = (double*)malloc((n+2)*sizeof(double));
46  fun_coef = (double*)malloc((n+2)*sizeof(double));
47
48  fun_type[0] = dp[3];
49        fun_coef[0] = fabs(dp[70]);
50        for (i =1; i<=n; i++){
51                sld[i] = dp[i+6];
52                thick_inter[i]= dp[i+16];
53                thick[i] = dp[i+26];
54                fun_type[i] = dp[i+36];
55                sld_im[i] = dp[i+46];
56                fun_coef[i] = fabs(dp[i+56]);
57                //printf("type_func2 =%g\n",fun_coef[i]);
58
59                total_thick += thick[i] + thick_inter[i];
60        }
61        sld[0] = sld_sub;
62        sld[n+1] = sld_super;
63        sld_im[0] = fabs(dp[1+66]);
64        sld_im[n+1] = fabs(dp[2+66]);
65        thick[0] = total_thick/5.0;
66        thick[n+1] = total_thick/5.0;
67        thick_inter[0] = thick_inter_sub;
68        thick_inter[n+1] = 0.0;
69        fun_coef[n+1] = 0.0;
70
71        nsl=npts;//21.0; //nsl = Num_sub_layer:  MUST ODD number in double //no other number works now
72
73        pi = 4.0*atan(1.0);
74    one = cassign(1.0,0.0);
75        Xn = cassign(0.0,0.0);
76        two = cassign(0.0,-2.0);
77
78        //Checking if floor is available.
79        //no imaginary sld inputs in this function yet
80    n_sub=cassign(1.0-sld_sub*pow(lamda,2.0)/(2.0*pi),pow(lamda,2.0)/(2.0*pi)*sld_im[0]);
81    n_sup=cassign(1.0-sld_super*pow(lamda,2.0)/(2.0*pi),pow(lamda,2.0)/(2.0*pi)*sld_im[n+1]);
82    ko2 = pow(2.0*pi/lamda,2.0);
83
84    phi = asin(lamda*q/(4.0*pi));
85    phi1 = cdiv(rcmult(phi,one),n_sup);
86    alpha = cmult(n_sup,ccos(phi1));
87        alpha2 = cmult(alpha,alpha);
88
89    nnp1=n_sub;
90    knp1=csqrt(rcmult(ko2,csub(cmult(nnp1,nnp1),alpha2)));  //nnp1*ko*sin(phinp1)
91    Xnp1=cassign(0.0,0.0);
92    dz = 0.0;
93    // iteration for # of layers +sub from the top
94    for (i=1;i<=n+1; i++){
95        //if (fun_coef[i]==0.0)
96        //      // this condition protects an error in numerical multiplication
97        //      fun_coef[i] = 1e-14;
98                //iteration for 9 sub-layers
99                for (j=0;j<2;j++){
100                        for (n_s=0;n_s<nsl; n_s++){
101                                // for flat layer
102                                if (j==1){
103                                        if (i==n+1)
104                                                break;
105                                        dz = thick[i];
106                                        sld_i = sld[i];
107                                        sldim_i = sld_im[i];
108                                }
109                                // for interface
110                                else{
111                                        dz = thick_inter[i-1]/nsl;
112                                        if (sld[i-1] == sld[i]){
113                                                sld_i = sld[i];
114                                        }
115                                        else{
116                                                sld_i = intersldfunc(fun_type[i-1],nsl, n_s+0.5, fun_coef[i-1], sld[i-1], sld[i]);
117                                        }
118                                        if (sld_im[i-1] == sld_im[i]){
119                                                sldim_i = sld_im[i];
120                                        }
121                                        else{
122                                                sldim_i = intersldfunc(fun_type[i-1],nsl, n_s+0.5, fun_coef[i-1], sld_im[i-1], sld_im[i]);
123                                        }
124                                }
125                                nn = cassign(1.0-sld_i*pow(lamda,2.0)/(2.0*pi),pow(lamda,2.0)/(2.0*pi)*sldim_i);
126                                nn2=cmult(nn,nn);
127
128                                kn=csqrt(rcmult(ko2,csub(nn2,alpha2)));        //nn*ko*sin(phin)
129                                an=cexp(rcmult(dz,cmult(two,kn)));
130
131                                fnm=csub(kn,knp1);
132                                fnp=cadd(kn,knp1);
133                                rn=cdiv(fnm,fnp);
134                                Xn=cmult(an,cdiv(cadd(rn,Xnp1),cadd(one,cmult(rn,Xnp1))));    //Xn=an*((rn+Xnp1*anp1)/(1+rn*Xnp1*anp1))
135
136                                Xnp1=Xn;
137                                knp1=kn;
138                                // no for-loop for flat layer
139                                if (j==1)
140                                        break;
141                        }
142                }
143    }
144    R=pow(Xn.re,2.0)+pow(Xn.im,2.0);
145    // This temperarily fixes the total reflection for Rfunction and linear.
146    // ToDo: Show why it happens that Xn.re=0 and Xn.im >1!
147    if (Xn.im == 0.0 || R > 1){
148        R=1.0;
149    }
150    R *= scale;
151    R += background;
152
153    free(fun_type);
154    free(sld);
155    free(sld_im);
156    free(thick_inter);
157    free(thick);
158    free(fun_coef);
159
160    return R;
161
162}
163
164/**
165 * Function to evaluate NR function
166 * @param pars: parameters of refl
167 * @param q: q-value
168 * @return: function value
169 */
170
171double refl_adv_analytical_1D(ReflAdvParameters *pars, double q) {
172        double dp[71];
173
174        dp[0] = pars->n_layers;
175        dp[1] = pars->scale;
176        dp[2] = pars->thick_inter0;
177        dp[3] = pars->func_inter0;
178        dp[4] = pars->sld_bottom0;
179        dp[5] = pars->sld_medium;
180        dp[6] = pars->background;
181
182        dp[7] = pars->sld_flat1;
183        dp[8] = pars->sld_flat2;
184        dp[9] = pars->sld_flat3;
185        dp[10] = pars->sld_flat4;
186        dp[11] = pars->sld_flat5;
187        dp[12] = pars->sld_flat6;
188        dp[13] = pars->sld_flat7;
189        dp[14] = pars->sld_flat8;
190        dp[15] = pars->sld_flat9;
191        dp[16] = pars->sld_flat10;
192
193        dp[17] = pars->thick_inter1;
194        dp[18] = pars->thick_inter2;
195        dp[19] = pars->thick_inter3;
196        dp[20] = pars->thick_inter4;
197        dp[21] = pars->thick_inter5;
198        dp[22] = pars->thick_inter6;
199        dp[23] = pars->thick_inter7;
200        dp[24] = pars->thick_inter8;
201        dp[25] = pars->thick_inter9;
202        dp[26] = pars->thick_inter10;
203
204        dp[27] = pars->thick_flat1;
205        dp[28] = pars->thick_flat2;
206        dp[29] = pars->thick_flat3;
207        dp[30] = pars->thick_flat4;
208        dp[31] = pars->thick_flat5;
209        dp[32] = pars->thick_flat6;
210        dp[33] = pars->thick_flat7;
211        dp[34] = pars->thick_flat8;
212        dp[35] = pars->thick_flat9;
213        dp[36] = pars->thick_flat10;
214
215        dp[37] = pars->func_inter1;
216        dp[38] = pars->func_inter2;
217        dp[39] = pars->func_inter3;
218        dp[40] = pars->func_inter4;
219        dp[41] = pars->func_inter5;
220        dp[42] = pars->func_inter6;
221        dp[43] = pars->func_inter7;
222        dp[44] = pars->func_inter8;
223        dp[45] = pars->func_inter9;
224        dp[46] = pars->func_inter10;
225
226        dp[47] = pars->sldIM_flat1;
227        dp[48] = pars->sldIM_flat2;
228        dp[49] = pars->sldIM_flat3;
229        dp[50] = pars->sldIM_flat4;
230        dp[51] = pars->sldIM_flat5;
231        dp[52] = pars->sldIM_flat6;
232        dp[53] = pars->sldIM_flat7;
233        dp[54] = pars->sldIM_flat8;
234        dp[55] = pars->sldIM_flat9;
235        dp[56] = pars->sldIM_flat10;
236
237        dp[57] = pars->nu_inter1;
238        dp[58] = pars->nu_inter2;
239        dp[59] = pars->nu_inter3;
240        dp[60] = pars->nu_inter4;
241        dp[61] = pars->nu_inter5;
242        dp[62] = pars->nu_inter6;
243        dp[63] = pars->nu_inter7;
244        dp[64] = pars->nu_inter8;
245        dp[65] = pars->nu_inter9;
246        dp[66] = pars->nu_inter10;
247
248        dp[67] = pars->sldIM_sub0;
249        dp[68] = pars->sldIM_medium;
250        dp[69] = pars->npts_inter;
251        dp[70] = pars->nu_inter0;
252
253        return re_adv_kernel(dp, q);
254}
255
256/**
257 * Function to evaluate NR function
258 * @param pars: parameters of NR
259 * @param q: q-value
260 * @return: function value
261 */
262double refl_adv_analytical_2D(ReflAdvParameters *pars, double q, double phi) {
263        return refl_adv_analytical_1D(pars,q);
264}
265
266double refl_adv_analytical_2DXY(ReflAdvParameters *pars, double qx, double qy){
267        return refl_adv_analytical_1D(pars,sqrt(qx*qx+qy*qy));
268}
Note: See TracBrowser for help on using the repository browser.