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

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

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