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

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

added point ranges to maxDistance functions to allow threading by breaking up the outer for loop.

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