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 |
---|
12 | HOST_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 | */ |
---|
30 | int |
---|
31 | DebyeSpheresX(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 | |
---|
139 | double 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 | |
---|
148 | double 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 | #pragma XOP_RESET_STRUCT_PACKING // All structures are 2-byte-aligned. |
---|
161 | // All structures are 2-byte-aligned. |
---|