1 | #include "simcylinder_fast.h" |
---|
2 | #include <stdio.h> |
---|
3 | #include <stdlib.h> |
---|
4 | #include <math.h> |
---|
5 | #include <time.h> |
---|
6 | #include <memory.h> |
---|
7 | |
---|
8 | |
---|
9 | double simcylinder_fast_analytical_1D_test(SimCylinderFParameters *pars, double q) { |
---|
10 | /*** |
---|
11 | * Things to keep here: |
---|
12 | * - volume calc |
---|
13 | * - point generation |
---|
14 | * |
---|
15 | */ |
---|
16 | // Check if Rho array is available |
---|
17 | int volume_points; |
---|
18 | int ptsGenerated, npts; |
---|
19 | double vol; |
---|
20 | |
---|
21 | int i, j, k, m; |
---|
22 | SpacePoint *p1; |
---|
23 | |
---|
24 | double sum; |
---|
25 | double theta_cyl, phi; |
---|
26 | double phase; |
---|
27 | |
---|
28 | double cos_term; |
---|
29 | double sin_term; |
---|
30 | double qx, qy; |
---|
31 | double pi_step; |
---|
32 | |
---|
33 | |
---|
34 | volume_points = (int)floor(pars->npoints); |
---|
35 | |
---|
36 | |
---|
37 | // Generate random points accross the volume |
---|
38 | if(pars->calcPars.isPointMemAllocated_2D==0) { |
---|
39 | pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint)); |
---|
40 | if(pars->calcPars.points_2D==NULL) { |
---|
41 | printf("Problem allocating memory for 2D volume points\n"); |
---|
42 | return -1.0; |
---|
43 | } |
---|
44 | |
---|
45 | ptsGenerated = simcylinder_fast_generatePoints(pars->calcPars.points_2D, |
---|
46 | volume_points, pars->radius, pars->length, (int)floor(pars->seed)); |
---|
47 | |
---|
48 | pars->calcPars.isPointMemAllocated_2D=1; |
---|
49 | } |
---|
50 | |
---|
51 | // Loop over theta_cyl |
---|
52 | |
---|
53 | |
---|
54 | sum = 0; |
---|
55 | npts = 21; |
---|
56 | // Allocate temporary memory |
---|
57 | p1 = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint)); |
---|
58 | |
---|
59 | pi_step = acos(-1.0)/npts; |
---|
60 | for(i=0; i<npts; i++) { |
---|
61 | theta_cyl = pi_step * i; |
---|
62 | // Rotate the points |
---|
63 | |
---|
64 | for(k=0; k<volume_points; k++) { |
---|
65 | p1[k] = fast_rotate(pars->calcPars.points_2D[k], theta_cyl, 0.0, 0.0); |
---|
66 | } |
---|
67 | |
---|
68 | // Loop over phi_cyl |
---|
69 | // Equivalent to looping over phi (detector) |
---|
70 | for(j=0; j<npts; j++) { |
---|
71 | phi = 2*pi_step * j; |
---|
72 | |
---|
73 | cos_term = 0.0; |
---|
74 | sin_term = 0.0; |
---|
75 | |
---|
76 | qx = q*cos(phi); |
---|
77 | qy = q*sin(phi); |
---|
78 | |
---|
79 | for(m=0;m<volume_points;m++) { |
---|
80 | phase = qx*p1[m].x + qy*p1[m].y; |
---|
81 | cos_term += cos(phase); |
---|
82 | sin_term += sin(phase); |
---|
83 | } |
---|
84 | |
---|
85 | sum += sin(theta_cyl)* (cos_term*cos_term + sin_term*sin_term); |
---|
86 | } |
---|
87 | |
---|
88 | } |
---|
89 | free(p1); |
---|
90 | // Calculate I(q,phi) and return that value |
---|
91 | vol = acos(-1.0)*pars->radius*pars->radius*pars->length; |
---|
92 | sum = 1.0e8/volume_points/volume_points*vol*sum*acos(-1.0)/2; |
---|
93 | return sum/npts/npts; |
---|
94 | } |
---|
95 | |
---|
96 | /// 1D scattering function |
---|
97 | double simcylinder_fast_analytical_1D(SimCylinderFParameters *pars, double q) { |
---|
98 | int i, j, npts, volume_points; |
---|
99 | double sum, vol, psi; |
---|
100 | |
---|
101 | //return simcylinder_fast_analytical_1D_test(pars, q); |
---|
102 | |
---|
103 | |
---|
104 | sum = 0; |
---|
105 | npts = 51; |
---|
106 | volume_points = (int)floor(pars->npoints); |
---|
107 | |
---|
108 | for(i=0; i<npts; i++) { |
---|
109 | psi = 2*acos(-1.0)/npts * i; |
---|
110 | pars->phi = psi; |
---|
111 | |
---|
112 | for(j=0; j<npts; j++) { |
---|
113 | pars->theta = acos(-1.0)/npts * j; |
---|
114 | sum += sin(pars->theta)*simcylinder_simple_analytical_2D(pars, q, 0.0); |
---|
115 | } |
---|
116 | } |
---|
117 | // Calculate I(q,phi) and return that value |
---|
118 | vol = acos(-1.0)*pars->radius*pars->radius*pars->length; |
---|
119 | sum = 1.0e8/volume_points/volume_points*vol*sum*acos(-1.0)/2; |
---|
120 | return sum/npts/npts; |
---|
121 | } |
---|
122 | |
---|
123 | |
---|
124 | /// 2D scattering function |
---|
125 | double simcylinder_fast_analytical_2D(SimCylinderFParameters *pars, double q, double phi) { |
---|
126 | int volume_points; |
---|
127 | int ptsGenerated; |
---|
128 | double vol; |
---|
129 | |
---|
130 | int i; |
---|
131 | SpacePoint p1; |
---|
132 | |
---|
133 | double cos_term; |
---|
134 | double sin_term; |
---|
135 | double qx, qy, qz; |
---|
136 | double phase; |
---|
137 | |
---|
138 | volume_points = (int)floor(pars->npoints); |
---|
139 | cos_term = 0.0; |
---|
140 | sin_term = 0.0; |
---|
141 | |
---|
142 | qx = q*cos(phi); |
---|
143 | qy = q*sin(phi); |
---|
144 | qz = 0.0; |
---|
145 | |
---|
146 | // Generate random points accross the volume |
---|
147 | if(pars->calcPars.isPointMemAllocated_2D==0) { |
---|
148 | pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint)); |
---|
149 | if(pars->calcPars.points_2D==NULL) { |
---|
150 | printf("Problem allocating memory for 2D volume points\n"); |
---|
151 | return -1.0; |
---|
152 | } |
---|
153 | |
---|
154 | ptsGenerated = simcylinder_fast_generatePoints(pars->calcPars.points_2D, |
---|
155 | volume_points, pars->radius, pars->length, (int)floor(pars->seed)); |
---|
156 | |
---|
157 | for(i=0;i<volume_points;i++) { |
---|
158 | p1 = fast_rotate(pars->calcPars.points_2D[i], pars->theta, pars->phi, 0.0); |
---|
159 | pars->calcPars.points_2D[i] = p1; |
---|
160 | } |
---|
161 | |
---|
162 | pars->calcPars.isPointMemAllocated_2D=1; |
---|
163 | } |
---|
164 | |
---|
165 | |
---|
166 | for(i=0;i<volume_points-1;i++) { |
---|
167 | p1 = pars->calcPars.points_2D[i]; |
---|
168 | phase = qx*p1.x + qy*p1.y; |
---|
169 | cos_term += cos(phase); |
---|
170 | sin_term += sin(phase); |
---|
171 | |
---|
172 | } |
---|
173 | |
---|
174 | // Calculate I(q,phi) and return that value |
---|
175 | vol = acos(-1.0)*pars->radius*pars->radius*pars->length; |
---|
176 | return 1.0e8/volume_points/volume_points*vol*acos(-1)/2*(cos_term*cos_term + sin_term*sin_term); |
---|
177 | } |
---|
178 | |
---|
179 | double simcylinder_simple_analytical_2D(SimCylinderFParameters *pars, double q, double phi) { |
---|
180 | int ptsGenerated; |
---|
181 | |
---|
182 | int i; |
---|
183 | SpacePoint p1; |
---|
184 | |
---|
185 | double cos_term; |
---|
186 | double sin_term; |
---|
187 | double qx, qy, qz; |
---|
188 | double phase; |
---|
189 | int volume_points; |
---|
190 | |
---|
191 | cos_term = 0.0; |
---|
192 | sin_term = 0.0; |
---|
193 | |
---|
194 | qx = q*cos(phi); |
---|
195 | qy = q*sin(phi); |
---|
196 | qz = 0.0; |
---|
197 | |
---|
198 | volume_points = (int)floor(pars->npoints); |
---|
199 | |
---|
200 | // Generate random points accross the volume |
---|
201 | if(pars->calcPars.isPointMemAllocated_2D==0) { |
---|
202 | pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint)); |
---|
203 | if(pars->calcPars.points_2D==NULL) { |
---|
204 | printf("Problem allocating memory for 2D volume points\n"); |
---|
205 | return -1.0; |
---|
206 | } |
---|
207 | |
---|
208 | ptsGenerated = simcylinder_fast_generatePoints(pars->calcPars.points_2D, |
---|
209 | volume_points, pars->radius, pars->length, (int)floor(pars->seed)); |
---|
210 | |
---|
211 | for(i=0;i<volume_points;i++) { |
---|
212 | p1 = fast_rotate(pars->calcPars.points_2D[i], pars->theta, pars->phi, 0.0); |
---|
213 | pars->calcPars.points_2D[i] = p1; |
---|
214 | } |
---|
215 | |
---|
216 | pars->calcPars.isPointMemAllocated_2D=1; |
---|
217 | } |
---|
218 | |
---|
219 | for(i=0;i<volume_points;i++) { |
---|
220 | p1 = fast_rotate(pars->calcPars.points_2D[i], pars->theta, pars->phi, 0.0); |
---|
221 | phase = qx*p1.x + qy*p1.y; |
---|
222 | cos_term += cos(phase); |
---|
223 | sin_term += sin(phase); |
---|
224 | } |
---|
225 | |
---|
226 | return (cos_term*cos_term + sin_term*sin_term); |
---|
227 | } |
---|
228 | |
---|
229 | /** |
---|
230 | * Rotation of pair correlation function |
---|
231 | */ |
---|
232 | SpacePoint fast_rotate(SpacePoint p, double theta, double phi, double omega) { |
---|
233 | SpacePoint new_point; |
---|
234 | double x_1, x_2; |
---|
235 | double y_1, y_2; |
---|
236 | double z_1, z_2; |
---|
237 | |
---|
238 | |
---|
239 | // Omega, around z-axis (doesn't change anything for cylindrical symmetry |
---|
240 | x_1 = p.x*cos(omega) - p.y*sin(omega); |
---|
241 | y_1 = p.x*sin(omega) + p.y*cos(omega); |
---|
242 | z_1 = p.z; |
---|
243 | |
---|
244 | // Theta, around y-axis |
---|
245 | x_2 = x_1*cos(theta) + z_1*sin(theta); |
---|
246 | y_2 = y_1; |
---|
247 | z_2 = -x_1*sin(theta) + z_1*cos(theta); |
---|
248 | |
---|
249 | // Phi, around z-axis |
---|
250 | new_point.x = x_2*cos(phi) - y_2*sin(phi); |
---|
251 | new_point.y = x_2*sin(phi) + y_2*cos(phi); |
---|
252 | new_point.z = z_2; |
---|
253 | |
---|
254 | return new_point; |
---|
255 | |
---|
256 | |
---|
257 | } |
---|
258 | |
---|
259 | |
---|
260 | |
---|
261 | /** |
---|
262 | * Generate points randomly accross the volume |
---|
263 | * @param points [SpacePoint*] Array of 3D points to be filled |
---|
264 | * @param n [int] Number of points to generat |
---|
265 | * @param radius [double] Radius of the sphere |
---|
266 | * @return Number of points generated |
---|
267 | */ |
---|
268 | int simcylinder_fast_generatePoints(SpacePoint * points, int n, double radius, double length, int seed) { |
---|
269 | int i; |
---|
270 | int testcounter; |
---|
271 | double x, y, z; |
---|
272 | |
---|
273 | // Create points |
---|
274 | // To have a uniform density, you want to generate |
---|
275 | // random points in a box and keep only those that are |
---|
276 | // within the volume. |
---|
277 | |
---|
278 | // Initialize random number generator |
---|
279 | //int seed; |
---|
280 | time_t now; |
---|
281 | |
---|
282 | time(&now); |
---|
283 | //seed = 10000; |
---|
284 | |
---|
285 | //seed = (int)floor(fmod(now,10000)); |
---|
286 | //seed = 10009; |
---|
287 | srand(seed); |
---|
288 | //printf("Seed = %i\n", seed); |
---|
289 | |
---|
290 | testcounter = 0; |
---|
291 | |
---|
292 | memset(points,0,n*sizeof(SpacePoint)); |
---|
293 | for(i=0;i<n;i++) { |
---|
294 | // Generate in a box centered around zero |
---|
295 | x = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius; |
---|
296 | y = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius; |
---|
297 | z = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * length/2.0; |
---|
298 | |
---|
299 | // reject those that are not within the volume |
---|
300 | if( sqrt(x*x+y*y) < radius && fabs(z)<length/2.0) { |
---|
301 | points[i].x = x; |
---|
302 | points[i].y = y; |
---|
303 | points[i].z = z; |
---|
304 | testcounter++; |
---|
305 | } else { |
---|
306 | i--; |
---|
307 | } |
---|
308 | } |
---|
309 | |
---|
310 | // Consistency check |
---|
311 | if(testcounter != n) { |
---|
312 | return -1; |
---|
313 | } |
---|
314 | |
---|
315 | return testcounter; |
---|
316 | } |
---|