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

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

Changes to the XOP code to upgrade to ToolKit? v6. Changes are the ones outlined in the Appendix A of the TK6 manual. SOme of the XOP support routines and the #pragma for 2-byte structures have changed. Per Howard Rodstein, there is no need to change the c files to cpp. c should work and compile just fine.

These changes work correctly on my mac. Next is to make sure that they work correctly on the two build machines.

File size: 9.4 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= (float*)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= (double*)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= (float*)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= (double*)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= (float*)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= (double*)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= (float*)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= (double*)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 = (double*)WaveData(p->SQHandle);
224        qw = (double*)WaveData(p->QHandle);                                     
225        cw = (double*)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 = (double*)WaveData(p->SQHandle);
310        qw = (double*)WaveData(p->QHandle);                                     
311        cw = (double*)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.