source: sasview/sansmodels/prototypes/src/simcylinder_fast.c @ 5548954

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 5548954 was 7df1a50, checked in by Jae Cho <jhjcho@…>, 13 years ago

moving a file
  • Property mode set to 100644
File size: 7.7 KB
Line 
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
9double 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
97double 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
125double 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
179double 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 */
232SpacePoint 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 */
268int 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}
Note: See TracBrowser for help on using the repository browser.