[35aface] | 1 | /** |
---|
| 2 | * Scattering model for a sphere |
---|
| 3 | */ |
---|
| 4 | |
---|
| 5 | #include <math.h> |
---|
| 6 | #include "refl.h" |
---|
| 7 | #include <stdio.h> |
---|
| 8 | #include <stdlib.h> |
---|
| 9 | |
---|
| 10 | #define lamda 4.62 |
---|
| 11 | |
---|
| 12 | typedef struct { |
---|
| 13 | double re; |
---|
| 14 | double im; |
---|
| 15 | } complex; |
---|
| 16 | |
---|
| 17 | typedef struct { |
---|
| 18 | complex a; |
---|
| 19 | complex b; |
---|
| 20 | complex c; |
---|
| 21 | complex d; |
---|
| 22 | } matrix; |
---|
| 23 | |
---|
| 24 | complex cassign(real, imag) |
---|
| 25 | double real, imag; |
---|
| 26 | { |
---|
| 27 | complex x; |
---|
| 28 | x.re = real; |
---|
| 29 | x.im = imag; |
---|
| 30 | return x; |
---|
| 31 | } |
---|
| 32 | |
---|
| 33 | |
---|
| 34 | complex cadd(x,y) |
---|
| 35 | complex x,y; |
---|
| 36 | { |
---|
| 37 | complex z; |
---|
| 38 | z.re = x.re + y.re; |
---|
| 39 | z.im = x.im + y.im; |
---|
| 40 | return z; |
---|
| 41 | } |
---|
| 42 | |
---|
| 43 | complex rcmult(x,y) |
---|
| 44 | double x; |
---|
| 45 | complex y; |
---|
| 46 | { |
---|
| 47 | complex z; |
---|
| 48 | z.re = x*y.re; |
---|
| 49 | z.im = x*y.im; |
---|
| 50 | return z; |
---|
| 51 | } |
---|
| 52 | |
---|
| 53 | complex csub(x,y) |
---|
| 54 | complex x,y; |
---|
| 55 | { |
---|
| 56 | complex z; |
---|
| 57 | z.re = x.re - y.re; |
---|
| 58 | z.im = x.im - y.im; |
---|
| 59 | return z; |
---|
| 60 | } |
---|
| 61 | |
---|
| 62 | |
---|
| 63 | complex cmult(x,y) |
---|
| 64 | complex x,y; |
---|
| 65 | { |
---|
| 66 | complex z; |
---|
| 67 | z.re = x.re*y.re - x.im*y.im; |
---|
| 68 | z.im = x.re*y.im + x.im*y.re; |
---|
| 69 | return z; |
---|
| 70 | } |
---|
| 71 | |
---|
| 72 | complex cdiv(x,y) |
---|
| 73 | complex x,y; |
---|
| 74 | { |
---|
| 75 | complex z; |
---|
| 76 | z.re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im); |
---|
| 77 | z.im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im); |
---|
| 78 | return z; |
---|
| 79 | } |
---|
| 80 | |
---|
| 81 | complex cexp(b) |
---|
| 82 | complex b; |
---|
| 83 | { |
---|
| 84 | complex z; |
---|
| 85 | double br,bi; |
---|
| 86 | br=b.re; |
---|
| 87 | bi=b.im; |
---|
| 88 | z.re = exp(br)*cos(bi); |
---|
| 89 | z.im = exp(br)*sin(bi); |
---|
| 90 | return z; |
---|
| 91 | } |
---|
| 92 | |
---|
| 93 | |
---|
| 94 | complex csqrt(z) /* see Schaum`s Math Handbook p. 22, 6.6 and 6.10 */ |
---|
| 95 | complex z; |
---|
| 96 | { |
---|
| 97 | complex c; |
---|
| 98 | double zr,zi,x,y,r,w; |
---|
| 99 | |
---|
| 100 | zr=z.re; |
---|
| 101 | zi=z.im; |
---|
| 102 | |
---|
| 103 | if (zr==0.0 && zi==0.0) |
---|
| 104 | { |
---|
| 105 | c.re=0.0; |
---|
| 106 | c.im=0.0; |
---|
| 107 | return c; |
---|
| 108 | } |
---|
| 109 | else |
---|
| 110 | { |
---|
| 111 | x=fabs(zr); |
---|
| 112 | y=fabs(zi); |
---|
| 113 | if (x>y) |
---|
| 114 | { |
---|
| 115 | r=y/x; |
---|
| 116 | w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r))); |
---|
| 117 | } |
---|
| 118 | else |
---|
| 119 | { |
---|
| 120 | r=x/y; |
---|
| 121 | w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r))); |
---|
| 122 | } |
---|
| 123 | if (zr >=0.0) |
---|
| 124 | { |
---|
| 125 | c.re=w; |
---|
| 126 | c.im=zi/(2.0*w); |
---|
| 127 | } |
---|
| 128 | else |
---|
| 129 | { |
---|
| 130 | c.im=(zi >= 0) ? w : -w; |
---|
| 131 | c.re=zi/(2.0*c.im); |
---|
| 132 | } |
---|
| 133 | return c; |
---|
| 134 | } |
---|
| 135 | } |
---|
| 136 | |
---|
| 137 | complex ccos(b) |
---|
| 138 | complex b; |
---|
| 139 | { |
---|
| 140 | complex neg,negb,zero,two,z,i,bi,negbi; |
---|
| 141 | zero = cassign(0.0,0.0); |
---|
| 142 | two = cassign(2.0,0.0); |
---|
| 143 | i = cassign(0.0,1.0); |
---|
| 144 | bi = cmult(b,i); |
---|
| 145 | negbi = csub(zero,bi); |
---|
| 146 | z = cdiv(cadd(cexp(bi),cexp(negbi)),two); |
---|
| 147 | return z; |
---|
| 148 | } |
---|
| 149 | |
---|
| 150 | double errfunc(n_sub, i) |
---|
| 151 | double n_sub; |
---|
| 152 | int i; |
---|
| 153 | { |
---|
| 154 | double bin_size, ind, func; |
---|
| 155 | ind = i; |
---|
| 156 | // i range = [ -4..4], x range = [ -2.5..2.5] |
---|
| 157 | bin_size = n_sub/2.0/2.5; //size of each sub-layer |
---|
| 158 | // rescale erf so that 0 < erf < 1 in -2.5 <= x <= 2.5 |
---|
| 159 | func = (erf(ind/bin_size)/2.0+0.5); |
---|
| 160 | return func; |
---|
| 161 | } |
---|
| 162 | |
---|
| 163 | double linefunc(n_sub, i) |
---|
| 164 | double n_sub; |
---|
| 165 | int i; |
---|
| 166 | { |
---|
| 167 | double bin_size, ind, func; |
---|
| 168 | ind = i + 0.5; |
---|
| 169 | // i range = [ -4..4], x range = [ -2.5..2.5] |
---|
| 170 | bin_size = 1.0/n_sub; //size of each sub-layer |
---|
| 171 | // rescale erf so that 0 < erf < 1 in -2.5 <= x <= 2.5 |
---|
| 172 | func = ((ind + floor(n_sub/2.0))*bin_size); |
---|
| 173 | return func; |
---|
| 174 | } |
---|
| 175 | |
---|
| 176 | double parabolic_r(n_sub, i) |
---|
| 177 | double n_sub; |
---|
| 178 | int i; |
---|
| 179 | { |
---|
| 180 | double bin_size, ind, func; |
---|
| 181 | ind = i + 0.5; |
---|
| 182 | // i range = [ -4..4], x range = [ 0..1] |
---|
| 183 | bin_size = 1.0/n_sub; //size of each sub-layer; n_sub = 0 is a singular point (error) |
---|
| 184 | func = ((ind + floor(n_sub/2.0))*bin_size)*((ind + floor(n_sub/2.0))*bin_size); |
---|
| 185 | return func; |
---|
| 186 | } |
---|
| 187 | |
---|
| 188 | double parabolic_l(n_sub, i) |
---|
| 189 | double n_sub; |
---|
| 190 | int i; |
---|
| 191 | { |
---|
| 192 | double bin_size,ind, func; |
---|
| 193 | ind = i + 0.5; |
---|
| 194 | bin_size = 1.0/n_sub; //size of each sub-layer; n_sub = 0 is a singular point (error) |
---|
| 195 | func =1.0-(((ind + floor(n_sub/2.0))*bin_size) - 1.0) *(((ind + floor(n_sub/2.0))*bin_size) - 1.0); |
---|
| 196 | return func; |
---|
| 197 | } |
---|
| 198 | |
---|
| 199 | double cubic_r(n_sub, i) |
---|
| 200 | double n_sub; |
---|
| 201 | int i; |
---|
| 202 | { |
---|
| 203 | double bin_size,ind, func; |
---|
| 204 | ind = i + 0.5; |
---|
| 205 | // i range = [ -4..4], x range = [ 0..1] |
---|
| 206 | bin_size = 1.0/n_sub; //size of each sub-layer; n_sub = 0 is a singular point (error) |
---|
| 207 | func = ((ind+ floor(n_sub/2.0))*bin_size)*((ind + floor(n_sub/2.0))*bin_size)*((ind + floor(n_sub/2.0))*bin_size); |
---|
| 208 | return func; |
---|
| 209 | } |
---|
| 210 | |
---|
| 211 | double cubic_l(n_sub, i) |
---|
| 212 | double n_sub; |
---|
| 213 | int i; |
---|
| 214 | { |
---|
| 215 | double bin_size,ind, func; |
---|
| 216 | ind = i + 0.5; |
---|
| 217 | bin_size = 1.0/n_sub; //size of each sub-layer; n_sub = 0 is a singular point (error) |
---|
| 218 | func = 1.0+(((ind + floor(n_sub/2.0)))*bin_size - 1.0)*(((ind + floor(n_sub/2.0)))*bin_size - 1.0)*(((ind + floor(n_sub/2.0)))*bin_size - 1.0); |
---|
| 219 | return func; |
---|
| 220 | } |
---|
| 221 | |
---|
| 222 | double interfunc(fun_type, n_sub, i, sld_l, sld_r) |
---|
| 223 | double n_sub, sld_l, sld_r; |
---|
| 224 | int fun_type, i; |
---|
| 225 | { |
---|
| 226 | double sld_i, func; |
---|
| 227 | switch(fun_type){ |
---|
| 228 | case 1 : |
---|
| 229 | func = linefunc(n_sub, i); |
---|
| 230 | break; |
---|
| 231 | case 2 : |
---|
| 232 | func = parabolic_r(n_sub, i); |
---|
| 233 | break; |
---|
| 234 | case 3 : |
---|
| 235 | func = parabolic_l(n_sub, i); |
---|
| 236 | break; |
---|
| 237 | case 4 : |
---|
| 238 | func = cubic_r(n_sub, i); |
---|
| 239 | break; |
---|
| 240 | case 5 : |
---|
| 241 | func = cubic_l(n_sub, i); |
---|
| 242 | break; |
---|
| 243 | default: |
---|
| 244 | func = errfunc(n_sub, i); |
---|
| 245 | break; |
---|
| 246 | } |
---|
| 247 | if (sld_r>sld_l){ |
---|
| 248 | sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick); |
---|
| 249 | } |
---|
| 250 | else if (sld_r<sld_l){ |
---|
| 251 | func = 1.0-func; |
---|
| 252 | sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick); |
---|
| 253 | } |
---|
| 254 | else{ |
---|
| 255 | sld_i = sld_r; |
---|
| 256 | } |
---|
| 257 | return sld_i; |
---|
| 258 | } |
---|
| 259 | |
---|
| 260 | double re_kernel(double dp[], double q) { |
---|
| 261 | int n = dp[0]; |
---|
| 262 | int i,j,fun_type[n+2]; |
---|
| 263 | |
---|
| 264 | double scale = dp[1]; |
---|
| 265 | double thick_inter_sub = dp[2]; |
---|
| 266 | double sld_sub = dp[4]; |
---|
| 267 | double sld_super = dp[5]; |
---|
| 268 | double background = dp[6]; |
---|
| 269 | |
---|
[339ce67] | 270 | double sld[n+2],sld_im[n+2],thick_inter[n+2],thick[n+2],total_thick; |
---|
[35aface] | 271 | fun_type[0] = dp[3]; |
---|
| 272 | for (i =1; i<=n; i++){ |
---|
| 273 | sld[i] = dp[i+6]; |
---|
| 274 | thick_inter[i]= dp[i+16]; |
---|
| 275 | thick[i] = dp[i+26]; |
---|
| 276 | fun_type[i] = dp[i+36]; |
---|
[339ce67] | 277 | sld_im[i] = dp[i+46]; |
---|
[35aface] | 278 | |
---|
| 279 | total_thick += thick[i] + thick_inter[i]; |
---|
| 280 | } |
---|
| 281 | sld[0] = sld_sub; |
---|
| 282 | sld[n+1] = sld_super; |
---|
[339ce67] | 283 | sld_im[0] = fabs(dp[0+56]); |
---|
| 284 | sld_im[n+1] = fabs(dp[1+56]); |
---|
[35aface] | 285 | thick[0] = total_thick/5.0; |
---|
| 286 | thick[n+1] = total_thick/5.0; |
---|
| 287 | thick_inter[0] = thick_inter_sub; |
---|
| 288 | thick_inter[n+1] = 0.0; |
---|
| 289 | |
---|
| 290 | double nsl=21.0; //nsl = Num_sub_layer: MUST ODD number in double //no other number works now |
---|
| 291 | int n_s, floor_nsl; |
---|
[339ce67] | 292 | double sld_i,sldim_i,dz,phi,R,ko2; |
---|
[35aface] | 293 | double sign,erfunc; |
---|
| 294 | double pi; |
---|
| 295 | complex inv_n,phi1,alpha,alpha2,kn,fnm,fnp,rn,Xn,nn,nn2,an,nnp1,one,zero,two,n_sub,n_sup,knp1,Xnp1; |
---|
| 296 | pi = 4.0*atan(1.0); |
---|
| 297 | one = cassign(1.0,0.0); |
---|
| 298 | //zero = cassign(0.0,0.0); |
---|
| 299 | two= cassign(0.0,-2.0); |
---|
| 300 | |
---|
| 301 | floor_nsl = floor(nsl/2.0); |
---|
| 302 | |
---|
| 303 | //Checking if floor is available. |
---|
| 304 | //no imaginary sld inputs in this function yet |
---|
[339ce67] | 305 | n_sub=cassign(1.0-sld_sub*pow(lamda,2.0)/(2.0*pi),pow(lamda,2.0)/(2.0*pi)*sld_im[0]); |
---|
| 306 | 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]); |
---|
[35aface] | 307 | ko2 = pow(2.0*pi/lamda,2.0); |
---|
| 308 | |
---|
| 309 | phi = asin(lamda*q/(4.0*pi)); |
---|
| 310 | phi1 = cdiv(rcmult(phi,one),n_sup); |
---|
| 311 | alpha = cmult(n_sup,ccos(phi1)); |
---|
| 312 | alpha2 = cmult(alpha,alpha); |
---|
| 313 | |
---|
| 314 | nnp1=n_sub; |
---|
| 315 | knp1=csqrt(rcmult(ko2,csub(cmult(nnp1,nnp1),alpha2))); //nnp1*ko*sin(phinp1) |
---|
| 316 | Xnp1=cassign(0.0,0.0); |
---|
| 317 | dz = 0.0; |
---|
| 318 | // iteration for # of layers +sub from the top |
---|
| 319 | for (i=1;i<=n+1; i++){ |
---|
| 320 | //iteration for 9 sub-layers |
---|
| 321 | for (j=0;j<2;j++){ |
---|
| 322 | for (n_s=-floor_nsl;n_s<=floor_nsl; n_s++){ |
---|
| 323 | if (j==1){ |
---|
| 324 | if (i==n+1) |
---|
| 325 | break; |
---|
| 326 | dz = thick[i]; |
---|
| 327 | sld_i = sld[i]; |
---|
[339ce67] | 328 | sldim_i = sld_im[i]; |
---|
[35aface] | 329 | } |
---|
| 330 | else{ |
---|
| 331 | dz = thick_inter[i-1]/nsl; |
---|
[339ce67] | 332 | if (sld[i-1] == sld[i]){ |
---|
| 333 | sld_i = sld[i]; |
---|
| 334 | } |
---|
| 335 | else{ |
---|
| 336 | sld_i = interfunc(fun_type[i-1],nsl, n_s, sld[i-1], sld[i]); |
---|
| 337 | } |
---|
| 338 | if (sld_im[i-1] == sld_im[i]){ |
---|
| 339 | sldim_i = sld_im[i]; |
---|
| 340 | } |
---|
| 341 | else{ |
---|
| 342 | sldim_i = interfunc(fun_type[i-1],nsl, n_s, sld_im[i-1], sld_im[i]); |
---|
| 343 | } |
---|
[35aface] | 344 | } |
---|
[339ce67] | 345 | nn = cassign(1.0-sld_i*pow(lamda,2.0)/(2.0*pi),pow(lamda,2.0)/(2.0*pi)*sldim_i); |
---|
[35aface] | 346 | nn2=cmult(nn,nn); |
---|
| 347 | |
---|
| 348 | kn=csqrt(rcmult(ko2,csub(nn2,alpha2))); //nn*ko*sin(phin) |
---|
| 349 | an=cexp(rcmult(dz,cmult(two,kn))); |
---|
| 350 | |
---|
| 351 | fnm=csub(kn,knp1); |
---|
| 352 | fnp=cadd(kn,knp1); |
---|
| 353 | rn=cdiv(fnm,fnp); |
---|
| 354 | Xn=cmult(an,cdiv(cadd(rn,Xnp1),cadd(one,cmult(rn,Xnp1)))); //Xn=an*((rn+Xnp1*anp1)/(1+rn*Xnp1*anp1)) |
---|
| 355 | |
---|
| 356 | Xnp1=Xn; |
---|
| 357 | knp1=kn; |
---|
| 358 | |
---|
| 359 | if (j==1) |
---|
| 360 | break; |
---|
| 361 | } |
---|
| 362 | } |
---|
| 363 | } |
---|
| 364 | R=pow(Xn.re,2.0)+pow(Xn.im,2.0); |
---|
| 365 | // This temperarily fixes the total reflection for Rfunction and linear. |
---|
| 366 | // ToDo: Show why it happens that Xn.re=0 and Xn.im >1! |
---|
| 367 | if (Xn.im == 0.0){ |
---|
| 368 | R=1.0; |
---|
| 369 | } |
---|
| 370 | R *= scale; |
---|
| 371 | R += background; |
---|
| 372 | |
---|
| 373 | return R; |
---|
| 374 | |
---|
| 375 | } |
---|
| 376 | /** |
---|
| 377 | * Function to evaluate 1D scattering function |
---|
| 378 | * @param pars: parameters of the sphere |
---|
| 379 | * @param q: q-value |
---|
| 380 | * @return: function value |
---|
| 381 | */ |
---|
| 382 | double refl_analytical_1D(ReflParameters *pars, double q) { |
---|
[339ce67] | 383 | double dp[59]; |
---|
[35aface] | 384 | |
---|
| 385 | dp[0] = pars->n_layers; |
---|
| 386 | dp[1] = pars->scale; |
---|
| 387 | dp[2] = pars->thick_inter0; |
---|
| 388 | dp[3] = pars->func_inter0; |
---|
| 389 | dp[4] = pars->sld_sub0; |
---|
| 390 | dp[5] = pars->sld_medium; |
---|
| 391 | dp[6] = pars->background; |
---|
| 392 | |
---|
| 393 | dp[7] = pars->sld_flat1; |
---|
| 394 | dp[8] = pars->sld_flat2; |
---|
| 395 | dp[9] = pars->sld_flat3; |
---|
| 396 | dp[10] = pars->sld_flat4; |
---|
| 397 | dp[11] = pars->sld_flat5; |
---|
| 398 | dp[12] = pars->sld_flat6; |
---|
| 399 | dp[13] = pars->sld_flat7; |
---|
| 400 | dp[14] = pars->sld_flat8; |
---|
| 401 | dp[15] = pars->sld_flat9; |
---|
| 402 | dp[16] = pars->sld_flat10; |
---|
| 403 | |
---|
| 404 | dp[17] = pars->thick_inter1; |
---|
| 405 | dp[18] = pars->thick_inter2; |
---|
| 406 | dp[19] = pars->thick_inter3; |
---|
| 407 | dp[20] = pars->thick_inter4; |
---|
| 408 | dp[21] = pars->thick_inter5; |
---|
| 409 | dp[22] = pars->thick_inter6; |
---|
| 410 | dp[23] = pars->thick_inter7; |
---|
| 411 | dp[24] = pars->thick_inter8; |
---|
| 412 | dp[25] = pars->thick_inter9; |
---|
| 413 | dp[26] = pars->thick_inter10; |
---|
| 414 | |
---|
| 415 | dp[27] = pars->thick_flat1; |
---|
| 416 | dp[28] = pars->thick_flat2; |
---|
| 417 | dp[29] = pars->thick_flat3; |
---|
| 418 | dp[30] = pars->thick_flat4; |
---|
| 419 | dp[31] = pars->thick_flat5; |
---|
| 420 | dp[32] = pars->thick_flat6; |
---|
| 421 | dp[33] = pars->thick_flat7; |
---|
| 422 | dp[34] = pars->thick_flat8; |
---|
| 423 | dp[35] = pars->thick_flat9; |
---|
| 424 | dp[36] = pars->thick_flat10; |
---|
| 425 | |
---|
| 426 | dp[37] = pars->func_inter1; |
---|
| 427 | dp[38] = pars->func_inter2; |
---|
| 428 | dp[39] = pars->func_inter3; |
---|
| 429 | dp[40] = pars->func_inter4; |
---|
| 430 | dp[41] = pars->func_inter5; |
---|
| 431 | dp[42] = pars->func_inter6; |
---|
| 432 | dp[43] = pars->func_inter7; |
---|
| 433 | dp[44] = pars->func_inter8; |
---|
| 434 | dp[45] = pars->func_inter9; |
---|
| 435 | dp[46] = pars->func_inter10; |
---|
| 436 | |
---|
[339ce67] | 437 | dp[47] = pars->sldIM_flat1; |
---|
| 438 | dp[48] = pars->sldIM_flat2; |
---|
| 439 | dp[49] = pars->sldIM_flat3; |
---|
| 440 | dp[50] = pars->sldIM_flat4; |
---|
| 441 | dp[51] = pars->sldIM_flat5; |
---|
| 442 | dp[52] = pars->sldIM_flat6; |
---|
| 443 | dp[53] = pars->sldIM_flat7; |
---|
| 444 | dp[54] = pars->sldIM_flat8; |
---|
| 445 | dp[55] = pars->sldIM_flat9; |
---|
| 446 | dp[56] = pars->sldIM_flat10; |
---|
| 447 | |
---|
| 448 | dp[57] = pars->sldIM_sub0; |
---|
| 449 | dp[58] = pars->sldIM_medium; |
---|
| 450 | |
---|
[35aface] | 451 | return re_kernel(dp, q); |
---|
| 452 | } |
---|
| 453 | |
---|
| 454 | /** |
---|
| 455 | * Function to evaluate 2D scattering function |
---|
| 456 | * @param pars: parameters of the sphere |
---|
| 457 | * @param q: q-value |
---|
| 458 | * @return: function value |
---|
| 459 | */ |
---|
| 460 | double refl_analytical_2D(ReflParameters *pars, double q, double phi) { |
---|
| 461 | return refl_analytical_1D(pars,q); |
---|
| 462 | } |
---|
| 463 | |
---|
| 464 | double refl_analytical_2DXY(ReflParameters *pars, double qx, double qy){ |
---|
| 465 | return refl_analytical_1D(pars,sqrt(qx*qx+qy*qy)); |
---|
| 466 | } |
---|