source: sans/XOP_Dev/SANSAnalysis/lib/libTwoPhase.c @ 453

Last change on this file since 453 was 453, checked in by srkline, 14 years ago

Additions to the library of the 2008 model functions. Direct proting of the Igor code, duplicated by these XOPs

File size: 10.4 KB
Line 
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...
9double
10TeubnerStreyModel(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
22double
23Power_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
38double
39Peak_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
54double
55Peak_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
71double
72Lorentz_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
87double
88Fractal(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
118double
119DAB_Model(double dp[], double q)
120{
121        double qval,inten;
122        double Izero, range, incoh;
123       
124        qval= q;
125        Izero = dp[0];
126        range = dp[1];
127        incoh = dp[2]; 
128       
129        inten = Izero/pow((1.0 + (qval*range)*(qval*range)),2) + incoh;
130       
131        return(inten);
132}
133
134// G. Beaucage's Unified Model (1-4 levels)
135//
136double
137OneLevel(double dp[], double q)
138{
139        double x,ans,erf1;
140        double G1,Rg1,B1,Pow1,bkg,scale;
141       
142        x=q;
143        scale = dp[0];
144        G1 = dp[1];
145        Rg1 = dp[2];
146        B1 = dp[3];
147        Pow1 = dp[4];
148        bkg = dp[5];
149       
150        erf1 = erf( (x*Rg1/sqrt(6.0)));
151       
152        ans = G1*exp(-x*x*Rg1*Rg1/3.0);
153        ans += B1*pow((erf1*erf1*erf1/x),Pow1);
154       
155        ans *= scale;
156        ans += bkg;
157        return(ans);
158}
159
160// G. Beaucage's Unified Model (1-4 levels)
161//
162double
163TwoLevel(double dp[], double q)
164{
165        double x;
166        double ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg;
167        double erf1,erf2,scale;
168       
169        x=q;
170       
171        scale = dp[0];
172        G1 = dp[1];     //equivalent to I(0)
173        Rg1 = dp[2];
174        B1 = dp[3];
175        Pow1 = dp[4];
176        G2 = dp[5];
177        Rg2 = dp[6];
178        B2 = dp[7];
179        Pow2 = dp[8];
180        bkg = dp[9];
181       
182        erf1 = erf( (x*Rg1/sqrt(6.0)) );
183        erf2 = erf( (x*Rg2/sqrt(6.0)) );
184        //Print erf1
185       
186        ans = G1*exp(-x*x*Rg1*Rg1/3.0);
187        ans += B1*exp(-x*x*Rg2*Rg2/3.0)*pow((erf1*erf1*erf1/x),Pow1);
188        ans += G2*exp(-x*x*Rg2*Rg2/3.0);
189        ans += B2*pow((erf2*erf2*erf2/x),Pow2);
190       
191        ans *= scale;
192        ans += bkg;
193       
194        return(ans);
195}
196
197// G. Beaucage's Unified Model (1-4 levels)
198//
199double
200ThreeLevel(double dp[], double q)
201{
202        double x;
203        double ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg;
204        double G3,Rg3,B3,Pow3,erf3;
205        double erf1,erf2,scale;
206       
207        x=q;
208       
209        scale = dp[0];
210        G1 = dp[1];     //equivalent to I(0)
211        Rg1 = dp[2];
212        B1 = dp[3];
213        Pow1 = dp[4];
214        G2 = dp[5];
215        Rg2 = dp[6];
216        B2 = dp[7];
217        Pow2 = dp[8];
218        G3 = dp[9];
219        Rg3 = dp[10];
220        B3 = dp[11];
221        Pow3 = dp[12];
222        bkg = dp[13];
223       
224        erf1 = erf( (x*Rg1/sqrt(6.0)) );
225        erf2 = erf( (x*Rg2/sqrt(6.0)) );
226        erf3 = erf( (x*Rg3/sqrt(6.0)) );
227        //Print erf1
228       
229        ans = G1*exp(-x*x*Rg1*Rg1/3.0) + B1*exp(-x*x*Rg2*Rg2/3.0)*pow((erf1*erf1*erf1/x),Pow1);
230        ans += G2*exp(-x*x*Rg2*Rg2/3.0) + B2*exp(-x*x*Rg3*Rg3/3.0)*pow((erf2*erf2*erf2/x),Pow2);
231        ans += G3*exp(-x*x*Rg3*Rg3/3.0) + B3*pow((erf3*erf3*erf3/x),Pow3);
232       
233        ans *= scale;
234        ans += bkg;
235       
236        return(ans);
237}
238
239// G. Beaucage's Unified Model (1-4 levels)
240//
241double
242FourLevel(double dp[], double q)
243{
244        double x;
245        double ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg;
246        double G3,Rg3,B3,Pow3,erf3;
247        double G4,Rg4,B4,Pow4,erf4;
248        double erf1,erf2,scale;
249       
250        x=q;
251       
252        scale = dp[0];
253        G1 = dp[1];     //equivalent to I(0)
254        Rg1 = dp[2];
255        B1 = dp[3];
256        Pow1 = dp[4];
257        G2 = dp[5];
258        Rg2 = dp[6];
259        B2 = dp[7];
260        Pow2 = dp[8];
261        G3 = dp[9];
262        Rg3 = dp[10];
263        B3 = dp[11];
264        Pow3 = dp[12];
265        G4 = dp[13];
266        Rg4 = dp[14];
267        B4 = dp[15];
268        Pow4 = dp[16];
269        bkg = dp[17];
270       
271        erf1 = erf( (x*Rg1/sqrt(6.0)) );
272        erf2 = erf( (x*Rg2/sqrt(6.0)) );
273        erf3 = erf( (x*Rg3/sqrt(6.0)) );
274        erf4 = erf( (x*Rg4/sqrt(6.0)) );
275       
276        ans = G1*exp(-x*x*Rg1*Rg1/3.0) + B1*exp(-x*x*Rg2*Rg2/3.0)*pow((erf1*erf1*erf1/x),Pow1);
277        ans += G2*exp(-x*x*Rg2*Rg2/3.0) + B2*exp(-x*x*Rg3*Rg3/3.0)*pow((erf2*erf2*erf2/x),Pow2);
278        ans += G3*exp(-x*x*Rg3*Rg3/3.0) + B3*exp(-x*x*Rg4*Rg4/3.0)*pow((erf3*erf3*erf3/x),Pow3);
279        ans += G4*exp(-x*x*Rg4*Rg4/3.0) + B4*pow((erf4*erf4*erf4/x),Pow4);
280       
281        ans *= scale;
282        ans += bkg;
283       
284    return(ans);
285}
286
287double
288BroadPeak(double dp[], double q)
289{
290        // variables are:                                                       
291        //[0] Porod term scaling
292        //[1] Porod exponent
293        //[2] Lorentzian term scaling
294        //[3] Lorentzian screening length [A]
295        //[4] peak location [1/A]
296        //[5] Lorentzian exponent
297        //[6] background
298       
299        double aa,nn,cc,LL,Qzero,mm,bgd,inten,qval;
300        qval= q;
301        aa = dp[0];
302        nn = dp[1];
303        cc = dp[2]; 
304        LL = dp[3]; 
305        Qzero = dp[4]; 
306        mm = dp[5]; 
307        bgd = dp[6]; 
308       
309        inten = aa/pow(qval,nn);
310        inten += cc/(1.0 + pow((fabs(qval-Qzero)*LL),mm) );
311        inten += bgd;
312       
313        return(inten);
314}
315
316double
317CorrLength(double dp[], double q)
318{
319        // variables are:                                                       
320        //[0] Porod term scaling
321        //[1] Porod exponent
322        //[2] Lorentzian term scaling
323        //[3] Lorentzian screening length [A]
324        //[4] Lorentzian exponent
325        //[5] background
326       
327        double aa,nn,cc,LL,mm,bgd,inten,qval;
328        qval= q;
329        aa = dp[0];
330        nn = dp[1];
331        cc = dp[2]; 
332        LL = dp[3]; 
333        mm = dp[4]; 
334        bgd = dp[5]; 
335       
336        inten = aa/pow(qval,nn);
337        inten += cc/(1.0 + pow((qval*LL),mm) );
338        inten += bgd;
339       
340        return(inten);
341}
342
343double
344TwoLorentzian(double dp[], double q)
345{
346        // variables are:                                                       
347        //[0] Lorentzian term scaling
348        //[1] Lorentzian screening length [A]
349        //[2] Lorentzian exponent
350        //[3] Lorentzian #2 term scaling
351        //[4] Lorentzian #2 screening length [A]
352        //[5] Lorentzian #2 exponent
353        //[6] background
354               
355        double aa,LL1,nn,cc,LL2,mm,bgd,inten,qval;
356        qval= q;
357        aa = dp[0];
358        LL1 = dp[1];
359        nn = dp[2]; 
360        cc = dp[3]; 
361        LL2 = dp[4]; 
362        mm = dp[5]; 
363        bgd= dp[6];
364       
365        inten = aa/(1.0 + pow((qval*LL1),nn) );
366        inten += cc/(1.0 + pow((qval*LL2),mm) );
367        inten += bgd;
368       
369        return(inten);
370}
371
372double
373TwoPowerLaw(double dp[], double q)
374{
375        //[0] Coefficient
376        //[1] (-) Power @ low Q
377        //[2] (-) Power @ high Q
378        //[3] crossover Q-value
379        //[4] incoherent background
380               
381        double A, m1,m2,qc,bgd,scale,inten,qval;
382        qval= q;
383        A = dp[0];
384        m1 = dp[1];
385        m2 = dp[2]; 
386        qc = dp[3]; 
387        bgd = dp[4]; 
388       
389        if(qval<=qc){
390                inten = A*pow(qval,-1.0*m1);
391        } else {
392                scale = A*pow(qc,-1.0*m1) / pow(qc,-1.0*m2);
393                inten = scale*pow(qval,-1.0*m2);
394        }
395       
396        inten += bgd;
397       
398        return(inten);
399}
400
401double
402PolyGaussCoil(double dp[], double x)
403{
404        //w[0] = scale
405        //w[1] = radius of gyration []
406        //w[2] = polydispersity, ratio of Mw/Mn
407        //w[3] = bkg [cm-1]
408               
409        double scale,bkg,Rg,uval,Mw_Mn,inten,xi;
410
411        scale = dp[0];
412        Rg = dp[1];
413        Mw_Mn = dp[2]; 
414        bkg = dp[3]; 
415       
416        uval = Mw_Mn - 1.0;
417        if(uval == 0) {
418                uval = 1e-6;            //avoid divide by zero error
419        }
420       
421        xi = Rg*Rg*x*x/(1.0+2.0*uval);
422       
423        if(xi < 1e-3) {
424                return(scale+bkg);              //limiting value
425        }
426       
427        inten = 2.0*(pow((1.0+uval*xi),(-1.0/uval))+xi-1.0);
428        inten /= (1.0+uval)*xi*xi;
429
430        inten *= scale;
431        //add in the background
432        inten += bkg;   
433        return(inten);
434}
435
436double
437GaussLorentzGel(double dp[], double x)
438{
439        //[0] Gaussian scale factor
440        //[1] Gaussian (static) screening length
441        //[2] Lorentzian (fluctuation) scale factor
442        //[3] Lorentzian screening length
443        //[4] incoherent background
444               
445        double Ig0,gg,Il0,ll,bgd,inten;
446
447        Ig0 = dp[0];
448        gg = dp[1];
449        Il0 = dp[2]; 
450        ll = dp[3];
451        bgd = dp[4]; 
452       
453        inten = Ig0*exp(-1.0*x*x*gg*gg/2.0) + Il0/(1 + (x*ll)*(x*ll)) + bgd;
454       
455        return(inten);
456}
457
458
459static double
460gammln(double xx) {
461       
462    double x,y,tmp,ser;
463    static double cof[6]={76.18009172947146,-86.50532032941677,
464                24.01409824083091,-1.231739572450155,
465                0.1208650973866179e-2,-0.5395239384953e-5};
466    int j;
467       
468    y=x=xx;
469    tmp=x+5.5;
470    tmp -= (x+0.5)*log(tmp);
471    ser=1.000000000190015;
472    for (j=0;j<=5;j++) ser += cof[j]/++y;
473    return -tmp+log(2.5066282746310005*ser/x);
474}
475
476double
477GaussianShell(double w[], double x)
478{
479        // variables are:                                                       
480        //[0] scale
481        //[1] radius ()
482        //[2] thick () (thickness parameter - this is the std. dev. of the Gaussian width of the shell)
483        //[3] polydispersity of the radius
484        //[4] sld shell (-2)
485        //[5] sld solvent
486        //[6] background (cm-1)
487       
488        double scale,rad,delrho,bkg,del,thick,pd,sig,pi;
489        double t1,t2,t3,t4,retval,exfact,vshell,vexcl,sldShell,sldSolvent;
490        scale = w[0];
491        rad = w[1];
492        thick = w[2];
493        pd = w[3];
494        sldShell = w[4];
495        sldSolvent = w[5];
496        bkg = w[6];
497       
498        delrho = w[4] - w[5];
499        sig = pd*rad;
500       
501        pi = 4.0*atan(1.0);
502       
503        ///APPROXIMATION (see eqn 4 - but not a bad approximation)
504        // del is the equivalent shell thickness with sharp boundaries, centered at mean radius
505        del = thick*sqrt(2.0*pi);
506       
507        // calculate the polydisperse shell volume and the excluded volume
508        vshell=4.0*pi/3.0*( pow((rad+del/2.0),3) - pow((rad-del/2.0),3) ) *(1.0+pd*pd);
509        vexcl=4.0*pi/3.0*( pow((rad+del/2.0),3) ) *(1.0+pd*pd);
510       
511        //intensity, eqn 9(a-d)
512        exfact = exp(-2.0*sig*sig*x*x);
513       
514        t1 = 0.5*x*x*thick*thick*thick*thick*(1.0+cos(2.0*x*rad)*exfact);
515        t2 = x*thick*thick*(rad*sin(2.0*x*rad) + 2.0*x*sig*sig*cos(2.0*x*rad))*exfact;
516        t3 = 0.5*rad*rad*(1.0-cos(2.0*x*rad)*exfact);
517        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);
518       
519        retval = t1+t2+t3+t4;
520        retval *= exp(-1.0*x*x*thick*thick);
521        retval *= (del*del/x/x);
522        retval *= 16.0*pi*pi*delrho*delrho*scale;
523        retval *= 1.0e8;
524       
525        //NORMALIZED by the AVERAGE shell volume, since scale is the volume fraction of material
526//      retval /= vshell
527        retval /= vexcl;
528        //re-normalize by polydisperse sphere volume, Gaussian distribution
529        retval /= (1.0+3.0*pd*pd);
530       
531        retval += bkg;
532       
533        return(retval);
534}
535
536
Note: See TracBrowser for help on using the repository browser.