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

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

Adding Threadsafe declarations to Debye Sphere functions to allow binning of distances to be threaded.

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