source: sans/Analysis/branches/ajj_23APR07/XOPs/SANSAnalysis/XOP/Func2D.c @ 235

Last change on this file since 235 was 235, checked in by srkline, 15 years ago

Changes to XOP library functions to convert models to take individual SLD's rather than contrast.

Individual SLD's are easier to work with - since they are experimental values, contrast is not.

This, unfortunately, means that what Mathieu has done may need to be tweaked, since our library has been tweaked.

File size: 14.3 KB
Line 
1/*
2 *  Func2D.c
3 *  SANSAnalysis
4 *
5 * Steve Kline Jan 2008
6 * these are the XOP calls to DANSE 2D functions
7 * from the library supplied by M. Doucet
8 *
9 */
10
11
12#include "XOPStandardHeaders.h"                 // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
13#include "SANSAnalysis.h"                               // structures, two-byte aligned
14#include "libSANSAnalysis.h"                    // functions from the libSANS   
15#include "danse.h"                                              // functions from the DANSE part of libSANS
16#include <math.h>
17#include "Func2D.h"                                             // declarations for the 2D functions in this file
18
19int
20Cylinder_2D(FitParams2DPtr p)
21{
22        double *dp;
23       
24        float *fp;                              // Pointer to single precision wave data.
25        double qx;
26        double qy;
27        double q, phi;
28        double pars[11];
29//      int i;
30//      char buf[256];
31       
32        if (p->waveHandle == NIL) {
33                SetNaN64(&p->result);
34                return NON_EXISTENT_WAVE;
35        }
36   
37        qx = p->qx;
38        qy = p->qy;
39       
40//      sprintf(buf, "Qx = %g, Qy = %g\r",qx, qy);
41//      XOPNotice(buf);
42
43        q = hypot(qx,qy);
44        phi = atan2(qy,qx);
45       
46        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
47                case NT_FP32:
48                        fp= WaveData(p->waveHandle);
49            SetNaN64(&p->result);
50                        return REQUIRES_SP_OR_DP_WAVE; //not quite true, but good enough for now AJJ 4/23/07           
51                case NT_FP64:
52                        dp= WaveData(p->waveHandle);
53                       
54//                      for(i=0; i<11; i++) {
55//                              pars[i] = dp[i];
56//                      }
57                        pars[0] = dp[0];
58                        pars[1] = dp[1];
59                        pars[2] = dp[2];
60                        pars[3] = dp[3] - dp[4];
61                        pars[4] = dp[5];
62                        pars[5] = dp[6];
63                        pars[6] = dp[7];
64                        pars[7] = dp[8];
65                        pars[8] = dp[9];
66                        pars[9] = dp[10];
67                        pars[10] = dp[11];
68
69                        p->result = disperse_cylinder_analytical_2D( pars, q, phi ); 
70                        return 0;
71                default:                                                                // We can't handle this wave data type.
72                        SetNaN64(&p->result);
73                        return REQUIRES_SP_OR_DP_WAVE;
74        }
75       
76        return 0;
77}
78
79int
80Cylinder_2D_Weight2D(FitParams2DWeightPtr p)
81{
82        double *dp;
83       
84        float *fp;                              // Pointer to single precision wave data.
85        double qx;
86        double qy;
87        double q, phi;
88        double *par_values;
89        double *weight_values;
90        double pars[13];
91        int i, i_theta;
92        double sum;
93        int n_slices;
94       
95        if (p->waveHandle == NIL) {
96                SetNaN64(&p->result);
97                return NON_EXISTENT_WAVE;
98        }
99   
100        qx = p->qx;
101        qy = p->qy;
102       
103        q = hypot(qx,qy);
104        phi = atan2(qy,qx);
105       
106        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
107                case NT_FP32:
108                        fp= WaveData(p->waveHandle);
109            SetNaN64(&p->result);
110                        return REQUIRES_SP_OR_DP_WAVE; //not quite true, but good enough for now AJJ 4/23/07           
111                case NT_FP64:
112                        dp= WaveData(p->waveHandle);
113                        par_values = WaveData(p->par_values);
114                        weight_values = WaveData(p->weight_values);
115                       
116                        for(i=0; i<13; i++) {
117                                pars[i] = dp[i];
118                        }
119
120                        sum = 0.0;
121                        n_slices = (int)floor(dp[13]);
122
123                        //XOPOKAlert("test","This is a test");
124
125                        if(n_slices == 0) {
126                                p->result =weight_dispersion( &disperse_cylinder_analytical_2D,
127                                                        par_values, weight_values, 
128                                                        (int)floor(pars[11]), (int)floor(pars[12]),
129                                                        pars, q, phi );                 
130                        } else {
131                                for(i_theta=0; i_theta<n_slices; i_theta++) {
132                                        SpinProcess();
133                                        // average over theta
134                                       
135                                        // For the cylinder model, theta_cyl=90 degrees
136                                        // will produce a NAN at phi_cyl=0 and 45 degrees
137                                        pars[5] = acos(-1.0)/n_slices * i_theta;
138                                        if( fabs(i_theta / n_slices) - 0.5 < 0.000001 ) {
139                                                //continue;
140                                                pars[5] += 0.00001;
141                                        }
142                                       
143                                        // Multiply by sin(theta) because we are integrating in
144                                        // spherical coordinates
145                                        sum += sin(pars[5])* weight_dispersion( &disperse_cylinder_analytical_2D,
146                                                                par_values, weight_values, 
147                                                                (int)floor(pars[11]), (int)floor(pars[12]),
148                                                                pars, q, phi ); 
149                                }
150                               
151                                p->result = sum/n_slices;
152                        }
153                       
154                        return 0;
155                default:                                                                // We can't handle this wave data type.
156                        SetNaN64(&p->result);
157                        return REQUIRES_SP_OR_DP_WAVE;
158        }
159       
160        return 0;
161}
162
163int
164CoreShellCylinder_2D(FitParams2DPtr p)
165{
166        double *dp;
167       
168        float *fp;                              // Pointer to single precision wave data.
169        double qx;
170        double qy;
171        double q, phi;
172        double pars[15];
173        int i;
174       
175        if (p->waveHandle == NIL) {
176                SetNaN64(&p->result);
177                return NON_EXISTENT_WAVE;
178        }
179   
180        qx = p->qx;
181        qy = p->qy;
182       
183        q = hypot(qx,qy);
184        phi = atan2(qy,qx);
185       
186        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
187                case NT_FP32:
188                        fp= WaveData(p->waveHandle);
189            SetNaN64(&p->result);
190                        return REQUIRES_SP_OR_DP_WAVE; //not quite true, but good enough for now AJJ 4/23/07           
191                case NT_FP64:
192                        dp= WaveData(p->waveHandle);
193                       
194                        for(i=0; i<15; i++) {
195                                pars[i] = dp[i];
196                        }
197
198                        p->result = disperse_core_shell_cylinder_analytical_2D( pars, q, phi ); 
199                        return 0;
200                default:                                                                // We can't handle this wave data type.
201                        SetNaN64(&p->result);
202                        return REQUIRES_SP_OR_DP_WAVE;
203        }
204       
205        return 0;
206}
207
208int
209CoreShellCylinder_2D_Weight2D(FitParams2DWeightPtr p)
210{
211        double *dp;
212       
213        float *fp;                              // Pointer to single precision wave data.
214        double qx;
215        double qy;
216        double q, phi;
217        double *par_values;
218        double *weight_values;
219        double pars[17];
220        int i, i_theta;
221        double sum;
222        int n_slices;
223       
224        if (p->waveHandle == NIL) {
225                SetNaN64(&p->result);
226                return NON_EXISTENT_WAVE;
227        }
228   
229        qx = p->qx;
230        qy = p->qy;
231       
232        q = hypot(qx,qy);
233        phi = atan2(qy,qx);
234       
235        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
236                case NT_FP32:
237                        fp= WaveData(p->waveHandle);
238            SetNaN64(&p->result);
239                        return REQUIRES_SP_OR_DP_WAVE; //not quite true, but good enough for now AJJ 4/23/07           
240                case NT_FP64:
241                        dp= WaveData(p->waveHandle);
242                        par_values = WaveData(p->par_values);
243                        weight_values = WaveData(p->weight_values);
244                       
245                        for(i=0; i<17; i++) {
246                                pars[i] = dp[i];
247                        }
248
249                        sum = 0.0;
250                        n_slices = (int)floor(dp[17]);
251
252                        if(n_slices == 0) {
253                                p->result =weight_dispersion( &disperse_core_shell_cylinder_analytical_2D,
254                                                        par_values, weight_values, 
255                                                        (int)floor(pars[15]), (int)floor(pars[16]),
256                                                        pars, q, phi );                 
257                        } else {
258                                for(i_theta=0; i_theta<n_slices; i_theta++) {
259                                        SpinProcess();
260                                        // average over theta
261                                       
262                                        // For the cylinder model, theta_cyl=90 degrees
263                                        // will produce a NAN at phi_cyl=0 and 45 degrees
264                                        // TODO: integrate from 0 to pi/2 instead of 0 to pi
265                                        pars[8] = acos(-1.0)/n_slices * i_theta;
266                                        if( fabs(i_theta / n_slices) - 0.5 < 0.000001 ) {
267                                                //continue;
268                                                pars[8] += 0.00001;
269                                        }
270                                       
271                                        // Multiply by sin(theta) because we are integrating in
272                                        // spherical coordinates
273                                        sum += sin(pars[8])* weight_dispersion( &disperse_core_shell_cylinder_analytical_2D,
274                                                                par_values, weight_values, 
275                                                                (int)floor(pars[15]), (int)floor(pars[16]),
276                                                                pars, q, phi ); 
277                                }
278                               
279                                p->result = sum/n_slices;
280                        }
281                       
282                        return 0;
283                default:                                                                // We can't handle this wave data type.
284                        SetNaN64(&p->result);
285                        return REQUIRES_SP_OR_DP_WAVE;
286        }
287       
288        return 0;
289}
290
291int
292Ellipsoid_2D(FitParams2DPtr p)
293{
294        double *dp;
295       
296        float *fp;                              // Pointer to single precision wave data.
297        double qx;
298        double qy;
299        double q, phi;
300        double pars[12];
301//      int i;
302       
303        if (p->waveHandle == NIL) {
304                SetNaN64(&p->result);
305                return NON_EXISTENT_WAVE;
306        }
307   
308        qx = p->qx;
309        qy = p->qy;
310       
311        q = hypot(qx,qy);
312        phi = atan2(qy,qx);
313       
314        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
315                case NT_FP32:
316                        fp= WaveData(p->waveHandle);
317            SetNaN64(&p->result);
318                        return REQUIRES_SP_OR_DP_WAVE; //not quite true, but good enough for now AJJ 4/23/07           
319                case NT_FP64:
320                        dp= WaveData(p->waveHandle);
321                       
322//                      for(i=0; i<12; i++) {
323//                              pars[i] = dp[i];
324//                      }
325                        pars[0] = dp[0];
326                        pars[1] = dp[1];
327                        pars[2] = dp[2];
328                        pars[3] = dp[3] - dp[4];
329                        pars[4] = dp[5];
330                        pars[5] = dp[6];
331                        pars[6] = dp[7];
332                        pars[7] = dp[8];
333                        pars[8] = dp[9];
334                        pars[9] = dp[10];
335                        pars[10] = dp[11];
336                        pars[11] = dp[12];
337
338                        //p->result = 1.0;
339                        p->result = disperse_ellipsoid_analytical_2D( pars, q, phi ); 
340                        return 0;
341                default:                                                                // We can't handle this wave data type.
342                        SetNaN64(&p->result);
343                        return REQUIRES_SP_OR_DP_WAVE;
344        }
345       
346        return 0;
347}
348
349int
350Ellipsoid_2D_Weight2D(FitParams2DWeightPtr p)
351{
352        double *dp;
353       
354        float *fp;                              // Pointer to single precision wave data.
355        double qx;
356        double qy;
357        double q, phi;
358        double *par_values;
359        double *weight_values;
360        double pars[14];
361        int i, i_theta;
362        double sum;
363        int n_slices;
364       
365        if (p->waveHandle == NIL) {
366                SetNaN64(&p->result);
367                return NON_EXISTENT_WAVE;
368        }
369   
370        qx = p->qx;
371        qy = p->qy;
372       
373        q = hypot(qx,qy);
374        phi = atan2(qy,qx);
375       
376        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
377                case NT_FP32:
378                        fp= WaveData(p->waveHandle);
379            SetNaN64(&p->result);
380                        return REQUIRES_SP_OR_DP_WAVE; //not quite true, but good enough for now AJJ 4/23/07           
381                case NT_FP64:
382                        dp= WaveData(p->waveHandle);
383                        par_values = WaveData(p->par_values);
384                        weight_values = WaveData(p->weight_values);
385                       
386                        for(i=0; i<14; i++) {
387                                pars[i] = dp[i];
388                        }
389
390                        sum = 0.0;
391                        n_slices = (int)floor(dp[14]);
392
393                        if(n_slices == 0) {
394                                p->result =weight_dispersion( &disperse_ellipsoid_analytical_2D,
395                                                        par_values, weight_values, 
396                                                        (int)floor(pars[12]), (int)floor(pars[13]),
397                                                        pars, q, phi );                 
398                        } else {
399                                for(i_theta=0; i_theta<n_slices; i_theta++) {
400                                        SpinProcess();
401                                        // average over theta
402                                       
403                                        // For the cylinder model, theta_cyl=90 degrees
404                                        // will produce a NAN at phi_cyl=0 and 45 degrees
405                                        pars[5] = acos(-1.0)/n_slices * i_theta;
406                                        if( fabs(i_theta / n_slices) - 0.5 < 0.000001 ) {
407                                                //continue;
408                                                pars[5] += 0.00001;
409                                        }
410                                       
411                                        // Multiply by sin(theta) because we are integrating in
412                                        // spherical coordinates
413                                        sum += sin(pars[5])* weight_dispersion( &disperse_ellipsoid_analytical_2D,
414                                                                par_values, weight_values, 
415                                                                (int)floor(pars[12]), (int)floor(pars[13]),
416                                                                pars, q, phi ); 
417                                }
418                               
419                                p->result = sum/n_slices;
420                        }
421                       
422                        return 0;
423                default:                                                                // We can't handle this wave data type.
424                        SetNaN64(&p->result);
425                        return REQUIRES_SP_OR_DP_WAVE;
426        }
427       
428        return 0;
429}
430
431
432int
433EllipticalCylinder_2D(FitParams2DPtr p)
434{
435        double *dp;
436       
437        float *fp;                              // Pointer to single precision wave data.
438        double qx;
439        double qy;
440        double q, phi;
441        double pars[14];
442//      int i;
443       
444        if (p->waveHandle == NIL) {
445                SetNaN64(&p->result);
446                return NON_EXISTENT_WAVE;
447        }
448   
449        qx = p->qx;
450        qy = p->qy;
451       
452        q = hypot(qx,qy);
453        phi = atan2(qy,qx);
454       
455        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
456                case NT_FP32:
457                        fp= WaveData(p->waveHandle);
458            SetNaN64(&p->result);
459                        return REQUIRES_SP_OR_DP_WAVE; //not quite true, but good enough for now AJJ 4/23/07           
460                case NT_FP64:
461                        dp= WaveData(p->waveHandle);
462                       
463//                      for(i=0; i<14; i++) {
464//                              pars[i] = dp[i];
465//                      }
466                        pars[0] = dp[0];
467                        pars[1] = dp[1];
468                        pars[2] = dp[2];
469                        pars[3] = dp[3];
470                        pars[4] = dp[4] - dp[5];
471                        pars[5] = dp[6];
472                        pars[6] = dp[7];
473                        pars[7] = dp[8];
474                        pars[8] = dp[9];
475                        pars[9] = dp[10];
476                        pars[10] = dp[11];
477                        pars[11] = dp[12];
478                        pars[12] = dp[13];
479                        pars[13] = dp[14];
480
481                       
482                        p->result = disperse_elliptical_cylinder_analytical_2D( pars, q, phi ); 
483                         
484                        return 0;
485                default:                                                                // We can't handle this wave data type.
486                        SetNaN64(&p->result);
487                        return REQUIRES_SP_OR_DP_WAVE;
488        }
489       
490        return 0;
491}
492
493int
494EllipticalCylinder_2D_Weight2D(FitParams2DWeightPtr p)
495{
496        double *dp;
497       
498        float *fp;                              // Pointer to single precision wave data.
499        double qx;
500        double qy;
501        double q, phi;
502        double *par_values;
503        double *weight_values;
504        double pars[16];
505        int i, i_theta;
506        double sum;
507        int n_slices;
508       
509        if (p->waveHandle == NIL) {
510                SetNaN64(&p->result);
511                return NON_EXISTENT_WAVE;
512        }
513   
514        qx = p->qx;
515        qy = p->qy;
516       
517        q = hypot(qx,qy);
518        phi = atan2(qy,qx);
519       
520        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
521                case NT_FP32:
522                        fp= WaveData(p->waveHandle);
523            SetNaN64(&p->result);
524                        return REQUIRES_SP_OR_DP_WAVE; //not quite true, but good enough for now AJJ 4/23/07           
525                case NT_FP64:
526                        dp= WaveData(p->waveHandle);
527                        par_values = WaveData(p->par_values);
528                        weight_values = WaveData(p->weight_values);
529                       
530                        for(i=0; i<16; i++) {
531                                pars[i] = dp[i];
532                        }
533
534                        sum = 0.0;
535                        n_slices = (int)floor(dp[16]);
536
537                        if(n_slices == 0) {
538                                p->result =weight_dispersion( &disperse_elliptical_cylinder_analytical_2D,
539                                                        par_values, weight_values, 
540                                                        (int)floor(pars[14]), (int)floor(pars[15]),
541                                                        pars, q, phi );                 
542                        } else {
543                                for(i_theta=0; i_theta<n_slices; i_theta++) {
544                                        SpinProcess();
545                                        // average over theta
546                                       
547                                        // For the cylinder model, theta_cyl=90 degrees
548                                        // will produce a NAN at phi_cyl=0 and 45 degrees
549                                        // TODO: integrate from 0 to pi/2 instead of 0 to pi
550                                       
551                                        pars[6] = acos(-1.0)/n_slices * i_theta;
552                                       
553                                         if( fabs(i_theta / n_slices) - 0.5 < 0.000001 ) {
554                                                //continue;
555                                                pars[6] += 0.00001;
556                                        }
557                                       
558                                        // Multiply by sin(theta) because we are integrating in
559                                        // spherical coordinates
560                                        sum += sin(pars[6])* weight_dispersion( &disperse_elliptical_cylinder_analytical_2D,
561                                                                par_values, weight_values, 
562                                                                (int)floor(pars[14]), (int)floor(pars[15]),
563                                                                pars, q, phi ); 
564                                }
565                               
566                                p->result = sum/n_slices;
567                        }
568                       
569                        return 0;
570                default:                                                                // We can't handle this wave data type.
571                        SetNaN64(&p->result);
572                        return REQUIRES_SP_OR_DP_WAVE;
573        }
574       
575        return 0;
576}
577
578int
579Sphere_2D(FitParams2DPtr p)
580{
581        double *dp;
582       
583        float *fp;                              // Pointer to single precision wave data.
584        double qx;
585        double qy;
586//      double q, phi;
587        double pars[5];
588        int i;
589//      char buf[256];
590       
591        if (p->waveHandle == NIL) {
592                SetNaN64(&p->result);
593                return NON_EXISTENT_WAVE;
594        }
595   
596        qx = p->qx;
597        qy = p->qy;
598       
599//      sprintf(buf, "Qx = %g, Qy = %g\r",qx, qy);
600//      XOPNotice(buf);
601       
602
603                //not needed for a symmetric scattering function like this
604//      q = hypot(qx,qy);
605//      phi = atan2(qy,qx);
606       
607        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
608                case NT_FP32:
609                        fp= WaveData(p->waveHandle);
610            SetNaN64(&p->result);
611                        return REQUIRES_SP_OR_DP_WAVE; //not quite true, but good enough for now AJJ 4/23/07           
612                case NT_FP64:
613                        dp= WaveData(p->waveHandle);
614                       
615                        for(i=0; i<5; i++) {
616                                pars[i] = dp[i];
617                        }
618
619                        p->result = SphereForm(pars, sqrt(qx*qx+qy*qy));                //kind of a trivial example...
620                        return 0;
621                default:                                                                // We can't handle this wave data type.
622                        SetNaN64(&p->result);
623                        return REQUIRES_SP_OR_DP_WAVE;
624        }
625       
626        return 0;
627}
Note: See TracBrowser for help on using the repository browser.