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

Last change on this file since 815 was 815, checked in by srkline, 11 years ago

fixed binning with different SLDs to make sure that it can handle SLDs that are negative.

  • Property svn:executable set to *
File size: 13.5 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        int 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        long indices[MAX_DIMENSIONS];
347//     
348        long 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 = (int) WavePoints(p->xwavH);   //wavePoints returns long, number of XYZ points
420    numBins = (int) WavePoints(p->bwavH);       //wavePoints returns long, number of points in bin wave
421       
422        xv = WaveData(p->xwavH);                //xyz locations
423        yv = WaveData(p->ywavH);
424        zv = WaveData(p->zwavH);
425        rho = WaveData(p->rhowavH);
426        SLDLook = WaveData(p->SLDLookH);
427        PSFid = WaveData(p->PSFidH);                    //this one is 2D
428       
429        p1 = (int) p->p1;
430        p2 = (int) p->p2;
431       
432        intSLD = (int) 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 = (int)(val/binWidth-0.5);
442                        if(binIndex > numBins -1 ) {
443                                //Print "bad index"
444                        } else {
445                                rhoi = (long) rho[i];                           //get the rho value at i and j
446                                rhoj = (long) rho[j];
447                                rii = (long) SLDLook[rhoi+intSLD];                      //rho i index
448                                rji = (long) SLDLook[rhoj+intSLD];                      //rho j index
449                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
450                                indices[0] = rii;
451                                indices[1] = rji;
452                                if (retVal = MDGetNumericWavePointValue(PSFwavH, indices, value))
453                                        return retVal;
454                                //PSFIndex = (long) PSFid[rii][rji];            //doesn't work
455                                PSFIndex = (long) 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] = binIndex;
462                                indices[1] = 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 = 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.