source: sans/XOP_Dev/MonteCarlo/DebyeSpheres.c @ 789

Last change on this file since 789 was 789, checked in by srkline, 12 years ago

Added XOP-ized functions for doing the distance-binned Debye Spheres calculation. Two versions are added, one that includes the SLDs and one that does not.

Also added a pseudo-random sequence generator, SobolX, that will generate pseudo-random 2D or 3D distributions. This fills space more uniformly than a random generation.

  • Property svn:executable set to *
File size: 13.2 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 void main(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        int i,j;
37    int 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 = (int) WavePoints(p->xwavH);   //wavePoints returns long, number of XYZ points
109        xv = WaveData(p->xwavH);                //xyz locations
110        yv = WaveData(p->ywavH);
111        zv = WaveData(p->zwavH);
112        rv = 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        int i,j;
173    int npt;
174       
175       
176       
177        // check for all of the required waves
178        if (p->zwavH == NIL) {
179                SetNaN64(&p->result);
180                return NON_EXISTENT_WAVE;
181        }
182        if (p->ywavH == NIL) {
183                SetNaN64(&p->result);
184                return NON_EXISTENT_WAVE;
185        }
186        if (p->xwavH == NIL) {
187                SetNaN64(&p->result);
188                return NON_EXISTENT_WAVE;
189        }
190
191    //check to see that all are double
192    if(WaveType(p->zwavH) != NT_FP64 ) {
193        SetNaN64(&p->result);
194        return kExpectedNT_FP64;
195    }
196    if(WaveType(p->ywavH) != NT_FP64 ) {
197        SetNaN64(&p->result);
198        return kExpectedNT_FP64;
199    }
200    if(WaveType(p->xwavH) != NT_FP64 ) {
201        SetNaN64(&p->result);
202        return kExpectedNT_FP64;
203    }
204       
205        //
206    npt = (int) WavePoints(p->xwavH);   //wavePoints returns long, number of XYZ points
207        xv = WaveData(p->xwavH);                //xyz locations
208        yv = WaveData(p->ywavH);
209        zv = WaveData(p->zwavH);
210       
211        dmax = 0;
212        //do the i!=j double loop, keeping the maximum distance
213
214        for(i=0;i<npt;i+=1) {
215                for(j=(i+1);j<npt;j+=1) {
216//                      dij=XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j]);
217                        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]);
218                        if(dij > dmax) {
219                                dmax = dij;
220                        }
221                }
222        }
223       
224    p->result = dmax;
225       
226    return 0;
227       
228}
229
230/*     
231 
232 given the distances XYZ as a triplet (on a unit grid)
233 return the binned histogram of distances
234 
235 */
236int
237binDistanceX(BinParamPtr p)
238{
239        double *xv,*yv,*zv,*bv;         //pointers to input xyz coordinates
240        int i,j;
241    int npt,numBins,binIndex;
242        double grid,binWidth,val;
243       
244       
245       
246        // check for all of the required waves
247        if (p->bwavH == NIL) {
248                SetNaN64(&p->result);
249                return NON_EXISTENT_WAVE;
250        }
251        if (p->zwavH == NIL) {
252                SetNaN64(&p->result);
253                return NON_EXISTENT_WAVE;
254        }
255        if (p->ywavH == NIL) {
256                SetNaN64(&p->result);
257                return NON_EXISTENT_WAVE;
258        }
259        if (p->xwavH == NIL) {
260                SetNaN64(&p->result);
261                return NON_EXISTENT_WAVE;
262        }
263       
264    //check to see that all are double
265        if(WaveType(p->bwavH) != NT_FP64 ) {
266        SetNaN64(&p->result);
267        return kExpectedNT_FP64;
268    }
269    if(WaveType(p->zwavH) != NT_FP64 ) {
270        SetNaN64(&p->result);
271        return kExpectedNT_FP64;
272    }
273    if(WaveType(p->ywavH) != NT_FP64 ) {
274        SetNaN64(&p->result);
275        return kExpectedNT_FP64;
276    }
277    if(WaveType(p->xwavH) != NT_FP64 ) {
278        SetNaN64(&p->result);
279        return kExpectedNT_FP64;
280    }
281       
282        //
283    npt = (int) WavePoints(p->xwavH);   //wavePoints returns long, number of XYZ points
284    numBins = (int) WavePoints(p->bwavH);       //wavePoints returns long, number of points in bin wave
285       
286        xv = WaveData(p->xwavH);                //xyz locations
287        yv = WaveData(p->ywavH);
288        zv = WaveData(p->zwavH);
289        bv = WaveData(p->bwavH);
290       
291        grid = p->grid;
292        binWidth = p->binWidth;
293       
294        //do the i!=j double loop,     
295        for(i=0;i<npt;i+=1) {
296                for(j=(i+1);j<npt;j+=1) {
297                        val = XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j])*grid;
298                        binIndex = (int)(val/binWidth-0.5);
299                        if(binIndex > numBins -1 ) {
300                                //Print "bad index"
301                        } else {
302                                bv[binIndex] += 1;
303                        }
304                       
305                }
306        }
307       
308    p->result = 0;
309       
310    return 0;
311       
312}
313
314
315/*     
316 
317 given the distances XYZ as a triplet (on a unit grid) and the SLD at each point,
318 return the binned histogram of distances for each of the parwise interactions
319
320 The returned binning is a matrix, and has to be assigned as such
321 
322 */
323int
324binSLDDistanceX(BinSLDParamPtr p)
325{
326        double *xv,*yv,*zv;             //pointers to input xyz coordinates
327        double *rho,*SLDLook,*PSFid;    // rho and the SLD lookup vector
328        int i,j;
329    int npt,numBins,binIndex;
330        double grid,binWidth,val,retVal;
331       
332// for accessing the 2D wave data to write the results 
333        waveHndl wavH,PSFwavH;
334//      long numDimensions;
335//      long dimensionSizes[MAX_DIMENSIONS+1];
336        double value[2];                                // Pointers used for double data.
337        long indices[MAX_DIMENSIONS];
338//     
339        long rhoi,rhoj,rii,rji,PSFIndex;
340       
341       
342        // check for all of the required waves
343        if (p->PSFidH == NIL) {
344                SetNaN64(&p->result);
345                return NON_EXISTENT_WAVE;
346        }
347        if (p->SLDLookH == NIL) {
348                SetNaN64(&p->result);
349                return NON_EXISTENT_WAVE;
350        }
351        if (p->bwavH == NIL) {
352                SetNaN64(&p->result);
353                return NON_EXISTENT_WAVE;
354        }
355        if (p->rhowavH == NIL) {
356                SetNaN64(&p->result);
357                return NON_EXISTENT_WAVE;
358        }
359        if (p->zwavH == NIL) {
360                SetNaN64(&p->result);
361                return NON_EXISTENT_WAVE;
362        }
363        if (p->ywavH == NIL) {
364                SetNaN64(&p->result);
365                return NON_EXISTENT_WAVE;
366        }
367        if (p->xwavH == NIL) {
368                SetNaN64(&p->result);
369                return NON_EXISTENT_WAVE;
370        }
371       
372    //check to see that all are double
373        if(WaveType(p->PSFidH) != NT_FP64 ) {
374        SetNaN64(&p->result);
375        return kExpectedNT_FP64;
376    }
377        if(WaveType(p->SLDLookH) != NT_FP64 ) {
378        SetNaN64(&p->result);
379        return kExpectedNT_FP64;
380    }
381        if(WaveType(p->bwavH) != NT_FP64 ) {
382        SetNaN64(&p->result);
383        return kExpectedNT_FP64;
384    }
385        if(WaveType(p->rhowavH) != NT_FP64 ) {
386        SetNaN64(&p->result);
387        return kExpectedNT_FP64;
388    }
389    if(WaveType(p->zwavH) != NT_FP64 ) {
390        SetNaN64(&p->result);
391        return kExpectedNT_FP64;
392    }
393    if(WaveType(p->ywavH) != NT_FP64 ) {
394        SetNaN64(&p->result);
395        return kExpectedNT_FP64;
396    }
397    if(WaveType(p->xwavH) != NT_FP64 ) {
398        SetNaN64(&p->result);
399        return kExpectedNT_FP64;
400    }
401       
402       
403        // access the 2D wave data for writing using the direct method
404        wavH = p->bwavH;
405        if (wavH == NIL)
406                return NOWAV;
407        //
408        PSFwavH = p->PSFidH;
409       
410    npt = (int) WavePoints(p->xwavH);   //wavePoints returns long, number of XYZ points
411    numBins = (int) WavePoints(p->bwavH);       //wavePoints returns long, number of points in bin wave
412       
413        xv = WaveData(p->xwavH);                //xyz locations
414        yv = WaveData(p->ywavH);
415        zv = WaveData(p->zwavH);
416        rho = WaveData(p->rhowavH);
417        SLDLook = WaveData(p->SLDLookH);
418        PSFid = WaveData(p->PSFidH);                    //this one is 2D
419       
420        grid = p->grid;
421        binWidth = p->binWidth;
422       
423        //do the i!=j double loop,     
424        for(i=0;i<npt;i+=1) {
425                for(j=(i+1);j<npt;j+=1) {
426                        val = XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j])*grid;
427                        binIndex = (int)(val/binWidth-0.5);
428                        if(binIndex > numBins -1 ) {
429                                //Print "bad index"
430                        } else {
431                                rhoi = (long) rho[i];                           //get the rho value at i and j
432                                rhoj = (long) rho[j];
433                                rii = (long) SLDLook[rhoi];                     //rho i index
434                                rji = (long) SLDLook[rhoj];                     //rho j index
435                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
436                                indices[0] = rii;
437                                indices[1] = rji;
438                                if (retVal = MDGetNumericWavePointValue(PSFwavH, indices, value))
439                                        return retVal;
440                                //PSFIndex = (long) PSFid[rii][rji];            //doesn't work
441                                PSFIndex = (long) value[0];
442                               
443                                //now do the assignment to the 2D
444                                // equivalent to binMatrix[binIndex][PSFIndex]
445                               
446                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
447                                indices[0] = binIndex;
448                                indices[1] = PSFIndex;
449                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
450                                        return retVal;
451                                value[0] += 1; // Real part
452                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
453                                        return retVal;
454                               
455                        }
456                       
457                }
458        }
459       
460    p->result = 0;
461       
462    return 0;
463       
464}
465
466
467///// this is directly from Numerical Recipes
468// -- I did change the float to double, since Igor treats all as double
469// and n is an int, not a pointer (seemed unnecessary)
470//
471#define MAXBIT 30
472#define MAXDIM 6
473static int iminarg1,iminarg2;
474#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ? (iminarg1) : (iminarg2))
475
476int
477SobolX(SobolParamPtr p)
478{
479        int j,k,l;
480        unsigned long i,im,ipp;
481        static double fac;
482        static unsigned long in,ix[MAXDIM+1],*iu[MAXBIT+1];
483        static unsigned long mdeg[MAXDIM+1]={0,1,2,3,3,4,4};
484        static unsigned long ip[MAXDIM+1]={0,0,1,1,2,1,4};
485        static unsigned long iv[MAXDIM*MAXBIT+1]={
486                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};
487       
488        static int initDone=0;
489        char buf[256];
490
491        int n=0;
492        double *x;              //output x vector
493       
494        // check for all of the required waves
495        if (p->bwavH == NIL) {
496                SetNaN64(&p->result);
497                return NON_EXISTENT_WAVE;
498        }
499       
500    //check to see that all are double
501        if(WaveType(p->bwavH) != NT_FP64 ) {
502        SetNaN64(&p->result);
503        return kExpectedNT_FP64;
504    }
505        x = WaveData(p->bwavH);
506        n = (int)(p->nIn);                      // not sure that the negative input will be properly cast to int
507       
508//      sprintf(buf, "input, recast n = %g  %d\r",p->nIn, n);
509//      XOPNotice(buf);
510       
511        if (n < 0) {
512               
513                if(initDone) {
514                        sprintf(buf, "Don't re-initialize\r");
515                        XOPNotice(buf);
516                        return 0;
517                }
518               
519                for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM) iu[j] = &iv[k];
520                for (k=1;k<=MAXDIM;k++) {
521                        for (j=1;j<=mdeg[k];j++) iu[j][k] <<= (MAXBIT-j);
522                        for (j=mdeg[k]+1;j<=MAXBIT;j++) {
523                                ipp=ip[k];
524                                i=iu[j-mdeg[k]][k];
525                                i ^= (i >> mdeg[k]);
526                                for (l=mdeg[k]-1;l>=1;l--) {
527                                        if (ipp & 1) i ^= iu[j-l][k];
528                                        ipp >>= 1;
529                                }
530                                iu[j][k]=i;
531                        }
532                }
533                fac=1.0/(1L << MAXBIT);
534                in=0;
535               
536                initDone=1;
537               
538                sprintf(buf, "Initialization loop done\r");
539                XOPNotice(buf);
540               
541        } else {
542                im=in;
543                for (j=1;j<=MAXBIT;j++) {
544                        if (!(im & 1)) break;
545                        im >>= 1;
546                }
547                if (j > MAXBIT) {
548                        sprintf(buf, "MAXBIT too small in sobseq\r");
549                        XOPNotice(buf);
550                }
551                im=(j-1)*MAXDIM;
552                for (k=1;k<=IMIN(n,MAXDIM);k++) {
553                        ix[k] ^= iv[im+k];
554                        x[k-1]=ix[k]*fac;                       /// this is a real array to send back, count this one from zero
555                        //sprintf(buf, "calculate x[%d] = %g\r",k,ix[k]*fac);
556                        //XOPNotice(buf);
557                }
558                in++;
559        }
560       
561//      sprintf(buf, "x[0],x[1] = %g  %g\r",x[0],x[1]);
562//      XOPNotice(buf);
563
564       
565        p->result = 0;
566       
567    return 0;
568}
569
570#undef MAXBIT
571#undef MAXDIM
572
573
574//#pragma XOP_RESET_STRUCT_PACKING                      // All structures are 2-byte-aligned.
Note: See TracBrowser for help on using the repository browser.