source: sans/XOP_Dev/MonteCarlo/DebyeSpheres.c

Last change on this file was 997, checked in by srkline, 7 years ago

Updates for Xcode 7, Toolkit 7, and 64-bit version

  • Property svn:executable set to *
File size: 13.8 KB
Line 
1/*      DebyeSpheres.c
2
3
4*/
5
6#include "XOPStandardHeaders.h"                 // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
7#include "DebyeSpheres.h"
8
9//#pragma XOP_SET_STRUCT_PACKING                        // All structures are 2-byte-aligned.
10
11// Prototypes
12HOST_IMPORT int XOPMain(IORecHandle ioRecHandle);
13
14// Custom error codes
15//#define REQUIRES_IGOR_200 1 + FIRST_XOP_ERR
16#define NON_EXISTENT_WAVE 2 + FIRST_XOP_ERR
17#define REQUIRES_SP_OR_DP_WAVE 3 + FIRST_XOP_ERR
18
19
20
21/*      Calculates the scattered intensity from a collection of spheres
22    using Debye's method.
23
24 everything was previously declared as float, some inputs were as double.
25 
26 -- to avoid confusion, calculate everything in double, pass everything in as double
27  -- this is more precision than is necessary, but avoids strange results from type casting
28     behind the scenes
29*/
30int
31DebyeSpheresX(AltiParamsPtr p)
32{
33        double qv;                              // input q-value
34        double ival;                            //output intensity value
35        double *xv,*yv,*zv,*rv;         //pointers to input xyz-rho coordinates
36        CountInt i,j;
37    CountInt npt;
38        double rval,grid,vol,fQR,dum,dij;
39
40       
41
42        // check for all of the required waves
43        if (p->rhowavH == NIL) {
44                SetNaN64(&p->result);
45                return NON_EXISTENT_WAVE;
46        }
47        if (p->zwavH == NIL) {
48                SetNaN64(&p->result);
49                return NON_EXISTENT_WAVE;
50        }
51        if (p->ywavH == NIL) {
52                SetNaN64(&p->result);
53                return NON_EXISTENT_WAVE;
54        }
55        if (p->xwavH == NIL) {
56                SetNaN64(&p->result);
57                return NON_EXISTENT_WAVE;
58        }
59
60
61    //check to see that all are float, not double
62        /*
63        if(WaveType(p->rhowavH) != NT_FP32 ) {
64        SetNaN64(&p->result);
65        return kExpectedNT_FP32;
66    }
67    if(WaveType(p->zwavH) != NT_FP32 ) {
68        SetNaN64(&p->result);
69        return kExpectedNT_FP32;
70    }
71    if(WaveType(p->ywavH) != NT_FP32 ) {
72        SetNaN64(&p->result);
73        return kExpectedNT_FP32;
74    }
75    if(WaveType(p->xwavH) != NT_FP32 ) {
76        SetNaN64(&p->result);
77        return kExpectedNT_FP32;
78    }
79        */
80       
81    //check to see that all are double
82        if(WaveType(p->rhowavH) != NT_FP64 ) {
83        SetNaN64(&p->result);
84        return kExpectedNT_FP64;
85    }
86    if(WaveType(p->zwavH) != NT_FP64 ) {
87        SetNaN64(&p->result);
88        return kExpectedNT_FP64;
89    }
90    if(WaveType(p->ywavH) != NT_FP64 ) {
91        SetNaN64(&p->result);
92        return kExpectedNT_FP64;
93    }
94    if(WaveType(p->xwavH) != NT_FP64 ) {
95        SetNaN64(&p->result);
96        return kExpectedNT_FP64;
97    }
98       
99       
100        // (NO) -- convert the input to float. Do all calculations as float.
101        // do everything in double. no reason to do in float, except for the previous Altivec incarnation
102       
103    rval = p->Rprimary;         // primary sphere radius
104    grid = p->grid;                     // calling program should set this to 0.62*Rprimary
105        qv = p->qval;
106
107//
108    npt = (CountInt) WavePoints(p->xwavH);      //wavePoints returns long, number of XYZ points
109        xv = (double*)WaveData(p->xwavH);               //xyz locations
110        yv = (double*)WaveData(p->ywavH);
111        zv = (double*)WaveData(p->zwavH);
112        rv = (double*)WaveData(p->rhowavH);
113   
114       
115        vol = 4.0*3.1415927/3.0*rval*rval*rval;
116        ival = 0.0;
117        fQR = PhiQR(qv,rval);
118        //do the i=j sum
119        for(i=0;i<npt;i+=1) {
120                dum = rv[i]*vol*fQR;
121                ival += dum*dum;
122        }
123        //do the i!=j double sum
124        dum = vol*vol*fQR*fQR;
125        for(i=0;i<npt;i+=1) {
126                for(j=(i+1);j<npt;j+=1) {
127                        dij=XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j]) * grid;
128                        ival += 2.0*rv[i]*rv[j]*dum*sin(dij*qv)/dij/qv;
129                }
130        }
131
132    p->result = ival;
133       
134    return 0;
135
136}
137
138
139double PhiQR(double qval, double rval)
140{
141        double retval,qr;
142       
143        qr = qval*rval;
144        retval = ( 3*(sin(qr) - qr*cos(qr))/qr/qr/qr );
145        return(retval);
146}
147
148double XYZDistance(double x1, double x2,double y1, double y2,double z1, double z2)
149{
150        double retval,dx,dy,dz;
151       
152        dx = (x1-x2);
153        dy = (y1-y2);
154        dz = (z1-z2);
155        retval = sqrt( dx*dx + dy*dy + dz*dz );
156        return(retval);
157}
158
159
160/*     
161 
162 given the distances XYZ as a triplet (on a unit grid)
163 return the maximum distance. The calling program must multiply by
164 the grid dimension to get real distance
165 
166 */
167int
168maxDistanceX(DistParamPtr p)
169{
170        double dmax,dij;                                //output dmax value, dij
171        double *xv,*yv,*zv;             //pointers to input xyz coordinates
172        CountInt i,j;
173    CountInt npt;
174        CountInt p1,p2;
175       
176        // check for all of the required waves
177        if (p->zwavH == NIL) {
178                SetNaN64(&p->result);
179                return NON_EXISTENT_WAVE;
180        }
181        if (p->ywavH == NIL) {
182                SetNaN64(&p->result);
183                return NON_EXISTENT_WAVE;
184        }
185        if (p->xwavH == NIL) {
186                SetNaN64(&p->result);
187                return NON_EXISTENT_WAVE;
188        }
189
190    //check to see that all are double
191    if(WaveType(p->zwavH) != NT_FP64 ) {
192        SetNaN64(&p->result);
193        return kExpectedNT_FP64;
194    }
195    if(WaveType(p->ywavH) != NT_FP64 ) {
196        SetNaN64(&p->result);
197        return kExpectedNT_FP64;
198    }
199    if(WaveType(p->xwavH) != NT_FP64 ) {
200        SetNaN64(&p->result);
201        return kExpectedNT_FP64;
202    }
203       
204        //
205    npt = (CountInt) WavePoints(p->xwavH);      //wavePoints returns long, number of XYZ points
206        xv = (double*)WaveData(p->xwavH);               //xyz locations
207        yv = (double*)WaveData(p->ywavH);
208        zv = (double*)WaveData(p->zwavH);
209       
210        p1 = (CountInt) p->p1;
211        p2 = (CountInt) p->p2;
212       
213        dmax = 0;
214        //do the i!=j double loop, keeping the maximum distance
215
216        for(i=p1;i<p2;i+=1) {
217                for(j=(i+1);j<npt;j+=1) {
218//                      dij=XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j]);
219                        dij = (xv[i]-xv[j])*(xv[i]-xv[j]) + (yv[i]-yv[j])*(yv[i]-yv[j]) + (zv[i]-zv[j])*(zv[i]-zv[j]);
220                        if(dij > dmax) {
221                                dmax = dij;
222                        }
223                }
224        }
225       
226    p->result = dmax;
227       
228    return 0;
229       
230}
231
232/*     
233 
234 given the distances XYZ as a triplet (on a unit grid)
235 return the binned histogram of distances
236 
237 */
238int
239binDistanceX(BinParamPtr p)
240{
241        double *xv,*yv,*zv,*bv;         //pointers to input xyz coordinates
242        CountInt i,j;
243    CountInt npt,numBins,binIndex;
244        double grid,binWidth,val;
245        CountInt p1,p2;
246       
247       
248       
249        // check for all of the required waves
250        if (p->bwavH == NIL) {
251                SetNaN64(&p->result);
252                return NON_EXISTENT_WAVE;
253        }
254        if (p->zwavH == NIL) {
255                SetNaN64(&p->result);
256                return NON_EXISTENT_WAVE;
257        }
258        if (p->ywavH == NIL) {
259                SetNaN64(&p->result);
260                return NON_EXISTENT_WAVE;
261        }
262        if (p->xwavH == NIL) {
263                SetNaN64(&p->result);
264                return NON_EXISTENT_WAVE;
265        }
266       
267    //check to see that all are double
268        if(WaveType(p->bwavH) != NT_FP64 ) {
269        SetNaN64(&p->result);
270        return kExpectedNT_FP64;
271    }
272    if(WaveType(p->zwavH) != NT_FP64 ) {
273        SetNaN64(&p->result);
274        return kExpectedNT_FP64;
275    }
276    if(WaveType(p->ywavH) != NT_FP64 ) {
277        SetNaN64(&p->result);
278        return kExpectedNT_FP64;
279    }
280    if(WaveType(p->xwavH) != NT_FP64 ) {
281        SetNaN64(&p->result);
282        return kExpectedNT_FP64;
283    }
284       
285        //
286    npt = (CountInt) WavePoints(p->xwavH);      //wavePoints returns long, number of XYZ points
287    numBins = (CountInt) WavePoints(p->bwavH);  //wavePoints returns long, number of points in bin wave
288        p1 = (CountInt) p->p1;
289        p2 = (CountInt) p->p2;
290       
291       
292        xv = (double*)WaveData(p->xwavH);               //xyz locations
293        yv = (double*)WaveData(p->ywavH);
294        zv = (double*)WaveData(p->zwavH);
295        bv = (double*)WaveData(p->bwavH);
296       
297        grid = p->grid;
298        binWidth = p->binWidth;
299       
300        //do the i!=j double loop,     
301        for(i=p1;i<p2;i+=1) {
302                for(j=(i+1);j<npt;j+=1) {
303                        val = XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j])*grid;
304                        binIndex = (CountInt)(val/binWidth-0.5);
305                        if(binIndex > numBins -1 ) {
306                                //Print "bad index"
307                        } else {
308                                bv[binIndex] += 1;
309                        }
310                       
311                }
312        }
313       
314    p->result = 0;
315       
316    return 0;
317       
318}
319
320
321/*     
322 
323 given the distances XYZ as a triplet (on a unit grid) and the SLD at each point,
324 return the binned histogram of distances for each of the parwise interactions
325
326 The returned binning is a matrix, and has to be assigned as such
327 
328 */
329int
330binSLDDistanceX(BinSLDParamPtr p)
331{
332        double *xv,*yv,*zv;             //pointers to input xyz coordinates
333        double *rho,*SLDLook,*PSFid;    // rho and the SLD lookup vector
334        CountInt i,j;
335    CountInt npt,numBins,binIndex;
336        double grid,binWidth,val,retVal;
337        CountInt p1,p2;
338        CountInt intSLD;
339
340       
341// for accessing the 2D wave data to write the results 
342        waveHndl wavH,PSFwavH;
343//      long numDimensions;
344//      long dimensionSizes[MAX_DIMENSIONS+1];
345        double value[2];                                // Pointers used for double data.
346        IndexInt indices[MAX_DIMENSIONS];
347//     
348        CountInt rhoi,rhoj,rii,rji,PSFIndex;
349       
350       
351        // check for all of the required waves
352        if (p->PSFidH == NIL) {
353                SetNaN64(&p->result);
354                return NON_EXISTENT_WAVE;
355        }
356        if (p->SLDLookH == NIL) {
357                SetNaN64(&p->result);
358                return NON_EXISTENT_WAVE;
359        }
360        if (p->bwavH == NIL) {
361                SetNaN64(&p->result);
362                return NON_EXISTENT_WAVE;
363        }
364        if (p->rhowavH == NIL) {
365                SetNaN64(&p->result);
366                return NON_EXISTENT_WAVE;
367        }
368        if (p->zwavH == NIL) {
369                SetNaN64(&p->result);
370                return NON_EXISTENT_WAVE;
371        }
372        if (p->ywavH == NIL) {
373                SetNaN64(&p->result);
374                return NON_EXISTENT_WAVE;
375        }
376        if (p->xwavH == NIL) {
377                SetNaN64(&p->result);
378                return NON_EXISTENT_WAVE;
379        }
380       
381    //check to see that all are double
382        if(WaveType(p->PSFidH) != NT_FP64 ) {
383        SetNaN64(&p->result);
384        return kExpectedNT_FP64;
385    }
386        if(WaveType(p->SLDLookH) != NT_FP64 ) {
387        SetNaN64(&p->result);
388        return kExpectedNT_FP64;
389    }
390        if(WaveType(p->bwavH) != NT_FP64 ) {
391        SetNaN64(&p->result);
392        return kExpectedNT_FP64;
393    }
394        if(WaveType(p->rhowavH) != NT_FP64 ) {
395        SetNaN64(&p->result);
396        return kExpectedNT_FP64;
397    }
398    if(WaveType(p->zwavH) != NT_FP64 ) {
399        SetNaN64(&p->result);
400        return kExpectedNT_FP64;
401    }
402    if(WaveType(p->ywavH) != NT_FP64 ) {
403        SetNaN64(&p->result);
404        return kExpectedNT_FP64;
405    }
406    if(WaveType(p->xwavH) != NT_FP64 ) {
407        SetNaN64(&p->result);
408        return kExpectedNT_FP64;
409    }
410       
411       
412        // access the 2D wave data for writing using the direct method
413        wavH = p->bwavH;
414        if (wavH == NIL)
415                return NOWAV;
416        //
417        PSFwavH = p->PSFidH;
418       
419    npt = (CountInt) WavePoints(p->xwavH);      //wavePoints returns long, number of XYZ points
420    numBins = (CountInt) WavePoints(p->bwavH);  //wavePoints returns long, number of points in bin wave
421       
422        xv = (double*)WaveData(p->xwavH);               //xyz locations
423        yv = (double*)WaveData(p->ywavH);
424        zv = (double*)WaveData(p->zwavH);
425        rho = (double*)WaveData(p->rhowavH);
426        SLDLook = (double*)WaveData(p->SLDLookH);
427        PSFid = (double*)WaveData(p->PSFidH);                   //this one is 2D
428       
429        p1 = (CountInt) p->p1;
430        p2 = (CountInt) p->p2;
431       
432        intSLD = (CountInt) p->minSLD;          //convert to int for use as index
433
434        grid = p->grid;
435        binWidth = p->binWidth;
436       
437        //do the i!=j double loop,     
438        for(i=p1;i<p2;i+=1) {
439                for(j=(i+1);j<npt;j+=1) {
440                        val = XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j])*grid;
441                        binIndex = (CountInt)(val/binWidth-0.5);
442                        if(binIndex > numBins -1 ) {
443                                //Print "bad index"
444                        } else {
445                                rhoi = (CountInt) rho[i];                               //get the rho value at i and j
446                                rhoj = (CountInt) rho[j];
447                                rii = (CountInt) SLDLook[rhoi+intSLD];                  //rho i index
448                                rji = (CountInt) SLDLook[rhoj+intSLD];                  //rho j index
449                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
450                                indices[0] = (IndexInt)rii;
451                                indices[1] = (IndexInt)rji;
452                                if (retVal = MDGetNumericWavePointValue(PSFwavH, indices, value))
453                                        return retVal;
454                                //PSFIndex = (long) PSFid[rii][rji];            //doesn't work
455                                PSFIndex = (CountInt) value[0];
456                               
457                                //now do the assignment to the 2D
458                                // equivalent to binMatrix[binIndex][PSFIndex]
459                               
460                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
461                                indices[0] = (IndexInt)binIndex;
462                                indices[1] = (IndexInt)PSFIndex;
463                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
464                                        return retVal;
465                                value[0] += 1; // Real part
466                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
467                                        return retVal;
468                               
469                        }
470                       
471                }
472        }
473       
474    p->result = 0;
475       
476    return 0;
477       
478}
479
480
481///// this is directly from Numerical Recipes
482// -- I did change the float to double, since Igor treats all as double
483// and n is an int, not a pointer (seemed unnecessary)
484//
485#define MAXBIT 30
486#define MAXDIM 6
487static int iminarg1,iminarg2;
488#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ? (iminarg1) : (iminarg2))
489
490int
491SobolX(SobolParamPtr p)
492{
493        int j,k,l;
494        unsigned long i,im,ipp;
495        static double fac;
496        static unsigned long in,ix[MAXDIM+1],*iu[MAXBIT+1];
497        static unsigned long mdeg[MAXDIM+1]={0,1,2,3,3,4,4};
498        static unsigned long ip[MAXDIM+1]={0,0,1,1,2,1,4};
499        static unsigned long iv[MAXDIM*MAXBIT+1]={
500                0,1,1,1,1,1,1,3,1,3,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9};
501       
502        static int initDone=0;
503        char buf[256];
504
505        int n=0;
506        double *x;              //output x vector
507       
508        // check for all of the required waves
509        if (p->bwavH == NIL) {
510                SetNaN64(&p->result);
511                return NON_EXISTENT_WAVE;
512        }
513       
514    //check to see that all are double
515        if(WaveType(p->bwavH) != NT_FP64 ) {
516        SetNaN64(&p->result);
517        return kExpectedNT_FP64;
518    }
519        x = (double*)WaveData(p->bwavH);
520        n = (int)(p->nIn);                      // not sure that the negative input will be properly cast to int
521       
522//      sprintf(buf, "input, recast n = %g  %d\r",p->nIn, n);
523//      XOPNotice(buf);
524       
525        if (n < 0) {
526               
527                if(initDone) {
528                        sprintf(buf, "Don't re-initialize\r");
529                        XOPNotice(buf);
530                        return 0;
531                }
532               
533                for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM) iu[j] = &iv[k];
534                for (k=1;k<=MAXDIM;k++) {
535                        for (j=1;j<=mdeg[k];j++) iu[j][k] <<= (MAXBIT-j);
536                        for (j=mdeg[k]+1;j<=MAXBIT;j++) {
537                                ipp=ip[k];
538                                i=iu[j-mdeg[k]][k];
539                                i ^= (i >> mdeg[k]);
540                                for (l=mdeg[k]-1;l>=1;l--) {
541                                        if (ipp & 1) i ^= iu[j-l][k];
542                                        ipp >>= 1;
543                                }
544                                iu[j][k]=i;
545                        }
546                }
547                fac=1.0/(1L << MAXBIT);
548                in=0;
549               
550                initDone=1;
551               
552                sprintf(buf, "Initialization loop done\r");
553                XOPNotice(buf);
554               
555        } else {
556                im=in;
557                for (j=1;j<=MAXBIT;j++) {
558                        if (!(im & 1)) break;
559                        im >>= 1;
560                }
561                if (j > MAXBIT) {
562                        sprintf(buf, "MAXBIT too small in sobseq\r");
563                        XOPNotice(buf);
564                }
565                im=(j-1)*MAXDIM;
566                for (k=1;k<=IMIN(n,MAXDIM);k++) {
567                        ix[k] ^= iv[im+k];
568                        x[k-1]=ix[k]*fac;                       /// this is a real array to send back, count this one from zero
569                        //sprintf(buf, "calculate x[%d] = %g\r",k,ix[k]*fac);
570                        //XOPNotice(buf);
571                }
572                in++;
573        }
574       
575//      sprintf(buf, "x[0],x[1] = %g  %g\r",x[0],x[1]);
576//      XOPNotice(buf);
577
578       
579        p->result = 0;
580       
581    return 0;
582}
583
584#undef MAXBIT
585#undef MAXDIM
586
587
588//#pragma XOP_RESET_STRUCT_PACKING                      // All structures are 2-byte-aligned.
Note: See TracBrowser for help on using the repository browser.