[230f479] | 1 | /* TwoPhaseFit.c |
---|
| 2 | |
---|
| 3 | */ |
---|
| 4 | |
---|
| 5 | #include "StandardHeaders.h" // Include ANSI headers, Mac headers |
---|
| 6 | #include "libTwoPhase.h" |
---|
| 7 | |
---|
| 8 | /* internal functions */ |
---|
| 9 | static double |
---|
| 10 | gammln(double xx) { |
---|
| 11 | |
---|
| 12 | double x,y,tmp,ser; |
---|
| 13 | static double cof[6]={76.18009172947146,-86.50532032941677, |
---|
| 14 | 24.01409824083091,-1.231739572450155, |
---|
| 15 | 0.1208650973866179e-2,-0.5395239384953e-5}; |
---|
| 16 | int j; |
---|
| 17 | |
---|
| 18 | y=x=xx; |
---|
| 19 | tmp=x+5.5; |
---|
| 20 | tmp -= (x+0.5)*log(tmp); |
---|
| 21 | ser=1.000000000190015; |
---|
| 22 | for (j=0;j<=5;j++) ser += cof[j]/++y; |
---|
| 23 | return -tmp+log(2.5066282746310005*ser/x); |
---|
| 24 | } |
---|
| 25 | |
---|
| 26 | // scattering from the Teubner-Strey model for microemulsions - hardly needs to be an XOP... |
---|
| 27 | double |
---|
| 28 | TeubnerStreyModel(double dp[], double q) |
---|
| 29 | { |
---|
| 30 | double inten,q2,q4; //my local names |
---|
| 31 | |
---|
| 32 | q2 = q*q; |
---|
| 33 | q4 = q2*q2; |
---|
| 34 | |
---|
| 35 | inten = 1.0/(dp[0]+dp[1]*q2+dp[2]*q4); |
---|
| 36 | inten += dp[3]; |
---|
| 37 | return(inten); |
---|
| 38 | } |
---|
| 39 | |
---|
| 40 | double |
---|
| 41 | Power_Law_Model(double dp[], double q) |
---|
| 42 | { |
---|
| 43 | double qval; |
---|
| 44 | double inten,A,m,bgd; //my local names |
---|
| 45 | |
---|
| 46 | qval= q; |
---|
| 47 | |
---|
| 48 | A = dp[0]; |
---|
| 49 | m = dp[1]; |
---|
| 50 | bgd = dp[2]; |
---|
| 51 | inten = A*pow(qval,-m) + bgd; |
---|
| 52 | return(inten); |
---|
| 53 | } |
---|
| 54 | |
---|
| 55 | |
---|
| 56 | double |
---|
| 57 | Peak_Lorentz_Model(double dp[], double q) |
---|
| 58 | { |
---|
| 59 | double qval; |
---|
| 60 | double inten,I0, qpk, dq,bgd; //my local names |
---|
| 61 | qval= q; |
---|
| 62 | |
---|
| 63 | I0 = dp[0]; |
---|
| 64 | qpk = dp[1]; |
---|
| 65 | dq = dp[2]; |
---|
| 66 | bgd = dp[3]; |
---|
| 67 | inten = I0/(1.0 + pow( (qval-qpk)/dq,2) ) + bgd; |
---|
| 68 | |
---|
| 69 | return(inten); |
---|
| 70 | } |
---|
| 71 | |
---|
| 72 | double |
---|
| 73 | Peak_Gauss_Model(double dp[], double q) |
---|
| 74 | { |
---|
| 75 | double qval; |
---|
| 76 | double inten,I0, qpk, dq,bgd; //my local names |
---|
| 77 | |
---|
| 78 | qval= q; |
---|
| 79 | |
---|
| 80 | I0 = dp[0]; |
---|
| 81 | qpk = dp[1]; |
---|
| 82 | dq = dp[2]; |
---|
| 83 | bgd = dp[3]; |
---|
| 84 | inten = I0*exp(-0.5*pow((qval-qpk)/dq,2))+ bgd; |
---|
| 85 | |
---|
| 86 | return(inten); |
---|
| 87 | } |
---|
| 88 | |
---|
| 89 | double |
---|
| 90 | Lorentz_Model(double dp[], double q) |
---|
| 91 | { |
---|
| 92 | double qval; |
---|
| 93 | double inten,I0, L,bgd; //my local names |
---|
| 94 | |
---|
| 95 | qval= q; |
---|
| 96 | |
---|
| 97 | I0 = dp[0]; |
---|
| 98 | L = dp[1]; |
---|
| 99 | bgd = dp[2]; |
---|
| 100 | inten = I0/(1.0 + (qval*L)*(qval*L)) + bgd; |
---|
| 101 | |
---|
| 102 | return(inten); |
---|
| 103 | } |
---|
| 104 | |
---|
| 105 | double |
---|
| 106 | Fractal(double dp[], double q) |
---|
| 107 | { |
---|
| 108 | double x,pi; |
---|
| 109 | double r0,Df,corr,phi,sldp,sldm,bkg; |
---|
| 110 | double pq,sq,ans; |
---|
| 111 | |
---|
| 112 | pi = 4.0*atan(1.0); |
---|
| 113 | x=q; |
---|
| 114 | |
---|
| 115 | phi = dp[0]; // volume fraction of building block spheres... |
---|
| 116 | r0 = dp[1]; // radius of building block |
---|
| 117 | Df = dp[2]; // fractal dimension |
---|
| 118 | corr = dp[3]; // correlation length of fractal-like aggregates |
---|
| 119 | sldp = dp[4]; // SLD of building block |
---|
| 120 | sldm = dp[5]; // SLD of matrix or solution |
---|
| 121 | bkg = dp[6]; // flat background |
---|
| 122 | |
---|
| 123 | //calculate P(q) for the spherical subunits, units cm-1 sr-1 |
---|
| 124 | pq = 1.0e8*phi*4.0/3.0*pi*r0*r0*r0*(sldp-sldm)*(sldp-sldm)*pow((3*(sin(x*r0) - x*r0*cos(x*r0))/pow((x*r0),3)),2); |
---|
| 125 | |
---|
| 126 | //calculate S(q) |
---|
| 127 | sq = Df*exp(gammln(Df-1.0))*sin((Df-1.0)*atan(x*corr)); |
---|
| 128 | sq /= pow((x*r0),Df) * pow((1.0 + 1.0/(x*corr)/(x*corr)),((Df-1.0)/2.0)); |
---|
| 129 | sq += 1.0; |
---|
| 130 | //combine and return |
---|
| 131 | ans = pq*sq + bkg; |
---|
| 132 | |
---|
| 133 | return(ans); |
---|
| 134 | } |
---|
| 135 | |
---|
| 136 | // 6 JUL 2009 SRK changed definition of Izero scale factor to be uncorrelated with range |
---|
| 137 | // |
---|
| 138 | double |
---|
| 139 | DAB_Model(double dp[], double q) |
---|
| 140 | { |
---|
| 141 | double qval,inten; |
---|
| 142 | double Izero, range, incoh; |
---|
| 143 | |
---|
| 144 | qval= q; |
---|
| 145 | Izero = dp[0]; |
---|
| 146 | range = dp[1]; |
---|
| 147 | incoh = dp[2]; |
---|
| 148 | |
---|
| 149 | inten = (Izero*range*range*range)/pow((1.0 + (qval*range)*(qval*range)),2) + incoh; |
---|
| 150 | |
---|
| 151 | return(inten); |
---|
| 152 | } |
---|
| 153 | |
---|
| 154 | // G. Beaucage's Unified Model (1-4 levels) |
---|
| 155 | // |
---|
| 156 | double |
---|
| 157 | OneLevel(double dp[], double q) |
---|
| 158 | { |
---|
| 159 | double x,ans,erf1; |
---|
| 160 | double G1,Rg1,B1,Pow1,bkg,scale; |
---|
| 161 | |
---|
| 162 | x=q; |
---|
| 163 | scale = dp[0]; |
---|
| 164 | G1 = dp[1]; |
---|
| 165 | Rg1 = dp[2]; |
---|
| 166 | B1 = dp[3]; |
---|
| 167 | Pow1 = dp[4]; |
---|
| 168 | bkg = dp[5]; |
---|
| 169 | |
---|
| 170 | erf1 = erf( (x*Rg1/sqrt(6.0))); |
---|
| 171 | |
---|
| 172 | ans = G1*exp(-x*x*Rg1*Rg1/3.0); |
---|
| 173 | ans += B1*pow((erf1*erf1*erf1/x),Pow1); |
---|
| 174 | |
---|
| 175 | if(x == 0) { |
---|
| 176 | ans = G1; |
---|
| 177 | } |
---|
| 178 | |
---|
| 179 | ans *= scale; |
---|
| 180 | ans += bkg; |
---|
| 181 | return(ans); |
---|
| 182 | } |
---|
| 183 | |
---|
| 184 | // G. Beaucage's Unified Model (1-4 levels) |
---|
| 185 | // |
---|
| 186 | double |
---|
| 187 | TwoLevel(double dp[], double q) |
---|
| 188 | { |
---|
| 189 | double x; |
---|
| 190 | double ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg; |
---|
| 191 | double erf1,erf2,scale; |
---|
| 192 | |
---|
| 193 | x=q; |
---|
| 194 | |
---|
| 195 | scale = dp[0]; |
---|
| 196 | G1 = dp[1]; //equivalent to I(0) |
---|
| 197 | Rg1 = dp[2]; |
---|
| 198 | B1 = dp[3]; |
---|
| 199 | Pow1 = dp[4]; |
---|
| 200 | G2 = dp[5]; |
---|
| 201 | Rg2 = dp[6]; |
---|
| 202 | B2 = dp[7]; |
---|
| 203 | Pow2 = dp[8]; |
---|
| 204 | bkg = dp[9]; |
---|
| 205 | |
---|
| 206 | erf1 = erf( (x*Rg1/sqrt(6.0)) ); |
---|
| 207 | erf2 = erf( (x*Rg2/sqrt(6.0)) ); |
---|
| 208 | //Print erf1 |
---|
| 209 | |
---|
| 210 | ans = G1*exp(-x*x*Rg1*Rg1/3.0); |
---|
| 211 | ans += B1*exp(-x*x*Rg2*Rg2/3.0)*pow((erf1*erf1*erf1/x),Pow1); |
---|
| 212 | ans += G2*exp(-x*x*Rg2*Rg2/3.0); |
---|
| 213 | ans += B2*pow((erf2*erf2*erf2/x),Pow2); |
---|
| 214 | |
---|
| 215 | if(x == 0) { |
---|
| 216 | ans = G1+G2; |
---|
| 217 | } |
---|
| 218 | |
---|
| 219 | ans *= scale; |
---|
| 220 | ans += bkg; |
---|
| 221 | |
---|
| 222 | return(ans); |
---|
| 223 | } |
---|
| 224 | |
---|
| 225 | // G. Beaucage's Unified Model (1-4 levels) |
---|
| 226 | // |
---|
| 227 | double |
---|
| 228 | ThreeLevel(double dp[], double q) |
---|
| 229 | { |
---|
| 230 | double x; |
---|
| 231 | double ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg; |
---|
| 232 | double G3,Rg3,B3,Pow3,erf3; |
---|
| 233 | double erf1,erf2,scale; |
---|
| 234 | |
---|
| 235 | x=q; |
---|
| 236 | |
---|
| 237 | scale = dp[0]; |
---|
| 238 | G1 = dp[1]; //equivalent to I(0) |
---|
| 239 | Rg1 = dp[2]; |
---|
| 240 | B1 = dp[3]; |
---|
| 241 | Pow1 = dp[4]; |
---|
| 242 | G2 = dp[5]; |
---|
| 243 | Rg2 = dp[6]; |
---|
| 244 | B2 = dp[7]; |
---|
| 245 | Pow2 = dp[8]; |
---|
| 246 | G3 = dp[9]; |
---|
| 247 | Rg3 = dp[10]; |
---|
| 248 | B3 = dp[11]; |
---|
| 249 | Pow3 = dp[12]; |
---|
| 250 | bkg = dp[13]; |
---|
| 251 | |
---|
| 252 | erf1 = erf( (x*Rg1/sqrt(6.0)) ); |
---|
| 253 | erf2 = erf( (x*Rg2/sqrt(6.0)) ); |
---|
| 254 | erf3 = erf( (x*Rg3/sqrt(6.0)) ); |
---|
| 255 | //Print erf1 |
---|
| 256 | |
---|
| 257 | ans = G1*exp(-x*x*Rg1*Rg1/3.0) + B1*exp(-x*x*Rg2*Rg2/3.0)*pow((erf1*erf1*erf1/x),Pow1); |
---|
| 258 | ans += G2*exp(-x*x*Rg2*Rg2/3.0) + B2*exp(-x*x*Rg3*Rg3/3.0)*pow((erf2*erf2*erf2/x),Pow2); |
---|
| 259 | ans += G3*exp(-x*x*Rg3*Rg3/3.0) + B3*pow((erf3*erf3*erf3/x),Pow3); |
---|
| 260 | |
---|
| 261 | if(x == 0) { |
---|
| 262 | ans = G1+G2+G3; |
---|
| 263 | } |
---|
| 264 | |
---|
| 265 | ans *= scale; |
---|
| 266 | ans += bkg; |
---|
| 267 | |
---|
| 268 | return(ans); |
---|
| 269 | } |
---|
| 270 | |
---|
| 271 | // G. Beaucage's Unified Model (1-4 levels) |
---|
| 272 | // |
---|
| 273 | double |
---|
| 274 | FourLevel(double dp[], double q) |
---|
| 275 | { |
---|
| 276 | double x; |
---|
| 277 | double ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg; |
---|
| 278 | double G3,Rg3,B3,Pow3,erf3; |
---|
| 279 | double G4,Rg4,B4,Pow4,erf4; |
---|
| 280 | double erf1,erf2,scale; |
---|
| 281 | |
---|
| 282 | x=q; |
---|
| 283 | |
---|
| 284 | scale = dp[0]; |
---|
| 285 | G1 = dp[1]; //equivalent to I(0) |
---|
| 286 | Rg1 = dp[2]; |
---|
| 287 | B1 = dp[3]; |
---|
| 288 | Pow1 = dp[4]; |
---|
| 289 | G2 = dp[5]; |
---|
| 290 | Rg2 = dp[6]; |
---|
| 291 | B2 = dp[7]; |
---|
| 292 | Pow2 = dp[8]; |
---|
| 293 | G3 = dp[9]; |
---|
| 294 | Rg3 = dp[10]; |
---|
| 295 | B3 = dp[11]; |
---|
| 296 | Pow3 = dp[12]; |
---|
| 297 | G4 = dp[13]; |
---|
| 298 | Rg4 = dp[14]; |
---|
| 299 | B4 = dp[15]; |
---|
| 300 | Pow4 = dp[16]; |
---|
| 301 | bkg = dp[17]; |
---|
| 302 | |
---|
| 303 | erf1 = erf( (x*Rg1/sqrt(6.0)) ); |
---|
| 304 | erf2 = erf( (x*Rg2/sqrt(6.0)) ); |
---|
| 305 | erf3 = erf( (x*Rg3/sqrt(6.0)) ); |
---|
| 306 | erf4 = erf( (x*Rg4/sqrt(6.0)) ); |
---|
| 307 | |
---|
| 308 | ans = G1*exp(-x*x*Rg1*Rg1/3.0) + B1*exp(-x*x*Rg2*Rg2/3.0)*pow((erf1*erf1*erf1/x),Pow1); |
---|
| 309 | ans += G2*exp(-x*x*Rg2*Rg2/3.0) + B2*exp(-x*x*Rg3*Rg3/3.0)*pow((erf2*erf2*erf2/x),Pow2); |
---|
| 310 | ans += G3*exp(-x*x*Rg3*Rg3/3.0) + B3*exp(-x*x*Rg4*Rg4/3.0)*pow((erf3*erf3*erf3/x),Pow3); |
---|
| 311 | ans += G4*exp(-x*x*Rg4*Rg4/3.0) + B4*pow((erf4*erf4*erf4/x),Pow4); |
---|
| 312 | |
---|
| 313 | if(x == 0) { |
---|
| 314 | ans = G1+G2+G3+G4; |
---|
| 315 | } |
---|
| 316 | |
---|
| 317 | ans *= scale; |
---|
| 318 | ans += bkg; |
---|
| 319 | |
---|
| 320 | return(ans); |
---|
| 321 | } |
---|
| 322 | |
---|
| 323 | double |
---|
| 324 | BroadPeak(double dp[], double q) |
---|
| 325 | { |
---|
| 326 | // variables are: |
---|
| 327 | //[0] Porod term scaling |
---|
| 328 | //[1] Porod exponent |
---|
| 329 | //[2] Lorentzian term scaling |
---|
| 330 | //[3] Lorentzian screening length [A] |
---|
| 331 | //[4] peak location [1/A] |
---|
| 332 | //[5] Lorentzian exponent |
---|
| 333 | //[6] background |
---|
| 334 | |
---|
| 335 | double aa,nn,cc,LL,Qzero,mm,bgd,inten,qval; |
---|
| 336 | qval= q; |
---|
| 337 | aa = dp[0]; |
---|
| 338 | nn = dp[1]; |
---|
| 339 | cc = dp[2]; |
---|
| 340 | LL = dp[3]; |
---|
| 341 | Qzero = dp[4]; |
---|
| 342 | mm = dp[5]; |
---|
| 343 | bgd = dp[6]; |
---|
| 344 | |
---|
| 345 | inten = aa/pow(qval,nn); |
---|
| 346 | inten += cc/(1.0 + pow((fabs(qval-Qzero)*LL),mm) ); |
---|
| 347 | inten += bgd; |
---|
| 348 | |
---|
| 349 | return(inten); |
---|
| 350 | } |
---|
| 351 | |
---|
| 352 | double |
---|
| 353 | CorrLength(double dp[], double q) |
---|
| 354 | { |
---|
| 355 | // variables are: |
---|
| 356 | //[0] Porod term scaling |
---|
| 357 | //[1] Porod exponent |
---|
| 358 | //[2] Lorentzian term scaling |
---|
| 359 | //[3] Lorentzian screening length [A] |
---|
| 360 | //[4] Lorentzian exponent |
---|
| 361 | //[5] background |
---|
| 362 | |
---|
| 363 | double aa,nn,cc,LL,mm,bgd,inten,qval; |
---|
| 364 | qval= q; |
---|
| 365 | aa = dp[0]; |
---|
| 366 | nn = dp[1]; |
---|
| 367 | cc = dp[2]; |
---|
| 368 | LL = dp[3]; |
---|
| 369 | mm = dp[4]; |
---|
| 370 | bgd = dp[5]; |
---|
| 371 | |
---|
| 372 | inten = aa/pow(qval,nn); |
---|
| 373 | inten += cc/(1.0 + pow((qval*LL),mm) ); |
---|
| 374 | inten += bgd; |
---|
| 375 | |
---|
| 376 | return(inten); |
---|
| 377 | } |
---|
| 378 | |
---|
| 379 | double |
---|
| 380 | TwoLorentzian(double dp[], double q) |
---|
| 381 | { |
---|
| 382 | // variables are: |
---|
| 383 | //[0] Lorentzian term scaling |
---|
| 384 | //[1] Lorentzian screening length [A] |
---|
| 385 | //[2] Lorentzian exponent |
---|
| 386 | //[3] Lorentzian #2 term scaling |
---|
| 387 | //[4] Lorentzian #2 screening length [A] |
---|
| 388 | //[5] Lorentzian #2 exponent |
---|
| 389 | //[6] background |
---|
| 390 | |
---|
| 391 | double aa,LL1,nn,cc,LL2,mm,bgd,inten,qval; |
---|
| 392 | qval= q; |
---|
| 393 | aa = dp[0]; |
---|
| 394 | LL1 = dp[1]; |
---|
| 395 | nn = dp[2]; |
---|
| 396 | cc = dp[3]; |
---|
| 397 | LL2 = dp[4]; |
---|
| 398 | mm = dp[5]; |
---|
| 399 | bgd= dp[6]; |
---|
| 400 | |
---|
| 401 | inten = aa/(1.0 + pow((qval*LL1),nn) ); |
---|
| 402 | inten += cc/(1.0 + pow((qval*LL2),mm) ); |
---|
| 403 | inten += bgd; |
---|
| 404 | |
---|
| 405 | return(inten); |
---|
| 406 | } |
---|
| 407 | |
---|
| 408 | double |
---|
| 409 | TwoPowerLaw(double dp[], double q) |
---|
| 410 | { |
---|
| 411 | //[0] Coefficient |
---|
| 412 | //[1] (-) Power @ low Q |
---|
| 413 | //[2] (-) Power @ high Q |
---|
| 414 | //[3] crossover Q-value |
---|
| 415 | //[4] incoherent background |
---|
| 416 | |
---|
| 417 | double A, m1,m2,qc,bgd,scale,inten,qval; |
---|
| 418 | qval= q; |
---|
| 419 | A = dp[0]; |
---|
| 420 | m1 = dp[1]; |
---|
| 421 | m2 = dp[2]; |
---|
| 422 | qc = dp[3]; |
---|
| 423 | bgd = dp[4]; |
---|
| 424 | |
---|
| 425 | if(qval<=qc){ |
---|
| 426 | inten = A*pow(qval,-1.0*m1); |
---|
| 427 | } else { |
---|
| 428 | scale = A*pow(qc,-1.0*m1) / pow(qc,-1.0*m2); |
---|
| 429 | inten = scale*pow(qval,-1.0*m2); |
---|
| 430 | } |
---|
| 431 | |
---|
| 432 | inten += bgd; |
---|
| 433 | |
---|
| 434 | return(inten); |
---|
| 435 | } |
---|
| 436 | |
---|
| 437 | double |
---|
| 438 | PolyGaussCoil(double dp[], double x) |
---|
| 439 | { |
---|
| 440 | //w[0] = scale |
---|
| 441 | //w[1] = radius of gyration [ᅵ] |
---|
| 442 | //w[2] = polydispersity, ratio of Mw/Mn |
---|
| 443 | //w[3] = bkg [cm-1] |
---|
| 444 | |
---|
| 445 | double scale,bkg,Rg,uval,Mw_Mn,inten,xi; |
---|
| 446 | |
---|
| 447 | scale = dp[0]; |
---|
| 448 | Rg = dp[1]; |
---|
| 449 | Mw_Mn = dp[2]; |
---|
| 450 | bkg = dp[3]; |
---|
| 451 | |
---|
| 452 | uval = Mw_Mn - 1.0; |
---|
| 453 | if(uval == 0.0) { |
---|
| 454 | uval = 1e-6; //avoid divide by zero error |
---|
| 455 | } |
---|
| 456 | |
---|
| 457 | xi = Rg*Rg*x*x/(1.0+2.0*uval); |
---|
| 458 | |
---|
| 459 | if(xi < 1e-3) { |
---|
| 460 | return(scale+bkg); //limiting value |
---|
| 461 | } |
---|
| 462 | |
---|
| 463 | inten = 2.0*(pow((1.0+uval*xi),(-1.0/uval))+xi-1.0); |
---|
| 464 | inten /= (1.0+uval)*xi*xi; |
---|
| 465 | |
---|
| 466 | inten *= scale; |
---|
| 467 | //add in the background |
---|
| 468 | inten += bkg; |
---|
| 469 | return(inten); |
---|
| 470 | } |
---|
| 471 | |
---|
| 472 | double |
---|
| 473 | GaussLorentzGel(double dp[], double x) |
---|
| 474 | { |
---|
| 475 | //[0] Gaussian scale factor |
---|
| 476 | //[1] Gaussian (static) screening length |
---|
| 477 | //[2] Lorentzian (fluctuation) scale factor |
---|
| 478 | //[3] Lorentzian screening length |
---|
| 479 | //[4] incoherent background |
---|
| 480 | |
---|
| 481 | double Ig0,gg,Il0,ll,bgd,inten; |
---|
| 482 | |
---|
| 483 | Ig0 = dp[0]; |
---|
| 484 | gg = dp[1]; |
---|
| 485 | Il0 = dp[2]; |
---|
| 486 | ll = dp[3]; |
---|
| 487 | bgd = dp[4]; |
---|
| 488 | |
---|
| 489 | inten = Ig0*exp(-1.0*x*x*gg*gg/2.0) + Il0/(1.0 + (x*ll)*(x*ll)) + bgd; |
---|
| 490 | |
---|
| 491 | return(inten); |
---|
| 492 | } |
---|
| 493 | |
---|
| 494 | |
---|
| 495 | double |
---|
| 496 | GaussianShell(double w[], double x) |
---|
| 497 | { |
---|
| 498 | // variables are: |
---|
| 499 | //[0] scale |
---|
| 500 | //[1] radius (ᅵ) |
---|
| 501 | //[2] thick (ᅵ) (thickness parameter - this is the std. dev. of the Gaussian width of the shell) |
---|
| 502 | //[3] polydispersity of the radius |
---|
| 503 | //[4] sld shell (ᅵ-2) |
---|
| 504 | //[5] sld solvent |
---|
| 505 | //[6] background (cm-1) |
---|
| 506 | |
---|
| 507 | double scale,rad,delrho,bkg,del,thick,pd,sig,pi; |
---|
| 508 | double t1,t2,t3,t4,retval,exfact,vshell,vexcl,sldShell,sldSolvent; |
---|
| 509 | scale = w[0]; |
---|
| 510 | rad = w[1]; |
---|
| 511 | thick = w[2]; |
---|
| 512 | pd = w[3]; |
---|
| 513 | sldShell = w[4]; |
---|
| 514 | sldSolvent = w[5]; |
---|
| 515 | bkg = w[6]; |
---|
| 516 | |
---|
| 517 | delrho = w[4] - w[5]; |
---|
| 518 | sig = pd*rad; |
---|
| 519 | |
---|
| 520 | pi = 4.0*atan(1.0); |
---|
| 521 | |
---|
| 522 | ///APPROXIMATION (see eqn 4 - but not a bad approximation) |
---|
| 523 | // del is the equivalent shell thickness with sharp boundaries, centered at mean radius |
---|
| 524 | del = thick*sqrt(2.0*pi); |
---|
| 525 | |
---|
| 526 | // calculate the polydisperse shell volume and the excluded volume |
---|
| 527 | vshell=4.0*pi/3.0*( pow((rad+del/2.0),3) - pow((rad-del/2.0),3) ) *(1.0+pd*pd); |
---|
| 528 | vexcl=4.0*pi/3.0*( pow((rad+del/2.0),3) ) *(1.0+pd*pd); |
---|
| 529 | |
---|
| 530 | //intensity, eqn 9(a-d) |
---|
| 531 | exfact = exp(-2.0*sig*sig*x*x); |
---|
| 532 | |
---|
| 533 | t1 = 0.5*x*x*thick*thick*thick*thick*(1.0+cos(2.0*x*rad)*exfact); |
---|
| 534 | t2 = x*thick*thick*(rad*sin(2.0*x*rad) + 2.0*x*sig*sig*cos(2.0*x*rad))*exfact; |
---|
| 535 | t3 = 0.5*rad*rad*(1.0-cos(2.0*x*rad)*exfact); |
---|
| 536 | t4 = 0.5*sig*sig*(1.0+4.0*x*rad*sin(2.0*x*rad)*exfact+cos(2.0*x*rad)*(4.0*sig*sig*x*x-1.0)*exfact); |
---|
| 537 | |
---|
| 538 | retval = t1+t2+t3+t4; |
---|
| 539 | retval *= exp(-1.0*x*x*thick*thick); |
---|
| 540 | retval *= (del*del/x/x); |
---|
| 541 | retval *= 16.0*pi*pi*delrho*delrho*scale; |
---|
| 542 | retval *= 1.0e8; |
---|
| 543 | |
---|
| 544 | //NORMALIZED by the AVERAGE shell volume, since scale is the volume fraction of material |
---|
| 545 | // retval /= vshell |
---|
| 546 | retval /= vexcl; |
---|
| 547 | //re-normalize by polydisperse sphere volume, Gaussian distribution |
---|
| 548 | retval /= (1.0+3.0*pd*pd); |
---|
| 549 | |
---|
| 550 | retval += bkg; |
---|
| 551 | |
---|
| 552 | return(retval); |
---|
| 553 | } |
---|
| 554 | |
---|
| 555 | |
---|