source: sasview/sansmodels/src/libigor/libTwoPhase.c @ 637631d

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 637631d was 67424cd, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Reorganizing models in preparation of cpp cleanup
  • Property mode set to 100644
File size: 10.7 KB
Line 
1/*      TwoPhaseFit.c
2
3*/
4
5#include "StandardHeaders.h"                    // Include ANSI headers, Mac headers
6#include "libTwoPhase.h"
7
8/* internal functions */
9static double
10gammln(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...
27double
28TeubnerStreyModel(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
40double
41Power_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
56double
57Peak_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
72double
73Peak_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
89double
90Lorentz_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
105double
106Fractal(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//
138double
139DAB_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//
156double
157OneLevel(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//
186double
187TwoLevel(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//
227double
228ThreeLevel(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//
273double
274FourLevel(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
323double
324BroadPeak(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
352double
353CorrLength(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
379double
380TwoLorentzian(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
408double
409TwoPowerLaw(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
437double
438PolyGaussCoil(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
472double
473GaussLorentzGel(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
495double
496GaussianShell(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
Note: See TracBrowser for help on using the repository browser.