[ae3ce4e] | 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 | /* CylinderForm : calculates the form factor of a cylinder at the give x-value p->x |
---|
| 11 | |
---|
| 12 | Warning: |
---|
| 13 | The call to WaveData() below returns a pointer to the middle |
---|
| 14 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 15 | calculations could cause memory to move, you should copy the coefficient |
---|
| 16 | values to local variables or an array before such operations. |
---|
| 17 | */ |
---|
| 18 | double |
---|
| 19 | CylinderForm(double dp[], double q) |
---|
| 20 | { |
---|
| 21 | int i; |
---|
| 22 | double Pi; |
---|
| 23 | double scale,radius,length,delrho,bkg,halfheight; //local variables of coefficient wave |
---|
| 24 | int nord=76; //order of integration |
---|
| 25 | double uplim,lolim; //upper and lower integration limits |
---|
| 26 | double summ,zi,yyy,answer,vcyl; //running tally of integration |
---|
[8e91f01] | 27 | |
---|
[ae3ce4e] | 28 | Pi = 4.0*atan(1.0); |
---|
| 29 | lolim = 0; |
---|
| 30 | uplim = Pi/2.0; |
---|
[8e91f01] | 31 | |
---|
[ae3ce4e] | 32 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 33 | |
---|
[ae3ce4e] | 34 | scale = dp[0]; //make local copies in case memory moves |
---|
| 35 | radius = dp[1]; |
---|
| 36 | length = dp[2]; |
---|
| 37 | delrho = dp[3]; |
---|
| 38 | bkg = dp[4]; |
---|
| 39 | halfheight = length/2.0; |
---|
| 40 | for(i=0;i<nord;i++) { |
---|
| 41 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 42 | yyy = Gauss76Wt[i] * CylKernel(q, radius, halfheight, zi); |
---|
| 43 | summ += yyy; |
---|
| 44 | } |
---|
[8e91f01] | 45 | |
---|
[ae3ce4e] | 46 | answer = (uplim-lolim)/2.0*summ; |
---|
| 47 | // Multiply by contrast^2 |
---|
| 48 | answer *= delrho*delrho; |
---|
| 49 | //normalize by cylinder volume |
---|
| 50 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
| 51 | vcyl=Pi*radius*radius*length; |
---|
| 52 | answer *= vcyl; |
---|
| 53 | //convert to [cm-1] |
---|
| 54 | answer *= 1.0e8; |
---|
| 55 | //Scale |
---|
| 56 | answer *= scale; |
---|
| 57 | // add in the background |
---|
| 58 | answer += bkg; |
---|
[8e91f01] | 59 | |
---|
[ae3ce4e] | 60 | return answer; |
---|
| 61 | } |
---|
| 62 | |
---|
| 63 | /* EllipCyl76X : calculates the form factor of a elliptical cylinder at the given x-value p->x |
---|
| 64 | |
---|
| 65 | Uses 76 pt Gaussian quadrature for both integrals |
---|
| 66 | |
---|
| 67 | Warning: |
---|
| 68 | The call to WaveData() below returns a pointer to the middle |
---|
| 69 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 70 | calculations could cause memory to move, you should copy the coefficient |
---|
| 71 | values to local variables or an array before such operations. |
---|
| 72 | */ |
---|
| 73 | double |
---|
| 74 | EllipCyl76(double dp[], double q) |
---|
| 75 | { |
---|
| 76 | int i,j; |
---|
| 77 | double Pi; |
---|
| 78 | double scale,ra,nu,length,delrho,bkg; //local variables of coefficient wave |
---|
| 79 | int nord=76; //order of integration |
---|
| 80 | double va,vb; //upper and lower integration limits |
---|
| 81 | double summ,zi,yyy,answer,vell; //running tally of integration |
---|
[975ec8e] | 82 | double summj,vaj,vbj,zij,arg, si; //for the inner integration |
---|
[8e91f01] | 83 | |
---|
[ae3ce4e] | 84 | Pi = 4.0*atan(1.0); |
---|
| 85 | va = 0; |
---|
| 86 | vb = 1; //orintational average, outer integral |
---|
| 87 | vaj=0; |
---|
| 88 | vbj=Pi; //endpoints of inner integral |
---|
[8e91f01] | 89 | |
---|
[ae3ce4e] | 90 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 91 | |
---|
[ae3ce4e] | 92 | scale = dp[0]; //make local copies in case memory moves |
---|
| 93 | ra = dp[1]; |
---|
| 94 | nu = dp[2]; |
---|
| 95 | length = dp[3]; |
---|
| 96 | delrho = dp[4]; |
---|
| 97 | bkg = dp[5]; |
---|
| 98 | for(i=0;i<nord;i++) { |
---|
| 99 | //setup inner integral over the ellipsoidal cross-section |
---|
| 100 | summj=0; |
---|
| 101 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy |
---|
| 102 | arg = ra*sqrt(1-zi*zi); |
---|
| 103 | for(j=0;j<nord;j++) { |
---|
| 104 | //76 gauss points for the inner integral as well |
---|
| 105 | zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy |
---|
| 106 | yyy = Gauss76Wt[j] * EllipCylKernel(q,arg,nu,zij); |
---|
| 107 | summj += yyy; |
---|
| 108 | } |
---|
| 109 | //now calculate the value of the inner integral |
---|
| 110 | answer = (vbj-vaj)/2.0*summj; |
---|
| 111 | //divide integral by Pi |
---|
| 112 | answer /=Pi; |
---|
[8e91f01] | 113 | |
---|
[ae3ce4e] | 114 | //now calculate outer integral |
---|
| 115 | arg = q*length*zi/2; |
---|
[975ec8e] | 116 | if (arg == 0){ |
---|
| 117 | si = 1; |
---|
| 118 | }else{ |
---|
| 119 | si = sin(arg) * sin(arg) / arg / arg; |
---|
| 120 | } |
---|
| 121 | yyy = Gauss76Wt[i] * answer * si; |
---|
[ae3ce4e] | 122 | summ += yyy; |
---|
| 123 | } |
---|
| 124 | answer = (vb-va)/2.0*summ; |
---|
| 125 | // Multiply by contrast^2 |
---|
| 126 | answer *= delrho*delrho; |
---|
| 127 | //normalize by cylinder volume |
---|
| 128 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
| 129 | vell = Pi*ra*(nu*ra)*length; |
---|
| 130 | answer *= vell; |
---|
| 131 | //convert to [cm-1] |
---|
| 132 | answer *= 1.0e8; |
---|
| 133 | //Scale |
---|
| 134 | answer *= scale; |
---|
| 135 | // add in the background |
---|
| 136 | answer += bkg; |
---|
[8e91f01] | 137 | |
---|
[ae3ce4e] | 138 | return answer; |
---|
| 139 | } |
---|
| 140 | |
---|
| 141 | /* EllipCyl20X : calculates the form factor of a elliptical cylinder at the given x-value p->x |
---|
| 142 | |
---|
| 143 | Uses 76 pt Gaussian quadrature for orientational integral |
---|
| 144 | Uses 20 pt quadrature for the inner integral over the elliptical cross-section |
---|
| 145 | |
---|
| 146 | Warning: |
---|
| 147 | The call to WaveData() below returns a pointer to the middle |
---|
| 148 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 149 | calculations could cause memory to move, you should copy the coefficient |
---|
| 150 | values to local variables or an array before such operations. |
---|
| 151 | */ |
---|
| 152 | double |
---|
| 153 | EllipCyl20(double dp[], double q) |
---|
| 154 | { |
---|
| 155 | int i,j; |
---|
| 156 | double Pi; |
---|
| 157 | double scale,ra,nu,length,delrho,bkg; //local variables of coefficient wave |
---|
| 158 | int nordi=76; //order of integration |
---|
| 159 | int nordj=20; |
---|
| 160 | double va,vb; //upper and lower integration limits |
---|
| 161 | double summ,zi,yyy,answer,vell; //running tally of integration |
---|
[975ec8e] | 162 | double summj,vaj,vbj,zij,arg,si; //for the inner integration |
---|
[8e91f01] | 163 | |
---|
[ae3ce4e] | 164 | Pi = 4.0*atan(1.0); |
---|
| 165 | va = 0; |
---|
| 166 | vb = 1; //orintational average, outer integral |
---|
| 167 | vaj=0; |
---|
| 168 | vbj=Pi; //endpoints of inner integral |
---|
[8e91f01] | 169 | |
---|
[ae3ce4e] | 170 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 171 | |
---|
[ae3ce4e] | 172 | scale = dp[0]; //make local copies in case memory moves |
---|
| 173 | ra = dp[1]; |
---|
| 174 | nu = dp[2]; |
---|
| 175 | length = dp[3]; |
---|
| 176 | delrho = dp[4]; |
---|
| 177 | bkg = dp[5]; |
---|
| 178 | for(i=0;i<nordi;i++) { |
---|
| 179 | //setup inner integral over the ellipsoidal cross-section |
---|
| 180 | summj=0; |
---|
| 181 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy |
---|
| 182 | arg = ra*sqrt(1-zi*zi); |
---|
| 183 | for(j=0;j<nordj;j++) { |
---|
| 184 | //20 gauss points for the inner integral |
---|
| 185 | zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy |
---|
| 186 | yyy = Gauss20Wt[j] * EllipCylKernel(q,arg,nu,zij); |
---|
| 187 | summj += yyy; |
---|
| 188 | } |
---|
| 189 | //now calculate the value of the inner integral |
---|
| 190 | answer = (vbj-vaj)/2.0*summj; |
---|
| 191 | //divide integral by Pi |
---|
| 192 | answer /=Pi; |
---|
[8e91f01] | 193 | |
---|
[ae3ce4e] | 194 | //now calculate outer integral |
---|
| 195 | arg = q*length*zi/2; |
---|
[975ec8e] | 196 | if (arg == 0){ |
---|
| 197 | si = 1; |
---|
| 198 | }else{ |
---|
| 199 | si = sin(arg) * sin(arg) / arg / arg; |
---|
| 200 | } |
---|
| 201 | yyy = Gauss76Wt[i] * answer * si; |
---|
[ae3ce4e] | 202 | summ += yyy; |
---|
| 203 | } |
---|
[8e91f01] | 204 | |
---|
[ae3ce4e] | 205 | answer = (vb-va)/2.0*summ; |
---|
| 206 | // Multiply by contrast^2 |
---|
| 207 | answer *= delrho*delrho; |
---|
| 208 | //normalize by cylinder volume |
---|
| 209 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
| 210 | vell = Pi*ra*(nu*ra)*length; |
---|
| 211 | answer *= vell; |
---|
| 212 | //convert to [cm-1] |
---|
| 213 | answer *= 1.0e8; |
---|
| 214 | //Scale |
---|
| 215 | answer *= scale; |
---|
| 216 | // add in the background |
---|
[8e91f01] | 217 | answer += bkg; |
---|
| 218 | |
---|
[ae3ce4e] | 219 | return answer; |
---|
| 220 | } |
---|
| 221 | |
---|
| 222 | /* TriaxialEllipsoidX : calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x |
---|
| 223 | |
---|
| 224 | Uses 76 pt Gaussian quadrature for both integrals |
---|
| 225 | |
---|
| 226 | Warning: |
---|
| 227 | The call to WaveData() below returns a pointer to the middle |
---|
| 228 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 229 | calculations could cause memory to move, you should copy the coefficient |
---|
| 230 | values to local variables or an array before such operations. |
---|
| 231 | */ |
---|
| 232 | double |
---|
| 233 | TriaxialEllipsoid(double dp[], double q) |
---|
| 234 | { |
---|
| 235 | int i,j; |
---|
| 236 | double Pi; |
---|
| 237 | double scale,aa,bb,cc,delrho,bkg; //local variables of coefficient wave |
---|
| 238 | int nordi=76; //order of integration |
---|
| 239 | int nordj=76; |
---|
| 240 | double va,vb; //upper and lower integration limits |
---|
| 241 | double summ,zi,yyy,answer; //running tally of integration |
---|
| 242 | double summj,vaj,vbj,zij; //for the inner integration |
---|
[8e91f01] | 243 | |
---|
[ae3ce4e] | 244 | Pi = 4.0*atan(1.0); |
---|
| 245 | va = 0; |
---|
| 246 | vb = 1; //orintational average, outer integral |
---|
| 247 | vaj = 0; |
---|
| 248 | vbj = 1; //endpoints of inner integral |
---|
[8e91f01] | 249 | |
---|
[ae3ce4e] | 250 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 251 | |
---|
[ae3ce4e] | 252 | scale = dp[0]; //make local copies in case memory moves |
---|
| 253 | aa = dp[1]; |
---|
| 254 | bb = dp[2]; |
---|
| 255 | cc = dp[3]; |
---|
| 256 | delrho = dp[4]; |
---|
| 257 | bkg = dp[5]; |
---|
| 258 | for(i=0;i<nordi;i++) { |
---|
| 259 | //setup inner integral over the ellipsoidal cross-section |
---|
| 260 | summj=0; |
---|
| 261 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy |
---|
| 262 | for(j=0;j<nordj;j++) { |
---|
| 263 | //20 gauss points for the inner integral |
---|
| 264 | zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy |
---|
| 265 | yyy = Gauss76Wt[j] * TriaxialKernel(q,aa,bb,cc,zi,zij); |
---|
| 266 | summj += yyy; |
---|
| 267 | } |
---|
| 268 | //now calculate the value of the inner integral |
---|
| 269 | answer = (vbj-vaj)/2.0*summj; |
---|
[8e91f01] | 270 | |
---|
[ae3ce4e] | 271 | //now calculate outer integral |
---|
| 272 | yyy = Gauss76Wt[i] * answer; |
---|
| 273 | summ += yyy; |
---|
| 274 | } //final scaling is done at the end of the function, after the NT_FP64 case |
---|
[8e91f01] | 275 | |
---|
[ae3ce4e] | 276 | answer = (vb-va)/2.0*summ; |
---|
| 277 | // Multiply by contrast^2 |
---|
| 278 | answer *= delrho*delrho; |
---|
| 279 | //normalize by ellipsoid volume |
---|
| 280 | answer *= 4*Pi/3*aa*bb*cc; |
---|
| 281 | //convert to [cm-1] |
---|
| 282 | answer *= 1.0e8; |
---|
| 283 | //Scale |
---|
| 284 | answer *= scale; |
---|
| 285 | // add in the background |
---|
| 286 | answer += bkg; |
---|
[8e91f01] | 287 | |
---|
[ae3ce4e] | 288 | return answer; |
---|
| 289 | } |
---|
| 290 | |
---|
| 291 | /* ParallelepipedX : calculates the form factor of a Parallelepiped (a rectangular solid) |
---|
| 292 | at the given x-value p->x |
---|
| 293 | |
---|
| 294 | Uses 76 pt Gaussian quadrature for both integrals |
---|
| 295 | |
---|
| 296 | Warning: |
---|
| 297 | The call to WaveData() below returns a pointer to the middle |
---|
| 298 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 299 | calculations could cause memory to move, you should copy the coefficient |
---|
| 300 | values to local variables or an array before such operations. |
---|
| 301 | */ |
---|
| 302 | double |
---|
| 303 | Parallelepiped(double dp[], double q) |
---|
| 304 | { |
---|
| 305 | int i,j; |
---|
| 306 | double scale,aa,bb,cc,delrho,bkg; //local variables of coefficient wave |
---|
| 307 | int nordi=76; //order of integration |
---|
| 308 | int nordj=76; |
---|
| 309 | double va,vb; //upper and lower integration limits |
---|
| 310 | double summ,yyy,answer; //running tally of integration |
---|
| 311 | double summj,vaj,vbj; //for the inner integration |
---|
| 312 | double mu,mudum,arg,sigma,uu,vol; |
---|
[8e91f01] | 313 | |
---|
| 314 | |
---|
[ae3ce4e] | 315 | // Pi = 4.0*atan(1.0); |
---|
| 316 | va = 0; |
---|
| 317 | vb = 1; //orintational average, outer integral |
---|
| 318 | vaj = 0; |
---|
| 319 | vbj = 1; //endpoints of inner integral |
---|
[8e91f01] | 320 | |
---|
[ae3ce4e] | 321 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 322 | |
---|
[ae3ce4e] | 323 | scale = dp[0]; //make local copies in case memory moves |
---|
| 324 | aa = dp[1]; |
---|
| 325 | bb = dp[2]; |
---|
| 326 | cc = dp[3]; |
---|
| 327 | delrho = dp[4]; |
---|
| 328 | bkg = dp[5]; |
---|
[8e91f01] | 329 | |
---|
[ae3ce4e] | 330 | mu = q*bb; |
---|
| 331 | vol = aa*bb*cc; |
---|
| 332 | // normalize all WRT bb |
---|
| 333 | aa = aa/bb; |
---|
| 334 | cc = cc/bb; |
---|
[8e91f01] | 335 | |
---|
[ae3ce4e] | 336 | for(i=0;i<nordi;i++) { |
---|
| 337 | //setup inner integral over the ellipsoidal cross-section |
---|
| 338 | summj=0; |
---|
| 339 | sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the outer dummy |
---|
[8e91f01] | 340 | |
---|
[ae3ce4e] | 341 | for(j=0;j<nordj;j++) { |
---|
| 342 | //76 gauss points for the inner integral |
---|
| 343 | uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the inner dummy |
---|
| 344 | mudum = mu*sqrt(1-sigma*sigma); |
---|
| 345 | yyy = Gauss76Wt[j] * PPKernel(aa,mudum,uu); |
---|
| 346 | summj += yyy; |
---|
| 347 | } |
---|
| 348 | //now calculate the value of the inner integral |
---|
| 349 | answer = (vbj-vaj)/2.0*summj; |
---|
[8e91f01] | 350 | |
---|
[ae3ce4e] | 351 | arg = mu*cc*sigma/2; |
---|
| 352 | if ( arg == 0 ) { |
---|
| 353 | answer *= 1; |
---|
| 354 | } else { |
---|
| 355 | answer *= sin(arg)*sin(arg)/arg/arg; |
---|
| 356 | } |
---|
[8e91f01] | 357 | |
---|
[ae3ce4e] | 358 | //now sum up the outer integral |
---|
| 359 | yyy = Gauss76Wt[i] * answer; |
---|
| 360 | summ += yyy; |
---|
| 361 | } //final scaling is done at the end of the function, after the NT_FP64 case |
---|
[8e91f01] | 362 | |
---|
[ae3ce4e] | 363 | answer = (vb-va)/2.0*summ; |
---|
| 364 | // Multiply by contrast^2 |
---|
| 365 | answer *= delrho*delrho; |
---|
| 366 | //normalize by volume |
---|
| 367 | answer *= vol; |
---|
| 368 | //convert to [cm-1] |
---|
| 369 | answer *= 1.0e8; |
---|
| 370 | //Scale |
---|
| 371 | answer *= scale; |
---|
| 372 | // add in the background |
---|
| 373 | answer += bkg; |
---|
[8e91f01] | 374 | |
---|
[ae3ce4e] | 375 | return answer; |
---|
| 376 | } |
---|
| 377 | |
---|
| 378 | /* HollowCylinderX : calculates the form factor of a Hollow Cylinder |
---|
| 379 | at the given x-value p->x |
---|
| 380 | |
---|
| 381 | Uses 76 pt Gaussian quadrature for the single integral |
---|
| 382 | |
---|
| 383 | Warning: |
---|
| 384 | The call to WaveData() below returns a pointer to the middle |
---|
| 385 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 386 | calculations could cause memory to move, you should copy the coefficient |
---|
| 387 | values to local variables or an array before such operations. |
---|
| 388 | */ |
---|
| 389 | double |
---|
| 390 | HollowCylinder(double dp[], double q) |
---|
| 391 | { |
---|
| 392 | int i; |
---|
| 393 | double scale,rcore,rshell,length,delrho,bkg; //local variables of coefficient wave |
---|
| 394 | int nord=76; //order of integration |
---|
| 395 | double va,vb,zi; //upper and lower integration limits |
---|
| 396 | double summ,answer,pi; //running tally of integration |
---|
[8e91f01] | 397 | |
---|
[ae3ce4e] | 398 | pi = 4.0*atan(1.0); |
---|
| 399 | va = 0; |
---|
| 400 | vb = 1; //limits of numerical integral |
---|
[8e91f01] | 401 | |
---|
[ae3ce4e] | 402 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 403 | |
---|
[ae3ce4e] | 404 | scale = dp[0]; //make local copies in case memory moves |
---|
| 405 | rcore = dp[1]; |
---|
| 406 | rshell = dp[2]; |
---|
| 407 | length = dp[3]; |
---|
| 408 | delrho = dp[4]; |
---|
| 409 | bkg = dp[5]; |
---|
[8e91f01] | 410 | |
---|
[ae3ce4e] | 411 | for(i=0;i<nord;i++) { |
---|
| 412 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; |
---|
| 413 | summ += Gauss76Wt[i] * HolCylKernel(q, rcore, rshell, length, zi); |
---|
| 414 | } |
---|
[8e91f01] | 415 | |
---|
[ae3ce4e] | 416 | answer = (vb-va)/2.0*summ; |
---|
| 417 | // Multiply by contrast^2 |
---|
| 418 | answer *= delrho*delrho; |
---|
| 419 | //normalize by volume |
---|
| 420 | answer *= pi*(rshell*rshell-rcore*rcore)*length; |
---|
| 421 | //convert to [cm-1] |
---|
| 422 | answer *= 1.0e8; |
---|
| 423 | //Scale |
---|
| 424 | answer *= scale; |
---|
| 425 | // add in the background |
---|
| 426 | answer += bkg; |
---|
[8e91f01] | 427 | |
---|
[ae3ce4e] | 428 | return answer; |
---|
| 429 | } |
---|
| 430 | |
---|
| 431 | /* EllipsoidFormX : calculates the form factor of an ellipsoid of revolution with semiaxes a:a:nua |
---|
| 432 | at the given x-value p->x |
---|
| 433 | |
---|
| 434 | Uses 76 pt Gaussian quadrature for the single integral |
---|
| 435 | |
---|
| 436 | Warning: |
---|
| 437 | The call to WaveData() below returns a pointer to the middle |
---|
| 438 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
| 439 | calculations could cause memory to move, you should copy the coefficient |
---|
| 440 | values to local variables or an array before such operations. |
---|
| 441 | */ |
---|
| 442 | double |
---|
| 443 | EllipsoidForm(double dp[], double q) |
---|
| 444 | { |
---|
| 445 | int i; |
---|
| 446 | double scale,a,nua,delrho,bkg; //local variables of coefficient wave |
---|
| 447 | int nord=76; //order of integration |
---|
| 448 | double va,vb,zi; //upper and lower integration limits |
---|
| 449 | double summ,answer,pi; //running tally of integration |
---|
[8e91f01] | 450 | |
---|
[ae3ce4e] | 451 | pi = 4.0*atan(1.0); |
---|
| 452 | va = 0; |
---|
| 453 | vb = 1; //limits of numerical integral |
---|
[8e91f01] | 454 | |
---|
[ae3ce4e] | 455 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 456 | |
---|
[ae3ce4e] | 457 | scale = dp[0]; //make local copies in case memory moves |
---|
| 458 | nua = dp[1]; |
---|
| 459 | a = dp[2]; |
---|
| 460 | delrho = dp[3]; |
---|
| 461 | bkg = dp[4]; |
---|
[8e91f01] | 462 | |
---|
[ae3ce4e] | 463 | for(i=0;i<nord;i++) { |
---|
| 464 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; |
---|
| 465 | summ += Gauss76Wt[i] * EllipsoidKernel(q, a, nua, zi); |
---|
| 466 | } |
---|
[8e91f01] | 467 | |
---|
[ae3ce4e] | 468 | answer = (vb-va)/2.0*summ; |
---|
| 469 | // Multiply by contrast^2 |
---|
| 470 | answer *= delrho*delrho; |
---|
| 471 | //normalize by volume |
---|
| 472 | answer *= 4*pi/3*a*a*nua; |
---|
| 473 | //convert to [cm-1] |
---|
| 474 | answer *= 1.0e8; |
---|
| 475 | //Scale |
---|
| 476 | answer *= scale; |
---|
| 477 | // add in the background |
---|
| 478 | answer += bkg; |
---|
[8e91f01] | 479 | |
---|
[ae3ce4e] | 480 | return answer; |
---|
| 481 | } |
---|
| 482 | |
---|
| 483 | |
---|
| 484 | /* Cyl_PolyRadiusX : calculates the form factor of a cylinder at the given x-value p->x |
---|
| 485 | the cylinder has a polydisperse cross section |
---|
| 486 | |
---|
| 487 | */ |
---|
| 488 | double |
---|
| 489 | Cyl_PolyRadius(double dp[], double q) |
---|
| 490 | { |
---|
| 491 | int i; |
---|
| 492 | double scale,radius,length,pd,delrho,bkg; //local variables of coefficient wave |
---|
| 493 | int nord=20; //order of integration |
---|
| 494 | double uplim,lolim; //upper and lower integration limits |
---|
| 495 | double summ,zi,yyy,answer,Vpoly; //running tally of integration |
---|
| 496 | double range,zz,Pi; |
---|
[8e91f01] | 497 | |
---|
[ae3ce4e] | 498 | Pi = 4.0*atan(1.0); |
---|
| 499 | range = 3.4; |
---|
[8e91f01] | 500 | |
---|
[ae3ce4e] | 501 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 502 | |
---|
[ae3ce4e] | 503 | scale = dp[0]; //make local copies in case memory moves |
---|
| 504 | radius = dp[1]; |
---|
| 505 | length = dp[2]; |
---|
| 506 | pd = dp[3]; |
---|
| 507 | delrho = dp[4]; |
---|
| 508 | bkg = dp[5]; |
---|
[8e91f01] | 509 | |
---|
[ae3ce4e] | 510 | zz = (1.0/pd)*(1.0/pd) - 1.0; |
---|
[8e91f01] | 511 | |
---|
[ae3ce4e] | 512 | lolim = radius*(1.0-range*pd); //set the upper/lower limits to cover the distribution |
---|
| 513 | if(lolim<0) { |
---|
| 514 | lolim = 0; |
---|
| 515 | } |
---|
| 516 | if(pd>0.3) { |
---|
| 517 | range = 3.4 + (pd-0.3)*18.0; |
---|
| 518 | } |
---|
| 519 | uplim = radius*(1.0+range*pd); |
---|
[8e91f01] | 520 | |
---|
[ae3ce4e] | 521 | for(i=0;i<nord;i++) { |
---|
| 522 | zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 523 | yyy = Gauss20Wt[i] * Cyl_PolyRadKernel(q, radius, length, zz, delrho, zi); |
---|
| 524 | summ += yyy; |
---|
| 525 | } |
---|
[8e91f01] | 526 | |
---|
[ae3ce4e] | 527 | answer = (uplim-lolim)/2.0*summ; |
---|
| 528 | //normalize by average cylinder volume |
---|
| 529 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
| 530 | Vpoly=Pi*radius*radius*length*(zz+2.0)/(zz+1.0); |
---|
| 531 | answer /= Vpoly; |
---|
| 532 | //convert to [cm-1] |
---|
| 533 | answer *= 1.0e8; |
---|
| 534 | //Scale |
---|
| 535 | answer *= scale; |
---|
| 536 | // add in the background |
---|
| 537 | answer += bkg; |
---|
[8e91f01] | 538 | |
---|
[ae3ce4e] | 539 | return answer; |
---|
| 540 | } |
---|
| 541 | |
---|
| 542 | /* Cyl_PolyLengthX : calculates the form factor of a cylinder at the given x-value p->x |
---|
| 543 | the cylinder has a polydisperse Length |
---|
| 544 | |
---|
| 545 | */ |
---|
| 546 | double |
---|
| 547 | Cyl_PolyLength(double dp[], double q) |
---|
| 548 | { |
---|
| 549 | int i; |
---|
| 550 | double scale,radius,length,pd,delrho,bkg; //local variables of coefficient wave |
---|
| 551 | int nord=20; //order of integration |
---|
| 552 | double uplim,lolim; //upper and lower integration limits |
---|
| 553 | double summ,zi,yyy,answer,Vpoly; //running tally of integration |
---|
| 554 | double range,zz,Pi; |
---|
[8e91f01] | 555 | |
---|
| 556 | |
---|
[ae3ce4e] | 557 | Pi = 4.0*atan(1.0); |
---|
| 558 | range = 3.4; |
---|
[8e91f01] | 559 | |
---|
[ae3ce4e] | 560 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 561 | |
---|
[ae3ce4e] | 562 | scale = dp[0]; //make local copies in case memory moves |
---|
| 563 | radius = dp[1]; |
---|
| 564 | length = dp[2]; |
---|
| 565 | pd = dp[3]; |
---|
| 566 | delrho = dp[4]; |
---|
| 567 | bkg = dp[5]; |
---|
[8e91f01] | 568 | |
---|
[ae3ce4e] | 569 | zz = (1.0/pd)*(1.0/pd) - 1.0; |
---|
[8e91f01] | 570 | |
---|
[ae3ce4e] | 571 | lolim = length*(1.0-range*pd); //set the upper/lower limits to cover the distribution |
---|
| 572 | if(lolim<0) { |
---|
| 573 | lolim = 0; |
---|
| 574 | } |
---|
| 575 | if(pd>0.3) { |
---|
| 576 | range = 3.4 + (pd-0.3)*18.0; |
---|
| 577 | } |
---|
| 578 | uplim = length*(1.0+range*pd); |
---|
[8e91f01] | 579 | |
---|
[ae3ce4e] | 580 | for(i=0;i<nord;i++) { |
---|
| 581 | zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 582 | yyy = Gauss20Wt[i] * Cyl_PolyLenKernel(q, radius, length, zz, delrho, zi); |
---|
| 583 | summ += yyy; |
---|
| 584 | } |
---|
[8e91f01] | 585 | |
---|
[ae3ce4e] | 586 | answer = (uplim-lolim)/2.0*summ; |
---|
| 587 | //normalize by average cylinder volume (first moment) |
---|
| 588 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
| 589 | Vpoly=Pi*radius*radius*length; |
---|
| 590 | answer /= Vpoly; |
---|
| 591 | //convert to [cm-1] |
---|
| 592 | answer *= 1.0e8; |
---|
| 593 | //Scale |
---|
| 594 | answer *= scale; |
---|
| 595 | // add in the background |
---|
| 596 | answer += bkg; |
---|
[8e91f01] | 597 | |
---|
[ae3ce4e] | 598 | return answer; |
---|
| 599 | } |
---|
| 600 | |
---|
| 601 | /* CoreShellCylinderX : calculates the form factor of a cylinder at the given x-value p->x |
---|
| 602 | the cylinder has a core-shell structure |
---|
| 603 | |
---|
| 604 | */ |
---|
| 605 | double |
---|
| 606 | CoreShellCylinder(double dp[], double q) |
---|
| 607 | { |
---|
| 608 | int i; |
---|
| 609 | double scale,rcore,length,bkg; //local variables of coefficient wave |
---|
| 610 | double thick,rhoc,rhos,rhosolv; |
---|
| 611 | int nord=76; //order of integration |
---|
| 612 | double uplim,lolim,halfheight; //upper and lower integration limits |
---|
| 613 | double summ,zi,yyy,answer,Vcyl; //running tally of integration |
---|
| 614 | double Pi; |
---|
[8e91f01] | 615 | |
---|
[ae3ce4e] | 616 | Pi = 4.0*atan(1.0); |
---|
[8e91f01] | 617 | |
---|
[ae3ce4e] | 618 | lolim = 0.0; |
---|
| 619 | uplim = Pi/2.0; |
---|
[8e91f01] | 620 | |
---|
[ae3ce4e] | 621 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 622 | |
---|
[ae3ce4e] | 623 | scale = dp[0]; //make local copies in case memory moves |
---|
| 624 | rcore = dp[1]; |
---|
| 625 | thick = dp[2]; |
---|
| 626 | length = dp[3]; |
---|
| 627 | rhoc = dp[4]; |
---|
| 628 | rhos = dp[5]; |
---|
| 629 | rhosolv = dp[6]; |
---|
| 630 | bkg = dp[7]; |
---|
[8e91f01] | 631 | |
---|
[ae3ce4e] | 632 | halfheight = length/2.0; |
---|
[8e91f01] | 633 | |
---|
[ae3ce4e] | 634 | for(i=0;i<nord;i++) { |
---|
| 635 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 636 | yyy = Gauss76Wt[i] * CoreShellCylKernel(q, rcore, thick, rhoc,rhos,rhosolv, halfheight, zi); |
---|
| 637 | summ += yyy; |
---|
| 638 | } |
---|
[8e91f01] | 639 | |
---|
[ae3ce4e] | 640 | answer = (uplim-lolim)/2.0*summ; |
---|
[8e91f01] | 641 | // length is the total core length |
---|
[ae3ce4e] | 642 | Vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2.0*thick); |
---|
| 643 | answer /= Vcyl; |
---|
| 644 | //convert to [cm-1] |
---|
| 645 | answer *= 1.0e8; |
---|
| 646 | //Scale |
---|
| 647 | answer *= scale; |
---|
| 648 | // add in the background |
---|
| 649 | answer += bkg; |
---|
[8e91f01] | 650 | |
---|
[ae3ce4e] | 651 | return answer; |
---|
| 652 | } |
---|
| 653 | |
---|
| 654 | |
---|
| 655 | /* PolyCoShCylinderX : calculates the form factor of a core-shell cylinder at the given x-value p->x |
---|
| 656 | the cylinder has a polydisperse CORE radius |
---|
| 657 | |
---|
| 658 | */ |
---|
| 659 | double |
---|
| 660 | PolyCoShCylinder(double dp[], double q) |
---|
| 661 | { |
---|
| 662 | int i; |
---|
| 663 | double scale,radius,length,sigma,bkg; //local variables of coefficient wave |
---|
| 664 | double rad,radthick,facthick,rhoc,rhos,rhosolv; |
---|
| 665 | int nord=20; //order of integration |
---|
| 666 | double uplim,lolim; //upper and lower integration limits |
---|
| 667 | double summ,yyy,answer,Vpoly; //running tally of integration |
---|
| 668 | double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr; |
---|
[8e91f01] | 669 | |
---|
[ae3ce4e] | 670 | Pi = 4.0*atan(1.0); |
---|
[8e91f01] | 671 | |
---|
[ae3ce4e] | 672 | summ = 0.0; //initialize intergral |
---|
| 673 | Rsqrsumm = 0.0; |
---|
[8e91f01] | 674 | |
---|
[ae3ce4e] | 675 | scale = dp[0]; |
---|
| 676 | radius = dp[1]; |
---|
| 677 | sigma = dp[2]; //sigma is the standard mean deviation |
---|
| 678 | length = dp[3]; |
---|
| 679 | radthick = dp[4]; |
---|
| 680 | facthick= dp[5]; |
---|
| 681 | rhoc = dp[6]; |
---|
| 682 | rhos = dp[7]; |
---|
| 683 | rhosolv = dp[8]; |
---|
| 684 | bkg = dp[9]; |
---|
[8e91f01] | 685 | |
---|
[ae3ce4e] | 686 | lolim = exp(log(radius)-(4.*sigma)); |
---|
| 687 | if (lolim<0) { |
---|
| 688 | lolim=0; //to avoid numerical error when va<0 (-ve r value) |
---|
| 689 | } |
---|
| 690 | uplim = exp(log(radius)+(4.*sigma)); |
---|
[8e91f01] | 691 | |
---|
[ae3ce4e] | 692 | for(i=0;i<nord;i++) { |
---|
| 693 | rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 694 | AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp(-(0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma))); |
---|
| 695 | yyy = AR* Gauss20Wt[i] * CSCylIntegration(q,rad,radthick,facthick,rhoc,rhos,rhosolv,length); |
---|
| 696 | Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick); //SRK normalize to total dimensions |
---|
| 697 | summ += yyy; |
---|
| 698 | Rsqrsumm += Rsqryyy; |
---|
| 699 | } |
---|
[8e91f01] | 700 | |
---|
[ae3ce4e] | 701 | answer = (uplim-lolim)/2.0*summ; |
---|
| 702 | Rsqr = (uplim-lolim)/2.0*Rsqrsumm; |
---|
| 703 | //normalize by average cylinder volume |
---|
| 704 | Vpoly = Pi*Rsqr*(length+2*facthick); |
---|
| 705 | answer /= Vpoly; |
---|
| 706 | //convert to [cm-1] |
---|
| 707 | answer *= 1.0e8; |
---|
| 708 | //Scale |
---|
| 709 | answer *= scale; |
---|
| 710 | // add in the background |
---|
| 711 | answer += bkg; |
---|
[8e91f01] | 712 | |
---|
[ae3ce4e] | 713 | return answer; |
---|
| 714 | } |
---|
| 715 | |
---|
| 716 | /* OblateFormX : calculates the form factor of a core-shell Oblate ellipsoid at the given x-value p->x |
---|
| 717 | the ellipsoid has a core-shell structure |
---|
| 718 | |
---|
| 719 | */ |
---|
| 720 | double |
---|
| 721 | OblateForm(double dp[], double q) |
---|
| 722 | { |
---|
| 723 | int i; |
---|
| 724 | double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg; |
---|
| 725 | int nord=76; //order of integration |
---|
| 726 | double uplim,lolim; //upper and lower integration limits |
---|
| 727 | double summ,zi,yyy,answer,oblatevol; //running tally of integration |
---|
| 728 | double Pi; |
---|
[8e91f01] | 729 | |
---|
[ae3ce4e] | 730 | Pi = 4.0*atan(1.0); |
---|
[8e91f01] | 731 | |
---|
[ae3ce4e] | 732 | lolim = 0.0; |
---|
| 733 | uplim = 1.0; |
---|
[8e91f01] | 734 | |
---|
[ae3ce4e] | 735 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 736 | |
---|
| 737 | |
---|
[ae3ce4e] | 738 | scale = dp[0]; //make local copies in case memory moves |
---|
| 739 | crmaj = dp[1]; |
---|
| 740 | crmin = dp[2]; |
---|
| 741 | trmaj = dp[3]; |
---|
| 742 | trmin = dp[4]; |
---|
| 743 | delpc = dp[5]; |
---|
| 744 | delps = dp[6]; |
---|
[8e91f01] | 745 | bkg = dp[7]; |
---|
| 746 | |
---|
[ae3ce4e] | 747 | for(i=0;i<nord;i++) { |
---|
| 748 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 749 | yyy = Gauss76Wt[i] * gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q); |
---|
| 750 | summ += yyy; |
---|
| 751 | } |
---|
[8e91f01] | 752 | |
---|
[ae3ce4e] | 753 | answer = (uplim-lolim)/2.0*summ; |
---|
| 754 | // normalize by particle volume |
---|
| 755 | oblatevol = 4*Pi/3*trmaj*trmaj*trmin; |
---|
| 756 | answer /= oblatevol; |
---|
[8e91f01] | 757 | |
---|
[ae3ce4e] | 758 | //convert to [cm-1] |
---|
| 759 | answer *= 1.0e8; |
---|
| 760 | //Scale |
---|
| 761 | answer *= scale; |
---|
| 762 | // add in the background |
---|
| 763 | answer += bkg; |
---|
[8e91f01] | 764 | |
---|
[ae3ce4e] | 765 | return answer; |
---|
| 766 | } |
---|
| 767 | |
---|
| 768 | /* ProlateFormX : calculates the form factor of a core-shell Prolate ellipsoid at the given x-value p->x |
---|
| 769 | the ellipsoid has a core-shell structure |
---|
| 770 | |
---|
| 771 | */ |
---|
| 772 | double |
---|
| 773 | ProlateForm(double dp[], double q) |
---|
| 774 | { |
---|
| 775 | int i; |
---|
| 776 | double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg; |
---|
| 777 | int nord=76; //order of integration |
---|
| 778 | double uplim,lolim; //upper and lower integration limits |
---|
| 779 | double summ,zi,yyy,answer,prolatevol; //running tally of integration |
---|
| 780 | double Pi; |
---|
[8e91f01] | 781 | |
---|
[ae3ce4e] | 782 | Pi = 4.0*atan(1.0); |
---|
[8e91f01] | 783 | |
---|
[ae3ce4e] | 784 | lolim = 0.0; |
---|
| 785 | uplim = 1.0; |
---|
[8e91f01] | 786 | |
---|
[ae3ce4e] | 787 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 788 | |
---|
[ae3ce4e] | 789 | scale = dp[0]; //make local copies in case memory moves |
---|
| 790 | crmaj = dp[1]; |
---|
| 791 | crmin = dp[2]; |
---|
| 792 | trmaj = dp[3]; |
---|
| 793 | trmin = dp[4]; |
---|
| 794 | delpc = dp[5]; |
---|
| 795 | delps = dp[6]; |
---|
[8e91f01] | 796 | bkg = dp[7]; |
---|
| 797 | |
---|
[ae3ce4e] | 798 | for(i=0;i<nord;i++) { |
---|
| 799 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 800 | yyy = Gauss76Wt[i] * gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q); |
---|
| 801 | summ += yyy; |
---|
| 802 | } |
---|
[8e91f01] | 803 | |
---|
[ae3ce4e] | 804 | answer = (uplim-lolim)/2.0*summ; |
---|
| 805 | // normalize by particle volume |
---|
| 806 | prolatevol = 4*Pi/3*trmaj*trmin*trmin; |
---|
| 807 | answer /= prolatevol; |
---|
[8e91f01] | 808 | |
---|
[ae3ce4e] | 809 | //convert to [cm-1] |
---|
| 810 | answer *= 1.0e8; |
---|
| 811 | //Scale |
---|
| 812 | answer *= scale; |
---|
| 813 | // add in the background |
---|
| 814 | answer += bkg; |
---|
[8e91f01] | 815 | |
---|
[ae3ce4e] | 816 | return answer; |
---|
| 817 | } |
---|
| 818 | |
---|
| 819 | |
---|
| 820 | /* StackedDiscsX : calculates the form factor of a stacked "tactoid" of core shell disks |
---|
| 821 | like clay platelets that are not exfoliated |
---|
| 822 | |
---|
| 823 | */ |
---|
| 824 | double |
---|
| 825 | StackedDiscs(double dp[], double q) |
---|
| 826 | { |
---|
| 827 | int i; |
---|
| 828 | double scale,length,bkg,rcore,thick,rhoc,rhol,rhosolv,N,gsd; //local variables of coefficient wave |
---|
| 829 | double va,vb,vcyl,summ,yyy,zi,halfheight,d,answer; |
---|
| 830 | int nord=76; //order of integration |
---|
| 831 | double Pi; |
---|
[8e91f01] | 832 | |
---|
| 833 | |
---|
[ae3ce4e] | 834 | Pi = 4.0*atan(1.0); |
---|
[8e91f01] | 835 | |
---|
[ae3ce4e] | 836 | va = 0.0; |
---|
| 837 | vb = Pi/2.0; |
---|
[8e91f01] | 838 | |
---|
[ae3ce4e] | 839 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 840 | |
---|
[ae3ce4e] | 841 | scale = dp[0]; |
---|
| 842 | rcore = dp[1]; |
---|
| 843 | length = dp[2]; |
---|
| 844 | thick = dp[3]; |
---|
| 845 | rhoc = dp[4]; |
---|
| 846 | rhol = dp[5]; |
---|
| 847 | rhosolv = dp[6]; |
---|
| 848 | N = dp[7]; |
---|
| 849 | gsd = dp[8]; |
---|
| 850 | bkg = dp[9]; |
---|
[8e91f01] | 851 | |
---|
[ae3ce4e] | 852 | d=2.0*thick+length; |
---|
| 853 | halfheight = length/2.0; |
---|
[8e91f01] | 854 | |
---|
[ae3ce4e] | 855 | for(i=0;i<nord;i++) { |
---|
| 856 | zi = ( Gauss76Z[i]*(vb-va) + vb + va )/2.0; |
---|
| 857 | yyy = Gauss76Wt[i] * Stackdisc_kern(q, rcore, rhoc,rhol,rhosolv, halfheight,thick,zi,gsd,d,N); |
---|
| 858 | summ += yyy; |
---|
| 859 | } |
---|
[8e91f01] | 860 | |
---|
[ae3ce4e] | 861 | answer = (vb-va)/2.0*summ; |
---|
[8e91f01] | 862 | // length is the total core length |
---|
[ae3ce4e] | 863 | vcyl=Pi*rcore*rcore*(2.0*thick+length)*N; |
---|
| 864 | answer /= vcyl; |
---|
| 865 | //Convert to [cm-1] |
---|
| 866 | answer *= 1.0e8; |
---|
| 867 | //Scale |
---|
| 868 | answer *= scale; |
---|
| 869 | // add in the background |
---|
| 870 | answer += bkg; |
---|
[8e91f01] | 871 | |
---|
[ae3ce4e] | 872 | return answer; |
---|
| 873 | } |
---|
| 874 | |
---|
| 875 | |
---|
| 876 | /* LamellarFFX : calculates the form factor of a lamellar structure - no S(q) effects included |
---|
| 877 | |
---|
| 878 | */ |
---|
| 879 | double |
---|
| 880 | LamellarFF(double dp[], double q) |
---|
| 881 | { |
---|
| 882 | double scale,del,sig,contr,bkg; //local variables of coefficient wave |
---|
| 883 | double inten, qval,Pq; |
---|
| 884 | double Pi; |
---|
[8e91f01] | 885 | |
---|
| 886 | |
---|
[ae3ce4e] | 887 | Pi = 4.0*atan(1.0); |
---|
| 888 | scale = dp[0]; |
---|
| 889 | del = dp[1]; |
---|
| 890 | sig = dp[2]*del; |
---|
| 891 | contr = dp[3]; |
---|
| 892 | bkg = dp[4]; |
---|
[34c3020] | 893 | qval = q; |
---|
[8e91f01] | 894 | |
---|
[ae3ce4e] | 895 | Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig)); |
---|
[8e91f01] | 896 | |
---|
[ae3ce4e] | 897 | inten = 2.0*Pi*scale*Pq/(qval*qval); //this is now dimensionless... |
---|
[8e91f01] | 898 | |
---|
[ae3ce4e] | 899 | inten /= del; //normalize by the thickness (in A) |
---|
[8e91f01] | 900 | |
---|
[ae3ce4e] | 901 | inten *= 1.0e8; // 1/A to 1/cm |
---|
[8e91f01] | 902 | |
---|
[ae3ce4e] | 903 | return(inten+bkg); |
---|
| 904 | } |
---|
[975ec8e] | 905 | /* LamellarFFX : calculates the form factor of a lamellar structure - no S(q) effects included |
---|
| 906 | -NO polydispersion included |
---|
| 907 | */ |
---|
| 908 | double lamellar_kernel(double dp[], double q){ |
---|
| 909 | double scale,del,sld_bi,sld_sol,contr,bkg; //local variables of coefficient wave |
---|
| 910 | double inten, qval,Pq; |
---|
| 911 | double Pi; |
---|
| 912 | |
---|
| 913 | |
---|
| 914 | Pi = 4.0*atan(1.0); |
---|
| 915 | scale = dp[0]; |
---|
| 916 | del = dp[1]; |
---|
| 917 | sld_bi = dp[2]; |
---|
| 918 | sld_sol = dp[3]; |
---|
| 919 | bkg = dp[4]; |
---|
| 920 | qval = q; |
---|
| 921 | contr = sld_bi -sld_sol; |
---|
[ae3ce4e] | 922 | |
---|
[975ec8e] | 923 | Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)); |
---|
| 924 | |
---|
| 925 | inten = 2.0*Pi*scale*Pq/(qval*qval); //this is now dimensionless... |
---|
| 926 | |
---|
| 927 | inten /= del; //normalize by the thickness (in A) |
---|
| 928 | |
---|
| 929 | inten *= 1.0e8; // 1/A to 1/cm |
---|
| 930 | |
---|
| 931 | return(inten+bkg); |
---|
| 932 | } |
---|
[ae3ce4e] | 933 | /* LamellarPSX : calculates the form factor of a lamellar structure - with S(q) effects included |
---|
| 934 | ------- |
---|
| 935 | ------- resolution effects ARE included, but only a CONSTANT default value, not the real q-dependent resolution!! |
---|
| 936 | |
---|
| 937 | */ |
---|
| 938 | double |
---|
| 939 | LamellarPS(double dp[], double q) |
---|
| 940 | { |
---|
| 941 | double scale,dd,del,sig,contr,NN,Cp,bkg; //local variables of coefficient wave |
---|
| 942 | double inten, qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ; |
---|
| 943 | double Pi,Euler,dQDefault,fii; |
---|
| 944 | int ii,NNint; |
---|
[8e91f01] | 945 | |
---|
[ae3ce4e] | 946 | Euler = 0.5772156649; // Euler's constant |
---|
[8e91f01] | 947 | dQDefault = 0.0;//0.0025; //[=] 1/A, q-resolution, default value |
---|
[ae3ce4e] | 948 | dQ = dQDefault; |
---|
[8e91f01] | 949 | |
---|
[ae3ce4e] | 950 | Pi = 4.0*atan(1.0); |
---|
| 951 | qval = q; |
---|
[8e91f01] | 952 | |
---|
[ae3ce4e] | 953 | scale = dp[0]; |
---|
| 954 | dd = dp[1]; |
---|
| 955 | del = dp[2]; |
---|
| 956 | sig = dp[3]*del; |
---|
| 957 | contr = dp[4]; |
---|
| 958 | NN = trunc(dp[5]); //be sure that NN is an integer |
---|
| 959 | Cp = dp[6]; |
---|
| 960 | bkg = dp[7]; |
---|
[8e91f01] | 961 | |
---|
[ae3ce4e] | 962 | Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig)); |
---|
[8e91f01] | 963 | |
---|
[ae3ce4e] | 964 | NNint = (int)NN; //cast to an integer for the loop |
---|
| 965 | ii=0; |
---|
| 966 | Sq = 0.0; |
---|
| 967 | for(ii=1;ii<(NNint-1);ii+=1) { |
---|
[8e91f01] | 968 | |
---|
[ae3ce4e] | 969 | fii = (double)ii; //do I really need to do this? |
---|
[8e91f01] | 970 | |
---|
[ae3ce4e] | 971 | temp = 0.0; |
---|
| 972 | alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler); |
---|
| 973 | t1 = 2.0*dQ*dQ*dd*dd*alpha; |
---|
| 974 | t2 = 2.0*qval*qval*dd*dd*alpha; |
---|
| 975 | t3 = dQ*dQ*dd*dd*ii*ii; |
---|
[8e91f01] | 976 | |
---|
[ae3ce4e] | 977 | temp = 1.0-ii/NN; |
---|
| 978 | temp *= cos(dd*qval*ii/(1.0+t1)); |
---|
| 979 | temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); |
---|
| 980 | temp /= sqrt(1.0+t1); |
---|
[8e91f01] | 981 | |
---|
[ae3ce4e] | 982 | Sq += temp; |
---|
| 983 | } |
---|
[8e91f01] | 984 | |
---|
[ae3ce4e] | 985 | Sq *= 2.0; |
---|
| 986 | Sq += 1.0; |
---|
[8e91f01] | 987 | |
---|
[ae3ce4e] | 988 | inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval); |
---|
[8e91f01] | 989 | |
---|
[ae3ce4e] | 990 | inten *= 1.0e8; // 1/A to 1/cm |
---|
[8e91f01] | 991 | |
---|
[ae3ce4e] | 992 | return(inten+bkg); |
---|
| 993 | } |
---|
| 994 | |
---|
| 995 | |
---|
| 996 | /* LamellarPS_HGX : calculates the form factor of a lamellar structure - with S(q) effects included |
---|
| 997 | ------- |
---|
| 998 | ------- resolution effects ARE included, but only a CONSTANT default value, not the real q-dependent resolution!! |
---|
| 999 | |
---|
| 1000 | */ |
---|
| 1001 | double |
---|
| 1002 | LamellarPS_HG(double dp[], double q) |
---|
| 1003 | { |
---|
| 1004 | double scale,dd,delT,delH,SLD_T,SLD_H,SLD_S,NN,Cp,bkg; //local variables of coefficient wave |
---|
| 1005 | double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ,drh,drt; |
---|
| 1006 | double Pi,Euler,dQDefault,fii; |
---|
| 1007 | int ii,NNint; |
---|
[8e91f01] | 1008 | |
---|
| 1009 | |
---|
[ae3ce4e] | 1010 | Euler = 0.5772156649; // Euler's constant |
---|
[96b59384] | 1011 | dQDefault = 0; //0.0025; //[=] 1/A, q-resolution, default value |
---|
[ae3ce4e] | 1012 | dQ = dQDefault; |
---|
[8e91f01] | 1013 | |
---|
[ae3ce4e] | 1014 | Pi = 4.0*atan(1.0); |
---|
| 1015 | qval= q; |
---|
[8e91f01] | 1016 | |
---|
[ae3ce4e] | 1017 | scale = dp[0]; |
---|
| 1018 | dd = dp[1]; |
---|
| 1019 | delT = dp[2]; |
---|
| 1020 | delH = dp[3]; |
---|
| 1021 | SLD_T = dp[4]; |
---|
| 1022 | SLD_H = dp[5]; |
---|
| 1023 | SLD_S = dp[6]; |
---|
| 1024 | NN = trunc(dp[7]); //be sure that NN is an integer |
---|
| 1025 | Cp = dp[8]; |
---|
| 1026 | bkg = dp[9]; |
---|
[8e91f01] | 1027 | |
---|
| 1028 | |
---|
[ae3ce4e] | 1029 | drh = SLD_H - SLD_S; |
---|
| 1030 | drt = SLD_T - SLD_S; //correction 13FEB06 by L.Porcar |
---|
[8e91f01] | 1031 | |
---|
[ae3ce4e] | 1032 | Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT); |
---|
| 1033 | Pq *= Pq; |
---|
| 1034 | Pq *= 4.0/(qval*qval); |
---|
[8e91f01] | 1035 | |
---|
[ae3ce4e] | 1036 | NNint = (int)NN; //cast to an integer for the loop |
---|
| 1037 | ii=0; |
---|
| 1038 | Sq = 0.0; |
---|
| 1039 | for(ii=1;ii<(NNint-1);ii+=1) { |
---|
[8e91f01] | 1040 | |
---|
[ae3ce4e] | 1041 | fii = (double)ii; //do I really need to do this? |
---|
[8e91f01] | 1042 | |
---|
[ae3ce4e] | 1043 | temp = 0.0; |
---|
| 1044 | alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler); |
---|
| 1045 | t1 = 2.0*dQ*dQ*dd*dd*alpha; |
---|
| 1046 | t2 = 2.0*qval*qval*dd*dd*alpha; |
---|
| 1047 | t3 = dQ*dQ*dd*dd*ii*ii; |
---|
[8e91f01] | 1048 | |
---|
[ae3ce4e] | 1049 | temp = 1.0-ii/NN; |
---|
| 1050 | temp *= cos(dd*qval*ii/(1.0+t1)); |
---|
| 1051 | temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); |
---|
| 1052 | temp /= sqrt(1.0+t1); |
---|
[8e91f01] | 1053 | |
---|
[ae3ce4e] | 1054 | Sq += temp; |
---|
| 1055 | } |
---|
[8e91f01] | 1056 | |
---|
[ae3ce4e] | 1057 | Sq *= 2.0; |
---|
| 1058 | Sq += 1.0; |
---|
[8e91f01] | 1059 | |
---|
[ae3ce4e] | 1060 | inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval); |
---|
[8e91f01] | 1061 | |
---|
[ae3ce4e] | 1062 | inten *= 1.0e8; // 1/A to 1/cm |
---|
[8e91f01] | 1063 | |
---|
[ae3ce4e] | 1064 | return(inten+bkg); |
---|
| 1065 | } |
---|
| 1066 | |
---|
| 1067 | /* LamellarFF_HGX : calculates the form factor of a lamellar structure - no S(q) effects included |
---|
| 1068 | but extra SLD for head groups is included |
---|
| 1069 | |
---|
| 1070 | */ |
---|
| 1071 | double |
---|
| 1072 | LamellarFF_HG(double dp[], double q) |
---|
| 1073 | { |
---|
| 1074 | double scale,delT,delH,slds,sldh,sldt,bkg; //local variables of coefficient wave |
---|
| 1075 | double inten, qval,Pq,drh,drt; |
---|
| 1076 | double Pi; |
---|
[8e91f01] | 1077 | |
---|
| 1078 | |
---|
[ae3ce4e] | 1079 | Pi = 4.0*atan(1.0); |
---|
| 1080 | qval= q; |
---|
| 1081 | scale = dp[0]; |
---|
| 1082 | delT = dp[1]; |
---|
| 1083 | delH = dp[2]; |
---|
| 1084 | sldt = dp[3]; |
---|
| 1085 | sldh = dp[4]; |
---|
| 1086 | slds = dp[5]; |
---|
| 1087 | bkg = dp[6]; |
---|
[8e91f01] | 1088 | |
---|
| 1089 | |
---|
[ae3ce4e] | 1090 | drh = sldh - slds; |
---|
| 1091 | drt = sldt - slds; //correction 13FEB06 by L.Porcar |
---|
[8e91f01] | 1092 | |
---|
[ae3ce4e] | 1093 | Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT); |
---|
| 1094 | Pq *= Pq; |
---|
| 1095 | Pq *= 4.0/(qval*qval); |
---|
[8e91f01] | 1096 | |
---|
[ae3ce4e] | 1097 | inten = 2.0*Pi*scale*Pq/(qval*qval); //dimensionless... |
---|
[8e91f01] | 1098 | |
---|
[ae3ce4e] | 1099 | inten /= 2.0*(delT+delH); //normalize by the bilayer thickness |
---|
[8e91f01] | 1100 | |
---|
[ae3ce4e] | 1101 | inten *= 1.0e8; // 1/A to 1/cm |
---|
[8e91f01] | 1102 | |
---|
[ae3ce4e] | 1103 | return(inten+bkg); |
---|
| 1104 | } |
---|
| 1105 | |
---|
| 1106 | /* FlexExclVolCylX : calculates the form factor of a flexible cylinder with a circular cross section |
---|
| 1107 | -- incorporates Wei-Ren Chen's fixes - 2006 |
---|
| 1108 | |
---|
| 1109 | */ |
---|
| 1110 | double |
---|
| 1111 | FlexExclVolCyl(double dp[], double q) |
---|
| 1112 | { |
---|
| 1113 | double scale,L,B,bkg,rad,qr,cont; |
---|
| 1114 | double Pi,flex,crossSect,answer; |
---|
[8e91f01] | 1115 | |
---|
| 1116 | |
---|
[ae3ce4e] | 1117 | Pi = 4.0*atan(1.0); |
---|
[8e91f01] | 1118 | |
---|
[ae3ce4e] | 1119 | scale = dp[0]; //make local copies in case memory moves |
---|
| 1120 | L = dp[1]; |
---|
| 1121 | B = dp[2]; |
---|
| 1122 | rad = dp[3]; |
---|
| 1123 | cont = dp[4]; |
---|
| 1124 | bkg = dp[5]; |
---|
[8e91f01] | 1125 | |
---|
| 1126 | |
---|
[ae3ce4e] | 1127 | qr = q*rad; |
---|
[8e91f01] | 1128 | |
---|
[ae3ce4e] | 1129 | flex = Sk_WR(q,L,B); |
---|
[8e91f01] | 1130 | |
---|
[ae3ce4e] | 1131 | crossSect = (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); |
---|
| 1132 | flex *= crossSect; |
---|
| 1133 | flex *= Pi*rad*rad*L; |
---|
| 1134 | flex *= cont*cont; |
---|
| 1135 | flex *= 1.0e8; |
---|
| 1136 | answer = scale*flex + bkg; |
---|
[8e91f01] | 1137 | |
---|
[ae3ce4e] | 1138 | return answer; |
---|
| 1139 | } |
---|
| 1140 | |
---|
| 1141 | /* FlexCyl_EllipX : calculates the form factor of a flexible cylinder with an elliptical cross section |
---|
| 1142 | -- incorporates Wei-Ren Chen's fixes - 2006 |
---|
| 1143 | |
---|
| 1144 | */ |
---|
| 1145 | double |
---|
| 1146 | FlexCyl_Ellip(double dp[], double q) |
---|
| 1147 | { |
---|
| 1148 | double scale,L,B,bkg,rad,qr,cont,ellRatio; |
---|
| 1149 | double Pi,flex,crossSect,answer; |
---|
[8e91f01] | 1150 | |
---|
| 1151 | |
---|
[ae3ce4e] | 1152 | Pi = 4.0*atan(1.0); |
---|
| 1153 | scale = dp[0]; //make local copies in case memory moves |
---|
| 1154 | L = dp[1]; |
---|
| 1155 | B = dp[2]; |
---|
| 1156 | rad = dp[3]; |
---|
| 1157 | ellRatio = dp[4]; |
---|
| 1158 | cont = dp[5]; |
---|
| 1159 | bkg = dp[6]; |
---|
[8e91f01] | 1160 | |
---|
[ae3ce4e] | 1161 | qr = q*rad; |
---|
[8e91f01] | 1162 | |
---|
[ae3ce4e] | 1163 | flex = Sk_WR(q,L,B); |
---|
[8e91f01] | 1164 | |
---|
[ae3ce4e] | 1165 | crossSect = EllipticalCross_fn(q,rad,(rad*ellRatio)); |
---|
| 1166 | flex *= crossSect; |
---|
| 1167 | flex *= Pi*rad*rad*ellRatio*L; |
---|
| 1168 | flex *= cont*cont; |
---|
| 1169 | flex *= 1.0e8; |
---|
| 1170 | answer = scale*flex + bkg; |
---|
[8e91f01] | 1171 | |
---|
[ae3ce4e] | 1172 | return answer; |
---|
| 1173 | } |
---|
| 1174 | |
---|
| 1175 | double |
---|
| 1176 | EllipticalCross_fn(double qq, double a, double b) |
---|
| 1177 | { |
---|
| 1178 | double uplim,lolim,Pi,summ,arg,zi,yyy,answer; |
---|
| 1179 | int i,nord=76; |
---|
[8e91f01] | 1180 | |
---|
[ae3ce4e] | 1181 | Pi = 4.0*atan(1.0); |
---|
| 1182 | lolim=0.0; |
---|
| 1183 | uplim=Pi/2.0; |
---|
| 1184 | summ=0.0; |
---|
[8e91f01] | 1185 | |
---|
[ae3ce4e] | 1186 | for(i=0;i<nord;i++) { |
---|
| 1187 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 1188 | arg = qq*sqrt(a*a*sin(zi)*sin(zi)+b*b*cos(zi)*cos(zi)); |
---|
| 1189 | yyy = pow((2.0 * NR_BessJ1(arg) / arg),2); |
---|
| 1190 | yyy *= Gauss76Wt[i]; |
---|
| 1191 | summ += yyy; |
---|
| 1192 | } |
---|
| 1193 | answer = (uplim-lolim)/2.0*summ; |
---|
| 1194 | answer *= 2.0/Pi; |
---|
| 1195 | return(answer); |
---|
[8e91f01] | 1196 | |
---|
[ae3ce4e] | 1197 | } |
---|
| 1198 | /* FlexCyl_PolyLenX : calculates the form factor of a flecible cylinder at the given x-value p->x |
---|
| 1199 | the cylinder has a polydisperse Length |
---|
| 1200 | |
---|
| 1201 | */ |
---|
| 1202 | double |
---|
| 1203 | FlexCyl_PolyLen(double dp[], double q) |
---|
| 1204 | { |
---|
| 1205 | int i; |
---|
| 1206 | double scale,radius,length,pd,delrho,bkg,lb; //local variables of coefficient wave |
---|
| 1207 | int nord=20; //order of integration |
---|
| 1208 | double uplim,lolim; //upper and lower integration limits |
---|
| 1209 | double summ,zi,yyy,answer,Vpoly; //running tally of integration |
---|
| 1210 | double range,zz,Pi; |
---|
[8e91f01] | 1211 | |
---|
[ae3ce4e] | 1212 | Pi = 4.0*atan(1.0); |
---|
| 1213 | range = 3.4; |
---|
[8e91f01] | 1214 | |
---|
[ae3ce4e] | 1215 | summ = 0.0; //initialize intergral |
---|
| 1216 | scale = dp[0]; //make local copies in case memory moves |
---|
| 1217 | length = dp[1]; //radius |
---|
| 1218 | pd = dp[2]; // average length |
---|
| 1219 | lb = dp[3]; |
---|
| 1220 | radius = dp[4]; |
---|
| 1221 | delrho = dp[5]; |
---|
| 1222 | bkg = dp[6]; |
---|
[8e91f01] | 1223 | |
---|
[ae3ce4e] | 1224 | zz = (1.0/pd)*(1.0/pd) - 1.0; |
---|
[8e91f01] | 1225 | |
---|
[ae3ce4e] | 1226 | lolim = length*(1.0-range*pd); //set the upper/lower limits to cover the distribution |
---|
| 1227 | if(lolim<0) { |
---|
| 1228 | lolim = 0; |
---|
| 1229 | } |
---|
| 1230 | if(pd>0.3) { |
---|
| 1231 | range = 3.4 + (pd-0.3)*18.0; |
---|
| 1232 | } |
---|
| 1233 | uplim = length*(1.0+range*pd); |
---|
[8e91f01] | 1234 | |
---|
[ae3ce4e] | 1235 | for(i=0;i<nord;i++) { |
---|
| 1236 | zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 1237 | yyy = Gauss20Wt[i] * FlePolyLen_kernel(q,radius,length,lb,zz,delrho,zi); |
---|
| 1238 | summ += yyy; |
---|
| 1239 | } |
---|
[8e91f01] | 1240 | |
---|
[ae3ce4e] | 1241 | answer = (uplim-lolim)/2.0*summ; |
---|
| 1242 | //normalize by average cylinder volume (first moment), using the average length |
---|
| 1243 | Vpoly=Pi*radius*radius*length; |
---|
| 1244 | answer /= Vpoly; |
---|
[8e91f01] | 1245 | |
---|
[ae3ce4e] | 1246 | answer *=delrho*delrho; |
---|
[8e91f01] | 1247 | |
---|
[ae3ce4e] | 1248 | //convert to [cm-1] |
---|
| 1249 | answer *= 1.0e8; |
---|
| 1250 | //Scale |
---|
| 1251 | answer *= scale; |
---|
| 1252 | // add in the background |
---|
| 1253 | answer += bkg; |
---|
[8e91f01] | 1254 | |
---|
[ae3ce4e] | 1255 | return answer; |
---|
| 1256 | } |
---|
| 1257 | |
---|
| 1258 | /* FlexCyl_PolyLenX : calculates the form factor of a flexible cylinder at the given x-value p->x |
---|
| 1259 | the cylinder has a polydisperse cross sectional radius |
---|
| 1260 | |
---|
| 1261 | */ |
---|
| 1262 | double |
---|
| 1263 | FlexCyl_PolyRad(double dp[], double q) |
---|
| 1264 | { |
---|
| 1265 | int i; |
---|
| 1266 | double scale,radius,length,pd,delrho,bkg,lb; //local variables of coefficient wave |
---|
| 1267 | int nord=76; //order of integration |
---|
| 1268 | double uplim,lolim; //upper and lower integration limits |
---|
| 1269 | double summ,zi,yyy,answer,Vpoly; //running tally of integration |
---|
| 1270 | double range,zz,Pi; |
---|
[8e91f01] | 1271 | |
---|
| 1272 | |
---|
[ae3ce4e] | 1273 | Pi = 4.0*atan(1.0); |
---|
| 1274 | range = 3.4; |
---|
[8e91f01] | 1275 | |
---|
[ae3ce4e] | 1276 | summ = 0.0; //initialize intergral |
---|
[8e91f01] | 1277 | |
---|
[ae3ce4e] | 1278 | scale = dp[0]; //make local copies in case memory moves |
---|
| 1279 | length = dp[1]; //radius |
---|
| 1280 | lb = dp[2]; // average length |
---|
| 1281 | radius = dp[3]; |
---|
| 1282 | pd = dp[4]; |
---|
| 1283 | delrho = dp[5]; |
---|
| 1284 | bkg = dp[6]; |
---|
[8e91f01] | 1285 | |
---|
[ae3ce4e] | 1286 | zz = (1.0/pd)*(1.0/pd) - 1.0; |
---|
[8e91f01] | 1287 | |
---|
[ae3ce4e] | 1288 | lolim = radius*(1.0-range*pd); //set the upper/lower limits to cover the distribution |
---|
| 1289 | if(lolim<0) { |
---|
| 1290 | lolim = 0; |
---|
| 1291 | } |
---|
| 1292 | if(pd>0.3) { |
---|
| 1293 | range = 3.4 + (pd-0.3)*18.0; |
---|
| 1294 | } |
---|
| 1295 | uplim = radius*(1.0+range*pd); |
---|
[8e91f01] | 1296 | |
---|
[ae3ce4e] | 1297 | for(i=0;i<nord;i++) { |
---|
| 1298 | //zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 1299 | //yyy = Gauss20Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi); |
---|
| 1300 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 1301 | yyy = Gauss76Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi); |
---|
| 1302 | summ += yyy; |
---|
| 1303 | } |
---|
[8e91f01] | 1304 | |
---|
[ae3ce4e] | 1305 | answer = (uplim-lolim)/2.0*summ; |
---|
| 1306 | //normalize by average cylinder volume (second moment), using the average radius |
---|
| 1307 | Vpoly = Pi*radius*radius*length*(zz+2.0)/(zz+1.0); |
---|
| 1308 | answer /= Vpoly; |
---|
[8e91f01] | 1309 | |
---|
[ae3ce4e] | 1310 | answer *=delrho*delrho; |
---|
[8e91f01] | 1311 | |
---|
[ae3ce4e] | 1312 | //convert to [cm-1] |
---|
| 1313 | answer *= 1.0e8; |
---|
| 1314 | //Scale |
---|
| 1315 | answer *= scale; |
---|
| 1316 | // add in the background |
---|
| 1317 | answer += bkg; |
---|
[8e91f01] | 1318 | |
---|
[ae3ce4e] | 1319 | return answer; |
---|
| 1320 | } |
---|
| 1321 | |
---|
| 1322 | /////////functions for WRC implementation of flexible cylinders |
---|
| 1323 | static double |
---|
| 1324 | Sk_WR(double q, double L, double b) |
---|
| 1325 | { |
---|
| 1326 | // |
---|
| 1327 | double p1,p2,p1short,p2short,q0,qconnect; |
---|
| 1328 | double C,epsilon,ans,q0short,Sexvmodify,pi; |
---|
[8e91f01] | 1329 | |
---|
[ae3ce4e] | 1330 | pi = 4.0*atan(1.0); |
---|
[8e91f01] | 1331 | |
---|
[ae3ce4e] | 1332 | p1 = 4.12; |
---|
| 1333 | p2 = 4.42; |
---|
| 1334 | p1short = 5.36; |
---|
| 1335 | p2short = 5.62; |
---|
| 1336 | q0 = 3.1; |
---|
| 1337 | qconnect = q0/b; |
---|
[8e91f01] | 1338 | // |
---|
[ae3ce4e] | 1339 | q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); |
---|
[8e91f01] | 1340 | |
---|
[ae3ce4e] | 1341 | // |
---|
| 1342 | if(L/b > 10.0) { |
---|
| 1343 | C = 3.06/pow((L/b),0.44); |
---|
| 1344 | epsilon = 0.176; |
---|
| 1345 | } else { |
---|
| 1346 | C = 1.0; |
---|
| 1347 | epsilon = 0.170; |
---|
| 1348 | } |
---|
| 1349 | // |
---|
[8e91f01] | 1350 | |
---|
[ae3ce4e] | 1351 | if( L > 4*b ) { // Longer Chains |
---|
| 1352 | if (q*b <= 3.1) { //Modified by Yun on Oct. 15, |
---|
| 1353 | Sexvmodify = Sexvnew(q, L, b); |
---|
| 1354 | 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); |
---|
| 1355 | } else { //q(i)*b > 3.1 |
---|
| 1356 | 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); |
---|
[8e91f01] | 1357 | } |
---|
[ae3ce4e] | 1358 | } else { //L <= 4*b Shorter Chains |
---|
| 1359 | if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { |
---|
| 1360 | if (q*b<=0.01) { |
---|
| 1361 | ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0; |
---|
| 1362 | } else { |
---|
| 1363 | ans = Sdebye1(q,L,b); |
---|
| 1364 | } |
---|
| 1365 | } else { //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) |
---|
| 1366 | 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); |
---|
| 1367 | } |
---|
| 1368 | } |
---|
[8e91f01] | 1369 | |
---|
[ae3ce4e] | 1370 | return(ans); |
---|
| 1371 | //return(a2long(q, L, b, p1, p2, q0)); |
---|
| 1372 | } |
---|
| 1373 | |
---|
| 1374 | //WR named this w (too generic) |
---|
| 1375 | static double |
---|
| 1376 | w_WR(double x) |
---|
| 1377 | { |
---|
| 1378 | double yy; |
---|
| 1379 | yy = 0.5*(1 + tanh((x - 1.523)/0.1477)); |
---|
[8e91f01] | 1380 | |
---|
[ae3ce4e] | 1381 | return (yy); |
---|
| 1382 | } |
---|
| 1383 | |
---|
| 1384 | // |
---|
| 1385 | static double |
---|
| 1386 | u1(double q, double L, double b) |
---|
| 1387 | { |
---|
| 1388 | double yy; |
---|
[8e91f01] | 1389 | |
---|
[ae3ce4e] | 1390 | yy = Rgsquareshort(q,L,b)*q*q; |
---|
[8e91f01] | 1391 | |
---|
[ae3ce4e] | 1392 | return (yy); |
---|
| 1393 | } |
---|
| 1394 | |
---|
| 1395 | // was named u |
---|
| 1396 | static double |
---|
| 1397 | u_WR(double q, double L, double b) |
---|
| 1398 | { |
---|
| 1399 | double yy; |
---|
| 1400 | yy = Rgsquare(q,L,b)*q*q; |
---|
| 1401 | return (yy); |
---|
| 1402 | } |
---|
| 1403 | |
---|
| 1404 | |
---|
| 1405 | |
---|
| 1406 | // |
---|
| 1407 | static double |
---|
| 1408 | Rgsquarezero(double q, double L, double b) |
---|
[8e91f01] | 1409 | { |
---|
[ae3ce4e] | 1410 | double yy; |
---|
| 1411 | 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)))); |
---|
[8e91f01] | 1412 | |
---|
[ae3ce4e] | 1413 | return (yy); |
---|
| 1414 | } |
---|
| 1415 | |
---|
| 1416 | // |
---|
| 1417 | static double |
---|
| 1418 | Rgsquareshort(double q, double L, double b) |
---|
[8e91f01] | 1419 | { |
---|
[ae3ce4e] | 1420 | double yy; |
---|
| 1421 | yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b); |
---|
[8e91f01] | 1422 | |
---|
[ae3ce4e] | 1423 | return (yy); |
---|
| 1424 | } |
---|
| 1425 | |
---|
| 1426 | // |
---|
| 1427 | static double |
---|
| 1428 | Rgsquare(double q, double L, double b) |
---|
| 1429 | { |
---|
| 1430 | double yy; |
---|
| 1431 | yy = AlphaSquare(L/b)*L*b/6.0; |
---|
[8e91f01] | 1432 | |
---|
[ae3ce4e] | 1433 | return (yy); |
---|
| 1434 | } |
---|
| 1435 | |
---|
| 1436 | // |
---|
| 1437 | static double |
---|
| 1438 | AlphaSquare(double x) |
---|
[8e91f01] | 1439 | { |
---|
[ae3ce4e] | 1440 | double yy; |
---|
| 1441 | yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); |
---|
[8e91f01] | 1442 | |
---|
[ae3ce4e] | 1443 | return (yy); |
---|
| 1444 | } |
---|
| 1445 | |
---|
| 1446 | // ?? funciton is not used - but should the log actually be log10??? |
---|
| 1447 | static double |
---|
| 1448 | miu(double x) |
---|
[8e91f01] | 1449 | { |
---|
[ae3ce4e] | 1450 | double yy; |
---|
| 1451 | 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))); |
---|
[8e91f01] | 1452 | |
---|
[ae3ce4e] | 1453 | return (yy); |
---|
| 1454 | } |
---|
| 1455 | |
---|
| 1456 | // |
---|
| 1457 | static double |
---|
| 1458 | Sdebye(double q, double L, double b) |
---|
[8e91f01] | 1459 | { |
---|
[ae3ce4e] | 1460 | double yy; |
---|
| 1461 | yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2)); |
---|
[8e91f01] | 1462 | |
---|
[ae3ce4e] | 1463 | return (yy); |
---|
| 1464 | } |
---|
| 1465 | |
---|
| 1466 | // |
---|
| 1467 | static double |
---|
| 1468 | Sdebye1(double q, double L, double b) |
---|
[8e91f01] | 1469 | { |
---|
[ae3ce4e] | 1470 | double yy; |
---|
| 1471 | yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) ); |
---|
[8e91f01] | 1472 | |
---|
[ae3ce4e] | 1473 | return (yy); |
---|
| 1474 | } |
---|
| 1475 | |
---|
| 1476 | // |
---|
| 1477 | static double |
---|
| 1478 | Sexv(double q, double L, double b) |
---|
[8e91f01] | 1479 | { |
---|
[ae3ce4e] | 1480 | double yy,C1,C2,C3,miu,Rg2; |
---|
| 1481 | C1=1.22; |
---|
| 1482 | C2=0.4288; |
---|
| 1483 | C3=-1.651; |
---|
| 1484 | miu = 0.585; |
---|
[8e91f01] | 1485 | |
---|
[ae3ce4e] | 1486 | Rg2 = Rgsquare(q,L,b); |
---|
[8e91f01] | 1487 | |
---|
[ae3ce4e] | 1488 | 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))); |
---|
[8e91f01] | 1489 | |
---|
[ae3ce4e] | 1490 | return (yy); |
---|
| 1491 | } |
---|
| 1492 | |
---|
| 1493 | // this must be WR modified version |
---|
| 1494 | static double |
---|
| 1495 | Sexvnew(double q, double L, double b) |
---|
[8e91f01] | 1496 | { |
---|
[ae3ce4e] | 1497 | double yy,C1,C2,C3,miu; |
---|
| 1498 | double del=1.05,C_star2,Rg2; |
---|
[8e91f01] | 1499 | |
---|
[ae3ce4e] | 1500 | C1=1.22; |
---|
| 1501 | C2=0.4288; |
---|
| 1502 | C3=-1.651; |
---|
| 1503 | miu = 0.585; |
---|
[8e91f01] | 1504 | |
---|
[ae3ce4e] | 1505 | //calculating the derivative to decide on the corection (cutoff) term? |
---|
| 1506 | // I have modified this from WRs original code |
---|
[8e91f01] | 1507 | |
---|
[ae3ce4e] | 1508 | if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) { |
---|
| 1509 | C_star2 = 0.0; |
---|
| 1510 | } else { |
---|
| 1511 | C_star2 = 1.0; |
---|
| 1512 | } |
---|
[8e91f01] | 1513 | |
---|
[ae3ce4e] | 1514 | Rg2 = Rgsquare(q,L,b); |
---|
[8e91f01] | 1515 | |
---|
[ae3ce4e] | 1516 | 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))); |
---|
[8e91f01] | 1517 | |
---|
[ae3ce4e] | 1518 | return (yy); |
---|
| 1519 | } |
---|
| 1520 | |
---|
| 1521 | // these are the messy ones |
---|
| 1522 | static double |
---|
| 1523 | a2short(double q, double L, double b, double p1short, double p2short, double q0) |
---|
[8e91f01] | 1524 | { |
---|
[ae3ce4e] | 1525 | double yy,Rg2_sh; |
---|
| 1526 | double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p; |
---|
| 1527 | double pi; |
---|
[8e91f01] | 1528 | |
---|
[ae3ce4e] | 1529 | E = 2.718281828459045091; |
---|
| 1530 | pi = 4.0*atan(1.0); |
---|
| 1531 | Rg2_sh = Rgsquareshort(q,L,b); |
---|
| 1532 | Rg2_sh2 = Rg2_sh*Rg2_sh; |
---|
| 1533 | t1 = ((q0*q0*Rg2_sh)/(b*b)); |
---|
| 1534 | Et1 = pow(E,t1); |
---|
[8e91f01] | 1535 | Emt1 =pow(E,-t1); |
---|
[ae3ce4e] | 1536 | q02 = q0*q0; |
---|
| 1537 | q0p = pow(q0,(-4.0 + p2short) ); |
---|
[8e91f01] | 1538 | |
---|
[ae3ce4e] | 1539 | //E is the number e |
---|
| 1540 | 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))))))); |
---|
[8e91f01] | 1541 | |
---|
[ae3ce4e] | 1542 | return (yy); |
---|
| 1543 | } |
---|
| 1544 | |
---|
| 1545 | // |
---|
| 1546 | static double |
---|
| 1547 | a1short(double q, double L, double b, double p1short, double p2short, double q0) |
---|
[8e91f01] | 1548 | { |
---|
[ae3ce4e] | 1549 | double yy,Rg2_sh; |
---|
| 1550 | double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3; |
---|
| 1551 | double pi; |
---|
[8e91f01] | 1552 | |
---|
[ae3ce4e] | 1553 | E = 2.718281828459045091; |
---|
| 1554 | pi = 4.0*atan(1.0); |
---|
| 1555 | Rg2_sh = Rgsquareshort(q,L,b); |
---|
| 1556 | Rg2_sh2 = Rg2_sh*Rg2_sh; |
---|
| 1557 | b3 = b*b*b; |
---|
| 1558 | t1 = ((q0*q0*Rg2_sh)/(b*b)); |
---|
| 1559 | Et1 = pow(E,t1); |
---|
[8e91f01] | 1560 | Emt1 =pow(E,-t1); |
---|
[ae3ce4e] | 1561 | q02 = q0*q0; |
---|
| 1562 | q0p = pow(q0,(-4.0 + p1short) ); |
---|
[8e91f01] | 1563 | |
---|
| 1564 | 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)))))); |
---|
| 1565 | |
---|
[ae3ce4e] | 1566 | return(yy); |
---|
| 1567 | } |
---|
| 1568 | |
---|
| 1569 | // this one will be lots of trouble |
---|
| 1570 | static double |
---|
| 1571 | a2long(double q, double L, double b, double p1, double p2, double q0) |
---|
| 1572 | { |
---|
| 1573 | double yy,C1,C2,C3,C4,C5,miu,C,Rg2; |
---|
| 1574 | double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi; |
---|
| 1575 | double E,b2,b3,b4,q02,q03,q04,q05,Rg22; |
---|
[8e91f01] | 1576 | |
---|
[ae3ce4e] | 1577 | pi = 4.0*atan(1.0); |
---|
| 1578 | E = 2.718281828459045091; |
---|
| 1579 | if( L/b > 10.0) { |
---|
| 1580 | C = 3.06/pow((L/b),0.44); |
---|
| 1581 | } else { |
---|
| 1582 | C = 1.0; |
---|
| 1583 | } |
---|
[8e91f01] | 1584 | |
---|
[ae3ce4e] | 1585 | C1 = 1.22; |
---|
| 1586 | C2 = 0.4288; |
---|
| 1587 | C3 = -1.651; |
---|
| 1588 | C4 = 1.523; |
---|
| 1589 | C5 = 0.1477; |
---|
| 1590 | miu = 0.585; |
---|
[8e91f01] | 1591 | |
---|
[ae3ce4e] | 1592 | Rg2 = Rgsquare(q,L,b); |
---|
| 1593 | Rg22 = Rg2*Rg2; |
---|
| 1594 | b2 = b*b; |
---|
| 1595 | b3 = b*b*b; |
---|
| 1596 | b4 = b3*b; |
---|
| 1597 | q02 = q0*q0; |
---|
| 1598 | q03 = q0*q0*q0; |
---|
| 1599 | q04 = q03*q0; |
---|
| 1600 | q05 = q04*q0; |
---|
[8e91f01] | 1601 | |
---|
[ae3ce4e] | 1602 | t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) )); |
---|
[8e91f01] | 1603 | |
---|
[ae3ce4e] | 1604 | t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14*b3*pow(E,(-((q02*Rg2)/b2))))/(15*q03*Rg2) + (2*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15*q02*Rg2)))*Rg2)/b)))/L; |
---|
[8e91f01] | 1605 | |
---|
[ae3ce4e] | 1606 | t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(2*C5); |
---|
[8e91f01] | 1607 | |
---|
[ae3ce4e] | 1608 | t4 = (b4*sqrt(Rg2)*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22); |
---|
[8e91f01] | 1609 | |
---|
[ae3ce4e] | 1610 | t5 = (2*b4*(((2*q0*Rg2)/b - (2*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); |
---|
[8e91f01] | 1611 | |
---|
[ae3ce4e] | 1612 | t6 = (8*b4*b*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22); |
---|
[8e91f01] | 1613 | |
---|
[ae3ce4e] | 1614 | t7 = (((-((3*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 3/miu)))/miu)) - (2*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 2/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 1/miu)))/miu)); |
---|
[8e91f01] | 1615 | |
---|
[ae3ce4e] | 1616 | t8 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); |
---|
[8e91f01] | 1617 | |
---|
[ae3ce4e] | 1618 | t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7*b2)/(15*q02*Rg2))) + (7*b2)/(15*q02*Rg2))))/L; |
---|
[8e91f01] | 1619 | |
---|
[ae3ce4e] | 1620 | t10 = (2*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); |
---|
[8e91f01] | 1621 | |
---|
| 1622 | |
---|
[ae3ce4e] | 1623 | yy = ((-1*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1)/miu))))*((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))); |
---|
[8e91f01] | 1624 | |
---|
[ae3ce4e] | 1625 | return (yy); |
---|
| 1626 | } |
---|
| 1627 | |
---|
| 1628 | //need to define this on my own |
---|
| 1629 | static double |
---|
| 1630 | sech_WR(double x) |
---|
[8e91f01] | 1631 | { |
---|
[ae3ce4e] | 1632 | return(1/cosh(x)); |
---|
| 1633 | } |
---|
| 1634 | |
---|
| 1635 | // |
---|
| 1636 | static double |
---|
| 1637 | a1long(double q, double L, double b, double p1, double p2, double q0) |
---|
[8e91f01] | 1638 | { |
---|
[ae3ce4e] | 1639 | double yy,C,C1,C2,C3,C4,C5,miu,Rg2; |
---|
| 1640 | double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15; |
---|
| 1641 | double E,pi; |
---|
| 1642 | double b2,b3,b4,q02,q03,q04,q05,Rg22; |
---|
[8e91f01] | 1643 | |
---|
[ae3ce4e] | 1644 | pi = 4.0*atan(1.0); |
---|
| 1645 | E = 2.718281828459045091; |
---|
[8e91f01] | 1646 | |
---|
[ae3ce4e] | 1647 | if( L/b > 10.0) { |
---|
| 1648 | C = 3.06/pow((L/b),0.44); |
---|
| 1649 | } else { |
---|
| 1650 | C = 1.0; |
---|
| 1651 | } |
---|
[8e91f01] | 1652 | |
---|
[ae3ce4e] | 1653 | C1 = 1.22; |
---|
| 1654 | C2 = 0.4288; |
---|
| 1655 | C3 = -1.651; |
---|
| 1656 | C4 = 1.523; |
---|
| 1657 | C5 = 0.1477; |
---|
| 1658 | miu = 0.585; |
---|
[8e91f01] | 1659 | |
---|
[ae3ce4e] | 1660 | Rg2 = Rgsquare(q,L,b); |
---|
| 1661 | Rg22 = Rg2*Rg2; |
---|
| 1662 | b2 = b*b; |
---|
| 1663 | b3 = b*b*b; |
---|
| 1664 | b4 = b3*b; |
---|
| 1665 | q02 = q0*q0; |
---|
| 1666 | q03 = q0*q0*q0; |
---|
| 1667 | q04 = q03*q0; |
---|
| 1668 | q05 = q04*q0; |
---|
[8e91f01] | 1669 | |
---|
[ae3ce4e] | 1670 | 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)))); |
---|
[8e91f01] | 1671 | |
---|
[ae3ce4e] | 1672 | 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)))))); |
---|
[8e91f01] | 1673 | |
---|
[ae3ce4e] | 1674 | 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)))); |
---|
[8e91f01] | 1675 | |
---|
[ae3ce4e] | 1676 | t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); |
---|
[8e91f01] | 1677 | |
---|
[ae3ce4e] | 1678 | t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)))); |
---|
[8e91f01] | 1679 | |
---|
[ae3ce4e] | 1680 | 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))); |
---|
[8e91f01] | 1681 | |
---|
[ae3ce4e] | 1682 | 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)); |
---|
[8e91f01] | 1683 | |
---|
[ae3ce4e] | 1684 | t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); |
---|
[8e91f01] | 1685 | |
---|
[ae3ce4e] | 1686 | 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)))))); |
---|
[8e91f01] | 1687 | |
---|
[ae3ce4e] | 1688 | 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)))))); |
---|
[8e91f01] | 1689 | |
---|
[ae3ce4e] | 1690 | 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)); |
---|
[8e91f01] | 1691 | |
---|
[ae3ce4e] | 1692 | t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); |
---|
[8e91f01] | 1693 | |
---|
[ae3ce4e] | 1694 | 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)))); |
---|
[8e91f01] | 1695 | |
---|
[ae3ce4e] | 1696 | 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)))))); |
---|
[8e91f01] | 1697 | |
---|
[ae3ce4e] | 1698 | 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)))); |
---|
[8e91f01] | 1699 | |
---|
| 1700 | |
---|
[ae3ce4e] | 1701 | 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))))))))))); |
---|
[8e91f01] | 1702 | |
---|
[ae3ce4e] | 1703 | return (yy); |
---|
| 1704 | } |
---|
| 1705 | |
---|
| 1706 | |
---|
| 1707 | |
---|
| 1708 | /////////////// |
---|
| 1709 | |
---|
| 1710 | // |
---|
| 1711 | // FUNCTION gfn2: CONTAINS F(Q,A,B,mu)**2 AS GIVEN |
---|
[8e91f01] | 1712 | // BY (53) AND (56,57) IN CHEN AND |
---|
[ae3ce4e] | 1713 | // KOTLARCHYK REFERENCE |
---|
| 1714 | // |
---|
| 1715 | // <PROLATE ELLIPSOIDS> |
---|
| 1716 | // |
---|
| 1717 | double |
---|
| 1718 | gfn2(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq) |
---|
| 1719 | { |
---|
| 1720 | // local variables |
---|
[975ec8e] | 1721 | double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,gfn2,pi43,gfn,Pi; |
---|
[8e91f01] | 1722 | |
---|
[ae3ce4e] | 1723 | Pi = 4.0*atan(1.0); |
---|
[8e91f01] | 1724 | |
---|
[ae3ce4e] | 1725 | pi43=4.0/3.0*Pi; |
---|
| 1726 | aa = crmaj; |
---|
| 1727 | bb = crmin; |
---|
| 1728 | u2 = (aa*aa*xx*xx + bb*bb*(1.0-xx*xx)); |
---|
| 1729 | ut2 = (trmaj*trmaj*xx*xx + trmin*trmin*(1.0-xx*xx)); |
---|
| 1730 | uq = sqrt(u2)*qq; |
---|
| 1731 | ut= sqrt(ut2)*qq; |
---|
| 1732 | vc = pi43*aa*bb*bb; |
---|
| 1733 | vt = pi43*trmaj*trmin*trmin; |
---|
[975ec8e] | 1734 | if (uq == 0){ |
---|
| 1735 | siq = 1/3; |
---|
| 1736 | }else{ |
---|
| 1737 | siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq; |
---|
| 1738 | } |
---|
| 1739 | if (ut == 0){ |
---|
| 1740 | sit = 1/3; |
---|
| 1741 | }else{ |
---|
| 1742 | sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut; |
---|
| 1743 | } |
---|
| 1744 | gfnc = 3.0*siq*vc*delpc; |
---|
| 1745 | gfnt = 3.0*sit*vt*delps; |
---|
[ae3ce4e] | 1746 | gfn = gfnc+gfnt; |
---|
| 1747 | gfn2 = gfn*gfn; |
---|
[8e91f01] | 1748 | |
---|
[ae3ce4e] | 1749 | return (gfn2); |
---|
| 1750 | } |
---|
| 1751 | |
---|
| 1752 | // |
---|
| 1753 | // FUNCTION gfn4: CONTAINS F(Q,A,B,MU)**2 AS GIVEN |
---|
| 1754 | // BY (53) & (58-59) IN CHEN AND |
---|
| 1755 | // KOTLARCHYK REFERENCE |
---|
| 1756 | // |
---|
| 1757 | // <OBLATE ELLIPSOID> |
---|
[8e91f01] | 1758 | // function gfn4 for oblate ellipsoids |
---|
[ae3ce4e] | 1759 | double |
---|
| 1760 | gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq) |
---|
| 1761 | { |
---|
| 1762 | // local variables |
---|
[975ec8e] | 1763 | double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi; |
---|
[8e91f01] | 1764 | |
---|
[ae3ce4e] | 1765 | Pi = 4.0*atan(1.0); |
---|
| 1766 | pi43=4.0/3.0*Pi; |
---|
| 1767 | aa = crmaj; |
---|
| 1768 | bb = crmin; |
---|
| 1769 | u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx)); |
---|
| 1770 | ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx)); |
---|
| 1771 | uq = sqrt(u2)*qq; |
---|
| 1772 | ut= sqrt(ut2)*qq; |
---|
| 1773 | vc = pi43*aa*aa*bb; |
---|
| 1774 | vt = pi43*trmaj*trmaj*trmin; |
---|
[975ec8e] | 1775 | if (uq == 0){ |
---|
| 1776 | siq = 1/3; |
---|
| 1777 | }else{ |
---|
| 1778 | siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq; |
---|
| 1779 | } |
---|
| 1780 | if (ut == 0){ |
---|
| 1781 | sit = 1/3; |
---|
| 1782 | }else{ |
---|
| 1783 | sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut; |
---|
| 1784 | } |
---|
| 1785 | gfnc = 3.0*siq*vc*delpc; |
---|
| 1786 | gfnt = 3.0*sit*vt*delps; |
---|
[ae3ce4e] | 1787 | tgfn = gfnc+gfnt; |
---|
| 1788 | gfn4 = tgfn*tgfn; |
---|
[8e91f01] | 1789 | |
---|
[ae3ce4e] | 1790 | return (gfn4); |
---|
| 1791 | } |
---|
| 1792 | |
---|
| 1793 | double |
---|
| 1794 | FlePolyLen_kernel(double q, double radius, double length, double lb, double zz, double delrho, double zi) |
---|
| 1795 | { |
---|
| 1796 | double Pq,vcyl,dl; |
---|
| 1797 | double Pi,qr; |
---|
[8e91f01] | 1798 | |
---|
[ae3ce4e] | 1799 | Pi = 4.0*atan(1.0); |
---|
| 1800 | qr = q*radius; |
---|
[8e91f01] | 1801 | |
---|
[ae3ce4e] | 1802 | Pq = Sk_WR(q,zi,lb); //does not have cross section term |
---|
[975ec8e] | 1803 | if (qr !=0){ |
---|
| 1804 | Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); |
---|
| 1805 | } //else Pk *=1; |
---|
[ae3ce4e] | 1806 | vcyl=Pi*radius*radius*zi; |
---|
| 1807 | Pq *= vcyl*vcyl; |
---|
[8e91f01] | 1808 | |
---|
[ae3ce4e] | 1809 | dl = SchulzPoint_cpr(zi,length,zz); |
---|
[8e91f01] | 1810 | return (Pq*dl); |
---|
| 1811 | |
---|
[ae3ce4e] | 1812 | } |
---|
| 1813 | |
---|
| 1814 | double |
---|
| 1815 | FlePolyRad_kernel(double q, double ravg, double Lc, double Lb, double zz, double delrho, double zi) |
---|
| 1816 | { |
---|
| 1817 | double Pq,vcyl,dr; |
---|
| 1818 | double Pi,qr; |
---|
[8e91f01] | 1819 | |
---|
[ae3ce4e] | 1820 | Pi = 4.0*atan(1.0); |
---|
| 1821 | qr = q*zi; |
---|
[8e91f01] | 1822 | |
---|
[ae3ce4e] | 1823 | Pq = Sk_WR(q,Lc,Lb); //does not have cross section term |
---|
[975ec8e] | 1824 | if (qr !=0){ |
---|
| 1825 | Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); |
---|
| 1826 | } |
---|
[8e91f01] | 1827 | |
---|
[ae3ce4e] | 1828 | vcyl=Pi*zi*zi*Lc; |
---|
| 1829 | Pq *= vcyl*vcyl; |
---|
[8e91f01] | 1830 | |
---|
[ae3ce4e] | 1831 | dr = SchulzPoint_cpr(zi,ravg,zz); |
---|
[8e91f01] | 1832 | return (Pq*dr); |
---|
| 1833 | |
---|
[ae3ce4e] | 1834 | } |
---|
| 1835 | |
---|
| 1836 | double |
---|
| 1837 | CSCylIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length) |
---|
| 1838 | { |
---|
| 1839 | double answer,halfheight,Pi; |
---|
| 1840 | double lolim,uplim,summ,yyy,zi; |
---|
| 1841 | int nord,i; |
---|
[8e91f01] | 1842 | |
---|
| 1843 | // set up the integration end points |
---|
[ae3ce4e] | 1844 | Pi = 4.0*atan(1.0); |
---|
| 1845 | nord = 76; |
---|
| 1846 | lolim = 0; |
---|
| 1847 | uplim = Pi/2; |
---|
| 1848 | halfheight = length/2.0; |
---|
[8e91f01] | 1849 | |
---|
[ae3ce4e] | 1850 | summ = 0.0; // initialize integral |
---|
| 1851 | i=0; |
---|
| 1852 | for(i=0;i<nord;i++) { |
---|
| 1853 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 1854 | yyy = Gauss76Wt[i] * CScyl(qq, rad, radthick, facthick, rhoc,rhos,rhosolv, halfheight, zi); |
---|
| 1855 | summ += yyy; |
---|
| 1856 | } |
---|
[8e91f01] | 1857 | |
---|
[ae3ce4e] | 1858 | // calculate value of integral to return |
---|
| 1859 | answer = (uplim-lolim)/2.0*summ; |
---|
| 1860 | return (answer); |
---|
| 1861 | } |
---|
| 1862 | |
---|
| 1863 | double |
---|
| 1864 | CScyl(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length, double dum) |
---|
[8e91f01] | 1865 | { |
---|
[ae3ce4e] | 1866 | // qq is the q-value for the calculation (1/A) |
---|
| 1867 | // radius is the core radius of the cylinder (A) |
---|
| 1868 | // radthick and facthick are the radial and face layer thicknesses |
---|
| 1869 | // rho(n) are the respective SLD's |
---|
[8e91f01] | 1870 | // length is the *Half* CORE-LENGTH of the cylinder |
---|
[ae3ce4e] | 1871 | // dum is the dummy variable for the integration (theta) |
---|
[8e91f01] | 1872 | |
---|
[975ec8e] | 1873 | double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval; |
---|
[ae3ce4e] | 1874 | double Pi; |
---|
[8e91f01] | 1875 | |
---|
| 1876 | Pi = 4.0*atan(1.0); |
---|
| 1877 | |
---|
[ae3ce4e] | 1878 | dr1 = rhoc-rhos; |
---|
| 1879 | dr2 = rhos-rhosolv; |
---|
| 1880 | vol1 = Pi*rad*rad*(2.0*length); |
---|
| 1881 | vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick); |
---|
[8e91f01] | 1882 | |
---|
[ae3ce4e] | 1883 | besarg1 = qq*rad*sin(dum); |
---|
| 1884 | besarg2 = qq*(rad+radthick)*sin(dum); |
---|
| 1885 | sinarg1 = qq*length*cos(dum); |
---|
| 1886 | sinarg2 = qq*(length+facthick)*cos(dum); |
---|
[975ec8e] | 1887 | if (besarg1 == 0){ |
---|
| 1888 | be1 = 0.5; |
---|
| 1889 | }else{ |
---|
| 1890 | be1 = NR_BessJ1(besarg1)/besarg1; |
---|
| 1891 | } |
---|
| 1892 | if (besarg2 == 0){ |
---|
| 1893 | be2 = 0.5; |
---|
| 1894 | }else{ |
---|
| 1895 | be2 = NR_BessJ1(besarg2)/besarg2; |
---|
| 1896 | } |
---|
| 1897 | if (sinarg1 == 0){ |
---|
| 1898 | si1 = 1; |
---|
| 1899 | }else{ |
---|
| 1900 | si1 = sin(sinarg1)/sinarg1; |
---|
| 1901 | } |
---|
| 1902 | if (besarg2 == 0){ |
---|
| 1903 | si2 = 1; |
---|
| 1904 | }else{ |
---|
| 1905 | si2 = sin(sinarg2)/sinarg2; |
---|
| 1906 | } |
---|
[8e91f01] | 1907 | |
---|
[975ec8e] | 1908 | t1 = 2.0*vol1*dr1*si1*be1; |
---|
| 1909 | t2 = 2.0*vol2*dr2*si2*be2; |
---|
[8e91f01] | 1910 | |
---|
[ae3ce4e] | 1911 | retval = ((t1+t2)*(t1+t2))*sin(dum); |
---|
| 1912 | return (retval); |
---|
[8e91f01] | 1913 | |
---|
[ae3ce4e] | 1914 | } |
---|
| 1915 | |
---|
| 1916 | |
---|
| 1917 | double |
---|
| 1918 | CoreShellCylKernel(double qq, double rcore, double thick, double rhoc, double rhos, double rhosolv, double length, double dum) |
---|
| 1919 | { |
---|
[8e91f01] | 1920 | |
---|
[975ec8e] | 1921 | double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval; |
---|
[ae3ce4e] | 1922 | double Pi; |
---|
[8e91f01] | 1923 | |
---|
[ae3ce4e] | 1924 | Pi = 4.0*atan(1.0); |
---|
[8e91f01] | 1925 | |
---|
[ae3ce4e] | 1926 | dr1 = rhoc-rhos; |
---|
| 1927 | dr2 = rhos-rhosolv; |
---|
| 1928 | vol1 = Pi*rcore*rcore*(2.0*length); |
---|
| 1929 | vol2 = Pi*(rcore+thick)*(rcore+thick)*(2.0*length+2.0*thick); |
---|
[8e91f01] | 1930 | |
---|
[ae3ce4e] | 1931 | besarg1 = qq*rcore*sin(dum); |
---|
| 1932 | besarg2 = qq*(rcore+thick)*sin(dum); |
---|
| 1933 | sinarg1 = qq*length*cos(dum); |
---|
| 1934 | sinarg2 = qq*(length+thick)*cos(dum); |
---|
[8e91f01] | 1935 | |
---|
[975ec8e] | 1936 | if (besarg1 == 0){ |
---|
| 1937 | be1 = 0.5; |
---|
| 1938 | }else{ |
---|
| 1939 | be1 = NR_BessJ1(besarg1)/besarg1; |
---|
| 1940 | } |
---|
| 1941 | if (besarg2 == 0){ |
---|
| 1942 | be2 = 0.5; |
---|
| 1943 | }else{ |
---|
| 1944 | be2 = NR_BessJ1(besarg2)/besarg2; |
---|
| 1945 | } |
---|
| 1946 | if (sinarg1 == 0){ |
---|
| 1947 | si1 = 1; |
---|
| 1948 | }else{ |
---|
| 1949 | si1 = sin(sinarg1)/sinarg1; |
---|
| 1950 | } |
---|
| 1951 | if (besarg2 == 0){ |
---|
| 1952 | si2 = 1; |
---|
| 1953 | }else{ |
---|
| 1954 | si2 = sin(sinarg2)/sinarg2; |
---|
| 1955 | } |
---|
| 1956 | |
---|
| 1957 | t1 = 2.0*vol1*dr1*si1*be1; |
---|
| 1958 | t2 = 2.0*vol2*dr2*si2*be2; |
---|
[8e91f01] | 1959 | |
---|
[ae3ce4e] | 1960 | retval = ((t1+t2)*(t1+t2))*sin(dum); |
---|
[8e91f01] | 1961 | |
---|
[ae3ce4e] | 1962 | return (retval); |
---|
| 1963 | } |
---|
| 1964 | |
---|
| 1965 | double |
---|
| 1966 | Cyl_PolyLenKernel(double q, double radius, double len_avg, double zz, double delrho, double dumLen) |
---|
| 1967 | { |
---|
[8e91f01] | 1968 | |
---|
[ae3ce4e] | 1969 | double halfheight,uplim,lolim,zi,summ,yyy,Pi; |
---|
| 1970 | double answer,dr,Vcyl; |
---|
| 1971 | int i,nord; |
---|
[8e91f01] | 1972 | |
---|
[ae3ce4e] | 1973 | Pi = 4.0*atan(1.0); |
---|
| 1974 | lolim = 0; |
---|
| 1975 | uplim = Pi/2.0; |
---|
| 1976 | halfheight = dumLen/2.0; |
---|
| 1977 | nord=20; |
---|
| 1978 | summ = 0.0; |
---|
[8e91f01] | 1979 | |
---|
[ae3ce4e] | 1980 | //do the cylinder orientational average |
---|
| 1981 | for(i=0;i<nord;i++) { |
---|
| 1982 | zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 1983 | yyy = Gauss20Wt[i] * CylKernel(q, radius, halfheight, zi); |
---|
| 1984 | summ += yyy; |
---|
| 1985 | } |
---|
| 1986 | answer = (uplim-lolim)/2.0*summ; |
---|
| 1987 | // Multiply by contrast^2 |
---|
| 1988 | answer *= delrho*delrho; |
---|
| 1989 | // don't do the normal scaling to volume here |
---|
| 1990 | // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder |
---|
| 1991 | Vcyl = Pi*radius*radius*dumLen; |
---|
| 1992 | answer *= Vcyl*Vcyl; |
---|
[8e91f01] | 1993 | |
---|
[ae3ce4e] | 1994 | dr = SchulzPoint_cpr(dumLen,len_avg,zz); |
---|
| 1995 | return(dr*answer); |
---|
| 1996 | } |
---|
| 1997 | |
---|
| 1998 | |
---|
| 1999 | double |
---|
| 2000 | Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N) |
---|
[8e91f01] | 2001 | { |
---|
[ae3ce4e] | 2002 | // qq is the q-value for the calculation (1/A) |
---|
| 2003 | // rcore is the core radius of the cylinder (A) |
---|
| 2004 | // rho(n) are the respective SLD's |
---|
| 2005 | // length is the *Half* CORE-LENGTH of the cylinder = L (A) |
---|
| 2006 | // dum is the dummy variable for the integration (x in Feigin's notation) |
---|
[8e91f01] | 2007 | |
---|
| 2008 | //Local variables |
---|
[975ec8e] | 2009 | double totald,dr1,dr2,besarg1,besarg2,be1,be2,si1,si2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt; |
---|
[ae3ce4e] | 2010 | double Pi; |
---|
| 2011 | int kk; |
---|
[8e91f01] | 2012 | |
---|
[ae3ce4e] | 2013 | Pi = 4.0*atan(1.0); |
---|
[8e91f01] | 2014 | |
---|
[ae3ce4e] | 2015 | dr1 = rhoc-rhosolv; |
---|
| 2016 | dr2 = rhol-rhosolv; |
---|
| 2017 | area = Pi*rcore*rcore; |
---|
| 2018 | totald=2.0*(thick+length); |
---|
[8e91f01] | 2019 | |
---|
[ae3ce4e] | 2020 | besarg1 = qq*rcore*sin(dum); |
---|
| 2021 | besarg2 = qq*rcore*sin(dum); |
---|
[8e91f01] | 2022 | |
---|
[ae3ce4e] | 2023 | sinarg1 = qq*length*cos(dum); |
---|
| 2024 | sinarg2 = qq*(length+thick)*cos(dum); |
---|
[8e91f01] | 2025 | |
---|
[975ec8e] | 2026 | if (besarg1 == 0){ |
---|
| 2027 | be1 = 0.5; |
---|
| 2028 | }else{ |
---|
| 2029 | be1 = NR_BessJ1(besarg1)/besarg1; |
---|
| 2030 | } |
---|
| 2031 | if (besarg2 == 0){ |
---|
| 2032 | be2 = 0.5; |
---|
| 2033 | }else{ |
---|
| 2034 | be2 = NR_BessJ1(besarg2)/besarg2; |
---|
| 2035 | } |
---|
| 2036 | if (sinarg1 == 0){ |
---|
| 2037 | si1 = 1; |
---|
| 2038 | }else{ |
---|
| 2039 | si1 = sin(sinarg1)/sinarg1; |
---|
| 2040 | } |
---|
| 2041 | if (besarg2 == 0){ |
---|
| 2042 | si2 = 1; |
---|
| 2043 | }else{ |
---|
| 2044 | si2 = sin(sinarg2)/sinarg2; |
---|
| 2045 | } |
---|
| 2046 | |
---|
| 2047 | t1 = 2*area*(2*length)*dr1*(si1)*(be1); |
---|
| 2048 | t2 = 2*area*dr2*(totald*si2-2*length*si1)*(be2); |
---|
[8e91f01] | 2049 | |
---|
[ae3ce4e] | 2050 | retval =((t1+t2)*(t1+t2))*sin(dum); |
---|
[8e91f01] | 2051 | |
---|
[ae3ce4e] | 2052 | // loop for the structure facture S(q) |
---|
| 2053 | sqq=0.0; |
---|
| 2054 | for(kk=1;kk<N;kk+=1) { |
---|
| 2055 | dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0; |
---|
| 2056 | sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt); |
---|
[8e91f01] | 2057 | } |
---|
| 2058 | |
---|
[ae3ce4e] | 2059 | // end of loop for S(q) |
---|
| 2060 | sqq=1.0+2.0*sqq/N; |
---|
[975ec8e] | 2061 | |
---|
[ae3ce4e] | 2062 | retval *= sqq; |
---|
[8e91f01] | 2063 | |
---|
[ae3ce4e] | 2064 | return(retval); |
---|
| 2065 | } |
---|
| 2066 | |
---|
| 2067 | |
---|
| 2068 | double |
---|
| 2069 | Cyl_PolyRadKernel(double q, double radius, double length, double zz, double delrho, double dumRad) |
---|
| 2070 | { |
---|
[8e91f01] | 2071 | |
---|
[ae3ce4e] | 2072 | double halfheight,uplim,lolim,zi,summ,yyy,Pi; |
---|
| 2073 | double answer,dr,Vcyl; |
---|
| 2074 | int i,nord; |
---|
[8e91f01] | 2075 | |
---|
[ae3ce4e] | 2076 | Pi = 4.0*atan(1.0); |
---|
| 2077 | lolim = 0; |
---|
| 2078 | uplim = Pi/2.0; |
---|
| 2079 | halfheight = length/2.0; |
---|
| 2080 | // nord=20; |
---|
| 2081 | nord=76; |
---|
| 2082 | summ = 0.0; |
---|
[8e91f01] | 2083 | |
---|
[ae3ce4e] | 2084 | //do the cylinder orientational average |
---|
| 2085 | // for(i=0;i<nord;i++) { |
---|
| 2086 | // zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 2087 | // yyy = Gauss20Wt[i] * CylKernel(q, dumRad, halfheight, zi); |
---|
| 2088 | // summ += yyy; |
---|
| 2089 | // } |
---|
| 2090 | for(i=0;i<nord;i++) { |
---|
| 2091 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
| 2092 | yyy = Gauss76Wt[i] * CylKernel(q, dumRad, halfheight, zi); |
---|
| 2093 | summ += yyy; |
---|
| 2094 | } |
---|
| 2095 | answer = (uplim-lolim)/2.0*summ; |
---|
| 2096 | // Multiply by contrast^2 |
---|
| 2097 | answer *= delrho*delrho; |
---|
| 2098 | // don't do the normal scaling to volume here |
---|
| 2099 | // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder |
---|
| 2100 | Vcyl = Pi*dumRad*dumRad*length; |
---|
| 2101 | answer *= Vcyl*Vcyl; |
---|
[8e91f01] | 2102 | |
---|
[ae3ce4e] | 2103 | dr = SchulzPoint_cpr(dumRad,radius,zz); |
---|
| 2104 | return(dr*answer); |
---|
| 2105 | } |
---|
| 2106 | |
---|
| 2107 | double |
---|
| 2108 | SchulzPoint_cpr(double dumRad, double radius, double zz) |
---|
| 2109 | { |
---|
| 2110 | double dr; |
---|
[8e91f01] | 2111 | |
---|
[ae3ce4e] | 2112 | dr = zz*log(dumRad) - gammaln(zz+1.0) + (zz+1.0)*log((zz+1.0)/radius)-(dumRad/radius*(zz+1.0)); |
---|
| 2113 | return(exp(dr)); |
---|
| 2114 | } |
---|
| 2115 | |
---|
| 2116 | static double |
---|
| 2117 | gammaln(double xx) |
---|
| 2118 | { |
---|
| 2119 | double x,y,tmp,ser; |
---|
| 2120 | static double cof[6]={76.18009172947146,-86.50532032941677, |
---|
| 2121 | 24.01409824083091,-1.231739572450155, |
---|
| 2122 | 0.1208650973866179e-2,-0.5395239384953e-5}; |
---|
| 2123 | int j; |
---|
[8e91f01] | 2124 | |
---|
[ae3ce4e] | 2125 | y=x=xx; |
---|
| 2126 | tmp=x+5.5; |
---|
| 2127 | tmp -= (x+0.5)*log(tmp); |
---|
| 2128 | ser=1.000000000190015; |
---|
| 2129 | for (j=0;j<=5;j++) ser += cof[j]/++y; |
---|
| 2130 | return -tmp+log(2.5066282746310005*ser/x); |
---|
| 2131 | } |
---|
| 2132 | |
---|
| 2133 | |
---|
| 2134 | double |
---|
| 2135 | EllipsoidKernel(double qq, double a, double nua, double dum) |
---|
| 2136 | { |
---|
| 2137 | double arg,nu,retval; //local variables |
---|
[8e91f01] | 2138 | |
---|
[ae3ce4e] | 2139 | nu = nua/a; |
---|
| 2140 | arg = qq*a*sqrt(1+dum*dum*(nu*nu-1)); |
---|
[975ec8e] | 2141 | if (arg == 0){ |
---|
| 2142 | retval =1/3; |
---|
| 2143 | }else{ |
---|
| 2144 | retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg); |
---|
| 2145 | } |
---|
[ae3ce4e] | 2146 | retval *= retval; |
---|
| 2147 | retval *= 9; |
---|
[8e91f01] | 2148 | |
---|
[ae3ce4e] | 2149 | return(retval); |
---|
| 2150 | }//Function EllipsoidKernel() |
---|
| 2151 | |
---|
| 2152 | double |
---|
| 2153 | HolCylKernel(double qq, double rcore, double rshell, double length, double dum) |
---|
| 2154 | { |
---|
| 2155 | double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval; //local variables |
---|
[8e91f01] | 2156 | |
---|
[ae3ce4e] | 2157 | gamma = rcore/rshell; |
---|
| 2158 | arg1 = qq*rshell*sqrt(1-dum*dum); //1=shell (outer radius) |
---|
| 2159 | arg2 = qq*rcore*sqrt(1-dum*dum); //2=core (inner radius) |
---|
[975ec8e] | 2160 | if (arg1 == 0){ |
---|
| 2161 | lam1 = 1; |
---|
| 2162 | }else{ |
---|
| 2163 | lam1 = 2*NR_BessJ1(arg1)/arg1; |
---|
| 2164 | } |
---|
| 2165 | if (arg2 == 0){ |
---|
| 2166 | lam2 = 1; |
---|
| 2167 | }else{ |
---|
| 2168 | lam2 = 2*NR_BessJ1(arg2)/arg2; |
---|
| 2169 | } |
---|
| 2170 | //Todo: Need to check psi behavior as gamma goes to 1. |
---|
[ae3ce4e] | 2171 | psi = 1/(1-gamma*gamma)*(lam1 - gamma*gamma*lam2); //SRK 10/19/00 |
---|
| 2172 | sinarg = qq*length*dum/2; |
---|
[975ec8e] | 2173 | if (sinarg == 0){ |
---|
| 2174 | t2 = 1; |
---|
| 2175 | }else{ |
---|
| 2176 | t2 = sin(sinarg)/sinarg; |
---|
| 2177 | } |
---|
[8e91f01] | 2178 | |
---|
[ae3ce4e] | 2179 | retval = psi*psi*t2*t2; |
---|
[8e91f01] | 2180 | |
---|
[ae3ce4e] | 2181 | return(retval); |
---|
| 2182 | }//Function HolCylKernel() |
---|
| 2183 | |
---|
| 2184 | double |
---|
| 2185 | PPKernel(double aa, double mu, double uu) |
---|
| 2186 | { |
---|
| 2187 | // mu passed in is really mu*sqrt(1-sig^2) |
---|
| 2188 | double arg1,arg2,Pi,tmp1,tmp2; //local variables |
---|
[8e91f01] | 2189 | |
---|
[ae3ce4e] | 2190 | Pi = 4.0*atan(1.0); |
---|
[8e91f01] | 2191 | |
---|
[ae3ce4e] | 2192 | //handle arg=0 separately, as sin(t)/t -> 1 as t->0 |
---|
| 2193 | arg1 = (mu/2)*cos(Pi*uu/2); |
---|
| 2194 | arg2 = (mu*aa/2)*sin(Pi*uu/2); |
---|
| 2195 | if(arg1==0) { |
---|
| 2196 | tmp1 = 1; |
---|
| 2197 | } else { |
---|
| 2198 | tmp1 = sin(arg1)*sin(arg1)/arg1/arg1; |
---|
| 2199 | } |
---|
[8e91f01] | 2200 | |
---|
[ae3ce4e] | 2201 | if (arg2==0) { |
---|
| 2202 | tmp2 = 1; |
---|
| 2203 | } else { |
---|
| 2204 | tmp2 = sin(arg2)*sin(arg2)/arg2/arg2; |
---|
| 2205 | } |
---|
[8e91f01] | 2206 | |
---|
[ae3ce4e] | 2207 | return (tmp1*tmp2); |
---|
[8e91f01] | 2208 | |
---|
[ae3ce4e] | 2209 | }//Function PPKernel() |
---|
| 2210 | |
---|
| 2211 | |
---|
| 2212 | double |
---|
| 2213 | TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy) |
---|
| 2214 | { |
---|
[8e91f01] | 2215 | |
---|
[ae3ce4e] | 2216 | double arg,val,pi; //local variables |
---|
[8e91f01] | 2217 | |
---|
[ae3ce4e] | 2218 | pi = 4.0*atan(1.0); |
---|
[8e91f01] | 2219 | |
---|
[ae3ce4e] | 2220 | arg = aa*aa*cos(pi*dx/2)*cos(pi*dx/2); |
---|
| 2221 | arg += bb*bb*sin(pi*dx/2)*sin(pi*dx/2)*(1-dy*dy); |
---|
| 2222 | arg += cc*cc*dy*dy; |
---|
| 2223 | arg = q*sqrt(arg); |
---|
[975ec8e] | 2224 | if (arg == 0){ |
---|
| 2225 | val = 1; // as arg --> 0, val should go to 1.0 |
---|
| 2226 | }else{ |
---|
| 2227 | val = 9 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ); |
---|
| 2228 | } |
---|
[ae3ce4e] | 2229 | return (val); |
---|
[8e91f01] | 2230 | |
---|
[ae3ce4e] | 2231 | }//Function TriaxialKernel() |
---|
| 2232 | |
---|
| 2233 | |
---|
| 2234 | double |
---|
| 2235 | CylKernel(double qq, double rr,double h, double theta) |
---|
| 2236 | { |
---|
[8e91f01] | 2237 | |
---|
[ae3ce4e] | 2238 | // qq is the q-value for the calculation (1/A) |
---|
| 2239 | // rr is the radius of the cylinder (A) |
---|
| 2240 | // h is the HALF-LENGTH of the cylinder = L/2 (A) |
---|
[8e91f01] | 2241 | |
---|
[975ec8e] | 2242 | double besarg,bj,retval,d1,t1,b1,t2,b2,siarg,be,si; //Local variables |
---|
[8e91f01] | 2243 | |
---|
| 2244 | |
---|
[ae3ce4e] | 2245 | besarg = qq*rr*sin(theta); |
---|
[975ec8e] | 2246 | siarg = qq * h * cos(theta); |
---|
[ae3ce4e] | 2247 | bj =NR_BessJ1(besarg); |
---|
[8e91f01] | 2248 | |
---|
[ae3ce4e] | 2249 | //* Computing 2nd power */ |
---|
[975ec8e] | 2250 | d1 = sin(siarg); |
---|
[ae3ce4e] | 2251 | t1 = d1 * d1; |
---|
| 2252 | //* Computing 2nd power */ |
---|
| 2253 | d1 = bj; |
---|
| 2254 | t2 = d1 * d1 * 4.0 * sin(theta); |
---|
| 2255 | //* Computing 2nd power */ |
---|
[975ec8e] | 2256 | d1 = siarg; |
---|
[ae3ce4e] | 2257 | b1 = d1 * d1; |
---|
| 2258 | //* Computing 2nd power */ |
---|
| 2259 | d1 = qq * rr * sin(theta); |
---|
| 2260 | b2 = d1 * d1; |
---|
[975ec8e] | 2261 | if (besarg == 0){ |
---|
| 2262 | be = sin(theta); |
---|
| 2263 | }else{ |
---|
| 2264 | be = t2 / b2; |
---|
| 2265 | } |
---|
| 2266 | if (siarg == 0){ |
---|
| 2267 | si = 1; |
---|
| 2268 | }else{ |
---|
| 2269 | si = t1 / b1; |
---|
| 2270 | } |
---|
| 2271 | retval = be * si; |
---|
[8e91f01] | 2272 | |
---|
[ae3ce4e] | 2273 | return (retval); |
---|
[8e91f01] | 2274 | |
---|
[ae3ce4e] | 2275 | }//Function CylKernel() |
---|
| 2276 | |
---|
| 2277 | double |
---|
| 2278 | EllipCylKernel(double qq, double ra,double nu, double theta) |
---|
| 2279 | { |
---|
[8e91f01] | 2280 | //this is the function LAMBDA1^2 in Feigin's notation |
---|
[ae3ce4e] | 2281 | // qq is the q-value for the calculation (1/A) |
---|
| 2282 | // ra is the transformed radius"a" in Feigin's notation |
---|
| 2283 | // nu is the ratio (major radius)/(minor radius) of the Ellipsoid [=] --- |
---|
| 2284 | // theta is the dummy variable of the integration |
---|
[8e91f01] | 2285 | |
---|
| 2286 | double retval,arg; //Local variables |
---|
| 2287 | |
---|
[ae3ce4e] | 2288 | arg = qq*ra*sqrt((1+nu*nu)/2+(1-nu*nu)*cos(theta)/2); |
---|
[975ec8e] | 2289 | if (arg == 0){ |
---|
| 2290 | retval = 1; |
---|
| 2291 | }else{ |
---|
| 2292 | retval = 2*NR_BessJ1(arg)/arg; |
---|
| 2293 | } |
---|
[8e91f01] | 2294 | |
---|
[ae3ce4e] | 2295 | //square it |
---|
| 2296 | retval *= retval; |
---|
[8e91f01] | 2297 | |
---|
[ae3ce4e] | 2298 | return(retval); |
---|
[8e91f01] | 2299 | |
---|
[ae3ce4e] | 2300 | }//Function EllipCylKernel() |
---|
| 2301 | |
---|
| 2302 | double NR_BessJ1(double x) |
---|
| 2303 | { |
---|
| 2304 | double ax,z; |
---|
| 2305 | double xx,y,ans,ans1,ans2; |
---|
[8e91f01] | 2306 | |
---|
[ae3ce4e] | 2307 | if ((ax=fabs(x)) < 8.0) { |
---|
| 2308 | y=x*x; |
---|
| 2309 | ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1 |
---|
| 2310 | +y*(-2972611.439+y*(15704.48260+y*(-30.16036606)))))); |
---|
| 2311 | ans2=144725228442.0+y*(2300535178.0+y*(18583304.74 |
---|
| 2312 | +y*(99447.43394+y*(376.9991397+y*1.0)))); |
---|
| 2313 | ans=ans1/ans2; |
---|
| 2314 | } else { |
---|
| 2315 | z=8.0/ax; |
---|
| 2316 | y=z*z; |
---|
| 2317 | xx=ax-2.356194491; |
---|
| 2318 | ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4 |
---|
| 2319 | +y*(0.2457520174e-5+y*(-0.240337019e-6)))); |
---|
| 2320 | ans2=0.04687499995+y*(-0.2002690873e-3 |
---|
| 2321 | +y*(0.8449199096e-5+y*(-0.88228987e-6 |
---|
| 2322 | +y*0.105787412e-6))); |
---|
| 2323 | ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2); |
---|
| 2324 | if (x < 0.0) ans = -ans; |
---|
| 2325 | } |
---|
[8e91f01] | 2326 | |
---|
[ae3ce4e] | 2327 | return(ans); |
---|
| 2328 | } |
---|