source: sasview/sansmodels/src/libigor/libTwoPhase.c @ f9a1279

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 f9a1279 was ae3ce4e, checked in by Mathieu Doucet <doucetm@…>, 17 years ago

Moving sansmodels to trunk
  • Property mode set to 100644
File size: 5.7 KB
RevLine 
[ae3ce4e]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
287
288static double
289gammln(double xx) {
290       
291    double x,y,tmp,ser;
292    static double cof[6]={76.18009172947146,-86.50532032941677,
293                24.01409824083091,-1.231739572450155,
294                0.1208650973866179e-2,-0.5395239384953e-5};
295    int j;
296       
297    y=x=xx;
298    tmp=x+5.5;
299    tmp -= (x+0.5)*log(tmp);
300    ser=1.000000000190015;
301    for (j=0;j<=5;j++) ser += cof[j]/++y;
302    return -tmp+log(2.5066282746310005*ser/x);
303}
304
305
Note: See TracBrowser for help on using the repository browser.