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

Last change on this file since 623 was 623, checked in by srkline, 13 years ago

Updated the MonteCarlo? code to allow 4 processors, but simply copying the function 4 times, and defining 4 different random number generators. Still can't figure out what the problem is with threading a single version, but not worth the effort. Copy/paste is way faster.

Also added some simple (non-optimized) calculations for using Debye's sphere method. These are largely undocumented at this point - so see the code. These are XOP versions of the old ipf code I've used in the past, and stripped of the now-obsolete AltiVec? code (I now lose the 4x speedup from the vectorization...)

  • Property svn:executable set to *
File size: 3.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        Warning:
25                The call to WaveData() below returns a pointer to the middle
26                of an unlocked Macintosh handle. In the unlikely event that your
27                calculations could cause memory to move, you should copy the coefficient
28                values to local variables or an array before such operations.
29*/
30int
31DebyeSpheresX(AltiParamsPtr p)
32{
33        float qv;                               // input q-value
34        float ival;                             //output intensity value
35        float *xv,*yv,*zv,*rv;          //pointers to input xyz-rho coordinates
36        int i,j;
37    int npt;
38        float 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        if(WaveType(p->rhowavH) != NT_FP32 ) {
63        SetNaN64(&p->result);
64        return kExpectedNT_FP32;
65    }
66    if(WaveType(p->zwavH) != NT_FP32 ) {
67        SetNaN64(&p->result);
68        return kExpectedNT_FP32;
69    }
70    if(WaveType(p->ywavH) != NT_FP32 ) {
71        SetNaN64(&p->result);
72        return kExpectedNT_FP32;
73    }
74    if(WaveType(p->xwavH) != NT_FP32 ) {
75        SetNaN64(&p->result);
76        return kExpectedNT_FP32;
77    }
78
79       
80        // convert the input to float. Do all calculations as float.
81
82    rval = (float) p->Rprimary;         // primary sphere radius
83    grid = (float) p->grid;                     // calling program should set this to 0.62*Rprimary
84        qv = (float) p->qval;
85
86//
87    npt = (int) WavePoints(p->xwavH);   //wavePoints returns long, number of XYZ points
88        xv = WaveData(p->xwavH);                //xyz locations
89        yv = WaveData(p->ywavH);
90        zv = WaveData(p->zwavH);
91        rv = WaveData(p->rhowavH);
92   
93       
94        vol = 4.0*3.1415927/3.0*rval*rval*rval;
95        ival = 0.0;
96        fQR = PhiQR(qv,rval);
97        //do the i=j sum
98        for(i=0;i<npt;i+=1) {
99                dum = rv[i]*vol*fQR;
100                ival += dum*dum;
101        }
102        //do the i!=j double sum
103        dum = vol*vol*fQR*fQR;
104        for(i=0;i<npt;i+=1) {
105                for(j=(i+1);j<npt;j+=1) {
106                        dij=XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j]) * grid;
107                        ival += 2.0*rv[i]*rv[j]*dum*sin(dij*qv)/dij/qv;
108                }
109        }
110
111    p->result= (DOUBLE) ival;
112       
113    return 0;
114
115}
116
117
118float PhiQR(float qval, float rval)
119{
120        float retval,qr;
121       
122        qr = qval*rval;
123        retval = (float) ( 3*(sin(qr) - qr*cos(qr))/qr/qr/qr );
124        return(retval);
125}
126
127float XYZDistance(float x1, float x2,float y1, float y2,float z1, float z2)
128{
129        float retval,dx,dy,dz;
130       
131        dx = (x1-x2);
132        dy = (y1-y2);
133        dz = (z1-z2);
134        retval = (float) sqrt( dx*dx + dy*dy + dz*dz );
135        return(retval);
136}
137
138
139#pragma XOP_RESET_STRUCT_PACKING                        // All structures are 2-byte-aligned.
140// All structures are 2-byte-aligned.
Note: See TracBrowser for help on using the repository browser.