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