source: sans/XOP_Dev/SANSAnalysis/XOP/StructureFactor.c @ 823

Last change on this file since 823 was 823, checked in by srkline, 11 years ago

changed 2-yukawa S(Q) calculation to return S(q)=1000 if the calculation fails, to be very obvious that the calculation failed. This change matches the Igof code.

File size: 9.2 KB
Line 
1/*      SimpleFit.c
2
3A simplified project designed to act as a template for your curve fitting function.
4The fitting function is a simple polynomial. It works but is of no practical use.
5*/
6
7
8#include "XOPStandardHeaders.h"                 // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
9#include "SANSAnalysis.h"
10#include "libSANSAnalysis.h"
11#include "StructureFactor.h"
12#include "2Y_OneYukawa.h"
13#include "2Y_TwoYukawa.h"
14
15
16//Hard Sphere Structure Factor
17//
18int
19HardSphereStructX(FitParamsPtr p)
20{
21        double *dp;                             // Pointer to double precision wave data.
22        float *fp;                              // Pointer to single precision wave data.
23        double q;
24       
25        if (p->waveHandle == NIL) {
26                SetNaN64(&p->result);
27                return NON_EXISTENT_WAVE;
28        }
29       
30        q= p->x;
31       
32        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
33                case NT_FP32:
34                        fp= WaveData(p->waveHandle);
35                        SetNaN64(&p->result);
36                        return REQUIRES_SP_OR_DP_WAVE; //not quite true, but good enough for now AJJ 4/23/07                   
37                case NT_FP64:
38                        dp= WaveData(p->waveHandle);
39                        p->result = HardSphereStruct(dp,q);
40                        return 0;
41                default:                                                                // We can't handle this wave data type.
42                        SetNaN64(&p->result);
43                        return REQUIRES_SP_OR_DP_WAVE;
44        }
45        return 0;
46}
47
48//Sticky Hard Sphere Structure Factor
49//
50int
51StickyHS_StructX(FitParamsPtr p)
52{
53        double *dp;                             // Pointer to double precision wave data.
54        float *fp;                              // Pointer to single precision wave data.
55        double q;
56       
57        if (p->waveHandle == NIL) {
58                SetNaN64(&p->result);
59                return NON_EXISTENT_WAVE;
60        }
61       
62        q= p->x;
63       
64        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
65                case NT_FP32:
66                        fp= WaveData(p->waveHandle);
67                        SetNaN64(&p->result);
68                        return REQUIRES_SP_OR_DP_WAVE; //not quite true, but good enough for now AJJ 4/23/07                   
69                case NT_FP64:
70                        dp= WaveData(p->waveHandle);
71                        p->result = StickyHS_Struct(dp,q);
72                        return 0;
73                default:                                                                // We can't handle this wave data type.
74                        SetNaN64(&p->result);
75                        return REQUIRES_SP_OR_DP_WAVE;
76        }
77        return 0;
78}
79
80
81
82//     SUBROUTINE SQWELL: CALCULATES THE STRUCTURE FACTOR FOR A
83//                        DISPERSION OF MONODISPERSE HARD SPHERES
84//     IN THE Mean Spherical APPROXIMATION ASSUMING THE SPHERES
85//     INTERACT THROUGH A SQUARE WELL POTENTIAL.
86//** not the best choice of closure ** see note below
87//     REFS:  SHARMA,SHARMA, PHYSICA 89A,(1977),212
88int
89SquareWellStructX(FitParamsPtr p)
90{
91        double *dp;                             // Pointer to double precision wave data.
92        float *fp;                              // Pointer to single precision wave data.
93        double q;
94       
95        if (p->waveHandle == NIL) {
96                SetNaN64(&p->result);
97                return NON_EXISTENT_WAVE;
98        }
99       
100        q= p->x;
101       
102        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
103                case NT_FP32:
104                        fp= WaveData(p->waveHandle);
105                        SetNaN64(&p->result);
106                        XOPNotice("I think it's single precision\r");
107                        return REQUIRES_SP_OR_DP_WAVE; //not quite true, but good enough for now AJJ 4/23/07                   
108                case NT_FP64:
109                        dp= WaveData(p->waveHandle);
110                        p->result = SquareWellStruct(dp,q);
111                        return 0;
112                default:                                                                // We can't handle this wave data type.
113                        SetNaN64(&p->result);
114                        XOPNotice("I don't know what it is\r");
115                        return REQUIRES_SP_OR_DP_WAVE;
116        }
117       
118        return 0;
119}
120
121// Hayter-Penfold (rescaled) MSA structure factor for screened Coulomb interactions
122//
123int
124HayterPenfoldMSAX(FitParamsPtr p)
125{
126        double *dp;                             // Pointer to double precision wave data.
127        float *fp;                              // Pointer to single precision wave data.
128        double q;
129       
130        if (p->waveHandle == NIL) {
131                SetNaN64(&p->result);
132                return NON_EXISTENT_WAVE;
133        }
134       
135        q= p->x;
136       
137        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
138                case NT_FP32:
139                        fp= WaveData(p->waveHandle);
140                        SetNaN64(&p->result);
141                        return REQUIRES_SP_OR_DP_WAVE; //not quite true, but good enough for now AJJ 4/23/07                   
142                case NT_FP64:
143                        dp= WaveData(p->waveHandle);
144                        p->result = HayterPenfoldMSA(dp,q);
145                        return 0;
146                default:                                                                // We can't handle this wave data type.
147                        SetNaN64(&p->result);
148                        return REQUIRES_SP_OR_DP_WAVE;
149        }
150       
151        return 0;
152}
153
154
155// called as DiamCylX(hcyl,rcyl)
156int
157DiamCylX(DiamParamsPtr p)
158{
159        double hcyl,rcyl;
160       
161        hcyl = p->p1;
162        rcyl = p->p2;
163   
164        p->result = DiamCyl(hcyl,rcyl);
165       
166        return(0);
167}
168
169//prolate OR oblate ellipsoids
170//aa is the axis of rotation
171//if aa>bb, then PROLATE
172//if aa<bb, then OBLATE
173// A. Isihara, J. Chem. Phys. 18, 1446 (1950)
174//returns DIAMETER
175// called as DiamEllipX(aa,bb)
176int
177DiamEllipX(DiamParamsPtr p)
178{
179       
180        double aa,bb;
181       
182        aa = p->p1;
183        bb = p->p2;
184       
185        p->result = DiamEllip(aa,bb);
186       
187        return(0);
188}
189
190
191//
192//
193int
194OneYukawaX(FitParamsPtr_Yuk p)
195{
196        double *sq;                     //pointer to sq wave
197        double *qw;                     //pointer to q wave
198        double *cw;                     //pointer to coef wave
199        long npnts,i;                   //number of points in waves
200       
201       
202        double a, b, c, d;     
203        int debug = 0,check,ok;                 //always keep this to zero...
204        //default parameters
205        double Z,K,phi,radius;
206        char buf[256];
207       
208        /* check that wave handles are all valid */
209        if (p->SQHandle == NIL) {
210                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
211                return(NON_EXISTENT_WAVE);
212        }
213        if (p->QHandle == NIL) {
214                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
215                return(NON_EXISTENT_WAVE);
216        }
217        if (p->CoefHandle == NIL) {
218                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
219                return(NON_EXISTENT_WAVE);
220        }       
221       
222        //get the wave data
223        sq = WaveData(p->SQHandle);
224        qw = WaveData(p->QHandle);                                     
225        cw = WaveData(p->CoefHandle);                                   
226        npnts = WavePoints(p->QHandle);                                         // Number of points in q wave.
227       
228        phi = cw[0];
229        radius = cw[1];
230        K = cw[2];
231        Z = cw[3];
232       
233        if(fabs(Z) < 0.1) Z = 0.1;                      // values near zero are very bad for the calculation
234       
235        //      sprintf(buf, "Input OK, phi,radius,K,Z = %g %g %g %g\r",phi,radius,K,Z);
236        //      XOPNotice(buf);
237       
238        debug = 0;
239        ok = Y_SolveEquations( Z, K, phi, &a, &b, &c, &d, debug );
240       
241        if( ok )
242        {
243                //              sprintf(buf, "Y_SolveEquations OK\r");
244                //              XOPNotice(buf);
245               
246                check = Y_CheckSolution( Z, K, phi, a, b, c, d );
247                if(debug) {
248                        sprintf(buf, "solution = (%g, %g, %g, %g) check = %d\r", a, b, c, d, check );
249                        XOPNotice(buf);
250                }
251        }
252       
253        //      sprintf(buf, "Y_CheckSolution OK\r");
254        //      XOPNotice(buf);
255       
256        // loop through and calculate the S(q), or return 1 if the solution failed
257        //      if(check) {             //from the converted c-code
258       
259        if(ok) {                //less restrictive, if a solution found, return it, even if the equations aren't quite satisfied
260               
261                for (i = 0; i < npnts; i++) {
262                        sq[i] = SqOneYukawa(qw[i]*radius*2.0, Z, K, phi, a, b, c, d);
263                }       
264        } else {
265                for (i = 0; i < npnts; i++) {
266                        sq[i] = 1000;           // return a really bogus value, as Yun suggests
267                }
268        }
269       
270        // mark the returned wave as updated so that the graph etc. automatically updates       
271        WaveHandleModified(p->SQHandle);
272       
273        p->retVal = 0;
274       
275        return 0;
276}
277
278
279int
280TwoYukawaX(FitParamsPtr_Yuk p)
281{
282        double *sq;                     //pointer to sq wave
283        double *qw;                     //pointer to q wave
284        double *cw;                     //pointer to coef wave
285        long npnts,i;                   //number of points in waves
286       
287       
288        double a, b, c1, c2, d1, d2;   
289        int debug = 0,check,ok;                 //always keep this to zero...
290        //default parameters
291        double Z1,K1,Z2,K2,phi,radius;
292        char buf[256];
293       
294        /* check that wave handles are all valid */
295        if (p->SQHandle == NIL) {
296                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
297                return(NON_EXISTENT_WAVE);
298        }
299        if (p->QHandle == NIL) {
300                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
301                return(NON_EXISTENT_WAVE);
302        }
303        if (p->CoefHandle == NIL) {
304                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
305                return(NON_EXISTENT_WAVE);
306        }       
307       
308        //get the wave data
309        sq = WaveData(p->SQHandle);
310        qw = WaveData(p->QHandle);                                     
311        cw = WaveData(p->CoefHandle);                                   
312        npnts = WavePoints(p->QHandle);                                         // Number of points in q wave.
313       
314        phi = cw[0];
315        radius = cw[1];
316        K1 = cw[2];
317        Z1 = cw[3];
318        K2 = cw[4];
319        Z2 = cw[5];
320       
321        if(fabs(Z1) < 0.001) Z1 = 0.001;                        // values near zero are very bad for the calculation
322        if(fabs(Z2) < 0.001) Z2 = 0.001;                        // values near zero are very bad for the calculation
323        if(fabs(K1) < 0.001) K1 = 0.001;                        // values near zero are very bad for the calculation
324        if(fabs(K2) < 0.001) K2 = 0.001;                        // values near zero are very bad for the calculation
325       
326       
327        //      sprintf(buf, "Input OK, phi,radius,K1,Z1,K2,Z2 = %g %g %g %g %g %g\r",phi,radius,K1,Z1,K2,Z2);
328        //      XOPNotice(buf);
329       
330        debug = 0;
331       
332        ok = TY_SolveEquations( Z1, Z2, K1, K2, phi, &a, &b, &c1, &c2, &d1, &d2, debug );
333        if( ok )
334        {
335                //              sprintf(buf, "TY_SolveEquations OK \r");
336                //              XOPNotice(buf);
337               
338               
339                check = TY_CheckSolution( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 );
340                if(debug) {
341                        sprintf(buf, "solution = (%g, %g, %g, %g, %g, %g) check = %d\r", a, b, c1, c2, d1, d2, check );
342                        XOPNotice(buf);
343                }
344        }
345        //      sprintf(buf, "TY_CheckSolution OK \r");
346        //      XOPNotice(buf);
347       
348        // loop through and calculate the S(q), or return 1 if the solution failed
349        //      if(check) {             //from the converted c-code
350       
351        if(ok) {                //less restrictive, if a solution found, return it, even if the equations aren't quite satisfied
352               
353               
354                for (i = 0; i < npnts; i++) {
355                        sq[i] = SqTwoYukawa(qw[i]*radius*2.0, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2);
356                }       
357        } else {
358                for (i = 0; i < npnts; i++) {
359                        sq[i] = 1000;           // return a really bogus value, as Yun suggests
360                }
361        }
362       
363        // mark the returned wave as updated so that the graph etc. automatically updates       
364        WaveHandleModified(p->SQHandle);
365       
366        p->retVal = 0;
367       
368        return 0;
369}
370
371
372
373
374
375///////////end of XOP
376
377
Note: See TracBrowser for help on using the repository browser.