source: sans/Analysis/branches/ajj_23APR07/XOPs/SANSAnalysis/TwoPhase.c @ 93

Last change on this file since 93 was 93, checked in by ajj, 15 years ago

Add combined XOP code. Currently only contains XCode project file to build Universal binary suitable for Igor 6.

File size: 19.3 KB
Line 
1/*      TwoPhaseFit.c
2
3*/
4
5#pragma XOP_SET_STRUCT_PACKING                  // All structures are 2-byte-aligned.
6
7#include "XOPStandardHeaders.h"                 // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
8#include "SANSAnalysis.h"
9#include "TwoPhase.h"
10
11// scattering from the Teubner-Strey model for microemulsions - hardly needs to be an XOP...
12int
13TeubnerStreyModelX(FitParamsPtr p)
14{
15        DOUBLE *dp;                             // Pointer to double precision wave data.
16        float *fp;                              // Pointer to single precision wave data.
17        DOUBLE q;
18        DOUBLE inten,q2,q4;             //my local names
19       
20        if (p->waveHandle == NIL) {
21                SetNaN64(&p->result);
22                return NON_EXISTENT_WAVE;
23        }
24       
25        q= p->x;
26        q2 = q*q;
27        q4 = q2*q2;
28       
29        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
30                case NT_FP32:
31                        fp= WaveData(p->waveHandle);
32//                        scale = fp[0];
33//                        radius = fp[1];
34//                        delrho = fp[2];
35//                        bkg = fp[3];
36               
37                        inten = 1.0/(fp[0]+fp[1]*q2+fp[2]*q4);
38                        inten += fp[3];
39                       
40                        break;
41                case NT_FP64:
42                        dp= WaveData(p->waveHandle);
43                       
44                        inten = 1.0/(dp[0]+dp[1]*q2+dp[2]*q4);
45                        inten += dp[3]; 
46                                                                     
47                        break;
48                default:                                                                // We can't handle this wave data type.
49                        SetNaN64(&p->result);
50                        return REQUIRES_SP_OR_DP_WAVE;
51        }
52               
53        p->result= (inten);     //scale, and add in the background
54        return 0;
55}
56
57int
58Power_Law_ModelX(FitParamsPtr p)
59{
60        DOUBLE *dp;                             // Pointer to double precision wave data.
61        float *fp;                              // Pointer to single precision wave data.
62        DOUBLE qval;
63        DOUBLE inten,A,m,bgd;           //my local names
64       
65        if (p->waveHandle == NIL) {
66                SetNaN64(&p->result);
67                return NON_EXISTENT_WAVE;
68        }
69       
70        qval= p->x;
71       
72        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
73                case NT_FP32:
74                        fp= WaveData(p->waveHandle);
75                        A = fp[0];
76                        m = fp[1];
77                        bgd = fp[2];
78                                       
79                        break;
80                case NT_FP64:
81                        dp= WaveData(p->waveHandle);
82                        A = dp[0];
83                        m = dp[1];
84                        bgd = dp[2];
85                                                                                             
86                        break;
87                default:                                                                // We can't handle this wave data type.
88                        SetNaN64(&p->result);
89                        return REQUIRES_SP_OR_DP_WAVE;
90        }
91       
92        inten = A*pow(qval,-m) + bgd;
93        p->result= (inten);     //scale, and add in the background
94        return 0;
95}
96
97
98int
99Peak_Lorentz_ModelX(FitParamsPtr p)
100{
101        DOUBLE *dp;                             // Pointer to double precision wave data.
102        float *fp;                              // Pointer to single precision wave data.
103        DOUBLE qval;
104        DOUBLE inten,I0, qpk, dq,bgd;           //my local names
105           
106        if (p->waveHandle == NIL) {
107                SetNaN64(&p->result);
108                return NON_EXISTENT_WAVE;
109        }
110       
111        qval= p->x;
112       
113        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
114                case NT_FP32:
115                        fp= WaveData(p->waveHandle);
116                        I0 = fp[0];
117                        qpk = fp[1];
118                        dq = fp[2];
119                        bgd = fp[3];
120                                       
121                        break;
122                case NT_FP64:
123                        dp= WaveData(p->waveHandle);
124                        I0 = dp[0];
125                        qpk = dp[1];
126                        dq = dp[2];
127                        bgd = dp[3];
128                                                                                             
129                        break;
130                default:                                                                // We can't handle this wave data type.
131                        SetNaN64(&p->result);
132                        return REQUIRES_SP_OR_DP_WAVE;
133        }
134       
135        inten = I0/(1.0 + pow( (qval-qpk)/dq,2) ) + bgd;
136
137        p->result= (inten);     //scale, and add in the background
138        return 0;
139}
140
141int
142Peak_Gauss_ModelX(FitParamsPtr p)
143{
144        DOUBLE *dp;                             // Pointer to double precision wave data.
145        float *fp;                              // Pointer to single precision wave data.
146        DOUBLE qval;
147        DOUBLE inten,I0, qpk, dq,bgd;           //my local names
148           
149        if (p->waveHandle == NIL) {
150                SetNaN64(&p->result);
151                return NON_EXISTENT_WAVE;
152        }
153       
154        qval= p->x;
155       
156        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
157                case NT_FP32:
158                        fp= WaveData(p->waveHandle);
159                        I0 = fp[0];
160                        qpk = fp[1];
161                        dq = fp[2];
162                        bgd = fp[3];
163                                       
164                        break;
165                case NT_FP64:
166                        dp= WaveData(p->waveHandle);
167                        I0 = dp[0];
168                        qpk = dp[1];
169                        dq = dp[2];
170                        bgd = dp[3];
171                                                                                             
172                        break;
173                default:                                                                // We can't handle this wave data type.
174                        SetNaN64(&p->result);
175                        return REQUIRES_SP_OR_DP_WAVE;
176        }
177       
178        inten = I0*exp(-0.5*pow((qval-qpk)/dq,2))+ bgd;
179
180        p->result= (inten);     //scale, and add in the background
181        return 0;
182}
183
184int
185Lorentz_ModelX(FitParamsPtr p)
186{
187        DOUBLE *dp;                             // Pointer to double precision wave data.
188        float *fp;                              // Pointer to single precision wave data.
189        DOUBLE qval;
190        DOUBLE inten,I0, L,bgd;         //my local names
191           
192        if (p->waveHandle == NIL) {
193                SetNaN64(&p->result);
194                return NON_EXISTENT_WAVE;
195        }
196           
197        qval= p->x;
198       
199        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
200                case NT_FP32:
201                        fp= WaveData(p->waveHandle);
202                        I0 = fp[0];
203                        L = fp[1];
204                        bgd = fp[2];
205                                       
206                        break;
207                case NT_FP64:
208                        dp= WaveData(p->waveHandle);
209                        I0 = dp[0];
210                        L = dp[1];
211                        bgd = dp[2];
212                                                                                             
213                        break;
214                default:                                                                // We can't handle this wave data type.
215                        SetNaN64(&p->result);
216                        return REQUIRES_SP_OR_DP_WAVE;
217        }
218       
219        inten = I0/(1.0 + (qval*L)*(qval*L)) + bgd;
220
221        p->result= (inten);     //scale, and add in the background
222        return 0;
223}
224
225int
226FractalX(FitParamsPtr p)
227{
228        DOUBLE *dp;                             // Pointer to double precision wave data.
229        float *fp;                              // Pointer to single precision wave data.
230        DOUBLE x,pi;
231        DOUBLE r0,Df,corr,phi,sldp,sldm,bkg;
232        DOUBLE pq,sq,ans;
233       
234        if (p->waveHandle == NIL) {
235                SetNaN64(&p->result);
236                return NON_EXISTENT_WAVE;
237        }
238       
239        pi = 4.0*atan(1.0);
240        x= p->x;
241       
242        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
243                case NT_FP32:
244                        fp= WaveData(p->waveHandle);
245                        phi = fp[0];            // volume fraction of building block spheres...
246                        r0 = fp[1];             //  radius of building block
247                        Df = fp[2];             //  fractal dimension
248                        corr = fp[3];           //  correlation length of fractal-like aggregates
249                        sldp = fp[4];           // SLD of building block
250                        sldm = fp[5];           // SLD of matrix or solution
251                        bkg = fp[6];            //  flat background
252                                                               
253                        break;
254                case NT_FP64:
255                        dp= WaveData(p->waveHandle);
256                        phi = dp[0];            // volume fraction of building block spheres...
257                        r0 = dp[1];             //  radius of building block
258                        Df = dp[2];             //  fractal dimension
259                        corr = dp[3];           //  correlation length of fractal-like aggregates
260                        sldp = dp[4];           // SLD of building block
261                        sldm = dp[5];           // SLD of matrix or solution
262                        bkg = dp[6];            //  flat background
263                                                                                             
264                        break;
265                default:                                                                // We can't handle this wave data type.
266                        SetNaN64(&p->result);
267                        return REQUIRES_SP_OR_DP_WAVE;
268        }
269       
270        //calculate P(q) for the spherical subunits, units cm-1 sr-1
271        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);
272
273        //calculate S(q)
274        sq = Df*exp(gammln(Df-1.0))*sin((Df-1.0)*atan(x*corr));
275        sq /= pow((x*r0),Df) * pow((1.0 + 1.0/(x*corr)/(x*corr)),((Df-1.0)/2.0));
276        sq += 1.0;
277        //combine and return
278        ans = pq*sq + bkg;
279       
280        p->result= (ans);       //scale, and add in the background
281        return 0;
282}
283
284int
285DAB_ModelX(FitParamsPtr p)
286{
287        DOUBLE *dp;                             // Pointer to double precision wave data.
288        float *fp;                              // Pointer to single precision wave data.
289        DOUBLE qval,inten;
290        DOUBLE Izero, range, incoh;
291       
292        if (p->waveHandle == NIL) {
293                SetNaN64(&p->result);
294                return NON_EXISTENT_WAVE;
295        }
296       
297        qval= p->x;
298       
299       
300        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
301                case NT_FP32:
302                        fp= WaveData(p->waveHandle);
303                        Izero = fp[0];
304                        range = fp[1];
305                        incoh = fp[2]; 
306                                                                                     
307                        break;
308                case NT_FP64:
309                        dp= WaveData(p->waveHandle);
310                        Izero = dp[0];
311                        range = dp[1];
312                        incoh = dp[2]; 
313                                                                                             
314                        break;
315                default:                                                                // We can't handle this wave data type.
316                        SetNaN64(&p->result);
317                        return REQUIRES_SP_OR_DP_WAVE;
318        }
319       
320        inten = Izero/pow((1.0 + (qval*range)*(qval*range)),2) + incoh;
321               
322        p->result= (inten);     //scale, and add in the background
323        return 0;
324}
325
326// G. Beaucage's Unified Model (1-4 levels)
327//
328int
329OneLevelX(FitParamsPtr p)
330{
331        DOUBLE *dp;                             // Pointer to double precision wave data.
332        float *fp;                              // Pointer to single precision wave data.
333        DOUBLE x,ans,erf1;
334        DOUBLE G1,Rg1,B1,Pow1,bkg,scale;
335       
336        if (p->waveHandle == NIL) {
337                SetNaN64(&p->result);
338                return NON_EXISTENT_WAVE;
339        }
340       
341        x= p->x;
342
343        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
344                case NT_FP32:
345                        fp= WaveData(p->waveHandle);
346                                                scale = fp[0];
347                        G1 = fp[1];
348                        Rg1 = fp[2];
349                        B1 = fp[3];
350                        Pow1 = fp[4];
351                        bkg = fp[5];
352                       
353                                                                                     
354                        break;
355                case NT_FP64:
356                        dp= WaveData(p->waveHandle);
357                                                scale = dp[0];
358                        G1 = dp[1];
359                        Rg1 = dp[2];
360                        B1 = dp[3];
361                        Pow1 = dp[4];
362                        bkg = dp[5];
363                       
364                                                                                             
365                        break;
366                default:                                                                // We can't handle this wave data type.
367                        SetNaN64(&p->result);
368                        return REQUIRES_SP_OR_DP_WAVE;
369        }
370       
371        erf1 = erf( (x*Rg1/sqrt(6.0)));
372       
373        ans = G1*exp(-x*x*Rg1*Rg1/3.0);
374        ans += B1*pow((erf1*erf1*erf1/x),Pow1);
375       
376        ans *= scale;
377        ans += bkg;
378               
379        p->result= ans; //scale, and add in the background
380        return 0;
381}
382
383// G. Beaucage's Unified Model (1-4 levels)
384//
385int
386TwoLevelX(FitParamsPtr p)
387{
388        DOUBLE *dp;                             // Pointer to double precision wave data.
389        float *fp;                              // Pointer to single precision wave data.
390        DOUBLE x;
391        DOUBLE ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg;
392        DOUBLE erf1,erf2,scale;
393       
394        if (p->waveHandle == NIL) {
395                SetNaN64(&p->result);
396                return NON_EXISTENT_WAVE;
397        }
398       
399        x= p->x;
400
401        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
402                case NT_FP32:
403                        fp= WaveData(p->waveHandle);
404                                                scale = fp[0];
405                        G1 = fp[1];     //equivalent to I(0)
406                        Rg1 = fp[2];
407                        B1 = fp[3];
408                        Pow1 = fp[4];
409                        G2 = fp[5];
410                        Rg2 = fp[6];
411                        B2 = fp[7];
412                        Pow2 = fp[8];
413                        bkg = fp[9];
414                                                                                                             
415                        break;
416                case NT_FP64:
417                        dp= WaveData(p->waveHandle);
418                                                scale = dp[0];
419                        G1 = dp[1];     //equivalent to I(0)
420                        Rg1 = dp[2];
421                        B1 = dp[3];
422                        Pow1 = dp[4];
423                        G2 = dp[5];
424                        Rg2 = dp[6];
425                        B2 = dp[7];
426                        Pow2 = dp[8];
427                        bkg = dp[9];
428                                                                                             
429                        break;
430                default:                                                                // We can't handle this wave data type.
431                        SetNaN64(&p->result);
432                        return REQUIRES_SP_OR_DP_WAVE;
433        }
434       
435        erf1 = erf( (x*Rg1/sqrt(6.0)) );
436        erf2 = erf( (x*Rg2/sqrt(6.0)) );
437        //Print erf1
438       
439        ans = G1*exp(-x*x*Rg1*Rg1/3.0);
440        ans += B1*exp(-x*x*Rg2*Rg2/3.0)*pow((erf1*erf1*erf1/x),Pow1);
441        ans += G2*exp(-x*x*Rg2*Rg2/3.0);
442        ans += B2*pow((erf2*erf2*erf2/x),Pow2);
443       
444        ans *= scale;
445        ans += bkg;
446               
447        p->result= ans; //scale, and add in the background
448        return 0;
449}
450
451// G. Beaucage's Unified Model (1-4 levels)
452//
453int
454ThreeLevelX(FitParamsPtr p)
455{
456        DOUBLE *dp;                             // Pointer to double precision wave data.
457        float *fp;                              // Pointer to single precision wave data.
458        DOUBLE x;
459        DOUBLE ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg;
460        DOUBLE G3,Rg3,B3,Pow3,erf3;
461        DOUBLE erf1,erf2,scale;
462       
463        if (p->waveHandle == NIL) {
464                SetNaN64(&p->result);
465                return NON_EXISTENT_WAVE;
466        }
467       
468        x= p->x;
469
470        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
471                case NT_FP32:
472                        fp= WaveData(p->waveHandle);
473                                                scale = fp[0];
474                        G1 = fp[1];     //equivalent to I(0)
475                        Rg1 = fp[2];
476                        B1 = fp[3];
477                        Pow1 = fp[4];
478                        G2 = fp[5];
479                        Rg2 = fp[6];
480                        B2 = fp[7];
481                        Pow2 = fp[8];
482                        G3 = fp[9];
483                        Rg3 = fp[10];
484                        B3 = fp[11];
485                        Pow3 = fp[12];
486                        bkg = fp[13];
487                                                                                     
488                        break;
489                case NT_FP64:
490                        dp= WaveData(p->waveHandle);
491                        scale = dp[0];
492                                                G1 = dp[1];     //equivalent to I(0)
493                        Rg1 = dp[2];
494                        B1 = dp[3];
495                        Pow1 = dp[4];
496                        G2 = dp[5];
497                        Rg2 = dp[6];
498                        B2 = dp[7];
499                        Pow2 = dp[8];
500                        G3 = dp[9];
501                        Rg3 = dp[10];
502                        B3 = dp[11];
503                        Pow3 = dp[12];
504                        bkg = dp[13];
505                                                                                                                     
506                        break;
507                default:                                                                // We can't handle this wave data type.
508                        SetNaN64(&p->result);
509                        return REQUIRES_SP_OR_DP_WAVE;
510        }
511       
512        erf1 = erf( (x*Rg1/sqrt(6.0)) );
513        erf2 = erf( (x*Rg2/sqrt(6.0)) );
514        erf3 = erf( (x*Rg3/sqrt(6.0)) );
515        //Print erf1
516       
517        ans = G1*exp(-x*x*Rg1*Rg1/3.0) + B1*exp(-x*x*Rg2*Rg2/3.0)*pow((erf1*erf1*erf1/x),Pow1);
518        ans += G2*exp(-x*x*Rg2*Rg2/3.0) + B2*exp(-x*x*Rg3*Rg3/3.0)*pow((erf2*erf2*erf2/x),Pow2);
519        ans += G3*exp(-x*x*Rg3*Rg3/3.0) + B3*pow((erf3*erf3*erf3/x),Pow3);
520       
521        ans *= scale;
522        ans += bkg;
523               
524        p->result= ans; //scale, and add in the background
525        return 0;
526}
527
528// G. Beaucage's Unified Model (1-4 levels)
529//
530int
531FourLevelX(FitParamsPtr p)
532{
533        DOUBLE *dp;                             // Pointer to double precision wave data.
534        float *fp;                              // Pointer to single precision wave data.
535        DOUBLE x;
536        DOUBLE ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg;
537        DOUBLE G3,Rg3,B3,Pow3,erf3;
538        DOUBLE G4,Rg4,B4,Pow4,erf4;
539        DOUBLE erf1,erf2,scale;
540               
541        if (p->waveHandle == NIL) {
542                SetNaN64(&p->result);
543                return NON_EXISTENT_WAVE;
544        }
545       
546        x= p->x;
547
548        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
549                case NT_FP32:
550                        fp= WaveData(p->waveHandle);
551                        scale = fp[0];
552                                                G1 = fp[1];     //equivalent to I(0)
553                        Rg1 = fp[2];
554                        B1 = fp[3];
555                        Pow1 = fp[4];
556                        G2 = fp[5];
557                        Rg2 = fp[6];
558                        B2 = fp[7];
559                        Pow2 = fp[8];
560                        G3 = fp[9];
561                        Rg3 = fp[10];
562                        B3 = fp[11];
563                        Pow3 = fp[12];
564                        G4 = fp[13];
565                        Rg4 = fp[14];
566                        B4 = fp[15];
567                        Pow4 = fp[16];
568                        bkg = fp[17];
569                                                                                     
570                        break;
571                case NT_FP64:
572                        dp= WaveData(p->waveHandle);
573                        scale = dp[0];
574                                                G1 = dp[1];     //equivalent to I(0)
575                        Rg1 = dp[2];
576                        B1 = dp[3];
577                        Pow1 = dp[4];
578                        G2 = dp[5];
579                        Rg2 = dp[6];
580                        B2 = dp[7];
581                        Pow2 = dp[8];
582                        G3 = dp[9];
583                        Rg3 = dp[10];
584                        B3 = dp[11];
585                        Pow3 = dp[12];
586                        G4 = dp[13];
587                        Rg4 = dp[14];
588                        B4 = dp[15];
589                        Pow4 = dp[16];
590                        bkg = dp[17];
591                                                                                                                     
592                        break;
593                default:                                                                // We can't handle this wave data type.
594                        SetNaN64(&p->result);
595                        return REQUIRES_SP_OR_DP_WAVE;
596        }
597       
598        erf1 = erf( (x*Rg1/sqrt(6.0)) );
599        erf2 = erf( (x*Rg2/sqrt(6.0)) );
600        erf3 = erf( (x*Rg3/sqrt(6.0)) );
601        erf4 = erf( (x*Rg4/sqrt(6.0)) );
602       
603        ans = G1*exp(-x*x*Rg1*Rg1/3.0) + B1*exp(-x*x*Rg2*Rg2/3.0)*pow((erf1*erf1*erf1/x),Pow1);
604        ans += G2*exp(-x*x*Rg2*Rg2/3.0) + B2*exp(-x*x*Rg3*Rg3/3.0)*pow((erf2*erf2*erf2/x),Pow2);
605        ans += G3*exp(-x*x*Rg3*Rg3/3.0) + B3*exp(-x*x*Rg4*Rg4/3.0)*pow((erf3*erf3*erf3/x),Pow3);
606        ans += G4*exp(-x*x*Rg4*Rg4/3.0) + B4*pow((erf4*erf4*erf4/x),Pow4);
607       
608        ans *= scale;
609        ans += bkg;
610               
611        p->result= ans; //scale, and add in the background
612        return 0;
613}
614
615
616static DOUBLE
617gammln(double xx) {
618
619    double x,y,tmp,ser;
620    static double cof[6]={76.18009172947146,-86.50532032941677,
621            24.01409824083091,-1.231739572450155,
622            0.1208650973866179e-2,-0.5395239384953e-5};
623    int j;
624
625    y=x=xx;
626    tmp=x+5.5;
627    tmp -= (x+0.5)*log(tmp);
628    ser=1.000000000190015;
629    for (j=0;j<=5;j++) ser += cof[j]/++y;
630    return -tmp+log(2.5066282746310005*ser/x);
631}
632
633
634
635#pragma XOP_RESET_STRUCT_PACKING                        // All structures are 2-byte-aligned.
636
637///////////end of XOP
638
639
Note: See TracBrowser for help on using the repository browser.