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

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 c09ace41 was 96656e3, checked in by Jae Cho <jhjcho@…>, 14 years ago

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