[97] | 1 | /* CylinderFit.c |
---|
| 2 | |
---|
| 3 | A simplified project designed to act as a template for your curve fitting function. |
---|
| 4 | The fitting function is a Cylinder form factor. No resolution effects are included (yet) |
---|
| 5 | */ |
---|
| 6 | |
---|
| 7 | #include "StandardHeaders.h" // Include ANSI headers, Mac headers |
---|
| 8 | #include "GaussWeights.h" |
---|
| 9 | #include "libCylinder.h" |
---|
| 10 | |
---|
| 11 | /* CylinderForm : calculates the form factor of a cylinder at the give x-value p->x |
---|
| 12 | |
---|
| 13 | Warning: |
---|
| 14 | The call to WaveData() below returns a pointer to the middle |
---|
| 15 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 16 | calculations could cause memory to move, you should copy the coefficient |
---|
| 17 | values to local variables or an array before such operations. |
---|
| 18 | */ |
---|
| 19 | double |
---|
| 20 | CylinderForm(double dp[], double q) |
---|
| 21 | { |
---|
| 22 | int i; |
---|
| 23 | double Pi; |
---|
[235] | 24 | double scale,radius,length,delrho,bkg,halfheight,sldCyl,sldSolv; //local variables of coefficient wave |
---|
[97] | 25 | int nord=76; //order of integration |
---|
| 26 | double uplim,lolim; //upper and lower integration limits |
---|
| 27 | double summ,zi,yyy,answer,vcyl; //running tally of integration |
---|
| 28 | |
---|
| 29 | Pi = 4.0*atan(1.0); |
---|
[632] | 30 | lolim = 0.0; |
---|
[97] | 31 | uplim = Pi/2.0; |
---|
| 32 | |
---|
| 33 | summ = 0.0; //initialize intergral |
---|
| 34 | |
---|
| 35 | scale = dp[0]; //make local copies in case memory moves |
---|
| 36 | radius = dp[1]; |
---|
| 37 | length = dp[2]; |
---|
[235] | 38 | sldCyl = dp[3]; |
---|
| 39 | sldSolv = dp[4]; |
---|
| 40 | bkg = dp[5]; |
---|
| 41 | |
---|
| 42 | delrho = sldCyl-sldSolv; |
---|
[97] | 43 | halfheight = length/2.0; |
---|
| 44 | for(i=0;i<nord;i++) { |
---|
| 45 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 46 | yyy = Gauss76Wt[i] * CylKernel(q, radius, halfheight, zi); |
---|
| 47 | summ += yyy; |
---|
| 48 | } |
---|
| 49 | |
---|
| 50 | answer = (uplim-lolim)/2.0*summ; |
---|
| 51 | // Multiply by contrast^2 |
---|
| 52 | answer *= delrho*delrho; |
---|
| 53 | //normalize by cylinder volume |
---|
| 54 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
| 55 | vcyl=Pi*radius*radius*length; |
---|
| 56 | answer *= vcyl; |
---|
| 57 | //convert to [cm-1] |
---|
| 58 | answer *= 1.0e8; |
---|
| 59 | //Scale |
---|
| 60 | answer *= scale; |
---|
| 61 | // add in the background |
---|
| 62 | answer += bkg; |
---|
| 63 | |
---|
| 64 | return answer; |
---|
| 65 | } |
---|
| 66 | |
---|
| 67 | /* EllipCyl76X : calculates the form factor of a elliptical cylinder at the given x-value p->x |
---|
| 68 | |
---|
| 69 | Uses 76 pt Gaussian quadrature for both integrals |
---|
| 70 | |
---|
| 71 | Warning: |
---|
| 72 | The call to WaveData() below returns a pointer to the middle |
---|
| 73 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 74 | calculations could cause memory to move, you should copy the coefficient |
---|
| 75 | values to local variables or an array before such operations. |
---|
| 76 | */ |
---|
| 77 | double |
---|
| 78 | EllipCyl76(double dp[], double q) |
---|
| 79 | { |
---|
| 80 | int i,j; |
---|
[235] | 81 | double Pi,slde,sld; |
---|
[97] | 82 | double scale,ra,nu,length,delrho,bkg; //local variables of coefficient wave |
---|
| 83 | int nord=76; //order of integration |
---|
| 84 | double va,vb; //upper and lower integration limits |
---|
| 85 | double summ,zi,yyy,answer,vell; //running tally of integration |
---|
[632] | 86 | double summj,vaj,vbj,zij,arg, si; //for the inner integration |
---|
| 87 | |
---|
[97] | 88 | Pi = 4.0*atan(1.0); |
---|
[632] | 89 | va = 0.0; |
---|
| 90 | vb = 1.0; //orintational average, outer integral |
---|
| 91 | vaj=0.0; |
---|
[97] | 92 | vbj=Pi; //endpoints of inner integral |
---|
| 93 | |
---|
| 94 | summ = 0.0; //initialize intergral |
---|
| 95 | |
---|
| 96 | scale = dp[0]; //make local copies in case memory moves |
---|
| 97 | ra = dp[1]; |
---|
| 98 | nu = dp[2]; |
---|
| 99 | length = dp[3]; |
---|
[235] | 100 | slde = dp[4]; |
---|
| 101 | sld = dp[5]; |
---|
| 102 | delrho = slde - sld; |
---|
| 103 | bkg = dp[6]; |
---|
| 104 | |
---|
[97] | 105 | for(i=0;i<nord;i++) { |
---|
| 106 | //setup inner integral over the ellipsoidal cross-section |
---|
| 107 | summj=0; |
---|
| 108 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy |
---|
[632] | 109 | arg = ra*sqrt(1.0-zi*zi); |
---|
[97] | 110 | for(j=0;j<nord;j++) { |
---|
| 111 | //76 gauss points for the inner integral as well |
---|
| 112 | zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy |
---|
| 113 | yyy = Gauss76Wt[j] * EllipCylKernel(q,arg,nu,zij); |
---|
| 114 | summj += yyy; |
---|
| 115 | } |
---|
| 116 | //now calculate the value of the inner integral |
---|
| 117 | answer = (vbj-vaj)/2.0*summj; |
---|
| 118 | //divide integral by Pi |
---|
| 119 | answer /=Pi; |
---|
| 120 | |
---|
| 121 | //now calculate outer integral |
---|
[632] | 122 | arg = q*length*zi/2.0; |
---|
| 123 | if (arg == 0.0){ |
---|
| 124 | si = 1.0; |
---|
| 125 | }else{ |
---|
| 126 | si = sin(arg) * sin(arg) / arg / arg; |
---|
| 127 | } |
---|
| 128 | yyy = Gauss76Wt[i] * answer * si; |
---|
[97] | 129 | summ += yyy; |
---|
| 130 | } |
---|
| 131 | answer = (vb-va)/2.0*summ; |
---|
| 132 | // Multiply by contrast^2 |
---|
| 133 | answer *= delrho*delrho; |
---|
| 134 | //normalize by cylinder volume |
---|
| 135 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
| 136 | vell = Pi*ra*(nu*ra)*length; |
---|
| 137 | answer *= vell; |
---|
| 138 | //convert to [cm-1] |
---|
| 139 | answer *= 1.0e8; |
---|
| 140 | //Scale |
---|
| 141 | answer *= scale; |
---|
| 142 | // add in the background |
---|
| 143 | answer += bkg; |
---|
| 144 | |
---|
| 145 | return answer; |
---|
| 146 | } |
---|
| 147 | |
---|
| 148 | /* EllipCyl20X : calculates the form factor of a elliptical cylinder at the given x-value p->x |
---|
| 149 | |
---|
| 150 | Uses 76 pt Gaussian quadrature for orientational integral |
---|
| 151 | Uses 20 pt quadrature for the inner integral over the elliptical cross-section |
---|
| 152 | |
---|
| 153 | Warning: |
---|
| 154 | The call to WaveData() below returns a pointer to the middle |
---|
| 155 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 156 | calculations could cause memory to move, you should copy the coefficient |
---|
| 157 | values to local variables or an array before such operations. |
---|
| 158 | */ |
---|
| 159 | double |
---|
| 160 | EllipCyl20(double dp[], double q) |
---|
| 161 | { |
---|
| 162 | int i,j; |
---|
[235] | 163 | double Pi,slde,sld; |
---|
[97] | 164 | double scale,ra,nu,length,delrho,bkg; //local variables of coefficient wave |
---|
| 165 | int nordi=76; //order of integration |
---|
| 166 | int nordj=20; |
---|
| 167 | double va,vb; //upper and lower integration limits |
---|
| 168 | double summ,zi,yyy,answer,vell; //running tally of integration |
---|
[632] | 169 | double summj,vaj,vbj,zij,arg,si; //for the inner integration |
---|
| 170 | |
---|
[97] | 171 | Pi = 4.0*atan(1.0); |
---|
[632] | 172 | va = 0.0; |
---|
| 173 | vb = 1.0; //orintational average, outer integral |
---|
| 174 | vaj=0.0; |
---|
[97] | 175 | vbj=Pi; //endpoints of inner integral |
---|
| 176 | |
---|
| 177 | summ = 0.0; //initialize intergral |
---|
| 178 | |
---|
| 179 | scale = dp[0]; //make local copies in case memory moves |
---|
| 180 | ra = dp[1]; |
---|
| 181 | nu = dp[2]; |
---|
| 182 | length = dp[3]; |
---|
[235] | 183 | slde = dp[4]; |
---|
| 184 | sld = dp[5]; |
---|
| 185 | delrho = slde - sld; |
---|
| 186 | bkg = dp[6]; |
---|
| 187 | |
---|
[97] | 188 | for(i=0;i<nordi;i++) { |
---|
| 189 | //setup inner integral over the ellipsoidal cross-section |
---|
| 190 | summj=0; |
---|
| 191 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy |
---|
[632] | 192 | arg = ra*sqrt(1.0-zi*zi); |
---|
[97] | 193 | for(j=0;j<nordj;j++) { |
---|
| 194 | //20 gauss points for the inner integral |
---|
| 195 | zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy |
---|
| 196 | yyy = Gauss20Wt[j] * EllipCylKernel(q,arg,nu,zij); |
---|
| 197 | summj += yyy; |
---|
| 198 | } |
---|
| 199 | //now calculate the value of the inner integral |
---|
| 200 | answer = (vbj-vaj)/2.0*summj; |
---|
| 201 | //divide integral by Pi |
---|
| 202 | answer /=Pi; |
---|
| 203 | |
---|
| 204 | //now calculate outer integral |
---|
[632] | 205 | arg = q*length*zi/2.0; |
---|
| 206 | if (arg == 0.0){ |
---|
| 207 | si = 1.0; |
---|
| 208 | }else{ |
---|
| 209 | si = sin(arg) * sin(arg) / arg / arg; |
---|
| 210 | } |
---|
| 211 | yyy = Gauss76Wt[i] * answer * si; |
---|
[97] | 212 | summ += yyy; |
---|
| 213 | } |
---|
| 214 | |
---|
| 215 | answer = (vb-va)/2.0*summ; |
---|
| 216 | // Multiply by contrast^2 |
---|
| 217 | answer *= delrho*delrho; |
---|
| 218 | //normalize by cylinder volume |
---|
| 219 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
| 220 | vell = Pi*ra*(nu*ra)*length; |
---|
| 221 | answer *= vell; |
---|
| 222 | //convert to [cm-1] |
---|
| 223 | answer *= 1.0e8; |
---|
| 224 | //Scale |
---|
| 225 | answer *= scale; |
---|
| 226 | // add in the background |
---|
[632] | 227 | answer += bkg; |
---|
| 228 | |
---|
[97] | 229 | return answer; |
---|
| 230 | } |
---|
| 231 | |
---|
| 232 | /* TriaxialEllipsoidX : calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x |
---|
| 233 | |
---|
| 234 | Uses 76 pt Gaussian quadrature for both integrals |
---|
| 235 | |
---|
| 236 | Warning: |
---|
| 237 | The call to WaveData() below returns a pointer to the middle |
---|
| 238 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 239 | calculations could cause memory to move, you should copy the coefficient |
---|
| 240 | values to local variables or an array before such operations. |
---|
| 241 | */ |
---|
| 242 | double |
---|
| 243 | TriaxialEllipsoid(double dp[], double q) |
---|
| 244 | { |
---|
| 245 | int i,j; |
---|
| 246 | double Pi; |
---|
| 247 | double scale,aa,bb,cc,delrho,bkg; //local variables of coefficient wave |
---|
| 248 | int nordi=76; //order of integration |
---|
| 249 | int nordj=76; |
---|
| 250 | double va,vb; //upper and lower integration limits |
---|
| 251 | double summ,zi,yyy,answer; //running tally of integration |
---|
[235] | 252 | double summj,vaj,vbj,zij,slde,sld; //for the inner integration |
---|
[97] | 253 | |
---|
| 254 | Pi = 4.0*atan(1.0); |
---|
[453] | 255 | va = 0.0; |
---|
| 256 | vb = 1.0; //orintational average, outer integral |
---|
| 257 | vaj = 0.0; |
---|
| 258 | vbj = 1.0; //endpoints of inner integral |
---|
[97] | 259 | |
---|
| 260 | summ = 0.0; //initialize intergral |
---|
| 261 | |
---|
| 262 | scale = dp[0]; //make local copies in case memory moves |
---|
| 263 | aa = dp[1]; |
---|
| 264 | bb = dp[2]; |
---|
| 265 | cc = dp[3]; |
---|
[235] | 266 | slde = dp[4]; |
---|
| 267 | sld = dp[5]; |
---|
| 268 | delrho = slde - sld; |
---|
| 269 | bkg = dp[6]; |
---|
[97] | 270 | for(i=0;i<nordi;i++) { |
---|
| 271 | //setup inner integral over the ellipsoidal cross-section |
---|
[453] | 272 | summj=0.0; |
---|
[97] | 273 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy |
---|
| 274 | for(j=0;j<nordj;j++) { |
---|
| 275 | //20 gauss points for the inner integral |
---|
| 276 | zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy |
---|
| 277 | yyy = Gauss76Wt[j] * TriaxialKernel(q,aa,bb,cc,zi,zij); |
---|
| 278 | summj += yyy; |
---|
| 279 | } |
---|
| 280 | //now calculate the value of the inner integral |
---|
| 281 | answer = (vbj-vaj)/2.0*summj; |
---|
| 282 | |
---|
| 283 | //now calculate outer integral |
---|
| 284 | yyy = Gauss76Wt[i] * answer; |
---|
| 285 | summ += yyy; |
---|
| 286 | } //final scaling is done at the end of the function, after the NT_FP64 case |
---|
| 287 | |
---|
| 288 | answer = (vb-va)/2.0*summ; |
---|
| 289 | // Multiply by contrast^2 |
---|
| 290 | answer *= delrho*delrho; |
---|
| 291 | //normalize by ellipsoid volume |
---|
[453] | 292 | answer *= 4.0*Pi/3.0*aa*bb*cc; |
---|
[97] | 293 | //convert to [cm-1] |
---|
| 294 | answer *= 1.0e8; |
---|
| 295 | //Scale |
---|
| 296 | answer *= scale; |
---|
| 297 | // add in the background |
---|
| 298 | answer += bkg; |
---|
| 299 | |
---|
| 300 | return answer; |
---|
| 301 | } |
---|
| 302 | |
---|
| 303 | /* ParallelepipedX : calculates the form factor of a Parallelepiped (a rectangular solid) |
---|
| 304 | at the given x-value p->x |
---|
| 305 | |
---|
| 306 | Uses 76 pt Gaussian quadrature for both integrals |
---|
| 307 | |
---|
| 308 | Warning: |
---|
| 309 | The call to WaveData() below returns a pointer to the middle |
---|
| 310 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 311 | calculations could cause memory to move, you should copy the coefficient |
---|
| 312 | values to local variables or an array before such operations. |
---|
| 313 | */ |
---|
| 314 | double |
---|
| 315 | Parallelepiped(double dp[], double q) |
---|
| 316 | { |
---|
| 317 | int i,j; |
---|
| 318 | double scale,aa,bb,cc,delrho,bkg; //local variables of coefficient wave |
---|
| 319 | int nordi=76; //order of integration |
---|
| 320 | int nordj=76; |
---|
| 321 | double va,vb; //upper and lower integration limits |
---|
| 322 | double summ,yyy,answer; //running tally of integration |
---|
| 323 | double summj,vaj,vbj; //for the inner integration |
---|
[235] | 324 | double mu,mudum,arg,sigma,uu,vol,sldp,sld; |
---|
[97] | 325 | |
---|
| 326 | |
---|
| 327 | // Pi = 4.0*atan(1.0); |
---|
[632] | 328 | va = 0.0; |
---|
| 329 | vb = 1.0; //orintational average, outer integral |
---|
| 330 | vaj = 0.0; |
---|
| 331 | vbj = 1.0; //endpoints of inner integral |
---|
| 332 | |
---|
[97] | 333 | summ = 0.0; //initialize intergral |
---|
| 334 | |
---|
| 335 | scale = dp[0]; //make local copies in case memory moves |
---|
| 336 | aa = dp[1]; |
---|
| 337 | bb = dp[2]; |
---|
| 338 | cc = dp[3]; |
---|
[235] | 339 | sldp = dp[4]; |
---|
| 340 | sld = dp[5]; |
---|
| 341 | delrho = sldp - sld; |
---|
| 342 | bkg = dp[6]; |
---|
[97] | 343 | |
---|
| 344 | mu = q*bb; |
---|
| 345 | vol = aa*bb*cc; |
---|
| 346 | // normalize all WRT bb |
---|
| 347 | aa = aa/bb; |
---|
| 348 | cc = cc/bb; |
---|
| 349 | |
---|
| 350 | for(i=0;i<nordi;i++) { |
---|
| 351 | //setup inner integral over the ellipsoidal cross-section |
---|
[632] | 352 | summj=0.0; |
---|
[97] | 353 | sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the outer dummy |
---|
| 354 | |
---|
| 355 | for(j=0;j<nordj;j++) { |
---|
| 356 | //76 gauss points for the inner integral |
---|
| 357 | uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the inner dummy |
---|
[632] | 358 | mudum = mu*sqrt(1.0-sigma*sigma); |
---|
[97] | 359 | yyy = Gauss76Wt[j] * PPKernel(aa,mudum,uu); |
---|
| 360 | summj += yyy; |
---|
| 361 | } |
---|
| 362 | //now calculate the value of the inner integral |
---|
| 363 | answer = (vbj-vaj)/2.0*summj; |
---|
[632] | 364 | |
---|
| 365 | arg = mu*cc*sigma/2.0; |
---|
| 366 | if ( arg == 0.0 ) { |
---|
| 367 | answer *= 1.0; |
---|
[97] | 368 | } else { |
---|
| 369 | answer *= sin(arg)*sin(arg)/arg/arg; |
---|
| 370 | } |
---|
| 371 | |
---|
| 372 | //now sum up the outer integral |
---|
| 373 | yyy = Gauss76Wt[i] * answer; |
---|
| 374 | summ += yyy; |
---|
| 375 | } //final scaling is done at the end of the function, after the NT_FP64 case |
---|
| 376 | |
---|
| 377 | answer = (vb-va)/2.0*summ; |
---|
| 378 | // Multiply by contrast^2 |
---|
| 379 | answer *= delrho*delrho; |
---|
| 380 | //normalize by volume |
---|
| 381 | answer *= vol; |
---|
| 382 | //convert to [cm-1] |
---|
| 383 | answer *= 1.0e8; |
---|
| 384 | //Scale |
---|
| 385 | answer *= scale; |
---|
| 386 | // add in the background |
---|
| 387 | answer += bkg; |
---|
| 388 | |
---|
| 389 | return answer; |
---|
| 390 | } |
---|
| 391 | |
---|
| 392 | /* HollowCylinderX : calculates the form factor of a Hollow Cylinder |
---|
| 393 | at the given x-value p->x |
---|
| 394 | |
---|
| 395 | Uses 76 pt Gaussian quadrature for the single integral |
---|
| 396 | |
---|
| 397 | Warning: |
---|
| 398 | The call to WaveData() below returns a pointer to the middle |
---|
| 399 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 400 | calculations could cause memory to move, you should copy the coefficient |
---|
| 401 | values to local variables or an array before such operations. |
---|
| 402 | */ |
---|
| 403 | double |
---|
| 404 | HollowCylinder(double dp[], double q) |
---|
| 405 | { |
---|
| 406 | int i; |
---|
| 407 | double scale,rcore,rshell,length,delrho,bkg; //local variables of coefficient wave |
---|
| 408 | int nord=76; //order of integration |
---|
| 409 | double va,vb,zi; //upper and lower integration limits |
---|
[235] | 410 | double summ,answer,pi,sldc,sld; //running tally of integration |
---|
[97] | 411 | |
---|
| 412 | pi = 4.0*atan(1.0); |
---|
[632] | 413 | va = 0.0; |
---|
| 414 | vb = 1.0; //limits of numerical integral |
---|
| 415 | |
---|
[97] | 416 | summ = 0.0; //initialize intergral |
---|
| 417 | |
---|
| 418 | scale = dp[0]; //make local copies in case memory moves |
---|
| 419 | rcore = dp[1]; |
---|
| 420 | rshell = dp[2]; |
---|
| 421 | length = dp[3]; |
---|
[235] | 422 | sldc = dp[4]; |
---|
| 423 | sld = dp[5]; |
---|
| 424 | delrho = sldc - sld; |
---|
| 425 | bkg = dp[6]; |
---|
[97] | 426 | |
---|
| 427 | for(i=0;i<nord;i++) { |
---|
| 428 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; |
---|
| 429 | summ += Gauss76Wt[i] * HolCylKernel(q, rcore, rshell, length, zi); |
---|
| 430 | } |
---|
| 431 | |
---|
| 432 | answer = (vb-va)/2.0*summ; |
---|
| 433 | // Multiply by contrast^2 |
---|
| 434 | answer *= delrho*delrho; |
---|
| 435 | //normalize by volume |
---|
| 436 | answer *= pi*(rshell*rshell-rcore*rcore)*length; |
---|
| 437 | //convert to [cm-1] |
---|
| 438 | answer *= 1.0e8; |
---|
| 439 | //Scale |
---|
| 440 | answer *= scale; |
---|
| 441 | // add in the background |
---|
| 442 | answer += bkg; |
---|
| 443 | |
---|
| 444 | return answer; |
---|
| 445 | } |
---|
| 446 | |
---|
| 447 | /* EllipsoidFormX : calculates the form factor of an ellipsoid of revolution with semiaxes a:a:nua |
---|
| 448 | at the given x-value p->x |
---|
| 449 | |
---|
| 450 | Uses 76 pt Gaussian quadrature for the single integral |
---|
| 451 | |
---|
| 452 | Warning: |
---|
| 453 | The call to WaveData() below returns a pointer to the middle |
---|
| 454 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 455 | calculations could cause memory to move, you should copy the coefficient |
---|
| 456 | values to local variables or an array before such operations. |
---|
| 457 | */ |
---|
| 458 | double |
---|
| 459 | EllipsoidForm(double dp[], double q) |
---|
| 460 | { |
---|
| 461 | int i; |
---|
| 462 | double scale,a,nua,delrho,bkg; //local variables of coefficient wave |
---|
| 463 | int nord=76; //order of integration |
---|
| 464 | double va,vb,zi; //upper and lower integration limits |
---|
[235] | 465 | double summ,answer,pi,slde,sld; //running tally of integration |
---|
[97] | 466 | |
---|
| 467 | pi = 4.0*atan(1.0); |
---|
[632] | 468 | va = 0.0; |
---|
| 469 | vb = 1.0; //limits of numerical integral |
---|
| 470 | |
---|
[97] | 471 | summ = 0.0; //initialize intergral |
---|
| 472 | |
---|
| 473 | scale = dp[0]; //make local copies in case memory moves |
---|
| 474 | nua = dp[1]; |
---|
| 475 | a = dp[2]; |
---|
[235] | 476 | slde = dp[3]; |
---|
| 477 | sld = dp[4]; |
---|
| 478 | delrho = slde - sld; |
---|
| 479 | bkg = dp[5]; |
---|
[97] | 480 | |
---|
| 481 | for(i=0;i<nord;i++) { |
---|
| 482 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; |
---|
| 483 | summ += Gauss76Wt[i] * EllipsoidKernel(q, a, nua, zi); |
---|
| 484 | } |
---|
| 485 | |
---|
| 486 | answer = (vb-va)/2.0*summ; |
---|
| 487 | // Multiply by contrast^2 |
---|
| 488 | answer *= delrho*delrho; |
---|
| 489 | //normalize by volume |
---|
[632] | 490 | answer *= 4.0*pi/3.0*a*a*nua; |
---|
[97] | 491 | //convert to [cm-1] |
---|
| 492 | answer *= 1.0e8; |
---|
| 493 | //Scale |
---|
| 494 | answer *= scale; |
---|
| 495 | // add in the background |
---|
| 496 | answer += bkg; |
---|
| 497 | |
---|
| 498 | return answer; |
---|
| 499 | } |
---|
| 500 | |
---|
| 501 | |
---|
| 502 | /* Cyl_PolyRadiusX : calculates the form factor of a cylinder at the given x-value p->x |
---|
| 503 | the cylinder has a polydisperse cross section |
---|
| 504 | |
---|
| 505 | */ |
---|
| 506 | double |
---|
| 507 | Cyl_PolyRadius(double dp[], double q) |
---|
| 508 | { |
---|
| 509 | int i; |
---|
| 510 | double scale,radius,length,pd,delrho,bkg; //local variables of coefficient wave |
---|
| 511 | int nord=20; //order of integration |
---|
| 512 | double uplim,lolim; //upper and lower integration limits |
---|
| 513 | double summ,zi,yyy,answer,Vpoly; //running tally of integration |
---|
[235] | 514 | double range,zz,Pi,sldc,sld; |
---|
[97] | 515 | |
---|
| 516 | Pi = 4.0*atan(1.0); |
---|
| 517 | range = 3.4; |
---|
| 518 | |
---|
| 519 | summ = 0.0; //initialize intergral |
---|
| 520 | |
---|
| 521 | scale = dp[0]; //make local copies in case memory moves |
---|
| 522 | radius = dp[1]; |
---|
| 523 | length = dp[2]; |
---|
| 524 | pd = dp[3]; |
---|
[235] | 525 | sldc = dp[4]; |
---|
| 526 | sld = dp[5]; |
---|
| 527 | delrho = sldc - sld; |
---|
| 528 | bkg = dp[6]; |
---|
[97] | 529 | |
---|
| 530 | zz = (1.0/pd)*(1.0/pd) - 1.0; |
---|
| 531 | |
---|
| 532 | lolim = radius*(1.0-range*pd); //set the upper/lower limits to cover the distribution |
---|
[632] | 533 | if(lolim<0.0) { |
---|
| 534 | lolim = 0.0; |
---|
[97] | 535 | } |
---|
| 536 | if(pd>0.3) { |
---|
| 537 | range = 3.4 + (pd-0.3)*18.0; |
---|
| 538 | } |
---|
| 539 | uplim = radius*(1.0+range*pd); |
---|
| 540 | |
---|
| 541 | for(i=0;i<nord;i++) { |
---|
| 542 | zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 543 | yyy = Gauss20Wt[i] * Cyl_PolyRadKernel(q, radius, length, zz, delrho, zi); |
---|
| 544 | summ += yyy; |
---|
| 545 | } |
---|
| 546 | |
---|
| 547 | answer = (uplim-lolim)/2.0*summ; |
---|
| 548 | //normalize by average cylinder volume |
---|
| 549 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
| 550 | Vpoly=Pi*radius*radius*length*(zz+2.0)/(zz+1.0); |
---|
| 551 | answer /= Vpoly; |
---|
| 552 | //convert to [cm-1] |
---|
| 553 | answer *= 1.0e8; |
---|
| 554 | //Scale |
---|
| 555 | answer *= scale; |
---|
| 556 | // add in the background |
---|
| 557 | answer += bkg; |
---|
| 558 | |
---|
| 559 | return answer; |
---|
| 560 | } |
---|
| 561 | |
---|
| 562 | /* Cyl_PolyLengthX : calculates the form factor of a cylinder at the given x-value p->x |
---|
| 563 | the cylinder has a polydisperse Length |
---|
| 564 | |
---|
| 565 | */ |
---|
| 566 | double |
---|
| 567 | Cyl_PolyLength(double dp[], double q) |
---|
| 568 | { |
---|
| 569 | int i; |
---|
| 570 | double scale,radius,length,pd,delrho,bkg; //local variables of coefficient wave |
---|
| 571 | int nord=20; //order of integration |
---|
| 572 | double uplim,lolim; //upper and lower integration limits |
---|
| 573 | double summ,zi,yyy,answer,Vpoly; //running tally of integration |
---|
[235] | 574 | double range,zz,Pi,sldc,sld; |
---|
[97] | 575 | |
---|
| 576 | |
---|
| 577 | Pi = 4.0*atan(1.0); |
---|
| 578 | range = 3.4; |
---|
| 579 | |
---|
| 580 | summ = 0.0; //initialize intergral |
---|
| 581 | |
---|
| 582 | scale = dp[0]; //make local copies in case memory moves |
---|
| 583 | radius = dp[1]; |
---|
| 584 | length = dp[2]; |
---|
| 585 | pd = dp[3]; |
---|
[235] | 586 | sldc = dp[4]; |
---|
| 587 | sld = dp[5]; |
---|
| 588 | delrho = sldc - sld; |
---|
| 589 | bkg = dp[6]; |
---|
[97] | 590 | |
---|
| 591 | zz = (1.0/pd)*(1.0/pd) - 1.0; |
---|
| 592 | |
---|
| 593 | lolim = length*(1.0-range*pd); //set the upper/lower limits to cover the distribution |
---|
[632] | 594 | if(lolim<0.0) { |
---|
| 595 | lolim = 0.0; |
---|
[97] | 596 | } |
---|
| 597 | if(pd>0.3) { |
---|
| 598 | range = 3.4 + (pd-0.3)*18.0; |
---|
| 599 | } |
---|
| 600 | uplim = length*(1.0+range*pd); |
---|
| 601 | |
---|
| 602 | for(i=0;i<nord;i++) { |
---|
| 603 | zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 604 | yyy = Gauss20Wt[i] * Cyl_PolyLenKernel(q, radius, length, zz, delrho, zi); |
---|
| 605 | summ += yyy; |
---|
| 606 | } |
---|
| 607 | |
---|
| 608 | answer = (uplim-lolim)/2.0*summ; |
---|
| 609 | //normalize by average cylinder volume (first moment) |
---|
| 610 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
| 611 | Vpoly=Pi*radius*radius*length; |
---|
| 612 | answer /= Vpoly; |
---|
| 613 | //convert to [cm-1] |
---|
| 614 | answer *= 1.0e8; |
---|
| 615 | //Scale |
---|
| 616 | answer *= scale; |
---|
| 617 | // add in the background |
---|
| 618 | answer += bkg; |
---|
| 619 | |
---|
| 620 | return answer; |
---|
| 621 | } |
---|
| 622 | |
---|
| 623 | /* CoreShellCylinderX : calculates the form factor of a cylinder at the given x-value p->x |
---|
| 624 | the cylinder has a core-shell structure |
---|
| 625 | |
---|
| 626 | */ |
---|
| 627 | double |
---|
| 628 | CoreShellCylinder(double dp[], double q) |
---|
| 629 | { |
---|
| 630 | int i; |
---|
| 631 | double scale,rcore,length,bkg; //local variables of coefficient wave |
---|
| 632 | double thick,rhoc,rhos,rhosolv; |
---|
| 633 | int nord=76; //order of integration |
---|
| 634 | double uplim,lolim,halfheight; //upper and lower integration limits |
---|
| 635 | double summ,zi,yyy,answer,Vcyl; //running tally of integration |
---|
| 636 | double Pi; |
---|
| 637 | |
---|
| 638 | Pi = 4.0*atan(1.0); |
---|
| 639 | |
---|
| 640 | lolim = 0.0; |
---|
| 641 | uplim = Pi/2.0; |
---|
| 642 | |
---|
| 643 | summ = 0.0; //initialize intergral |
---|
| 644 | |
---|
| 645 | scale = dp[0]; //make local copies in case memory moves |
---|
| 646 | rcore = dp[1]; |
---|
| 647 | thick = dp[2]; |
---|
| 648 | length = dp[3]; |
---|
| 649 | rhoc = dp[4]; |
---|
| 650 | rhos = dp[5]; |
---|
| 651 | rhosolv = dp[6]; |
---|
| 652 | bkg = dp[7]; |
---|
| 653 | |
---|
| 654 | halfheight = length/2.0; |
---|
| 655 | |
---|
| 656 | for(i=0;i<nord;i++) { |
---|
| 657 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 658 | yyy = Gauss76Wt[i] * CoreShellCylKernel(q, rcore, thick, rhoc,rhos,rhosolv, halfheight, zi); |
---|
| 659 | summ += yyy; |
---|
| 660 | } |
---|
| 661 | |
---|
| 662 | answer = (uplim-lolim)/2.0*summ; |
---|
| 663 | // length is the total core length |
---|
| 664 | Vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2.0*thick); |
---|
| 665 | answer /= Vcyl; |
---|
| 666 | //convert to [cm-1] |
---|
| 667 | answer *= 1.0e8; |
---|
| 668 | //Scale |
---|
| 669 | answer *= scale; |
---|
| 670 | // add in the background |
---|
| 671 | answer += bkg; |
---|
| 672 | |
---|
| 673 | return answer; |
---|
| 674 | } |
---|
| 675 | |
---|
| 676 | |
---|
| 677 | /* PolyCoShCylinderX : calculates the form factor of a core-shell cylinder at the given x-value p->x |
---|
| 678 | the cylinder has a polydisperse CORE radius |
---|
| 679 | |
---|
| 680 | */ |
---|
| 681 | double |
---|
| 682 | PolyCoShCylinder(double dp[], double q) |
---|
| 683 | { |
---|
| 684 | int i; |
---|
| 685 | double scale,radius,length,sigma,bkg; //local variables of coefficient wave |
---|
| 686 | double rad,radthick,facthick,rhoc,rhos,rhosolv; |
---|
| 687 | int nord=20; //order of integration |
---|
| 688 | double uplim,lolim; //upper and lower integration limits |
---|
| 689 | double summ,yyy,answer,Vpoly; //running tally of integration |
---|
| 690 | double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr; |
---|
[501] | 691 | |
---|
[97] | 692 | Pi = 4.0*atan(1.0); |
---|
| 693 | |
---|
| 694 | summ = 0.0; //initialize intergral |
---|
| 695 | Rsqrsumm = 0.0; |
---|
| 696 | |
---|
| 697 | scale = dp[0]; |
---|
| 698 | radius = dp[1]; |
---|
| 699 | sigma = dp[2]; //sigma is the standard mean deviation |
---|
| 700 | length = dp[3]; |
---|
| 701 | radthick = dp[4]; |
---|
| 702 | facthick= dp[5]; |
---|
| 703 | rhoc = dp[6]; |
---|
| 704 | rhos = dp[7]; |
---|
| 705 | rhosolv = dp[8]; |
---|
| 706 | bkg = dp[9]; |
---|
| 707 | |
---|
| 708 | lolim = exp(log(radius)-(4.*sigma)); |
---|
[632] | 709 | if (lolim<0.0) { |
---|
| 710 | lolim=0.0; //to avoid numerical error when va<0 (-ve r value) |
---|
[97] | 711 | } |
---|
| 712 | uplim = exp(log(radius)+(4.*sigma)); |
---|
| 713 | |
---|
| 714 | for(i=0;i<nord;i++) { |
---|
| 715 | rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 716 | AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp(-(0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma))); |
---|
| 717 | yyy = AR* Gauss20Wt[i] * CSCylIntegration(q,rad,radthick,facthick,rhoc,rhos,rhosolv,length); |
---|
| 718 | Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick); //SRK normalize to total dimensions |
---|
| 719 | summ += yyy; |
---|
| 720 | Rsqrsumm += Rsqryyy; |
---|
| 721 | } |
---|
| 722 | |
---|
| 723 | answer = (uplim-lolim)/2.0*summ; |
---|
| 724 | Rsqr = (uplim-lolim)/2.0*Rsqrsumm; |
---|
| 725 | //normalize by average cylinder volume |
---|
| 726 | Vpoly = Pi*Rsqr*(length+2*facthick); |
---|
| 727 | answer /= Vpoly; |
---|
| 728 | //convert to [cm-1] |
---|
| 729 | answer *= 1.0e8; |
---|
| 730 | //Scale |
---|
| 731 | answer *= scale; |
---|
| 732 | // add in the background |
---|
| 733 | answer += bkg; |
---|
| 734 | |
---|
| 735 | return answer; |
---|
| 736 | } |
---|
| 737 | |
---|
| 738 | /* OblateFormX : calculates the form factor of a core-shell Oblate ellipsoid at the given x-value p->x |
---|
| 739 | the ellipsoid has a core-shell structure |
---|
| 740 | |
---|
| 741 | */ |
---|
| 742 | double |
---|
| 743 | OblateForm(double dp[], double q) |
---|
| 744 | { |
---|
| 745 | int i; |
---|
| 746 | double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg; |
---|
| 747 | int nord=76; //order of integration |
---|
| 748 | double uplim,lolim; //upper and lower integration limits |
---|
| 749 | double summ,zi,yyy,answer,oblatevol; //running tally of integration |
---|
[235] | 750 | double Pi,sldc,slds,sld; |
---|
[97] | 751 | |
---|
| 752 | Pi = 4.0*atan(1.0); |
---|
| 753 | |
---|
| 754 | lolim = 0.0; |
---|
| 755 | uplim = 1.0; |
---|
| 756 | |
---|
| 757 | summ = 0.0; //initialize intergral |
---|
| 758 | |
---|
| 759 | |
---|
| 760 | scale = dp[0]; //make local copies in case memory moves |
---|
| 761 | crmaj = dp[1]; |
---|
| 762 | crmin = dp[2]; |
---|
| 763 | trmaj = dp[3]; |
---|
| 764 | trmin = dp[4]; |
---|
[235] | 765 | sldc = dp[5]; |
---|
| 766 | slds = dp[6]; |
---|
| 767 | sld = dp[7]; |
---|
| 768 | delpc = sldc - slds; //core - shell |
---|
| 769 | delps = slds - sld; //shell - solvent |
---|
[632] | 770 | bkg = dp[8]; |
---|
| 771 | |
---|
[97] | 772 | for(i=0;i<nord;i++) { |
---|
| 773 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 774 | yyy = Gauss76Wt[i] * gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q); |
---|
| 775 | summ += yyy; |
---|
| 776 | } |
---|
| 777 | |
---|
| 778 | answer = (uplim-lolim)/2.0*summ; |
---|
| 779 | // normalize by particle volume |
---|
[632] | 780 | oblatevol = 4.0*Pi/3.0*trmaj*trmaj*trmin; |
---|
[97] | 781 | answer /= oblatevol; |
---|
| 782 | |
---|
| 783 | //convert to [cm-1] |
---|
| 784 | answer *= 1.0e8; |
---|
| 785 | //Scale |
---|
| 786 | answer *= scale; |
---|
| 787 | // add in the background |
---|
| 788 | answer += bkg; |
---|
| 789 | |
---|
| 790 | return answer; |
---|
| 791 | } |
---|
| 792 | |
---|
| 793 | /* ProlateFormX : calculates the form factor of a core-shell Prolate ellipsoid at the given x-value p->x |
---|
| 794 | the ellipsoid has a core-shell structure |
---|
| 795 | |
---|
| 796 | */ |
---|
| 797 | double |
---|
| 798 | ProlateForm(double dp[], double q) |
---|
| 799 | { |
---|
| 800 | int i; |
---|
| 801 | double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg; |
---|
| 802 | int nord=76; //order of integration |
---|
| 803 | double uplim,lolim; //upper and lower integration limits |
---|
| 804 | double summ,zi,yyy,answer,prolatevol; //running tally of integration |
---|
[235] | 805 | double Pi,sldc,slds,sld; |
---|
[97] | 806 | |
---|
| 807 | Pi = 4.0*atan(1.0); |
---|
| 808 | |
---|
| 809 | lolim = 0.0; |
---|
| 810 | uplim = 1.0; |
---|
| 811 | |
---|
| 812 | summ = 0.0; //initialize intergral |
---|
| 813 | |
---|
| 814 | scale = dp[0]; //make local copies in case memory moves |
---|
| 815 | crmaj = dp[1]; |
---|
| 816 | crmin = dp[2]; |
---|
| 817 | trmaj = dp[3]; |
---|
| 818 | trmin = dp[4]; |
---|
[235] | 819 | sldc = dp[5]; |
---|
| 820 | slds = dp[6]; |
---|
| 821 | sld = dp[7]; |
---|
| 822 | delpc = sldc - slds; //core - shell |
---|
| 823 | delps = slds - sld; //shell - sovent |
---|
[632] | 824 | bkg = dp[8]; |
---|
| 825 | |
---|
[97] | 826 | for(i=0;i<nord;i++) { |
---|
| 827 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 828 | yyy = Gauss76Wt[i] * gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q); |
---|
| 829 | summ += yyy; |
---|
| 830 | } |
---|
| 831 | |
---|
| 832 | answer = (uplim-lolim)/2.0*summ; |
---|
| 833 | // normalize by particle volume |
---|
[632] | 834 | prolatevol = 4.0*Pi/3.0*trmaj*trmin*trmin; |
---|
[97] | 835 | answer /= prolatevol; |
---|
| 836 | |
---|
| 837 | //convert to [cm-1] |
---|
| 838 | answer *= 1.0e8; |
---|
| 839 | //Scale |
---|
| 840 | answer *= scale; |
---|
| 841 | // add in the background |
---|
| 842 | answer += bkg; |
---|
| 843 | |
---|
| 844 | return answer; |
---|
| 845 | } |
---|
| 846 | |
---|
| 847 | |
---|
| 848 | /* StackedDiscsX : calculates the form factor of a stacked "tactoid" of core shell disks |
---|
| 849 | like clay platelets that are not exfoliated |
---|
| 850 | |
---|
| 851 | */ |
---|
| 852 | double |
---|
| 853 | StackedDiscs(double dp[], double q) |
---|
| 854 | { |
---|
| 855 | int i; |
---|
| 856 | double scale,length,bkg,rcore,thick,rhoc,rhol,rhosolv,N,gsd; //local variables of coefficient wave |
---|
| 857 | double va,vb,vcyl,summ,yyy,zi,halfheight,d,answer; |
---|
| 858 | int nord=76; //order of integration |
---|
| 859 | double Pi; |
---|
| 860 | |
---|
| 861 | |
---|
| 862 | Pi = 4.0*atan(1.0); |
---|
| 863 | |
---|
| 864 | va = 0.0; |
---|
| 865 | vb = Pi/2.0; |
---|
| 866 | |
---|
| 867 | summ = 0.0; //initialize intergral |
---|
| 868 | |
---|
| 869 | scale = dp[0]; |
---|
| 870 | rcore = dp[1]; |
---|
| 871 | length = dp[2]; |
---|
| 872 | thick = dp[3]; |
---|
| 873 | rhoc = dp[4]; |
---|
| 874 | rhol = dp[5]; |
---|
| 875 | rhosolv = dp[6]; |
---|
| 876 | N = dp[7]; |
---|
| 877 | gsd = dp[8]; |
---|
| 878 | bkg = dp[9]; |
---|
| 879 | |
---|
| 880 | d=2.0*thick+length; |
---|
| 881 | halfheight = length/2.0; |
---|
| 882 | |
---|
| 883 | for(i=0;i<nord;i++) { |
---|
| 884 | zi = ( Gauss76Z[i]*(vb-va) + vb + va )/2.0; |
---|
| 885 | yyy = Gauss76Wt[i] * Stackdisc_kern(q, rcore, rhoc,rhol,rhosolv, halfheight,thick,zi,gsd,d,N); |
---|
| 886 | summ += yyy; |
---|
| 887 | } |
---|
| 888 | |
---|
| 889 | answer = (vb-va)/2.0*summ; |
---|
| 890 | // length is the total core length |
---|
| 891 | vcyl=Pi*rcore*rcore*(2.0*thick+length)*N; |
---|
| 892 | answer /= vcyl; |
---|
| 893 | //Convert to [cm-1] |
---|
| 894 | answer *= 1.0e8; |
---|
| 895 | //Scale |
---|
| 896 | answer *= scale; |
---|
| 897 | // add in the background |
---|
| 898 | answer += bkg; |
---|
| 899 | |
---|
| 900 | return answer; |
---|
| 901 | } |
---|
| 902 | |
---|
| 903 | |
---|
| 904 | /* LamellarFFX : calculates the form factor of a lamellar structure - no S(q) effects included |
---|
| 905 | |
---|
| 906 | */ |
---|
| 907 | double |
---|
| 908 | LamellarFF(double dp[], double q) |
---|
| 909 | { |
---|
| 910 | double scale,del,sig,contr,bkg; //local variables of coefficient wave |
---|
| 911 | double inten, qval,Pq; |
---|
[235] | 912 | double Pi,sldb,sld; |
---|
[97] | 913 | |
---|
| 914 | |
---|
| 915 | Pi = 4.0*atan(1.0); |
---|
| 916 | scale = dp[0]; |
---|
| 917 | del = dp[1]; |
---|
| 918 | sig = dp[2]*del; |
---|
[235] | 919 | sldb = dp[3]; |
---|
| 920 | sld = dp[4]; |
---|
| 921 | contr = sldb - sld; |
---|
| 922 | bkg = dp[5]; |
---|
[156] | 923 | qval=q; |
---|
[97] | 924 | |
---|
| 925 | Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig)); |
---|
| 926 | |
---|
| 927 | inten = 2.0*Pi*scale*Pq/(qval*qval); //this is now dimensionless... |
---|
| 928 | |
---|
| 929 | inten /= del; //normalize by the thickness (in A) |
---|
| 930 | |
---|
| 931 | inten *= 1.0e8; // 1/A to 1/cm |
---|
| 932 | |
---|
| 933 | return(inten+bkg); |
---|
| 934 | } |
---|
| 935 | |
---|
| 936 | /* LamellarPSX : calculates the form factor of a lamellar structure - with S(q) effects included |
---|
[356] | 937 | --- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the |
---|
| 938 | model is smeared just like any other function |
---|
| 939 | */ |
---|
[97] | 940 | double |
---|
| 941 | LamellarPS(double dp[], double q) |
---|
| 942 | { |
---|
| 943 | double scale,dd,del,sig,contr,NN,Cp,bkg; //local variables of coefficient wave |
---|
[235] | 944 | double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ; |
---|
| 945 | double Pi,Euler,dQDefault,fii,sldb,sld; |
---|
[97] | 946 | int ii,NNint; |
---|
[235] | 947 | // char buf[256]; |
---|
| 948 | |
---|
[97] | 949 | |
---|
| 950 | Euler = 0.5772156649; // Euler's constant |
---|
[356] | 951 | // dQDefault = 0.0025; //[=] 1/A, q-resolution, default value |
---|
| 952 | dQDefault = 0.0; |
---|
[97] | 953 | dQ = dQDefault; |
---|
| 954 | |
---|
| 955 | Pi = 4.0*atan(1.0); |
---|
| 956 | qval = q; |
---|
| 957 | |
---|
| 958 | scale = dp[0]; |
---|
| 959 | dd = dp[1]; |
---|
| 960 | del = dp[2]; |
---|
| 961 | sig = dp[3]*del; |
---|
[235] | 962 | sldb = dp[4]; |
---|
| 963 | sld = dp[5]; |
---|
| 964 | contr = sldb - sld; |
---|
| 965 | NN = trunc(dp[6]); //be sure that NN is an integer |
---|
| 966 | Cp = dp[7]; |
---|
| 967 | bkg = dp[8]; |
---|
[97] | 968 | |
---|
| 969 | Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig)); |
---|
| 970 | |
---|
| 971 | NNint = (int)NN; //cast to an integer for the loop |
---|
[235] | 972 | |
---|
| 973 | // sprintf(buf, "qval = %g\r", qval); |
---|
| 974 | // XOPNotice(buf); |
---|
| 975 | |
---|
[97] | 976 | ii=0; |
---|
| 977 | Sq = 0.0; |
---|
| 978 | for(ii=1;ii<(NNint-1);ii+=1) { |
---|
| 979 | |
---|
| 980 | fii = (double)ii; //do I really need to do this? |
---|
| 981 | |
---|
| 982 | temp = 0.0; |
---|
[235] | 983 | alpha = Cp/4.0/Pi/Pi*(log(Pi*fii) + Euler); |
---|
[97] | 984 | t1 = 2.0*dQ*dQ*dd*dd*alpha; |
---|
| 985 | t2 = 2.0*qval*qval*dd*dd*alpha; |
---|
[235] | 986 | t3 = dQ*dQ*dd*dd*fii*fii; |
---|
[97] | 987 | |
---|
[235] | 988 | temp = 1.0-fii/NN; |
---|
| 989 | temp *= cos(dd*qval*fii/(1.0+t1)); |
---|
[97] | 990 | temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); |
---|
| 991 | temp /= sqrt(1.0+t1); |
---|
| 992 | |
---|
| 993 | Sq += temp; |
---|
| 994 | } |
---|
| 995 | |
---|
| 996 | Sq *= 2.0; |
---|
| 997 | Sq += 1.0; |
---|
| 998 | |
---|
| 999 | inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval); |
---|
| 1000 | |
---|
| 1001 | inten *= 1.0e8; // 1/A to 1/cm |
---|
| 1002 | |
---|
| 1003 | return(inten+bkg); |
---|
| 1004 | } |
---|
| 1005 | |
---|
| 1006 | |
---|
| 1007 | /* LamellarPS_HGX : calculates the form factor of a lamellar structure - with S(q) effects included |
---|
[356] | 1008 | --- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the |
---|
| 1009 | model is smeared just like any other function |
---|
| 1010 | */ |
---|
[97] | 1011 | double |
---|
| 1012 | LamellarPS_HG(double dp[], double q) |
---|
| 1013 | { |
---|
| 1014 | double scale,dd,delT,delH,SLD_T,SLD_H,SLD_S,NN,Cp,bkg; //local variables of coefficient wave |
---|
| 1015 | double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ,drh,drt; |
---|
| 1016 | double Pi,Euler,dQDefault,fii; |
---|
| 1017 | int ii,NNint; |
---|
| 1018 | |
---|
| 1019 | |
---|
| 1020 | Euler = 0.5772156649; // Euler's constant |
---|
[356] | 1021 | // dQDefault = 0.0025; //[=] 1/A, q-resolution, default value |
---|
| 1022 | dQDefault = 0.0; |
---|
[97] | 1023 | dQ = dQDefault; |
---|
| 1024 | |
---|
| 1025 | Pi = 4.0*atan(1.0); |
---|
| 1026 | qval= q; |
---|
| 1027 | |
---|
| 1028 | scale = dp[0]; |
---|
| 1029 | dd = dp[1]; |
---|
| 1030 | delT = dp[2]; |
---|
| 1031 | delH = dp[3]; |
---|
| 1032 | SLD_T = dp[4]; |
---|
| 1033 | SLD_H = dp[5]; |
---|
| 1034 | SLD_S = dp[6]; |
---|
| 1035 | NN = trunc(dp[7]); //be sure that NN is an integer |
---|
| 1036 | Cp = dp[8]; |
---|
| 1037 | bkg = dp[9]; |
---|
| 1038 | |
---|
| 1039 | |
---|
| 1040 | drh = SLD_H - SLD_S; |
---|
| 1041 | drt = SLD_T - SLD_S; //correction 13FEB06 by L.Porcar |
---|
| 1042 | |
---|
| 1043 | Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT); |
---|
| 1044 | Pq *= Pq; |
---|
| 1045 | Pq *= 4.0/(qval*qval); |
---|
| 1046 | |
---|
| 1047 | NNint = (int)NN; //cast to an integer for the loop |
---|
| 1048 | ii=0; |
---|
| 1049 | Sq = 0.0; |
---|
| 1050 | for(ii=1;ii<(NNint-1);ii+=1) { |
---|
| 1051 | |
---|
| 1052 | fii = (double)ii; //do I really need to do this? |
---|
| 1053 | |
---|
| 1054 | temp = 0.0; |
---|
| 1055 | alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler); |
---|
| 1056 | t1 = 2.0*dQ*dQ*dd*dd*alpha; |
---|
| 1057 | t2 = 2.0*qval*qval*dd*dd*alpha; |
---|
| 1058 | t3 = dQ*dQ*dd*dd*ii*ii; |
---|
| 1059 | |
---|
| 1060 | temp = 1.0-ii/NN; |
---|
| 1061 | temp *= cos(dd*qval*ii/(1.0+t1)); |
---|
| 1062 | temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); |
---|
| 1063 | temp /= sqrt(1.0+t1); |
---|
| 1064 | |
---|
| 1065 | Sq += temp; |
---|
| 1066 | } |
---|
| 1067 | |
---|
| 1068 | Sq *= 2.0; |
---|
| 1069 | Sq += 1.0; |
---|
| 1070 | |
---|
| 1071 | inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval); |
---|
| 1072 | |
---|
| 1073 | inten *= 1.0e8; // 1/A to 1/cm |
---|
| 1074 | |
---|
| 1075 | return(inten+bkg); |
---|
| 1076 | } |
---|
| 1077 | |
---|
| 1078 | /* LamellarFF_HGX : calculates the form factor of a lamellar structure - no S(q) effects included |
---|
| 1079 | but extra SLD for head groups is included |
---|
| 1080 | |
---|
| 1081 | */ |
---|
| 1082 | double |
---|
| 1083 | LamellarFF_HG(double dp[], double q) |
---|
| 1084 | { |
---|
| 1085 | double scale,delT,delH,slds,sldh,sldt,bkg; //local variables of coefficient wave |
---|
| 1086 | double inten, qval,Pq,drh,drt; |
---|
| 1087 | double Pi; |
---|
| 1088 | |
---|
| 1089 | |
---|
| 1090 | Pi = 4.0*atan(1.0); |
---|
| 1091 | qval= q; |
---|
| 1092 | scale = dp[0]; |
---|
| 1093 | delT = dp[1]; |
---|
| 1094 | delH = dp[2]; |
---|
| 1095 | sldt = dp[3]; |
---|
| 1096 | sldh = dp[4]; |
---|
| 1097 | slds = dp[5]; |
---|
| 1098 | bkg = dp[6]; |
---|
| 1099 | |
---|
| 1100 | |
---|
| 1101 | drh = sldh - slds; |
---|
| 1102 | drt = sldt - slds; //correction 13FEB06 by L.Porcar |
---|
| 1103 | |
---|
| 1104 | Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT); |
---|
| 1105 | Pq *= Pq; |
---|
| 1106 | Pq *= 4.0/(qval*qval); |
---|
| 1107 | |
---|
| 1108 | inten = 2.0*Pi*scale*Pq/(qval*qval); //dimensionless... |
---|
| 1109 | |
---|
| 1110 | inten /= 2.0*(delT+delH); //normalize by the bilayer thickness |
---|
| 1111 | |
---|
| 1112 | inten *= 1.0e8; // 1/A to 1/cm |
---|
| 1113 | |
---|
| 1114 | return(inten+bkg); |
---|
| 1115 | } |
---|
| 1116 | |
---|
| 1117 | /* FlexExclVolCylX : calculates the form factor of a flexible cylinder with a circular cross section |
---|
| 1118 | -- incorporates Wei-Ren Chen's fixes - 2006 |
---|
| 1119 | |
---|
| 1120 | */ |
---|
| 1121 | double |
---|
| 1122 | FlexExclVolCyl(double dp[], double q) |
---|
| 1123 | { |
---|
[235] | 1124 | double scale,L,B,bkg,rad,qr,cont,sldc,slds; |
---|
[97] | 1125 | double Pi,flex,crossSect,answer; |
---|
| 1126 | |
---|
| 1127 | |
---|
| 1128 | Pi = 4.0*atan(1.0); |
---|
| 1129 | |
---|
| 1130 | scale = dp[0]; //make local copies in case memory moves |
---|
| 1131 | L = dp[1]; |
---|
| 1132 | B = dp[2]; |
---|
| 1133 | rad = dp[3]; |
---|
[235] | 1134 | sldc = dp[4]; |
---|
| 1135 | slds = dp[5]; |
---|
| 1136 | cont = sldc-slds; |
---|
| 1137 | bkg = dp[6]; |
---|
[97] | 1138 | |
---|
| 1139 | |
---|
| 1140 | qr = q*rad; |
---|
| 1141 | |
---|
| 1142 | flex = Sk_WR(q,L,B); |
---|
| 1143 | |
---|
| 1144 | crossSect = (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); |
---|
| 1145 | flex *= crossSect; |
---|
| 1146 | flex *= Pi*rad*rad*L; |
---|
| 1147 | flex *= cont*cont; |
---|
| 1148 | flex *= 1.0e8; |
---|
| 1149 | answer = scale*flex + bkg; |
---|
| 1150 | |
---|
| 1151 | return answer; |
---|
| 1152 | } |
---|
| 1153 | |
---|
| 1154 | /* FlexCyl_EllipX : calculates the form factor of a flexible cylinder with an elliptical cross section |
---|
| 1155 | -- incorporates Wei-Ren Chen's fixes - 2006 |
---|
| 1156 | |
---|
| 1157 | */ |
---|
| 1158 | double |
---|
| 1159 | FlexCyl_Ellip(double dp[], double q) |
---|
| 1160 | { |
---|
[235] | 1161 | double scale,L,B,bkg,rad,qr,cont,ellRatio,slds,sldc; |
---|
[97] | 1162 | double Pi,flex,crossSect,answer; |
---|
| 1163 | |
---|
| 1164 | |
---|
| 1165 | Pi = 4.0*atan(1.0); |
---|
| 1166 | scale = dp[0]; //make local copies in case memory moves |
---|
| 1167 | L = dp[1]; |
---|
| 1168 | B = dp[2]; |
---|
| 1169 | rad = dp[3]; |
---|
| 1170 | ellRatio = dp[4]; |
---|
[235] | 1171 | sldc = dp[5]; |
---|
| 1172 | slds = dp[6]; |
---|
| 1173 | bkg = dp[7]; |
---|
[97] | 1174 | |
---|
[235] | 1175 | cont = sldc - slds; |
---|
[97] | 1176 | qr = q*rad; |
---|
| 1177 | |
---|
| 1178 | flex = Sk_WR(q,L,B); |
---|
| 1179 | |
---|
| 1180 | crossSect = EllipticalCross_fn(q,rad,(rad*ellRatio)); |
---|
| 1181 | flex *= crossSect; |
---|
| 1182 | flex *= Pi*rad*rad*ellRatio*L; |
---|
| 1183 | flex *= cont*cont; |
---|
| 1184 | flex *= 1.0e8; |
---|
| 1185 | answer = scale*flex + bkg; |
---|
| 1186 | |
---|
| 1187 | return answer; |
---|
| 1188 | } |
---|
| 1189 | |
---|
| 1190 | double |
---|
| 1191 | EllipticalCross_fn(double qq, double a, double b) |
---|
| 1192 | { |
---|
| 1193 | double uplim,lolim,Pi,summ,arg,zi,yyy,answer; |
---|
| 1194 | int i,nord=76; |
---|
| 1195 | |
---|
| 1196 | Pi = 4.0*atan(1.0); |
---|
| 1197 | lolim=0.0; |
---|
| 1198 | uplim=Pi/2.0; |
---|
| 1199 | summ=0.0; |
---|
| 1200 | |
---|
| 1201 | for(i=0;i<nord;i++) { |
---|
| 1202 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 1203 | arg = qq*sqrt(a*a*sin(zi)*sin(zi)+b*b*cos(zi)*cos(zi)); |
---|
| 1204 | yyy = pow((2.0 * NR_BessJ1(arg) / arg),2); |
---|
| 1205 | yyy *= Gauss76Wt[i]; |
---|
| 1206 | summ += yyy; |
---|
| 1207 | } |
---|
| 1208 | answer = (uplim-lolim)/2.0*summ; |
---|
| 1209 | answer *= 2.0/Pi; |
---|
| 1210 | return(answer); |
---|
| 1211 | |
---|
| 1212 | } |
---|
| 1213 | /* FlexCyl_PolyLenX : calculates the form factor of a flecible cylinder at the given x-value p->x |
---|
| 1214 | the cylinder has a polydisperse Length |
---|
| 1215 | |
---|
| 1216 | */ |
---|
| 1217 | double |
---|
| 1218 | FlexCyl_PolyLen(double dp[], double q) |
---|
| 1219 | { |
---|
| 1220 | int i; |
---|
[235] | 1221 | double scale,radius,length,pd,bkg,lb,delrho,sldc,slds; //local variables of coefficient wave |
---|
[97] | 1222 | int nord=20; //order of integration |
---|
| 1223 | double uplim,lolim; //upper and lower integration limits |
---|
| 1224 | double summ,zi,yyy,answer,Vpoly; //running tally of integration |
---|
| 1225 | double range,zz,Pi; |
---|
| 1226 | |
---|
| 1227 | Pi = 4.0*atan(1.0); |
---|
| 1228 | range = 3.4; |
---|
| 1229 | |
---|
| 1230 | summ = 0.0; //initialize intergral |
---|
| 1231 | scale = dp[0]; //make local copies in case memory moves |
---|
| 1232 | length = dp[1]; //radius |
---|
| 1233 | pd = dp[2]; // average length |
---|
| 1234 | lb = dp[3]; |
---|
| 1235 | radius = dp[4]; |
---|
[235] | 1236 | sldc = dp[5]; |
---|
| 1237 | slds = dp[6]; |
---|
| 1238 | bkg = dp[7]; |
---|
[97] | 1239 | |
---|
[235] | 1240 | delrho = sldc - slds; |
---|
[97] | 1241 | zz = (1.0/pd)*(1.0/pd) - 1.0; |
---|
| 1242 | |
---|
| 1243 | lolim = length*(1.0-range*pd); //set the upper/lower limits to cover the distribution |
---|
[632] | 1244 | if(lolim<0.0) { |
---|
| 1245 | lolim = 0.0; |
---|
[97] | 1246 | } |
---|
| 1247 | if(pd>0.3) { |
---|
| 1248 | range = 3.4 + (pd-0.3)*18.0; |
---|
| 1249 | } |
---|
| 1250 | uplim = length*(1.0+range*pd); |
---|
| 1251 | |
---|
| 1252 | for(i=0;i<nord;i++) { |
---|
| 1253 | zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 1254 | yyy = Gauss20Wt[i] * FlePolyLen_kernel(q,radius,length,lb,zz,delrho,zi); |
---|
| 1255 | summ += yyy; |
---|
| 1256 | } |
---|
| 1257 | |
---|
| 1258 | answer = (uplim-lolim)/2.0*summ; |
---|
| 1259 | //normalize by average cylinder volume (first moment), using the average length |
---|
| 1260 | Vpoly=Pi*radius*radius*length; |
---|
| 1261 | answer /= Vpoly; |
---|
| 1262 | |
---|
| 1263 | answer *=delrho*delrho; |
---|
| 1264 | |
---|
| 1265 | //convert to [cm-1] |
---|
| 1266 | answer *= 1.0e8; |
---|
| 1267 | //Scale |
---|
| 1268 | answer *= scale; |
---|
| 1269 | // add in the background |
---|
| 1270 | answer += bkg; |
---|
| 1271 | |
---|
| 1272 | return answer; |
---|
| 1273 | } |
---|
| 1274 | |
---|
| 1275 | /* FlexCyl_PolyLenX : calculates the form factor of a flexible cylinder at the given x-value p->x |
---|
| 1276 | the cylinder has a polydisperse cross sectional radius |
---|
| 1277 | |
---|
| 1278 | */ |
---|
| 1279 | double |
---|
| 1280 | FlexCyl_PolyRad(double dp[], double q) |
---|
| 1281 | { |
---|
| 1282 | int i; |
---|
[235] | 1283 | double scale,radius,length,pd,delrho,bkg,lb,sldc,slds; //local variables of coefficient wave |
---|
[97] | 1284 | int nord=76; //order of integration |
---|
| 1285 | double uplim,lolim; //upper and lower integration limits |
---|
| 1286 | double summ,zi,yyy,answer,Vpoly; //running tally of integration |
---|
| 1287 | double range,zz,Pi; |
---|
| 1288 | |
---|
| 1289 | |
---|
| 1290 | Pi = 4.0*atan(1.0); |
---|
| 1291 | range = 3.4; |
---|
| 1292 | |
---|
| 1293 | summ = 0.0; //initialize intergral |
---|
| 1294 | |
---|
| 1295 | scale = dp[0]; //make local copies in case memory moves |
---|
| 1296 | length = dp[1]; //radius |
---|
| 1297 | lb = dp[2]; // average length |
---|
| 1298 | radius = dp[3]; |
---|
| 1299 | pd = dp[4]; |
---|
[235] | 1300 | sldc = dp[5]; |
---|
| 1301 | slds = dp[6]; |
---|
| 1302 | bkg = dp[7]; |
---|
[97] | 1303 | |
---|
[235] | 1304 | delrho = sldc-slds; |
---|
[97] | 1305 | zz = (1.0/pd)*(1.0/pd) - 1.0; |
---|
| 1306 | |
---|
| 1307 | lolim = radius*(1.0-range*pd); //set the upper/lower limits to cover the distribution |
---|
[632] | 1308 | if(lolim<0.0) { |
---|
| 1309 | lolim = 0.0; |
---|
[97] | 1310 | } |
---|
| 1311 | if(pd>0.3) { |
---|
| 1312 | range = 3.4 + (pd-0.3)*18.0; |
---|
| 1313 | } |
---|
| 1314 | uplim = radius*(1.0+range*pd); |
---|
| 1315 | |
---|
| 1316 | for(i=0;i<nord;i++) { |
---|
| 1317 | //zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 1318 | //yyy = Gauss20Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi); |
---|
| 1319 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 1320 | yyy = Gauss76Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi); |
---|
| 1321 | summ += yyy; |
---|
| 1322 | } |
---|
| 1323 | |
---|
| 1324 | answer = (uplim-lolim)/2.0*summ; |
---|
| 1325 | //normalize by average cylinder volume (second moment), using the average radius |
---|
| 1326 | Vpoly = Pi*radius*radius*length*(zz+2.0)/(zz+1.0); |
---|
| 1327 | answer /= Vpoly; |
---|
| 1328 | |
---|
| 1329 | answer *=delrho*delrho; |
---|
| 1330 | |
---|
| 1331 | //convert to [cm-1] |
---|
| 1332 | answer *= 1.0e8; |
---|
| 1333 | //Scale |
---|
| 1334 | answer *= scale; |
---|
| 1335 | // add in the background |
---|
| 1336 | answer += bkg; |
---|
| 1337 | |
---|
| 1338 | return answer; |
---|
| 1339 | } |
---|
| 1340 | |
---|
| 1341 | /////////functions for WRC implementation of flexible cylinders |
---|
| 1342 | static double |
---|
| 1343 | Sk_WR(double q, double L, double b) |
---|
| 1344 | { |
---|
| 1345 | // |
---|
| 1346 | double p1,p2,p1short,p2short,q0,qconnect; |
---|
| 1347 | double C,epsilon,ans,q0short,Sexvmodify,pi; |
---|
| 1348 | |
---|
| 1349 | pi = 4.0*atan(1.0); |
---|
| 1350 | |
---|
| 1351 | p1 = 4.12; |
---|
| 1352 | p2 = 4.42; |
---|
| 1353 | p1short = 5.36; |
---|
| 1354 | p2short = 5.62; |
---|
| 1355 | q0 = 3.1; |
---|
| 1356 | qconnect = q0/b; |
---|
| 1357 | // |
---|
| 1358 | q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); |
---|
| 1359 | |
---|
| 1360 | // |
---|
| 1361 | if(L/b > 10.0) { |
---|
| 1362 | C = 3.06/pow((L/b),0.44); |
---|
| 1363 | epsilon = 0.176; |
---|
| 1364 | } else { |
---|
| 1365 | C = 1.0; |
---|
| 1366 | epsilon = 0.170; |
---|
| 1367 | } |
---|
| 1368 | // |
---|
| 1369 | |
---|
| 1370 | if( L > 4*b ) { // Longer Chains |
---|
| 1371 | if (q*b <= 3.1) { //Modified by Yun on Oct. 15, |
---|
| 1372 | Sexvmodify = Sexvnew(q, L, b); |
---|
| 1373 | ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L); |
---|
| 1374 | } else { //q(i)*b > 3.1 |
---|
| 1375 | ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L); |
---|
| 1376 | } |
---|
| 1377 | } else { //L <= 4*b Shorter Chains |
---|
| 1378 | if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { |
---|
| 1379 | if (q*b<=0.01) { |
---|
| 1380 | ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0; |
---|
| 1381 | } else { |
---|
| 1382 | ans = Sdebye1(q,L,b); |
---|
| 1383 | } |
---|
| 1384 | } else { //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) |
---|
| 1385 | ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + pi/(q*L); |
---|
| 1386 | } |
---|
| 1387 | } |
---|
| 1388 | |
---|
| 1389 | return(ans); |
---|
| 1390 | //return(a2long(q, L, b, p1, p2, q0)); |
---|
| 1391 | } |
---|
| 1392 | |
---|
| 1393 | //WR named this w (too generic) |
---|
| 1394 | static double |
---|
| 1395 | w_WR(double x) |
---|
| 1396 | { |
---|
| 1397 | double yy; |
---|
| 1398 | yy = 0.5*(1 + tanh((x - 1.523)/0.1477)); |
---|
| 1399 | |
---|
| 1400 | return (yy); |
---|
| 1401 | } |
---|
| 1402 | |
---|
| 1403 | // |
---|
| 1404 | static double |
---|
| 1405 | u1(double q, double L, double b) |
---|
| 1406 | { |
---|
| 1407 | double yy; |
---|
| 1408 | |
---|
| 1409 | yy = Rgsquareshort(q,L,b)*q*q; |
---|
| 1410 | |
---|
| 1411 | return (yy); |
---|
| 1412 | } |
---|
| 1413 | |
---|
| 1414 | // was named u |
---|
| 1415 | static double |
---|
| 1416 | u_WR(double q, double L, double b) |
---|
| 1417 | { |
---|
| 1418 | double yy; |
---|
| 1419 | yy = Rgsquare(q,L,b)*q*q; |
---|
| 1420 | return (yy); |
---|
| 1421 | } |
---|
| 1422 | |
---|
| 1423 | |
---|
| 1424 | |
---|
| 1425 | // |
---|
| 1426 | static double |
---|
| 1427 | Rgsquarezero(double q, double L, double b) |
---|
| 1428 | { |
---|
| 1429 | double yy; |
---|
| 1430 | yy = (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) - 0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b)))); |
---|
| 1431 | |
---|
| 1432 | return (yy); |
---|
| 1433 | } |
---|
| 1434 | |
---|
| 1435 | // |
---|
| 1436 | static double |
---|
| 1437 | Rgsquareshort(double q, double L, double b) |
---|
| 1438 | { |
---|
| 1439 | double yy; |
---|
| 1440 | yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b); |
---|
| 1441 | |
---|
| 1442 | return (yy); |
---|
| 1443 | } |
---|
| 1444 | |
---|
| 1445 | // |
---|
| 1446 | static double |
---|
| 1447 | Rgsquare(double q, double L, double b) |
---|
| 1448 | { |
---|
| 1449 | double yy; |
---|
| 1450 | yy = AlphaSquare(L/b)*L*b/6.0; |
---|
| 1451 | |
---|
| 1452 | return (yy); |
---|
| 1453 | } |
---|
| 1454 | |
---|
| 1455 | // |
---|
| 1456 | static double |
---|
| 1457 | AlphaSquare(double x) |
---|
| 1458 | { |
---|
| 1459 | double yy; |
---|
| 1460 | yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); |
---|
| 1461 | |
---|
| 1462 | return (yy); |
---|
| 1463 | } |
---|
| 1464 | |
---|
| 1465 | // ?? funciton is not used - but should the log actually be log10??? |
---|
| 1466 | static double |
---|
| 1467 | miu(double x) |
---|
| 1468 | { |
---|
| 1469 | double yy; |
---|
| 1470 | yy = (1.0/8.0)*(9.0*x - 2.0 + 2.0*log(1.0 + x)/x)*exp(1.0/2.565*(1.0/x + (1.0 - 1.0/(x*x))*log(1.0 + x))); |
---|
| 1471 | |
---|
| 1472 | return (yy); |
---|
| 1473 | } |
---|
| 1474 | |
---|
| 1475 | // |
---|
| 1476 | static double |
---|
| 1477 | Sdebye(double q, double L, double b) |
---|
| 1478 | { |
---|
| 1479 | double yy; |
---|
| 1480 | yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2)); |
---|
| 1481 | |
---|
| 1482 | return (yy); |
---|
| 1483 | } |
---|
| 1484 | |
---|
| 1485 | // |
---|
| 1486 | static double |
---|
| 1487 | Sdebye1(double q, double L, double b) |
---|
| 1488 | { |
---|
| 1489 | double yy; |
---|
| 1490 | yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) ); |
---|
| 1491 | |
---|
| 1492 | return (yy); |
---|
| 1493 | } |
---|
| 1494 | |
---|
| 1495 | // |
---|
| 1496 | static double |
---|
| 1497 | Sexv(double q, double L, double b) |
---|
| 1498 | { |
---|
| 1499 | double yy,C1,C2,C3,miu,Rg2; |
---|
| 1500 | C1=1.22; |
---|
| 1501 | C2=0.4288; |
---|
| 1502 | C3=-1.651; |
---|
| 1503 | miu = 0.585; |
---|
| 1504 | |
---|
| 1505 | Rg2 = Rgsquare(q,L,b); |
---|
| 1506 | |
---|
| 1507 | yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu))); |
---|
| 1508 | |
---|
| 1509 | return (yy); |
---|
| 1510 | } |
---|
| 1511 | |
---|
| 1512 | // this must be WR modified version |
---|
| 1513 | static double |
---|
| 1514 | Sexvnew(double q, double L, double b) |
---|
| 1515 | { |
---|
| 1516 | double yy,C1,C2,C3,miu; |
---|
| 1517 | double del=1.05,C_star2,Rg2; |
---|
| 1518 | |
---|
| 1519 | C1=1.22; |
---|
| 1520 | C2=0.4288; |
---|
| 1521 | C3=-1.651; |
---|
| 1522 | miu = 0.585; |
---|
| 1523 | |
---|
| 1524 | //calculating the derivative to decide on the corection (cutoff) term? |
---|
| 1525 | // I have modified this from WRs original code |
---|
| 1526 | |
---|
| 1527 | if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) { |
---|
| 1528 | C_star2 = 0.0; |
---|
| 1529 | } else { |
---|
| 1530 | C_star2 = 1.0; |
---|
| 1531 | } |
---|
| 1532 | |
---|
| 1533 | Rg2 = Rgsquare(q,L,b); |
---|
| 1534 | |
---|
| 1535 | yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu))); |
---|
| 1536 | |
---|
| 1537 | return (yy); |
---|
| 1538 | } |
---|
| 1539 | |
---|
| 1540 | // these are the messy ones |
---|
| 1541 | static double |
---|
| 1542 | a2short(double q, double L, double b, double p1short, double p2short, double q0) |
---|
| 1543 | { |
---|
| 1544 | double yy,Rg2_sh; |
---|
| 1545 | double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p; |
---|
| 1546 | double pi; |
---|
| 1547 | |
---|
| 1548 | E = 2.718281828459045091; |
---|
| 1549 | pi = 4.0*atan(1.0); |
---|
| 1550 | Rg2_sh = Rgsquareshort(q,L,b); |
---|
| 1551 | Rg2_sh2 = Rg2_sh*Rg2_sh; |
---|
| 1552 | t1 = ((q0*q0*Rg2_sh)/(b*b)); |
---|
| 1553 | Et1 = pow(E,t1); |
---|
| 1554 | Emt1 =pow(E,-t1); |
---|
| 1555 | q02 = q0*q0; |
---|
| 1556 | q0p = pow(q0,(-4.0 + p2short) ); |
---|
| 1557 | |
---|
| 1558 | //E is the number e |
---|
| 1559 | yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short + 2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p1short*pi*q02*q0*Rg2_sh2))))))); |
---|
| 1560 | |
---|
| 1561 | return (yy); |
---|
| 1562 | } |
---|
| 1563 | |
---|
| 1564 | // |
---|
| 1565 | static double |
---|
| 1566 | a1short(double q, double L, double b, double p1short, double p2short, double q0) |
---|
| 1567 | { |
---|
| 1568 | double yy,Rg2_sh; |
---|
| 1569 | double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3; |
---|
| 1570 | double pi; |
---|
| 1571 | |
---|
| 1572 | E = 2.718281828459045091; |
---|
| 1573 | pi = 4.0*atan(1.0); |
---|
| 1574 | Rg2_sh = Rgsquareshort(q,L,b); |
---|
| 1575 | Rg2_sh2 = Rg2_sh*Rg2_sh; |
---|
| 1576 | b3 = b*b*b; |
---|
| 1577 | t1 = ((q0*q0*Rg2_sh)/(b*b)); |
---|
| 1578 | Et1 = pow(E,t1); |
---|
| 1579 | Emt1 =pow(E,-t1); |
---|
| 1580 | q02 = q0*q0; |
---|
| 1581 | q0p = pow(q0,(-4.0 + p1short) ); |
---|
[632] | 1582 | |
---|
| 1583 | yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p2short*pi*q02*q0*Rg2_sh2)))))); |
---|
| 1584 | |
---|
[97] | 1585 | return(yy); |
---|
| 1586 | } |
---|
| 1587 | |
---|
| 1588 | // this one will be lots of trouble |
---|
| 1589 | static double |
---|
| 1590 | a2long(double q, double L, double b, double p1, double p2, double q0) |
---|
| 1591 | { |
---|
| 1592 | double yy,C1,C2,C3,C4,C5,miu,C,Rg2; |
---|
| 1593 | double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi; |
---|
| 1594 | double E,b2,b3,b4,q02,q03,q04,q05,Rg22; |
---|
| 1595 | |
---|
| 1596 | pi = 4.0*atan(1.0); |
---|
| 1597 | E = 2.718281828459045091; |
---|
| 1598 | if( L/b > 10.0) { |
---|
| 1599 | C = 3.06/pow((L/b),0.44); |
---|
| 1600 | } else { |
---|
| 1601 | C = 1.0; |
---|
| 1602 | } |
---|
| 1603 | |
---|
| 1604 | C1 = 1.22; |
---|
| 1605 | C2 = 0.4288; |
---|
| 1606 | C3 = -1.651; |
---|
| 1607 | C4 = 1.523; |
---|
| 1608 | C5 = 0.1477; |
---|
| 1609 | miu = 0.585; |
---|
| 1610 | |
---|
| 1611 | Rg2 = Rgsquare(q,L,b); |
---|
| 1612 | Rg22 = Rg2*Rg2; |
---|
| 1613 | b2 = b*b; |
---|
| 1614 | b3 = b*b*b; |
---|
| 1615 | b4 = b3*b; |
---|
| 1616 | q02 = q0*q0; |
---|
| 1617 | q03 = q0*q0*q0; |
---|
| 1618 | q04 = q03*q0; |
---|
| 1619 | q05 = q04*q0; |
---|
| 1620 | |
---|
| 1621 | t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) )); |
---|
[632] | 1622 | |
---|
| 1623 | t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L; |
---|
| 1624 | |
---|
| 1625 | t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2.0))/(2.0*C5); |
---|
| 1626 | |
---|
| 1627 | t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22); |
---|
| 1628 | |
---|
| 1629 | t5 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); |
---|
| 1630 | |
---|
| 1631 | t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22); |
---|
| 1632 | |
---|
| 1633 | t7 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu)); |
---|
| 1634 | |
---|
| 1635 | t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); |
---|
| 1636 | |
---|
| 1637 | t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L; |
---|
| 1638 | |
---|
| 1639 | t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); |
---|
| 1640 | |
---|
| 1641 | |
---|
| 1642 | yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))); |
---|
| 1643 | |
---|
[97] | 1644 | return (yy); |
---|
| 1645 | } |
---|
| 1646 | |
---|
| 1647 | //need to define this on my own |
---|
| 1648 | static double |
---|
| 1649 | sech_WR(double x) |
---|
| 1650 | { |
---|
| 1651 | return(1/cosh(x)); |
---|
| 1652 | } |
---|
| 1653 | |
---|
| 1654 | // |
---|
| 1655 | static double |
---|
| 1656 | a1long(double q, double L, double b, double p1, double p2, double q0) |
---|
| 1657 | { |
---|
| 1658 | double yy,C,C1,C2,C3,C4,C5,miu,Rg2; |
---|
| 1659 | double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15; |
---|
| 1660 | double E,pi; |
---|
| 1661 | double b2,b3,b4,q02,q03,q04,q05,Rg22; |
---|
| 1662 | |
---|
| 1663 | pi = 4.0*atan(1.0); |
---|
| 1664 | E = 2.718281828459045091; |
---|
| 1665 | |
---|
| 1666 | if( L/b > 10.0) { |
---|
| 1667 | C = 3.06/pow((L/b),0.44); |
---|
| 1668 | } else { |
---|
| 1669 | C = 1.0; |
---|
| 1670 | } |
---|
| 1671 | |
---|
| 1672 | C1 = 1.22; |
---|
| 1673 | C2 = 0.4288; |
---|
| 1674 | C3 = -1.651; |
---|
| 1675 | C4 = 1.523; |
---|
| 1676 | C5 = 0.1477; |
---|
| 1677 | miu = 0.585; |
---|
| 1678 | |
---|
| 1679 | Rg2 = Rgsquare(q,L,b); |
---|
| 1680 | Rg22 = Rg2*Rg2; |
---|
| 1681 | b2 = b*b; |
---|
| 1682 | b3 = b*b*b; |
---|
| 1683 | b4 = b3*b; |
---|
| 1684 | q02 = q0*q0; |
---|
| 1685 | q03 = q0*q0*q0; |
---|
| 1686 | q04 = q03*q0; |
---|
| 1687 | q05 = q04*q0; |
---|
| 1688 | |
---|
| 1689 | t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2)))); |
---|
| 1690 | |
---|
| 1691 | t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); |
---|
| 1692 | |
---|
| 1693 | t3 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu)))); |
---|
| 1694 | |
---|
| 1695 | t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); |
---|
| 1696 | |
---|
| 1697 | t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)))); |
---|
| 1698 | |
---|
| 1699 | t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2)))*Rg2)/b))); |
---|
| 1700 | |
---|
| 1701 | t7 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); |
---|
| 1702 | |
---|
| 1703 | t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); |
---|
| 1704 | |
---|
| 1705 | t9 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); |
---|
| 1706 | |
---|
| 1707 | t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); |
---|
| 1708 | |
---|
| 1709 | t11 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu)); |
---|
| 1710 | |
---|
| 1711 | t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); |
---|
| 1712 | |
---|
| 1713 | t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02* Rg2))) + (7.0*b2)/(15.0*q02*Rg2)))); |
---|
| 1714 | |
---|
| 1715 | t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); |
---|
| 1716 | |
---|
| 1717 | t15 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu)))); |
---|
| 1718 | |
---|
| 1719 | |
---|
| 1720 | yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) + 1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*(((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -t8/(C5*q04*Rg22) +t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) +t13/L +t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))))); |
---|
| 1721 | |
---|
| 1722 | return (yy); |
---|
| 1723 | } |
---|
| 1724 | |
---|
| 1725 | |
---|
| 1726 | |
---|
| 1727 | /////////////// |
---|
| 1728 | |
---|
| 1729 | // |
---|
| 1730 | // FUNCTION gfn2: CONTAINS F(Q,A,B,mu)**2 AS GIVEN |
---|
| 1731 | // BY (53) AND (56,57) IN CHEN AND |
---|
| 1732 | // KOTLARCHYK REFERENCE |
---|
| 1733 | // |
---|
| 1734 | // <PROLATE ELLIPSOIDS> |
---|
| 1735 | // |
---|
| 1736 | double |
---|
| 1737 | gfn2(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq) |
---|
| 1738 | { |
---|
| 1739 | // local variables |
---|
[632] | 1740 | double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,gfn2,pi43,gfn,Pi; |
---|
| 1741 | |
---|
[97] | 1742 | Pi = 4.0*atan(1.0); |
---|
| 1743 | |
---|
| 1744 | pi43=4.0/3.0*Pi; |
---|
| 1745 | aa = crmaj; |
---|
| 1746 | bb = crmin; |
---|
| 1747 | u2 = (aa*aa*xx*xx + bb*bb*(1.0-xx*xx)); |
---|
| 1748 | ut2 = (trmaj*trmaj*xx*xx + trmin*trmin*(1.0-xx*xx)); |
---|
| 1749 | uq = sqrt(u2)*qq; |
---|
| 1750 | ut= sqrt(ut2)*qq; |
---|
| 1751 | vc = pi43*aa*bb*bb; |
---|
| 1752 | vt = pi43*trmaj*trmin*trmin; |
---|
[632] | 1753 | if (uq == 0.0){ |
---|
| 1754 | siq = 1.0/3.0; |
---|
| 1755 | }else{ |
---|
| 1756 | siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq; |
---|
| 1757 | } |
---|
| 1758 | if (ut == 0.0){ |
---|
| 1759 | sit = 1.0/3.0; |
---|
| 1760 | }else{ |
---|
| 1761 | sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut; |
---|
| 1762 | } |
---|
| 1763 | gfnc = 3.0*siq*vc*delpc; |
---|
| 1764 | gfnt = 3.0*sit*vt*delps; |
---|
[97] | 1765 | gfn = gfnc+gfnt; |
---|
| 1766 | gfn2 = gfn*gfn; |
---|
| 1767 | |
---|
| 1768 | return (gfn2); |
---|
| 1769 | } |
---|
| 1770 | |
---|
| 1771 | // |
---|
| 1772 | // FUNCTION gfn4: CONTAINS F(Q,A,B,MU)**2 AS GIVEN |
---|
| 1773 | // BY (53) & (58-59) IN CHEN AND |
---|
| 1774 | // KOTLARCHYK REFERENCE |
---|
| 1775 | // |
---|
| 1776 | // <OBLATE ELLIPSOID> |
---|
| 1777 | // function gfn4 for oblate ellipsoids |
---|
| 1778 | double |
---|
| 1779 | gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq) |
---|
| 1780 | { |
---|
| 1781 | // local variables |
---|
[632] | 1782 | double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi; |
---|
| 1783 | |
---|
[97] | 1784 | Pi = 4.0*atan(1.0); |
---|
| 1785 | pi43=4.0/3.0*Pi; |
---|
| 1786 | aa = crmaj; |
---|
| 1787 | bb = crmin; |
---|
| 1788 | u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx)); |
---|
| 1789 | ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx)); |
---|
| 1790 | uq = sqrt(u2)*qq; |
---|
| 1791 | ut= sqrt(ut2)*qq; |
---|
| 1792 | vc = pi43*aa*aa*bb; |
---|
| 1793 | vt = pi43*trmaj*trmaj*trmin; |
---|
[632] | 1794 | if (uq == 0.0){ |
---|
| 1795 | siq = 1.0/3.0; |
---|
| 1796 | }else{ |
---|
| 1797 | siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq; |
---|
| 1798 | } |
---|
| 1799 | if (ut == 0.0){ |
---|
| 1800 | sit = 1.0/3.0; |
---|
| 1801 | }else{ |
---|
| 1802 | sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut; |
---|
| 1803 | } |
---|
| 1804 | gfnc = 3.0*siq*vc*delpc; |
---|
| 1805 | gfnt = 3.0*sit*vt*delps; |
---|
[97] | 1806 | tgfn = gfnc+gfnt; |
---|
| 1807 | gfn4 = tgfn*tgfn; |
---|
| 1808 | |
---|
| 1809 | return (gfn4); |
---|
| 1810 | } |
---|
| 1811 | |
---|
| 1812 | double |
---|
| 1813 | FlePolyLen_kernel(double q, double radius, double length, double lb, double zz, double delrho, double zi) |
---|
| 1814 | { |
---|
| 1815 | double Pq,vcyl,dl; |
---|
| 1816 | double Pi,qr; |
---|
| 1817 | |
---|
| 1818 | Pi = 4.0*atan(1.0); |
---|
| 1819 | qr = q*radius; |
---|
| 1820 | |
---|
| 1821 | Pq = Sk_WR(q,zi,lb); //does not have cross section term |
---|
[632] | 1822 | if (qr !=0){ |
---|
| 1823 | Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); |
---|
| 1824 | } |
---|
[97] | 1825 | vcyl=Pi*radius*radius*zi; |
---|
| 1826 | Pq *= vcyl*vcyl; |
---|
| 1827 | |
---|
| 1828 | dl = SchulzPoint_cpr(zi,length,zz); |
---|
| 1829 | return (Pq*dl); |
---|
| 1830 | |
---|
| 1831 | } |
---|
| 1832 | |
---|
| 1833 | double |
---|
| 1834 | FlePolyRad_kernel(double q, double ravg, double Lc, double Lb, double zz, double delrho, double zi) |
---|
| 1835 | { |
---|
| 1836 | double Pq,vcyl,dr; |
---|
| 1837 | double Pi,qr; |
---|
| 1838 | |
---|
| 1839 | Pi = 4.0*atan(1.0); |
---|
| 1840 | qr = q*zi; |
---|
| 1841 | |
---|
| 1842 | Pq = Sk_WR(q,Lc,Lb); //does not have cross section term |
---|
[632] | 1843 | if (qr !=0){ |
---|
| 1844 | Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); |
---|
| 1845 | } |
---|
| 1846 | |
---|
[97] | 1847 | vcyl=Pi*zi*zi*Lc; |
---|
| 1848 | Pq *= vcyl*vcyl; |
---|
| 1849 | |
---|
| 1850 | dr = SchulzPoint_cpr(zi,ravg,zz); |
---|
| 1851 | return (Pq*dr); |
---|
| 1852 | |
---|
| 1853 | } |
---|
| 1854 | |
---|
| 1855 | double |
---|
| 1856 | CSCylIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length) |
---|
| 1857 | { |
---|
| 1858 | double answer,halfheight,Pi; |
---|
| 1859 | double lolim,uplim,summ,yyy,zi; |
---|
| 1860 | int nord,i; |
---|
| 1861 | |
---|
| 1862 | // set up the integration end points |
---|
| 1863 | Pi = 4.0*atan(1.0); |
---|
| 1864 | nord = 76; |
---|
[632] | 1865 | lolim = 0.0; |
---|
| 1866 | uplim = Pi/2.0; |
---|
[97] | 1867 | halfheight = length/2.0; |
---|
| 1868 | |
---|
| 1869 | summ = 0.0; // initialize integral |
---|
| 1870 | i=0; |
---|
| 1871 | for(i=0;i<nord;i++) { |
---|
| 1872 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 1873 | yyy = Gauss76Wt[i] * CScyl(qq, rad, radthick, facthick, rhoc,rhos,rhosolv, halfheight, zi); |
---|
| 1874 | summ += yyy; |
---|
| 1875 | } |
---|
| 1876 | |
---|
| 1877 | // calculate value of integral to return |
---|
| 1878 | answer = (uplim-lolim)/2.0*summ; |
---|
| 1879 | return (answer); |
---|
| 1880 | } |
---|
| 1881 | |
---|
| 1882 | double |
---|
| 1883 | CScyl(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length, double dum) |
---|
| 1884 | { |
---|
| 1885 | // qq is the q-value for the calculation (1/A) |
---|
| 1886 | // radius is the core radius of the cylinder (A) |
---|
| 1887 | // radthick and facthick are the radial and face layer thicknesses |
---|
| 1888 | // rho(n) are the respective SLD's |
---|
| 1889 | // length is the *Half* CORE-LENGTH of the cylinder |
---|
| 1890 | // dum is the dummy variable for the integration (theta) |
---|
[632] | 1891 | |
---|
| 1892 | double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval; |
---|
[97] | 1893 | double Pi; |
---|
| 1894 | |
---|
| 1895 | Pi = 4.0*atan(1.0); |
---|
| 1896 | |
---|
| 1897 | dr1 = rhoc-rhos; |
---|
| 1898 | dr2 = rhos-rhosolv; |
---|
| 1899 | vol1 = Pi*rad*rad*(2.0*length); |
---|
| 1900 | vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick); |
---|
| 1901 | |
---|
| 1902 | besarg1 = qq*rad*sin(dum); |
---|
| 1903 | besarg2 = qq*(rad+radthick)*sin(dum); |
---|
| 1904 | sinarg1 = qq*length*cos(dum); |
---|
| 1905 | sinarg2 = qq*(length+facthick)*cos(dum); |
---|
[632] | 1906 | if (besarg1 == 0.0){ |
---|
| 1907 | be1 = 0.5; |
---|
| 1908 | }else{ |
---|
| 1909 | be1 = NR_BessJ1(besarg1)/besarg1; |
---|
| 1910 | } |
---|
| 1911 | if (besarg2 == 0.0){ |
---|
| 1912 | be2 = 0.5; |
---|
| 1913 | }else{ |
---|
| 1914 | be2 = NR_BessJ1(besarg2)/besarg2; |
---|
| 1915 | } |
---|
| 1916 | if (sinarg1 == 0.0){ |
---|
| 1917 | si1 = 1.0; |
---|
| 1918 | }else{ |
---|
| 1919 | si1 = sin(sinarg1)/sinarg1; |
---|
| 1920 | } |
---|
| 1921 | if (besarg2 == 0.0){ |
---|
| 1922 | si2 = 1.0; |
---|
| 1923 | }else{ |
---|
| 1924 | si2 = sin(sinarg2)/sinarg2; |
---|
| 1925 | } |
---|
| 1926 | |
---|
| 1927 | t1 = 2.0*vol1*dr1*si1*be1; |
---|
| 1928 | t2 = 2.0*vol2*dr2*si2*be2; |
---|
| 1929 | |
---|
[97] | 1930 | retval = ((t1+t2)*(t1+t2))*sin(dum); |
---|
| 1931 | return (retval); |
---|
| 1932 | |
---|
| 1933 | } |
---|
| 1934 | |
---|
| 1935 | |
---|
| 1936 | double |
---|
| 1937 | CoreShellCylKernel(double qq, double rcore, double thick, double rhoc, double rhos, double rhosolv, double length, double dum) |
---|
| 1938 | { |
---|
[632] | 1939 | |
---|
| 1940 | double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval; |
---|
[97] | 1941 | double Pi; |
---|
| 1942 | |
---|
| 1943 | Pi = 4.0*atan(1.0); |
---|
| 1944 | |
---|
| 1945 | dr1 = rhoc-rhos; |
---|
| 1946 | dr2 = rhos-rhosolv; |
---|
| 1947 | vol1 = Pi*rcore*rcore*(2.0*length); |
---|
| 1948 | vol2 = Pi*(rcore+thick)*(rcore+thick)*(2.0*length+2.0*thick); |
---|
| 1949 | |
---|
| 1950 | besarg1 = qq*rcore*sin(dum); |
---|
| 1951 | besarg2 = qq*(rcore+thick)*sin(dum); |
---|
| 1952 | sinarg1 = qq*length*cos(dum); |
---|
| 1953 | sinarg2 = qq*(length+thick)*cos(dum); |
---|
[632] | 1954 | |
---|
| 1955 | if (besarg1 == 0.0){ |
---|
| 1956 | be1 = 0.5; |
---|
| 1957 | }else{ |
---|
| 1958 | be1 = NR_BessJ1(besarg1)/besarg1; |
---|
| 1959 | } |
---|
| 1960 | if (besarg2 == 0.0){ |
---|
| 1961 | be2 = 0.5; |
---|
| 1962 | }else{ |
---|
| 1963 | be2 = NR_BessJ1(besarg2)/besarg2; |
---|
| 1964 | } |
---|
| 1965 | if (sinarg1 == 0.0){ |
---|
| 1966 | si1 = 1.0; |
---|
| 1967 | }else{ |
---|
| 1968 | si1 = sin(sinarg1)/sinarg1; |
---|
| 1969 | } |
---|
| 1970 | if (besarg2 == 0.0){ |
---|
| 1971 | si2 = 1.0; |
---|
| 1972 | }else{ |
---|
| 1973 | si2 = sin(sinarg2)/sinarg2; |
---|
| 1974 | } |
---|
| 1975 | |
---|
| 1976 | t1 = 2.0*vol1*dr1*si1*be1; |
---|
| 1977 | t2 = 2.0*vol2*dr2*si2*be2; |
---|
| 1978 | |
---|
[97] | 1979 | retval = ((t1+t2)*(t1+t2))*sin(dum); |
---|
| 1980 | |
---|
| 1981 | return (retval); |
---|
| 1982 | } |
---|
| 1983 | |
---|
| 1984 | double |
---|
| 1985 | Cyl_PolyLenKernel(double q, double radius, double len_avg, double zz, double delrho, double dumLen) |
---|
| 1986 | { |
---|
| 1987 | |
---|
| 1988 | double halfheight,uplim,lolim,zi,summ,yyy,Pi; |
---|
| 1989 | double answer,dr,Vcyl; |
---|
| 1990 | int i,nord; |
---|
| 1991 | |
---|
| 1992 | Pi = 4.0*atan(1.0); |
---|
[632] | 1993 | lolim = 0.0; |
---|
[97] | 1994 | uplim = Pi/2.0; |
---|
| 1995 | halfheight = dumLen/2.0; |
---|
| 1996 | nord=20; |
---|
| 1997 | summ = 0.0; |
---|
| 1998 | |
---|
| 1999 | //do the cylinder orientational average |
---|
| 2000 | for(i=0;i<nord;i++) { |
---|
| 2001 | zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 2002 | yyy = Gauss20Wt[i] * CylKernel(q, radius, halfheight, zi); |
---|
| 2003 | summ += yyy; |
---|
| 2004 | } |
---|
| 2005 | answer = (uplim-lolim)/2.0*summ; |
---|
| 2006 | // Multiply by contrast^2 |
---|
| 2007 | answer *= delrho*delrho; |
---|
| 2008 | // don't do the normal scaling to volume here |
---|
| 2009 | // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder |
---|
| 2010 | Vcyl = Pi*radius*radius*dumLen; |
---|
| 2011 | answer *= Vcyl*Vcyl; |
---|
| 2012 | |
---|
| 2013 | dr = SchulzPoint_cpr(dumLen,len_avg,zz); |
---|
| 2014 | return(dr*answer); |
---|
| 2015 | } |
---|
| 2016 | |
---|
| 2017 | |
---|
| 2018 | double |
---|
| 2019 | Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N) |
---|
| 2020 | { |
---|
| 2021 | // qq is the q-value for the calculation (1/A) |
---|
| 2022 | // rcore is the core radius of the cylinder (A) |
---|
| 2023 | // rho(n) are the respective SLD's |
---|
| 2024 | // length is the *Half* CORE-LENGTH of the cylinder = L (A) |
---|
| 2025 | // dum is the dummy variable for the integration (x in Feigin's notation) |
---|
[632] | 2026 | |
---|
| 2027 | //Local variables |
---|
| 2028 | double totald,dr1,dr2,besarg1,besarg2,be1,be2,si1,si2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt; |
---|
[97] | 2029 | double Pi; |
---|
| 2030 | int kk; |
---|
| 2031 | |
---|
| 2032 | Pi = 4.0*atan(1.0); |
---|
| 2033 | |
---|
| 2034 | dr1 = rhoc-rhosolv; |
---|
| 2035 | dr2 = rhol-rhosolv; |
---|
| 2036 | area = Pi*rcore*rcore; |
---|
| 2037 | totald=2.0*(thick+length); |
---|
| 2038 | |
---|
| 2039 | besarg1 = qq*rcore*sin(dum); |
---|
| 2040 | besarg2 = qq*rcore*sin(dum); |
---|
| 2041 | |
---|
| 2042 | sinarg1 = qq*length*cos(dum); |
---|
| 2043 | sinarg2 = qq*(length+thick)*cos(dum); |
---|
[632] | 2044 | |
---|
| 2045 | if (besarg1 == 0.0){ |
---|
| 2046 | be1 = 0.5; |
---|
| 2047 | }else{ |
---|
| 2048 | be1 = NR_BessJ1(besarg1)/besarg1; |
---|
| 2049 | } |
---|
| 2050 | if (besarg2 == 0.0){ |
---|
| 2051 | be2 = 0.5; |
---|
| 2052 | }else{ |
---|
| 2053 | be2 = NR_BessJ1(besarg2)/besarg2; |
---|
| 2054 | } |
---|
| 2055 | if (sinarg1 == 0.0){ |
---|
| 2056 | si1 = 1.0; |
---|
| 2057 | }else{ |
---|
| 2058 | si1 = sin(sinarg1)/sinarg1; |
---|
| 2059 | } |
---|
| 2060 | if (besarg2 == 0.0){ |
---|
| 2061 | si2 = 1.0; |
---|
| 2062 | }else{ |
---|
| 2063 | si2 = sin(sinarg2)/sinarg2; |
---|
| 2064 | } |
---|
| 2065 | |
---|
| 2066 | t1 = 2.0*area*(2.0*length)*dr1*(si1)*(be1); |
---|
| 2067 | t2 = 2.0*area*dr2*(totald*si2-2.0*length*si1)*(be2); |
---|
| 2068 | |
---|
[97] | 2069 | retval =((t1+t2)*(t1+t2))*sin(dum); |
---|
| 2070 | |
---|
| 2071 | // loop for the structure facture S(q) |
---|
| 2072 | sqq=0.0; |
---|
| 2073 | for(kk=1;kk<N;kk+=1) { |
---|
| 2074 | dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0; |
---|
| 2075 | sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt); |
---|
| 2076 | } |
---|
| 2077 | |
---|
| 2078 | // end of loop for S(q) |
---|
| 2079 | sqq=1.0+2.0*sqq/N; |
---|
| 2080 | retval *= sqq; |
---|
| 2081 | |
---|
| 2082 | return(retval); |
---|
| 2083 | } |
---|
| 2084 | |
---|
| 2085 | |
---|
| 2086 | double |
---|
| 2087 | Cyl_PolyRadKernel(double q, double radius, double length, double zz, double delrho, double dumRad) |
---|
| 2088 | { |
---|
| 2089 | |
---|
| 2090 | double halfheight,uplim,lolim,zi,summ,yyy,Pi; |
---|
| 2091 | double answer,dr,Vcyl; |
---|
| 2092 | int i,nord; |
---|
| 2093 | |
---|
| 2094 | Pi = 4.0*atan(1.0); |
---|
[632] | 2095 | lolim = 0.0; |
---|
[97] | 2096 | uplim = Pi/2.0; |
---|
| 2097 | halfheight = length/2.0; |
---|
| 2098 | // nord=20; |
---|
| 2099 | nord=76; |
---|
| 2100 | summ = 0.0; |
---|
| 2101 | |
---|
| 2102 | //do the cylinder orientational average |
---|
| 2103 | // for(i=0;i<nord;i++) { |
---|
| 2104 | // zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 2105 | // yyy = Gauss20Wt[i] * CylKernel(q, dumRad, halfheight, zi); |
---|
| 2106 | // summ += yyy; |
---|
| 2107 | // } |
---|
| 2108 | for(i=0;i<nord;i++) { |
---|
| 2109 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 2110 | yyy = Gauss76Wt[i] * CylKernel(q, dumRad, halfheight, zi); |
---|
| 2111 | summ += yyy; |
---|
| 2112 | } |
---|
| 2113 | answer = (uplim-lolim)/2.0*summ; |
---|
| 2114 | // Multiply by contrast^2 |
---|
| 2115 | answer *= delrho*delrho; |
---|
| 2116 | // don't do the normal scaling to volume here |
---|
| 2117 | // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder |
---|
| 2118 | Vcyl = Pi*dumRad*dumRad*length; |
---|
| 2119 | answer *= Vcyl*Vcyl; |
---|
| 2120 | |
---|
| 2121 | dr = SchulzPoint_cpr(dumRad,radius,zz); |
---|
| 2122 | return(dr*answer); |
---|
| 2123 | } |
---|
| 2124 | |
---|
| 2125 | double |
---|
| 2126 | SchulzPoint_cpr(double dumRad, double radius, double zz) |
---|
| 2127 | { |
---|
| 2128 | double dr; |
---|
| 2129 | |
---|
| 2130 | dr = zz*log(dumRad) - gammaln(zz+1.0) + (zz+1.0)*log((zz+1.0)/radius)-(dumRad/radius*(zz+1.0)); |
---|
| 2131 | return(exp(dr)); |
---|
| 2132 | } |
---|
| 2133 | |
---|
| 2134 | static double |
---|
| 2135 | gammaln(double xx) |
---|
| 2136 | { |
---|
| 2137 | double x,y,tmp,ser; |
---|
| 2138 | static double cof[6]={76.18009172947146,-86.50532032941677, |
---|
| 2139 | 24.01409824083091,-1.231739572450155, |
---|
| 2140 | 0.1208650973866179e-2,-0.5395239384953e-5}; |
---|
| 2141 | int j; |
---|
| 2142 | |
---|
| 2143 | y=x=xx; |
---|
| 2144 | tmp=x+5.5; |
---|
| 2145 | tmp -= (x+0.5)*log(tmp); |
---|
| 2146 | ser=1.000000000190015; |
---|
| 2147 | for (j=0;j<=5;j++) ser += cof[j]/++y; |
---|
| 2148 | return -tmp+log(2.5066282746310005*ser/x); |
---|
| 2149 | } |
---|
| 2150 | |
---|
| 2151 | |
---|
| 2152 | double |
---|
| 2153 | EllipsoidKernel(double qq, double a, double nua, double dum) |
---|
| 2154 | { |
---|
| 2155 | double arg,nu,retval; //local variables |
---|
| 2156 | |
---|
| 2157 | nu = nua/a; |
---|
[632] | 2158 | arg = qq*a*sqrt(1.0+dum*dum*(nu*nu-1.0)); |
---|
| 2159 | if (arg == 0.0){ |
---|
| 2160 | retval =1.0/3.0; |
---|
| 2161 | }else{ |
---|
| 2162 | retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg); |
---|
| 2163 | } |
---|
[97] | 2164 | retval *= retval; |
---|
[632] | 2165 | retval *= 9.0; |
---|
| 2166 | |
---|
[97] | 2167 | return(retval); |
---|
| 2168 | }//Function EllipsoidKernel() |
---|
| 2169 | |
---|
| 2170 | double |
---|
| 2171 | HolCylKernel(double qq, double rcore, double rshell, double length, double dum) |
---|
| 2172 | { |
---|
| 2173 | double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval; //local variables |
---|
| 2174 | |
---|
| 2175 | gamma = rcore/rshell; |
---|
[632] | 2176 | arg1 = qq*rshell*sqrt(1.0-dum*dum); //1=shell (outer radius) |
---|
| 2177 | arg2 = qq*rcore*sqrt(1.0-dum*dum); //2=core (inner radius) |
---|
| 2178 | if (arg1 == 0.0){ |
---|
| 2179 | lam1 = 1.0; |
---|
| 2180 | }else{ |
---|
| 2181 | lam1 = 2.0*NR_BessJ1(arg1)/arg1; |
---|
| 2182 | } |
---|
| 2183 | if (arg2 == 0.0){ |
---|
| 2184 | lam2 = 1.0; |
---|
| 2185 | }else{ |
---|
| 2186 | lam2 = 2.0*NR_BessJ1(arg2)/arg2; |
---|
| 2187 | } |
---|
| 2188 | //Todo: Need to check psi behavior as gamma goes to 1. |
---|
| 2189 | psi = 1.0/(1.0-gamma*gamma)*(lam1 - gamma*gamma*lam2); //SRK 10/19/00 |
---|
| 2190 | sinarg = qq*length*dum/2.0; |
---|
| 2191 | if (sinarg == 0.0){ |
---|
| 2192 | t2 = 1.0; |
---|
| 2193 | }else{ |
---|
| 2194 | t2 = sin(sinarg)/sinarg; |
---|
| 2195 | } |
---|
| 2196 | |
---|
[97] | 2197 | retval = psi*psi*t2*t2; |
---|
| 2198 | |
---|
| 2199 | return(retval); |
---|
| 2200 | }//Function HolCylKernel() |
---|
| 2201 | |
---|
| 2202 | double |
---|
| 2203 | PPKernel(double aa, double mu, double uu) |
---|
| 2204 | { |
---|
| 2205 | // mu passed in is really mu*sqrt(1-sig^2) |
---|
| 2206 | double arg1,arg2,Pi,tmp1,tmp2; //local variables |
---|
| 2207 | |
---|
| 2208 | Pi = 4.0*atan(1.0); |
---|
| 2209 | |
---|
| 2210 | //handle arg=0 separately, as sin(t)/t -> 1 as t->0 |
---|
[632] | 2211 | arg1 = (mu/2.0)*cos(Pi*uu/2.0); |
---|
| 2212 | arg2 = (mu*aa/2.0)*sin(Pi*uu/2.0); |
---|
| 2213 | if(arg1==0.0) { |
---|
| 2214 | tmp1 = 1.0; |
---|
[97] | 2215 | } else { |
---|
| 2216 | tmp1 = sin(arg1)*sin(arg1)/arg1/arg1; |
---|
| 2217 | } |
---|
[632] | 2218 | |
---|
| 2219 | if (arg2==0.0) { |
---|
| 2220 | tmp2 = 1.0; |
---|
[97] | 2221 | } else { |
---|
| 2222 | tmp2 = sin(arg2)*sin(arg2)/arg2/arg2; |
---|
| 2223 | } |
---|
| 2224 | |
---|
| 2225 | return (tmp1*tmp2); |
---|
| 2226 | |
---|
| 2227 | }//Function PPKernel() |
---|
| 2228 | |
---|
| 2229 | |
---|
| 2230 | double |
---|
| 2231 | TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy) |
---|
| 2232 | { |
---|
| 2233 | |
---|
| 2234 | double arg,val,pi; //local variables |
---|
| 2235 | |
---|
| 2236 | pi = 4.0*atan(1.0); |
---|
| 2237 | |
---|
| 2238 | arg = aa*aa*cos(pi*dx/2)*cos(pi*dx/2); |
---|
| 2239 | arg += bb*bb*sin(pi*dx/2)*sin(pi*dx/2)*(1-dy*dy); |
---|
| 2240 | arg += cc*cc*dy*dy; |
---|
| 2241 | arg = q*sqrt(arg); |
---|
[632] | 2242 | if (arg == 0.0){ |
---|
| 2243 | val = 1.0; // as arg --> 0, val should go to 1.0 |
---|
| 2244 | }else{ |
---|
| 2245 | val = 9.0 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ); |
---|
| 2246 | } |
---|
[97] | 2247 | return (val); |
---|
| 2248 | |
---|
| 2249 | }//Function TriaxialKernel() |
---|
| 2250 | |
---|
| 2251 | |
---|
| 2252 | double |
---|
| 2253 | CylKernel(double qq, double rr,double h, double theta) |
---|
| 2254 | { |
---|
| 2255 | |
---|
| 2256 | // qq is the q-value for the calculation (1/A) |
---|
| 2257 | // rr is the radius of the cylinder (A) |
---|
| 2258 | // h is the HALF-LENGTH of the cylinder = L/2 (A) |
---|
[632] | 2259 | |
---|
| 2260 | double besarg,bj,retval,d1,t1,b1,t2,b2,siarg,be,si; //Local variables |
---|
| 2261 | |
---|
| 2262 | |
---|
[97] | 2263 | besarg = qq*rr*sin(theta); |
---|
[632] | 2264 | siarg = qq * h * cos(theta); |
---|
[97] | 2265 | bj =NR_BessJ1(besarg); |
---|
| 2266 | |
---|
| 2267 | //* Computing 2nd power */ |
---|
[632] | 2268 | d1 = sin(siarg); |
---|
[97] | 2269 | t1 = d1 * d1; |
---|
| 2270 | //* Computing 2nd power */ |
---|
| 2271 | d1 = bj; |
---|
| 2272 | t2 = d1 * d1 * 4.0 * sin(theta); |
---|
| 2273 | //* Computing 2nd power */ |
---|
[632] | 2274 | d1 = siarg; |
---|
[97] | 2275 | b1 = d1 * d1; |
---|
| 2276 | //* Computing 2nd power */ |
---|
| 2277 | d1 = qq * rr * sin(theta); |
---|
| 2278 | b2 = d1 * d1; |
---|
[632] | 2279 | if (besarg == 0.0){ |
---|
| 2280 | be = sin(theta); |
---|
| 2281 | }else{ |
---|
| 2282 | be = t2 / b2; |
---|
| 2283 | } |
---|
| 2284 | if (siarg == 0.0){ |
---|
| 2285 | si = 1.0; |
---|
| 2286 | }else{ |
---|
| 2287 | si = t1 / b1; |
---|
| 2288 | } |
---|
| 2289 | retval = be * si; |
---|
| 2290 | |
---|
[97] | 2291 | return (retval); |
---|
| 2292 | |
---|
| 2293 | }//Function CylKernel() |
---|
| 2294 | |
---|
| 2295 | double |
---|
| 2296 | EllipCylKernel(double qq, double ra,double nu, double theta) |
---|
| 2297 | { |
---|
| 2298 | //this is the function LAMBDA1^2 in Feigin's notation |
---|
| 2299 | // qq is the q-value for the calculation (1/A) |
---|
| 2300 | // ra is the transformed radius"a" in Feigin's notation |
---|
| 2301 | // nu is the ratio (major radius)/(minor radius) of the Ellipsoid [=] --- |
---|
| 2302 | // theta is the dummy variable of the integration |
---|
| 2303 | |
---|
| 2304 | double retval,arg; //Local variables |
---|
| 2305 | |
---|
[632] | 2306 | arg = qq*ra*sqrt((1.0+nu*nu)/2+(1.0-nu*nu)*cos(theta)/2); |
---|
| 2307 | if (arg == 0.0){ |
---|
| 2308 | retval = 1.0; |
---|
| 2309 | }else{ |
---|
| 2310 | retval = 2.0*NR_BessJ1(arg)/arg; |
---|
| 2311 | } |
---|
| 2312 | |
---|
[97] | 2313 | //square it |
---|
| 2314 | retval *= retval; |
---|
| 2315 | |
---|
| 2316 | return(retval); |
---|
| 2317 | |
---|
| 2318 | }//Function EllipCylKernel() |
---|
| 2319 | |
---|
| 2320 | double NR_BessJ1(double x) |
---|
| 2321 | { |
---|
| 2322 | double ax,z; |
---|
| 2323 | double xx,y,ans,ans1,ans2; |
---|
| 2324 | |
---|
| 2325 | if ((ax=fabs(x)) < 8.0) { |
---|
| 2326 | y=x*x; |
---|
| 2327 | ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1 |
---|
| 2328 | +y*(-2972611.439+y*(15704.48260+y*(-30.16036606)))))); |
---|
| 2329 | ans2=144725228442.0+y*(2300535178.0+y*(18583304.74 |
---|
| 2330 | +y*(99447.43394+y*(376.9991397+y*1.0)))); |
---|
| 2331 | ans=ans1/ans2; |
---|
| 2332 | } else { |
---|
| 2333 | z=8.0/ax; |
---|
| 2334 | y=z*z; |
---|
| 2335 | xx=ax-2.356194491; |
---|
| 2336 | ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4 |
---|
| 2337 | +y*(0.2457520174e-5+y*(-0.240337019e-6)))); |
---|
| 2338 | ans2=0.04687499995+y*(-0.2002690873e-3 |
---|
| 2339 | +y*(0.8449199096e-5+y*(-0.88228987e-6 |
---|
| 2340 | +y*0.105787412e-6))); |
---|
| 2341 | ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2); |
---|
| 2342 | if (x < 0.0) ans = -ans; |
---|
| 2343 | } |
---|
| 2344 | |
---|
| 2345 | return(ans); |
---|
| 2346 | } |
---|
[453] | 2347 | |
---|
| 2348 | /* Lamellar_ParaCrystal - Pedersen's model |
---|
| 2349 | |
---|
| 2350 | */ |
---|
| 2351 | double |
---|
| 2352 | Lamellar_ParaCrystal(double w[], double q) |
---|
| 2353 | { |
---|
| 2354 | // Input (fitting) variables are: |
---|
| 2355 | //[0] scale factor |
---|
| 2356 | //[1] thickness |
---|
| 2357 | //[2] number of layers |
---|
| 2358 | //[3] spacing between layers |
---|
| 2359 | //[4] polydispersity of spacing |
---|
| 2360 | //[5] SLD lamellar |
---|
| 2361 | //[6] SLD solvent |
---|
| 2362 | //[7] incoherent background |
---|
| 2363 | // give them nice names |
---|
| 2364 | double inten,qval,scale,th,nl,davg,pd,contr,bkg,xn; |
---|
| 2365 | double xi,ww,Pbil,Znq,Snq,an,sldLayer,sldSolvent,pi; |
---|
| 2366 | long n1,n2; |
---|
| 2367 | |
---|
| 2368 | pi = 4.0*atan(1.0); |
---|
| 2369 | scale = w[0]; |
---|
| 2370 | th = w[1]; |
---|
| 2371 | nl = w[2]; |
---|
| 2372 | davg = w[3]; |
---|
| 2373 | pd = w[4]; |
---|
| 2374 | sldLayer = w[5]; |
---|
| 2375 | sldSolvent = w[6]; |
---|
| 2376 | bkg = w[7]; |
---|
| 2377 | |
---|
| 2378 | contr = w[5] - w[6]; |
---|
| 2379 | qval = q; |
---|
| 2380 | |
---|
| 2381 | //get the fractional part of nl, to determine the "mixing" of N's |
---|
| 2382 | |
---|
| 2383 | n1 = trunc(nl); //rounds towards zero |
---|
| 2384 | n2 = n1 + 1; |
---|
| 2385 | xn = (double)n2 - nl; //fractional contribution of n1 |
---|
| 2386 | |
---|
| 2387 | ww = exp(-qval*qval*pd*pd*davg*davg/2.0); |
---|
| 2388 | |
---|
| 2389 | //calculate the n1 contribution |
---|
| 2390 | an = paraCryst_an(ww,qval,davg,n1); |
---|
| 2391 | Snq = paraCryst_sn(ww,qval,davg,n1,an); |
---|
| 2392 | |
---|
| 2393 | Znq = xn*Snq; |
---|
| 2394 | |
---|
| 2395 | //calculate the n2 contribution |
---|
| 2396 | an = paraCryst_an(ww,qval,davg,n2); |
---|
| 2397 | Snq = paraCryst_sn(ww,qval,davg,n2,an); |
---|
| 2398 | |
---|
| 2399 | Znq += (1.0-xn)*Snq; |
---|
| 2400 | |
---|
| 2401 | //and the independent contribution |
---|
| 2402 | Znq += (1.0-ww*ww)/(1.0+ww*ww-2.0*ww*cos(qval*davg)); |
---|
| 2403 | |
---|
| 2404 | //the limit when NL approaches infinity |
---|
| 2405 | // Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg)) |
---|
| 2406 | |
---|
| 2407 | xi = th/2.0; //use 1/2 the bilayer thickness |
---|
| 2408 | Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi)); |
---|
| 2409 | |
---|
| 2410 | inten = 2.0*pi*contr*contr*Pbil*Znq/(qval*qval); |
---|
| 2411 | inten *= 1.0e8; |
---|
| 2412 | |
---|
| 2413 | return(scale*inten+bkg); |
---|
| 2414 | } |
---|
| 2415 | |
---|
| 2416 | // functions for the lamellar paracrystal model |
---|
| 2417 | double |
---|
| 2418 | paraCryst_sn(double ww, double qval, double davg, long nl, double an) { |
---|
| 2419 | |
---|
| 2420 | double Snq; |
---|
| 2421 | |
---|
| 2422 | Snq = an/( (double)nl*pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) ); |
---|
| 2423 | |
---|
| 2424 | return(Snq); |
---|
| 2425 | } |
---|
| 2426 | |
---|
| 2427 | |
---|
| 2428 | double |
---|
| 2429 | paraCryst_an(double ww, double qval, double davg, long nl) { |
---|
| 2430 | |
---|
| 2431 | double an; |
---|
| 2432 | |
---|
| 2433 | an = 4.0*ww*ww - 2.0*(ww*ww*ww+ww)*cos(qval*davg); |
---|
| 2434 | an -= 4.0*pow(ww,(nl+2))*cos((double)nl*qval*davg); |
---|
| 2435 | an += 2.0*pow(ww,(nl+3))*cos((double)(nl-1)*qval*davg); |
---|
| 2436 | an += 2.0*pow(ww,(nl+1))*cos((double)(nl+1)*qval*davg); |
---|
| 2437 | |
---|
| 2438 | return(an); |
---|
| 2439 | } |
---|
| 2440 | |
---|
| 2441 | |
---|
| 2442 | /* Spherocylinder : |
---|
| 2443 | |
---|
| 2444 | Uses 76 pt Gaussian quadrature for both integrals |
---|
| 2445 | */ |
---|
| 2446 | double |
---|
| 2447 | Spherocylinder(double w[], double x) |
---|
| 2448 | { |
---|
| 2449 | int i,j; |
---|
| 2450 | double Pi; |
---|
| 2451 | double scale,contr,bkg,sldc,slds; |
---|
| 2452 | double len,rad,hDist,endRad; |
---|
| 2453 | int nordi=76; //order of integration |
---|
| 2454 | int nordj=76; |
---|
| 2455 | double va,vb; //upper and lower integration limits |
---|
| 2456 | double summ,zi,yyy,answer; //running tally of integration |
---|
| 2457 | double summj,vaj,vbj,zij; //for the inner integration |
---|
| 2458 | double SphCyl_tmp[7],arg1,arg2,inner; |
---|
| 2459 | |
---|
| 2460 | |
---|
| 2461 | scale = w[0]; |
---|
| 2462 | rad = w[1]; |
---|
| 2463 | len = w[2]; |
---|
| 2464 | sldc = w[3]; |
---|
| 2465 | slds = w[4]; |
---|
| 2466 | bkg = w[5]; |
---|
| 2467 | |
---|
| 2468 | SphCyl_tmp[0] = w[0]; |
---|
| 2469 | SphCyl_tmp[1] = w[1]; |
---|
| 2470 | SphCyl_tmp[2] = w[2]; |
---|
| 2471 | SphCyl_tmp[3] = w[1]; //end radius is same as cylinder radius |
---|
| 2472 | SphCyl_tmp[4] = w[3]; |
---|
| 2473 | SphCyl_tmp[5] = w[4]; |
---|
| 2474 | SphCyl_tmp[6] = w[5]; |
---|
| 2475 | |
---|
| 2476 | hDist = 0; //by definition for this model |
---|
| 2477 | endRad = rad; |
---|
| 2478 | |
---|
| 2479 | contr = sldc-slds; |
---|
| 2480 | |
---|
| 2481 | Pi = 4.0*atan(1.0); |
---|
| 2482 | va = 0.0; |
---|
| 2483 | vb = Pi/2.0; //orintational average, outer integral |
---|
| 2484 | vaj = -1.0*hDist/endRad; |
---|
| 2485 | vbj = 1.0; //endpoints of inner integral |
---|
| 2486 | |
---|
| 2487 | summ = 0.0; //initialize intergral |
---|
| 2488 | |
---|
| 2489 | for(i=0;i<nordi;i++) { |
---|
| 2490 | //setup inner integral over the ellipsoidal cross-section |
---|
| 2491 | summj=0.0; |
---|
| 2492 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy |
---|
| 2493 | |
---|
| 2494 | for(j=0;j<nordj;j++) { |
---|
| 2495 | //20 gauss points for the inner integral |
---|
| 2496 | zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy |
---|
| 2497 | yyy = Gauss76Wt[j] * SphCyl_kernel(SphCyl_tmp,x,zij,zi); |
---|
| 2498 | summj += yyy; |
---|
| 2499 | } |
---|
| 2500 | //now calculate the value of the inner integral |
---|
| 2501 | inner = (vbj-vaj)/2.0*summj; |
---|
| 2502 | inner *= 4.0*Pi*endRad*endRad*endRad; |
---|
| 2503 | |
---|
| 2504 | //now calculate outer integrand |
---|
| 2505 | arg1 = x*len/2.0*cos(zi); |
---|
| 2506 | arg2 = x*rad*sin(zi); |
---|
| 2507 | yyy = inner; |
---|
| 2508 | |
---|
[632] | 2509 | if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h |
---|
[453] | 2510 | yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; |
---|
| 2511 | } else { |
---|
| 2512 | yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; |
---|
| 2513 | } |
---|
| 2514 | yyy *= yyy; |
---|
| 2515 | yyy *= sin(zi); // = |A(q)|^2*sin(theta) |
---|
| 2516 | yyy *= Gauss76Wt[i]; |
---|
| 2517 | summ += yyy; |
---|
| 2518 | } //final scaling is done at the end of the function, after the NT_FP64 case |
---|
| 2519 | |
---|
| 2520 | answer = (vb-va)/2.0*summ; |
---|
| 2521 | |
---|
| 2522 | answer /= Pi*rad*rad*len + Pi*4.0*endRad*endRad*endRad/3.0; //divide by volume |
---|
| 2523 | answer *= 1.0e8; //convert to cm^-1 |
---|
| 2524 | answer *= contr*contr; |
---|
| 2525 | answer *= scale; |
---|
| 2526 | answer += bkg; |
---|
| 2527 | |
---|
| 2528 | return answer; |
---|
| 2529 | } |
---|
| 2530 | |
---|
| 2531 | |
---|
| 2532 | // inner integral of the sphereocylinder model, special case of hDist = 0 |
---|
| 2533 | // |
---|
| 2534 | double |
---|
| 2535 | SphCyl_kernel(double w[], double x, double tt, double theta) { |
---|
| 2536 | |
---|
| 2537 | double val,arg1,arg2; |
---|
| 2538 | double scale,bkg,sldc,slds; |
---|
| 2539 | double len,rad,hDist,endRad; |
---|
| 2540 | scale = w[0]; |
---|
| 2541 | rad = w[1]; |
---|
| 2542 | len = w[2]; |
---|
| 2543 | endRad = w[3]; |
---|
| 2544 | sldc = w[4]; |
---|
| 2545 | slds = w[5]; |
---|
| 2546 | bkg = w[6]; |
---|
| 2547 | |
---|
| 2548 | hDist = 0.0; |
---|
| 2549 | |
---|
| 2550 | arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0); |
---|
| 2551 | arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt); |
---|
| 2552 | |
---|
| 2553 | val = cos(arg1)*(1.0-tt*tt)*NR_BessJ1(arg2)/arg2; |
---|
| 2554 | |
---|
| 2555 | return(val); |
---|
| 2556 | } |
---|
| 2557 | |
---|
| 2558 | |
---|
| 2559 | /* Convex Lens : special case where L ~ 0 and hDist < 0 |
---|
| 2560 | |
---|
| 2561 | Uses 76 pt Gaussian quadrature for both integrals |
---|
| 2562 | */ |
---|
| 2563 | double |
---|
| 2564 | ConvexLens(double w[], double x) |
---|
| 2565 | { |
---|
| 2566 | int i,j; |
---|
| 2567 | double Pi; |
---|
| 2568 | double scale,contr,bkg,sldc,slds; |
---|
| 2569 | double len,rad,hDist,endRad; |
---|
| 2570 | int nordi=76; //order of integration |
---|
| 2571 | int nordj=76; |
---|
| 2572 | double va,vb; //upper and lower integration limits |
---|
| 2573 | double summ,zi,yyy,answer; //running tally of integration |
---|
| 2574 | double summj,vaj,vbj,zij; //for the inner integration |
---|
| 2575 | double CLens_tmp[7],arg1,arg2,inner,hh; |
---|
| 2576 | |
---|
| 2577 | |
---|
| 2578 | scale = w[0]; |
---|
| 2579 | rad = w[1]; |
---|
| 2580 | // len = w[2] |
---|
| 2581 | endRad = w[2]; |
---|
| 2582 | sldc = w[3]; |
---|
| 2583 | slds = w[4]; |
---|
| 2584 | bkg = w[5]; |
---|
| 2585 | |
---|
| 2586 | len = 0.01; |
---|
| 2587 | |
---|
| 2588 | CLens_tmp[0] = w[0]; |
---|
| 2589 | CLens_tmp[1] = w[1]; |
---|
| 2590 | CLens_tmp[2] = len; //length is some small number, essentially zero |
---|
| 2591 | CLens_tmp[3] = w[2]; |
---|
| 2592 | CLens_tmp[4] = w[3]; |
---|
| 2593 | CLens_tmp[5] = w[4]; |
---|
| 2594 | CLens_tmp[6] = w[5]; |
---|
| 2595 | |
---|
| 2596 | hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad)); //by definition for this model |
---|
| 2597 | |
---|
| 2598 | contr = sldc-slds; |
---|
| 2599 | |
---|
| 2600 | Pi = 4.0*atan(1.0); |
---|
| 2601 | va = 0.0; |
---|
| 2602 | vb = Pi/2.0; //orintational average, outer integral |
---|
| 2603 | vaj = -1.0*hDist/endRad; |
---|
| 2604 | vbj = 1.0; //endpoints of inner integral |
---|
| 2605 | |
---|
| 2606 | summ = 0.0; //initialize intergral |
---|
| 2607 | |
---|
| 2608 | for(i=0;i<nordi;i++) { |
---|
| 2609 | //setup inner integral over the ellipsoidal cross-section |
---|
| 2610 | summj=0.0; |
---|
| 2611 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy |
---|
| 2612 | |
---|
| 2613 | for(j=0;j<nordj;j++) { |
---|
| 2614 | //20 gauss points for the inner integral |
---|
| 2615 | zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy |
---|
| 2616 | yyy = Gauss76Wt[j] * ConvLens_kernel(CLens_tmp,x,zij,zi); |
---|
| 2617 | summj += yyy; |
---|
| 2618 | } |
---|
| 2619 | //now calculate the value of the inner integral |
---|
| 2620 | inner = (vbj-vaj)/2.0*summj; |
---|
| 2621 | inner *= 4.0*Pi*endRad*endRad*endRad; |
---|
| 2622 | |
---|
| 2623 | //now calculate outer integrand |
---|
| 2624 | arg1 = x*len/2.0*cos(zi); |
---|
| 2625 | arg2 = x*rad*sin(zi); |
---|
| 2626 | yyy = inner; |
---|
| 2627 | |
---|
[632] | 2628 | if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h |
---|
[453] | 2629 | yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; |
---|
| 2630 | } else { |
---|
| 2631 | yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; |
---|
| 2632 | } |
---|
| 2633 | yyy *= yyy; |
---|
| 2634 | yyy *= sin(zi); // = |A(q)|^2*sin(theta) |
---|
| 2635 | yyy *= Gauss76Wt[i]; |
---|
| 2636 | summ += yyy; |
---|
| 2637 | } //final scaling is done at the end of the function, after the NT_FP64 case |
---|
| 2638 | |
---|
| 2639 | answer = (vb-va)/2.0*summ; |
---|
| 2640 | |
---|
| 2641 | hh = fabs(hDist); //need positive value for spherical cap volume |
---|
| 2642 | answer /= 2.0*(1.0/3.0*Pi*(endRad-hh)*(endRad-hh)*(2.0*endRad+hh)); //divide by volume |
---|
| 2643 | answer *= 1.0e8; //convert to cm^-1 |
---|
| 2644 | answer *= contr*contr; |
---|
| 2645 | answer *= scale; |
---|
| 2646 | answer += bkg; |
---|
| 2647 | |
---|
| 2648 | return answer; |
---|
| 2649 | } |
---|
| 2650 | |
---|
| 2651 | /* Capped Cylinder : special case where L is nonzero and hDist < 0 |
---|
| 2652 | |
---|
| 2653 | -- uses the same Kernel as the Convex Lens |
---|
| 2654 | |
---|
| 2655 | Uses 76 pt Gaussian quadrature for both integrals |
---|
| 2656 | */ |
---|
| 2657 | double |
---|
| 2658 | CappedCylinder(double w[], double x) |
---|
| 2659 | { |
---|
| 2660 | int i,j; |
---|
| 2661 | double Pi; |
---|
| 2662 | double scale,contr,bkg,sldc,slds; |
---|
| 2663 | double len,rad,hDist,endRad; |
---|
| 2664 | int nordi=76; //order of integration |
---|
| 2665 | int nordj=76; |
---|
| 2666 | double va,vb; //upper and lower integration limits |
---|
| 2667 | double summ,zi,yyy,answer; //running tally of integration |
---|
| 2668 | double summj,vaj,vbj,zij; //for the inner integration |
---|
| 2669 | double arg1,arg2,inner,hh; |
---|
| 2670 | |
---|
| 2671 | |
---|
| 2672 | scale = w[0]; |
---|
| 2673 | rad = w[1]; |
---|
| 2674 | len = w[2]; |
---|
| 2675 | endRad = w[3]; |
---|
| 2676 | sldc = w[4]; |
---|
| 2677 | slds = w[5]; |
---|
| 2678 | bkg = w[6]; |
---|
| 2679 | |
---|
| 2680 | hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad)); //by definition for this model |
---|
| 2681 | |
---|
| 2682 | contr = sldc-slds; |
---|
| 2683 | |
---|
| 2684 | Pi = 4.0*atan(1.0); |
---|
| 2685 | va = 0.0; |
---|
| 2686 | vb = Pi/2.0; //orintational average, outer integral |
---|
| 2687 | vaj = -1.0*hDist/endRad; |
---|
| 2688 | vbj = 1.0; //endpoints of inner integral |
---|
| 2689 | |
---|
| 2690 | summ = 0.0; //initialize intergral |
---|
| 2691 | |
---|
| 2692 | for(i=0;i<nordi;i++) { |
---|
| 2693 | //setup inner integral over the ellipsoidal cross-section |
---|
| 2694 | summj=0.0; |
---|
| 2695 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy |
---|
| 2696 | |
---|
| 2697 | for(j=0;j<nordj;j++) { |
---|
| 2698 | //20 gauss points for the inner integral |
---|
| 2699 | zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy |
---|
| 2700 | yyy = Gauss76Wt[j] * ConvLens_kernel(w,x,zij,zi); //uses the same kernel as ConvexLens, except here L != 0 |
---|
| 2701 | summj += yyy; |
---|
| 2702 | } |
---|
| 2703 | //now calculate the value of the inner integral |
---|
| 2704 | inner = (vbj-vaj)/2.0*summj; |
---|
| 2705 | inner *= 4.0*Pi*endRad*endRad*endRad; |
---|
| 2706 | |
---|
| 2707 | //now calculate outer integrand |
---|
| 2708 | arg1 = x*len/2.0*cos(zi); |
---|
| 2709 | arg2 = x*rad*sin(zi); |
---|
| 2710 | yyy = inner; |
---|
| 2711 | |
---|
[632] | 2712 | if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h |
---|
[453] | 2713 | yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; |
---|
| 2714 | } else { |
---|
| 2715 | yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; |
---|
| 2716 | } |
---|
| 2717 | yyy *= yyy; |
---|
| 2718 | yyy *= sin(zi); // = |A(q)|^2*sin(theta) |
---|
| 2719 | yyy *= Gauss76Wt[i]; |
---|
| 2720 | summ += yyy; |
---|
| 2721 | } //final scaling is done at the end of the function, after the NT_FP64 case |
---|
| 2722 | |
---|
| 2723 | answer = (vb-va)/2.0*summ; |
---|
| 2724 | |
---|
| 2725 | hh = fabs(hDist); //need positive value for spherical cap volume |
---|
| 2726 | answer /= Pi*rad*rad*len + 2.0*(1.0/3.0*Pi*(endRad-hh)*(endRad-hh)*(2.0*endRad+hh)); //divide by volume |
---|
| 2727 | answer *= 1.0e8; //convert to cm^-1 |
---|
| 2728 | answer *= contr*contr; |
---|
| 2729 | answer *= scale; |
---|
| 2730 | answer += bkg; |
---|
| 2731 | |
---|
| 2732 | return answer; |
---|
| 2733 | } |
---|
| 2734 | |
---|
| 2735 | |
---|
| 2736 | |
---|
| 2737 | // inner integral of the ConvexLens model, special case where L ~ 0 and hDist < 0 |
---|
| 2738 | // |
---|
| 2739 | double |
---|
| 2740 | ConvLens_kernel(double w[], double x, double tt, double theta) { |
---|
| 2741 | |
---|
| 2742 | double val,arg1,arg2; |
---|
| 2743 | double scale,bkg,sldc,slds; |
---|
| 2744 | double len,rad,hDist,endRad; |
---|
| 2745 | scale = w[0]; |
---|
| 2746 | rad = w[1]; |
---|
| 2747 | len = w[2]; |
---|
| 2748 | endRad = w[3]; |
---|
| 2749 | sldc = w[4]; |
---|
| 2750 | slds = w[5]; |
---|
| 2751 | bkg = w[6]; |
---|
| 2752 | |
---|
| 2753 | hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad)); |
---|
| 2754 | |
---|
| 2755 | arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0); |
---|
| 2756 | arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt); |
---|
| 2757 | |
---|
| 2758 | val = cos(arg1)*(1.0-tt*tt)*NR_BessJ1(arg2)/arg2; |
---|
| 2759 | |
---|
| 2760 | return(val); |
---|
| 2761 | } |
---|
| 2762 | |
---|
| 2763 | |
---|
| 2764 | /* Dumbbell : special case where L ~ 0 and hDist > 0 |
---|
| 2765 | |
---|
| 2766 | Uses 76 pt Gaussian quadrature for both integrals |
---|
| 2767 | */ |
---|
| 2768 | double |
---|
| 2769 | Dumbbell(double w[], double x) |
---|
| 2770 | { |
---|
| 2771 | int i,j; |
---|
| 2772 | double Pi; |
---|
| 2773 | double scale,contr,bkg,sldc,slds; |
---|
| 2774 | double len,rad,hDist,endRad; |
---|
| 2775 | int nordi=76; //order of integration |
---|
| 2776 | int nordj=76; |
---|
| 2777 | double va,vb; //upper and lower integration limits |
---|
| 2778 | double summ,zi,yyy,answer; //running tally of integration |
---|
| 2779 | double summj,vaj,vbj,zij; //for the inner integration |
---|
| 2780 | double Dumb_tmp[7],arg1,arg2,inner; |
---|
| 2781 | |
---|
| 2782 | |
---|
| 2783 | scale = w[0]; |
---|
| 2784 | rad = w[1]; |
---|
| 2785 | // len = w[2] |
---|
| 2786 | endRad = w[2]; |
---|
| 2787 | sldc = w[3]; |
---|
| 2788 | slds = w[4]; |
---|
| 2789 | bkg = w[5]; |
---|
| 2790 | |
---|
| 2791 | len = 0.01; |
---|
| 2792 | |
---|
| 2793 | Dumb_tmp[0] = w[0]; |
---|
| 2794 | Dumb_tmp[1] = w[1]; |
---|
| 2795 | Dumb_tmp[2] = len; //length is some small number, essentially zero |
---|
| 2796 | Dumb_tmp[3] = w[2]; |
---|
| 2797 | Dumb_tmp[4] = w[3]; |
---|
| 2798 | Dumb_tmp[5] = w[4]; |
---|
| 2799 | Dumb_tmp[6] = w[5]; |
---|
| 2800 | |
---|
| 2801 | hDist = sqrt(fabs(endRad*endRad-rad*rad)); //by definition for this model |
---|
| 2802 | |
---|
| 2803 | contr = sldc-slds; |
---|
| 2804 | |
---|
| 2805 | Pi = 4.0*atan(1.0); |
---|
| 2806 | va = 0.0; |
---|
| 2807 | vb = Pi/2.0; //orintational average, outer integral |
---|
| 2808 | vaj = -1.0*hDist/endRad; |
---|
| 2809 | vbj = 1.0; //endpoints of inner integral |
---|
| 2810 | |
---|
| 2811 | summ = 0.0; //initialize intergral |
---|
| 2812 | |
---|
| 2813 | for(i=0;i<nordi;i++) { |
---|
| 2814 | //setup inner integral over the ellipsoidal cross-section |
---|
| 2815 | summj=0.0; |
---|
| 2816 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy |
---|
| 2817 | |
---|
| 2818 | for(j=0;j<nordj;j++) { |
---|
| 2819 | //20 gauss points for the inner integral |
---|
| 2820 | zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy |
---|
| 2821 | yyy = Gauss76Wt[j] * Dumb_kernel(Dumb_tmp,x,zij,zi); |
---|
| 2822 | summj += yyy; |
---|
| 2823 | } |
---|
| 2824 | //now calculate the value of the inner integral |
---|
| 2825 | inner = (vbj-vaj)/2.0*summj; |
---|
| 2826 | inner *= 4.0*Pi*endRad*endRad*endRad; |
---|
| 2827 | |
---|
| 2828 | //now calculate outer integrand |
---|
| 2829 | arg1 = x*len/2.0*cos(zi); |
---|
| 2830 | arg2 = x*rad*sin(zi); |
---|
| 2831 | yyy = inner; |
---|
| 2832 | |
---|
[632] | 2833 | if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h |
---|
[453] | 2834 | yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; |
---|
| 2835 | } else { |
---|
| 2836 | yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; |
---|
| 2837 | } |
---|
| 2838 | yyy *= yyy; |
---|
| 2839 | yyy *= sin(zi); // = |A(q)|^2*sin(theta) |
---|
| 2840 | yyy *= Gauss76Wt[i]; |
---|
| 2841 | summ += yyy; |
---|
| 2842 | } //final scaling is done at the end of the function, after the NT_FP64 case |
---|
| 2843 | |
---|
| 2844 | answer = (vb-va)/2.0*summ; |
---|
| 2845 | |
---|
| 2846 | answer /= 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0); //divide by volume |
---|
| 2847 | answer *= 1.0e8; //convert to cm^-1 |
---|
| 2848 | answer *= contr*contr; |
---|
| 2849 | answer *= scale; |
---|
| 2850 | answer += bkg; |
---|
| 2851 | |
---|
| 2852 | return answer; |
---|
| 2853 | } |
---|
| 2854 | |
---|
| 2855 | |
---|
| 2856 | /* Barbell : "normal" case where L is nonzero 0 and hDist > 0 |
---|
| 2857 | |
---|
| 2858 | -- uses the same kernel as the Dumbbell case |
---|
| 2859 | |
---|
| 2860 | Uses 76 pt Gaussian quadrature for both integrals |
---|
| 2861 | */ |
---|
| 2862 | double |
---|
| 2863 | Barbell(double w[], double x) |
---|
| 2864 | { |
---|
| 2865 | int i,j; |
---|
| 2866 | double Pi; |
---|
| 2867 | double scale,contr,bkg,sldc,slds; |
---|
| 2868 | double len,rad,hDist,endRad; |
---|
| 2869 | int nordi=76; //order of integration |
---|
| 2870 | int nordj=76; |
---|
| 2871 | double va,vb; //upper and lower integration limits |
---|
| 2872 | double summ,zi,yyy,answer; //running tally of integration |
---|
| 2873 | double summj,vaj,vbj,zij; //for the inner integration |
---|
| 2874 | double arg1,arg2,inner; |
---|
| 2875 | |
---|
| 2876 | |
---|
| 2877 | scale = w[0]; |
---|
| 2878 | rad = w[1]; |
---|
| 2879 | len = w[2]; |
---|
| 2880 | endRad = w[3]; |
---|
| 2881 | sldc = w[4]; |
---|
| 2882 | slds = w[5]; |
---|
| 2883 | bkg = w[6]; |
---|
| 2884 | |
---|
| 2885 | hDist = sqrt(fabs(endRad*endRad-rad*rad)); //by definition for this model |
---|
| 2886 | |
---|
| 2887 | contr = sldc-slds; |
---|
| 2888 | |
---|
| 2889 | Pi = 4.0*atan(1.0); |
---|
| 2890 | va = 0.0; |
---|
| 2891 | vb = Pi/2.0; //orintational average, outer integral |
---|
| 2892 | vaj = -1.0*hDist/endRad; |
---|
| 2893 | vbj = 1.0; //endpoints of inner integral |
---|
| 2894 | |
---|
| 2895 | summ = 0.0; //initialize intergral |
---|
| 2896 | |
---|
| 2897 | for(i=0;i<nordi;i++) { |
---|
| 2898 | //setup inner integral over the ellipsoidal cross-section |
---|
| 2899 | summj=0.0; |
---|
| 2900 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "theta" dummy |
---|
| 2901 | |
---|
| 2902 | for(j=0;j<nordj;j++) { |
---|
| 2903 | //20 gauss points for the inner integral |
---|
| 2904 | zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy |
---|
| 2905 | yyy = Gauss76Wt[j] * Dumb_kernel(w,x,zij,zi); //uses the same Kernel as the Dumbbell, here L>0 |
---|
| 2906 | summj += yyy; |
---|
| 2907 | } |
---|
| 2908 | //now calculate the value of the inner integral |
---|
| 2909 | inner = (vbj-vaj)/2.0*summj; |
---|
| 2910 | inner *= 4.0*Pi*endRad*endRad*endRad; |
---|
| 2911 | |
---|
| 2912 | //now calculate outer integrand |
---|
| 2913 | arg1 = x*len/2.0*cos(zi); |
---|
| 2914 | arg2 = x*rad*sin(zi); |
---|
| 2915 | yyy = inner; |
---|
| 2916 | |
---|
[632] | 2917 | if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h |
---|
[453] | 2918 | yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; |
---|
| 2919 | } else { |
---|
| 2920 | yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; |
---|
| 2921 | } |
---|
| 2922 | yyy *= yyy; |
---|
| 2923 | yyy *= sin(zi); // = |A(q)|^2*sin(theta) |
---|
| 2924 | yyy *= Gauss76Wt[i]; |
---|
| 2925 | summ += yyy; |
---|
| 2926 | } //final scaling is done at the end of the function, after the NT_FP64 case |
---|
| 2927 | |
---|
| 2928 | answer = (vb-va)/2.0*summ; |
---|
| 2929 | |
---|
| 2930 | answer /= Pi*rad*rad*len + 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0); //divide by volume |
---|
| 2931 | answer *= 1.0e8; //convert to cm^-1 |
---|
| 2932 | answer *= contr*contr; |
---|
| 2933 | answer *= scale; |
---|
| 2934 | answer += bkg; |
---|
| 2935 | |
---|
| 2936 | return answer; |
---|
| 2937 | } |
---|
| 2938 | |
---|
| 2939 | |
---|
| 2940 | |
---|
| 2941 | // inner integral of the Dumbbell model, special case where L ~ 0 and hDist > 0 |
---|
| 2942 | // |
---|
| 2943 | // inner integral of the Barbell model if L is nonzero |
---|
| 2944 | // |
---|
| 2945 | double |
---|
| 2946 | Dumb_kernel(double w[], double x, double tt, double theta) { |
---|
| 2947 | |
---|
| 2948 | double val,arg1,arg2; |
---|
| 2949 | double scale,bkg,sldc,slds; |
---|
| 2950 | double len,rad,hDist,endRad; |
---|
| 2951 | scale = w[0]; |
---|
| 2952 | rad = w[1]; |
---|
| 2953 | len = w[2]; |
---|
| 2954 | endRad = w[3]; |
---|
| 2955 | sldc = w[4]; |
---|
| 2956 | slds = w[5]; |
---|
| 2957 | bkg = w[6]; |
---|
| 2958 | |
---|
| 2959 | hDist = sqrt(fabs(endRad*endRad-rad*rad)); |
---|
| 2960 | |
---|
| 2961 | arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0); |
---|
| 2962 | arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt); |
---|
| 2963 | |
---|
| 2964 | val = cos(arg1)*(1.0-tt*tt)*NR_BessJ1(arg2)/arg2; |
---|
| 2965 | |
---|
| 2966 | return(val); |
---|
| 2967 | } |
---|
[500] | 2968 | |
---|
[501] | 2969 | double PolyCoreBicelle(double dp[], double q) |
---|
| 2970 | { |
---|
[500] | 2971 | int i; |
---|
| 2972 | int nord = 20; |
---|
| 2973 | double scale, length, sigma, bkg, radius, radthick, facthick; |
---|
| 2974 | double rhoc, rhoh, rhor, rhosolv; |
---|
| 2975 | double answer, Vpoly; |
---|
[501] | 2976 | double Pi,lolim,uplim,summ,yyy,rad,AR,Rsqr,Rsqrsumm,Rsqryyy; |
---|
[500] | 2977 | |
---|
[501] | 2978 | scale = dp[0]; |
---|
| 2979 | radius = dp[1]; |
---|
| 2980 | sigma = dp[2]; //sigma is the standard mean deviation |
---|
| 2981 | length = dp[3]; |
---|
| 2982 | radthick = dp[4]; |
---|
| 2983 | facthick= dp[5]; |
---|
| 2984 | rhoc = dp[6]; |
---|
| 2985 | rhoh = dp[7]; |
---|
| 2986 | rhor=dp[8]; |
---|
| 2987 | rhosolv = dp[9]; |
---|
| 2988 | bkg = dp[10]; |
---|
[500] | 2989 | |
---|
[501] | 2990 | Pi = 4.0*atan(1.0); |
---|
| 2991 | |
---|
[500] | 2992 | lolim = exp(log(radius)-(4.*sigma)); |
---|
[632] | 2993 | if (lolim<0.0) { |
---|
| 2994 | lolim=0.0; //to avoid numerical error when va<0 (-ve r value) |
---|
[500] | 2995 | } |
---|
| 2996 | uplim = exp(log(radius)+(4.*sigma)); |
---|
| 2997 | |
---|
[501] | 2998 | summ = 0.0; |
---|
| 2999 | Rsqrsumm = 0.0; |
---|
| 3000 | |
---|
[500] | 3001 | for(i=0;i<nord;i++) { |
---|
| 3002 | rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 3003 | AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp(-(0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma))); |
---|
| 3004 | yyy = AR* Gauss20Wt[i] * BicelleIntegration(q,rad,radthick,facthick,rhoc,rhoh,rhor,rhosolv,length); |
---|
| 3005 | Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick); //SRK normalize to total dimensions |
---|
| 3006 | summ += yyy; |
---|
| 3007 | Rsqrsumm += Rsqryyy; |
---|
| 3008 | } |
---|
| 3009 | |
---|
| 3010 | answer = (uplim-lolim)/2.0*summ; |
---|
| 3011 | Rsqr = (uplim-lolim)/2.0*Rsqrsumm; |
---|
| 3012 | //normalize by average cylinder volume |
---|
| 3013 | Vpoly = Pi*Rsqr*(length+2*facthick); |
---|
| 3014 | answer /= Vpoly; |
---|
| 3015 | //convert to [cm-1] |
---|
| 3016 | answer *= 1.0e8; |
---|
| 3017 | //Scale |
---|
| 3018 | answer *= scale; |
---|
| 3019 | // add in the background |
---|
| 3020 | answer += bkg; |
---|
[501] | 3021 | |
---|
| 3022 | return answer; |
---|
[500] | 3023 | |
---|
| 3024 | } |
---|
| 3025 | |
---|
| 3026 | double |
---|
| 3027 | BicelleIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length){ |
---|
| 3028 | |
---|
| 3029 | double answer,halfheight,Pi; |
---|
| 3030 | double lolim,uplim,summ,yyy,zi; |
---|
| 3031 | int nord,i; |
---|
| 3032 | |
---|
| 3033 | // set up the integration end points |
---|
| 3034 | Pi = 4.0*atan(1.0); |
---|
| 3035 | nord = 76; |
---|
[632] | 3036 | lolim = 0.0; |
---|
[500] | 3037 | uplim = Pi/2; |
---|
| 3038 | halfheight = length/2.0; |
---|
| 3039 | |
---|
| 3040 | summ = 0.0; // initialize integral |
---|
| 3041 | i=0; |
---|
| 3042 | for(i=0;i<nord;i++) { |
---|
| 3043 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 3044 | yyy = Gauss76Wt[i] * BicelleKernel(qq, rad, radthick, facthick, rhoc, rhoh, rhor,rhosolv, halfheight, zi); |
---|
| 3045 | summ += yyy; |
---|
| 3046 | } |
---|
| 3047 | |
---|
| 3048 | // calculate value of integral to return |
---|
| 3049 | answer = (uplim-lolim)/2.0*summ; |
---|
| 3050 | return(answer); |
---|
| 3051 | } |
---|
| 3052 | |
---|
| 3053 | double |
---|
[594] | 3054 | BicelleKernel(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length, double dum) |
---|
| 3055 | { |
---|
[500] | 3056 | double dr1,dr2,dr3; |
---|
| 3057 | double besarg1,besarg2; |
---|
| 3058 | double vol1,vol2,vol3; |
---|
| 3059 | double sinarg1,sinarg2; |
---|
| 3060 | double t1,t2,t3; |
---|
| 3061 | double retval; |
---|
| 3062 | |
---|
| 3063 | double Pi = 4.0*atan(1.0); |
---|
| 3064 | |
---|
| 3065 | dr1 = rhoc-rhoh; |
---|
| 3066 | dr2 = rhor-rhosolv; |
---|
| 3067 | dr3= rhoh-rhor; |
---|
[632] | 3068 | vol1 = Pi*rad*rad*(2.0*length); |
---|
| 3069 | vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick); |
---|
| 3070 | vol3= Pi*(rad)*(rad)*(2.0*length+2.0*facthick); |
---|
[500] | 3071 | besarg1 = qq*rad*sin(dum); |
---|
| 3072 | besarg2 = qq*(rad+radthick)*sin(dum); |
---|
| 3073 | sinarg1 = qq*length*cos(dum); |
---|
| 3074 | sinarg2 = qq*(length+facthick)*cos(dum); |
---|
| 3075 | |
---|
[632] | 3076 | t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1; |
---|
| 3077 | t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2; |
---|
| 3078 | t3 = 2.0*vol3*dr3*sin(sinarg2)/sinarg2*NR_BessJ1(besarg1)/besarg1; |
---|
[500] | 3079 | |
---|
| 3080 | retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum); |
---|
[501] | 3081 | return(retval); |
---|
[500] | 3082 | |
---|
| 3083 | } |
---|
[501] | 3084 | |
---|
| 3085 | |
---|
| 3086 | double |
---|
[594] | 3087 | CSPPKernel(double dp[], double mu, double uu) |
---|
| 3088 | { |
---|
[501] | 3089 | double aa,bb,cc, ta,tb,tc; |
---|
| 3090 | double Vin,Vot,V1,V2; |
---|
| 3091 | double rhoA,rhoB,rhoC, rhoP, rhosolv; |
---|
[594] | 3092 | double dr0, drA,drB, drC; |
---|
[501] | 3093 | double arg1,arg2,arg3,arg4,t1,t2, t3, t4; |
---|
[594] | 3094 | double Pi,retVal; |
---|
[501] | 3095 | |
---|
| 3096 | aa = dp[1]; |
---|
| 3097 | bb = dp[2]; |
---|
| 3098 | cc = dp[3]; |
---|
| 3099 | ta = dp[4]; |
---|
| 3100 | tb = dp[5]; |
---|
| 3101 | tc = dp[6]; |
---|
| 3102 | rhoA=dp[7]; |
---|
| 3103 | rhoB=dp[8]; |
---|
| 3104 | rhoC=dp[9]; |
---|
| 3105 | rhoP=dp[10]; |
---|
| 3106 | rhosolv=dp[11]; |
---|
| 3107 | dr0=rhoP-rhosolv; |
---|
| 3108 | drA=rhoA-rhosolv; |
---|
| 3109 | drB=rhoB-rhosolv; |
---|
| 3110 | drC=rhoC-rhosolv; |
---|
| 3111 | Vin=(aa*bb*cc); |
---|
[594] | 3112 | Vot=(aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc); |
---|
| 3113 | V1=(2.0*ta*bb*cc); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) |
---|
| 3114 | V2=(2.0*aa*tb*cc); // incorrect V2(aa*bb*cc+2*aa*tb*cc) |
---|
[501] | 3115 | aa = aa/bb; |
---|
[594] | 3116 | ta=(aa+2.0*ta)/bb; |
---|
| 3117 | tb=(aa+2.0*tb)/bb; |
---|
[501] | 3118 | |
---|
| 3119 | Pi = 4.0*atan(1.0); |
---|
| 3120 | |
---|
[594] | 3121 | arg1 = (mu*aa/2.0)*sin(Pi*uu/2.0); |
---|
| 3122 | arg2 = (mu/2.0)*cos(Pi*uu/2.0); |
---|
| 3123 | arg3= (mu*ta/2.0)*sin(Pi*uu/2.0); |
---|
| 3124 | arg4= (mu*tb/2.0)*cos(Pi*uu/2.0); |
---|
[501] | 3125 | |
---|
[632] | 3126 | if(arg1==0.0){ |
---|
[594] | 3127 | t1 = 1.0; |
---|
[501] | 3128 | } else { |
---|
| 3129 | t1 = (sin(arg1)/arg1); //defn for CSPP model sin(arg1)/arg1 test: (sin(arg1)/arg1)*(sin(arg1)/arg1) |
---|
| 3130 | } |
---|
[632] | 3131 | if(arg2==0.0){ |
---|
[594] | 3132 | t2 = 1.0; |
---|
[501] | 3133 | } else { |
---|
| 3134 | t2 = (sin(arg2)/arg2); //defn for CSPP model sin(arg2)/arg2 test: (sin(arg2)/arg2)*(sin(arg2)/arg2) |
---|
| 3135 | } |
---|
[632] | 3136 | if(arg3==0.0){ |
---|
[594] | 3137 | t3 = 1.0; |
---|
[501] | 3138 | } else { |
---|
| 3139 | t3 = sin(arg3)/arg3; |
---|
| 3140 | } |
---|
[632] | 3141 | if(arg4==0.0){ |
---|
[594] | 3142 | t4 = 1.0; |
---|
[501] | 3143 | } else { |
---|
| 3144 | t4 = sin(arg4)/arg4; |
---|
| 3145 | } |
---|
[594] | 3146 | retVal =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 ); // correct FF : square of sum of phase factors |
---|
[501] | 3147 | return(retVal); |
---|
[594] | 3148 | |
---|
| 3149 | } |
---|
| 3150 | |
---|
| 3151 | /* CSParallelepiped : calculates the form factor of a Parallelepiped with a core-shell structure |
---|
| 3152 | -- different SLDs can be used for the face and rim |
---|
| 3153 | |
---|
| 3154 | Uses 76 pt Gaussian quadrature for both integrals |
---|
| 3155 | */ |
---|
| 3156 | double |
---|
| 3157 | CSParallelepiped(double dp[], double q) |
---|
| 3158 | { |
---|
| 3159 | int i,j; |
---|
| 3160 | double scale,aa,bb,cc,ta,tb,tc,rhoA,rhoB,rhoC,rhoP,rhosolv,bkg; //local variables of coefficient wave |
---|
| 3161 | int nordi=76; //order of integration |
---|
| 3162 | int nordj=76; |
---|
| 3163 | double va,vb; //upper and lower integration limits |
---|
| 3164 | double summ,yyy,answer; //running tally of integration |
---|
| 3165 | double summj,vaj,vbj; //for the inner integration |
---|
| 3166 | double mu,mudum,arg,sigma,uu,vol; |
---|
[501] | 3167 | |
---|
[594] | 3168 | |
---|
| 3169 | // Pi = 4.0*atan(1.0); |
---|
| 3170 | va = 0.0; |
---|
| 3171 | vb = 1.0; //orintational average, outer integral |
---|
| 3172 | vaj = 0.0; |
---|
| 3173 | vbj = 1.0; //endpoints of inner integral |
---|
| 3174 | |
---|
| 3175 | summ = 0.0; //initialize intergral |
---|
| 3176 | |
---|
| 3177 | scale = dp[0]; |
---|
| 3178 | aa = dp[1]; |
---|
| 3179 | bb = dp[2]; |
---|
| 3180 | cc = dp[3]; |
---|
| 3181 | ta = dp[4]; |
---|
| 3182 | tb = dp[5]; |
---|
| 3183 | tc = dp[6]; // is 0 at the moment |
---|
| 3184 | rhoA=dp[7]; //rim A SLD |
---|
| 3185 | rhoB=dp[8]; //rim B SLD |
---|
| 3186 | rhoC=dp[9]; //rim C SLD |
---|
| 3187 | rhoP = dp[10]; //Parallelpiped core SLD |
---|
| 3188 | rhosolv=dp[11]; // Solvent SLD |
---|
| 3189 | bkg = dp[12]; |
---|
| 3190 | |
---|
| 3191 | mu = q*bb; |
---|
| 3192 | vol = aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc; //calculate volume before rescaling |
---|
| 3193 | |
---|
| 3194 | // do the rescaling here, not in the kernel |
---|
| 3195 | // normalize all WRT bb |
---|
| 3196 | aa = aa/bb; |
---|
| 3197 | cc = cc/bb; |
---|
| 3198 | |
---|
| 3199 | for(i=0;i<nordi;i++) { |
---|
| 3200 | //setup inner integral over the ellipsoidal cross-section |
---|
[632] | 3201 | summj=0.0; |
---|
[594] | 3202 | sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the outer dummy |
---|
| 3203 | |
---|
| 3204 | for(j=0;j<nordj;j++) { |
---|
| 3205 | //76 gauss points for the inner integral |
---|
| 3206 | uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the inner dummy |
---|
| 3207 | mudum = mu*sqrt(1.0-sigma*sigma); |
---|
| 3208 | yyy = Gauss76Wt[j] * CSPPKernel(dp,mudum,uu); |
---|
| 3209 | summj += yyy; |
---|
| 3210 | } |
---|
| 3211 | //now calculate the value of the inner integral |
---|
| 3212 | answer = (vbj-vaj)/2.0*summj; |
---|
| 3213 | |
---|
| 3214 | //finish the outer integral cc already scaled |
---|
| 3215 | arg = mu*cc*sigma/2.0; |
---|
[632] | 3216 | if ( arg == 0.0 ) { |
---|
[594] | 3217 | answer *= 1.0; |
---|
| 3218 | } else { |
---|
| 3219 | answer *= sin(arg)*sin(arg)/arg/arg; |
---|
| 3220 | } |
---|
| 3221 | |
---|
| 3222 | //now sum up the outer integral |
---|
| 3223 | yyy = Gauss76Wt[i] * answer; |
---|
| 3224 | summ += yyy; |
---|
| 3225 | } //final scaling is done at the end of the function, after the NT_FP64 case |
---|
| 3226 | |
---|
| 3227 | answer = (vb-va)/2.0*summ; |
---|
[501] | 3228 | |
---|
[594] | 3229 | //normalize by volume |
---|
| 3230 | answer /= vol; |
---|
| 3231 | //convert to [cm-1] |
---|
| 3232 | answer *= 1.0e8; |
---|
| 3233 | //Scale |
---|
| 3234 | answer *= scale; |
---|
| 3235 | // add in the background |
---|
| 3236 | answer += bkg; |
---|
| 3237 | |
---|
| 3238 | return answer; |
---|
| 3239 | } |
---|
[501] | 3240 | |
---|