source: sasview/sansmodels/src/sans/models/libigor/libTwoPhase.c @ 34c2649

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 34c2649 was 34c2649, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Re #4 Fixed warnings
  • Property mode set to 100644
File size: 10.2 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
118// 6 JUL 2009 SRK changed definition of Izero scale factor to be uncorrelated with range
119//
120double
121DAB_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//
138double
139OneLevel(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//
168double
169TwoLevel(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//
209double
210ThreeLevel(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//
255double
256FourLevel(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
305double
306BroadPeak(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
334double
335CorrLength(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
361double
362TwoLorentzian(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
390double
391TwoPowerLaw(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
419double
420PolyGaussCoil(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
454double
455GaussLorentzGel(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
477double
478GaussianShell(double w[], double x)
479{
480        // variables are:                                                       
481        //[0] scale
482        //[1] radius (ï¿œ)
483        //[2] thick (ï¿œ) (thickness parameter - this is the std. dev. of the Gaussian width of the shell)
484        //[3] polydispersity of the radius
485        //[4] sld shell (ï¿œ-2)
486        //[5] sld solvent
487        //[6] background (cm-1)
488       
489        double scale,rad,delrho,bkg,del,thick,pd,sig,pi;
490        double t1,t2,t3,t4,retval,exfact,vshell,vexcl,sldShell,sldSolvent;
491        scale = w[0];
492        rad = w[1];
493        thick = w[2];
494        pd = w[3];
495        sldShell = w[4];
496        sldSolvent = w[5];
497        bkg = w[6];
498       
499        delrho = w[4] - w[5];
500        sig = pd*rad;
501       
502        pi = 4.0*atan(1.0);
503       
504        ///APPROXIMATION (see eqn 4 - but not a bad approximation)
505        // del is the equivalent shell thickness with sharp boundaries, centered at mean radius
506        del = thick*sqrt(2.0*pi);
507       
508        // calculate the polydisperse shell volume and the excluded volume
509        vshell=4.0*pi/3.0*( pow((rad+del/2.0),3) - pow((rad-del/2.0),3) ) *(1.0+pd*pd);
510        vexcl=4.0*pi/3.0*( pow((rad+del/2.0),3) ) *(1.0+pd*pd);
511       
512        //intensity, eqn 9(a-d)
513        exfact = exp(-2.0*sig*sig*x*x);
514       
515        t1 = 0.5*x*x*thick*thick*thick*thick*(1.0+cos(2.0*x*rad)*exfact);
516        t2 = x*thick*thick*(rad*sin(2.0*x*rad) + 2.0*x*sig*sig*cos(2.0*x*rad))*exfact;
517        t3 = 0.5*rad*rad*(1.0-cos(2.0*x*rad)*exfact);
518        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);
519       
520        retval = t1+t2+t3+t4;
521        retval *= exp(-1.0*x*x*thick*thick);
522        retval *= (del*del/x/x);
523        retval *= 16.0*pi*pi*delrho*delrho*scale;
524        retval *= 1.0e8;
525       
526        //NORMALIZED by the AVERAGE shell volume, since scale is the volume fraction of material
527//      retval /= vshell
528        retval /= vexcl;
529        //re-normalize by polydisperse sphere volume, Gaussian distribution
530        retval /= (1.0+3.0*pd*pd);
531       
532        retval += bkg;
533       
534        return(retval);
535}
536
537
Note: See TracBrowser for help on using the repository browser.