source: sasview/sansmodels/prototypes/src/simcylinder.c @ ea1012b

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 ea1012b was 7df1a50, checked in by Jae Cho <jhjcho@…>, 13 years ago

moving a file
  • Property mode set to 100644
File size: 9.6 KB
Line 
1#include "simcylinder.h"
2#include <stdio.h>
3#include <stdlib.h>
4#include <math.h>
5#include <time.h>
6#include <memory.h>
7#include "modelCalculations.h"
8#include "libCylinder.h"
9
10double test_analytical_2D(SimCylinderParameters *pars, double q, double phi) {
11        double cyl_x, cyl_y, cyl_z;
12        double q_x, q_y, q_z, lenq;
13        double theta, alpha, f, vol, sin_val, cos_val;
14        double answer;
15       
16    // Cylinder orientation
17    cyl_x = sin(pars->theta) * cos(pars->phi);
18    cyl_y = sin(pars->theta) * sin(pars->phi);
19    cyl_z = cos(pars->theta);
20     
21    // q orientation
22   
23    q_x = cos(phi);
24    q_y = sin(phi);
25    q_z = 0;
26   
27    // Length of q vector
28    //lenq = sqrt(q_x*q_x + q_y*q_y + q_z*q_z);
29   
30    // Normalize unit vector q
31    //q_x = q_x/lenq;
32    //q_y = q_y/lenq;
33    //q_z = q_z/lenq;
34   
35    // Compute the angle btw vector q and the
36    // axis of the cylinder
37    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
38   
39    // The following test should always pass
40    if (fabs(cos_val)>1.0) {
41        printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
42        return 0;
43    }
44   
45        alpha = acos( cos_val );
46       
47        // Call the IGOR library function to get the kernel
48        answer = acos(-1.0)/2.0 * CylKernel(q, pars->radius, pars->length/2.0, alpha)/sin(alpha);
49       
50        //normalize by cylinder volume
51        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
52    vol = acos(-1.0) * pars->radius * pars->radius * pars->length;
53        answer *= vol;
54       
55        //Scale
56        answer *= pars->scale;
57       
58        //convert to [cm-1]
59        //answer *= 1.0e8;
60       
61        return answer;
62}
63
64
65/// 1D scattering function
66double simcylinder_analytical_1D(SimCylinderParameters *pars, double q) {
67        /***
68         * Things to keep here:
69         *      - volume calc
70         *  - point generation
71         * 
72         */
73       
74        // Check if Rho array is available
75        int volume_points;
76        int r_points;
77        int ptsGenerated;
78       
79        double r_step;
80        double vol;
81        double bin_width;
82        double twopiq;
83        double tmp;
84               
85        //printf("returning %g %g\n",test_analytical_2D(pars, q, 0.0), pars->length);
86        return test_analytical_2D(pars, q, 0.0);
87       
88       
89        // These should be parameters
90        vol = 1.0*acos(-1.0)*pars->radius*pars->radius*pars->length;
91       
92       
93        r_points      = pars->calcPars.r_points;
94        volume_points = pars->calcPars.volume_points;
95       
96        if(pars->calcPars.isPointMemAllocated==0) {
97               
98                // Call modelCalc function here to init_volume
99               
100                //twopiq        = (2*acos(-1.0)*pars->qmax);
101                //tmp           = twopiq*twopiq*twopiq*vol;
102                //volume_points = (int) floor(tmp);
103                //r_points = (int) floor(10.0*pow(tmp,0.3333));
104                //r_points      = (int) floor(pow(tmp,0.3333));
105                //printf("v, r = %d, %d\n",volume_points, r_points);
106               
107                // TEST
108                volume_points = 1000;
109                r_points = 1000;
110               
111                pars->calcPars.volume_points = volume_points;
112                pars->calcPars.r_points      = r_points;
113               
114                // Memory allocation
115                pars->calcPars.points = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
116                if(pars->calcPars.points==NULL) {
117                        printf("Problem allocating memory for 1D volume points\n");
118                        return -1.0;
119                }
120                pars->calcPars.isPointMemAllocated=1;
121        }
122
123    r_step = sqrt(pars->radius*pars->radius*4.0+pars->length*pars->length)/(double)r_points;
124       
125        //printf("step=%g  r=%g  l=%g\n", r_step, pars->radius, pars->length);
126        if(pars->calcPars.isRhoAvailable==0) {
127                // Generate random points accross the volume
128                ptsGenerated = simcylinder_generatePoints(pars->calcPars.points, volume_points, pars->radius, pars->length);
129               
130                // Consistency check
131                if(ptsGenerated <= 0) {
132                        // Set error code here
133                        return 0;
134                }
135               
136                // Allocate memory
137               
138                pars->calcPars.rho = (double*) malloc(r_points*sizeof(double));
139                if(pars->calcPars.rho==NULL){
140                        printf("Problem allocating memory for 1D correlation points\n");
141                        return 0;
142                }
143               
144                               
145                if(modelcalculations_calculatePairCorrelation_1D(pars->calcPars.points, volume_points, pars->calcPars.rho, r_points, r_step)<=0){
146                        printf("Error occured!\n");
147                        return 0;
148                };
149                pars->calcPars.isRhoAvailable=1;
150        } 
151       
152        // Calculate I(q,phi) and return that value
153       
154        return acos(-1.0)/2.0*modelcalculations_calculateIq_1D(pars->calcPars.rho, r_points, r_step, q) * vol;
155
156}
157
158/// 1D scattering function
159double simcylinder_analytical_2D(SimCylinderParameters *pars, double q, double phi) {
160        // Check if Rho array is available
161        int volume_points;
162        int r_points;
163        int ptsGenerated;
164        double bin_width;
165        double r_step;
166        double vol;
167        double retval;
168       
169        // These should be parameters
170        //r_points      = pars->calcPars.r_points;
171        //volume_points = pars->calcPars.volume_points;
172        volume_points = 50000;
173        r_points = 100;
174        if(pars->calcPars.isPointMemAllocated_2D==0) {
175                //volume_points = 1000;
176                //r_points = 100;
177                pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
178                if(pars->calcPars.points_2D==NULL) {
179                        printf("Problem allocating memory for 2D volume points\n");
180                        return -1.0;
181                }
182                pars->calcPars.isPointMemAllocated_2D=1;
183        }
184                 
185    r_step = sqrt(pars->radius*pars->radius*4.0+pars->length*pars->length)/(double)r_points;
186        if(pars->calcPars.isRhoAvailable_2D==0) {
187            // Initialize random number generator
188                int seed;   
189                seed = 10000;
190                srand(seed);
191               
192                // Generate random points accross the volume
193                ptsGenerated = simcylinder_generatePoints(pars->calcPars.points_2D, volume_points, pars->radius, pars->length);
194               
195                // Calculate correlation function
196                pars->calcPars.rho_2D = (float*) malloc(r_points*r_points*sizeof(float));
197                if(pars->calcPars.rho_2D==NULL){
198                        printf("Problem allocating memory for 2D correlations points\n");
199                        return -1.0;
200                }
201                //if(modelcalculations_calculatePairCorrelation_2D(pars->calcPars.points_2D, volume_points, pars->calcPars.rho_2D, r_points, r_step)==NULL){
202                if(modelcalculations_calculatePairCorrelation_2D_vector(pars->calcPars.points_2D, 
203                        volume_points, pars->calcPars.rho_2D, r_points, r_step,
204                        pars->theta, pars->phi,0.0)==0){
205                        return 0;
206                };
207                pars->calcPars.isRhoAvailable_2D=1;
208        } 
209        // Calculate I(q,phi) and return that value
210        vol = acos(-1.0)*pars->radius*pars->radius*pars->length;
211       
212        //printf("in ana_2D %f %f\n",q, phi);
213        retval = modelcalculations_calculateIq_2D(pars->calcPars.rho_2D, r_points, r_step, q, phi)*vol; 
214        //printf("I=%g %f %f\n",retval, q, pars->theta);
215        return acos(-1.0)/2.0*retval;
216}
217
218/// 1D scattering function
219double simcylinder_analytical_2D_3Darray(SimCylinderParameters *pars, double q, double phi) {
220        // Check if Rho array is available
221        int volume_points;
222        int r_points;
223        int ptsGenerated;
224        double bin_width;
225        double r_step;
226        double vol;
227       
228        // These should be parameters
229        volume_points = 5000;
230        r_points = 100;
231       
232        if(pars->calcPars.isPointMemAllocated_2D==0) {
233                pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
234                if(pars->calcPars.points_2D==NULL) {
235                        printf("Problem allocating memory for 2D volume points\n");
236                        return -1.0;
237                }
238                pars->calcPars.isPointMemAllocated_2D=1;
239        }
240                 
241    r_step = sqrt(pars->radius*pars->radius*4.0+pars->length*pars->length)/(double)r_points;
242       
243        if(pars->calcPars.isRhoAvailable_2D==0) {
244            // Initialize random number generator
245                int seed;   
246                seed = 10000;
247                srand(seed);
248                 
249                // Generate random points accross the volume
250                ptsGenerated = simcylinder_generatePoints(pars->calcPars.points_2D, volume_points, pars->radius, pars->length);
251               
252                // Calculate correlation function
253                // TODO: Check memory leaks
254                printf("A\n");
255                pars->calcPars.rho_2D = (float*) malloc(r_points*r_points*r_points*sizeof(float));
256                if(pars->calcPars.rho_2D==NULL){
257                        printf("Problem allocating memory for 2D correlations points\n");
258                        return -1.0;
259                }
260                printf("B\n");
261                if(modelcalculations_calculatePairCorrelation_2D_3Darray(pars->calcPars.points_2D, 
262                        volume_points, pars->calcPars.rho_2D, r_points, r_step)<0){
263                        return 0;
264                };
265                printf("C\n");
266                pars->calcPars.isRhoAvailable_2D=1;
267        } 
268        // Calculate I(q,phi) and return that value
269        vol = acos(-1.0)*pars->radius*pars->radius*pars->length;
270       
271        printf("in ana_2D %f %f\n",q, phi);
272        return modelcalculations_calculateIq_2D_3Darray(pars->calcPars.rho_2D, r_points, r_step, q, phi)*vol;
273       
274       
275}
276
277
278/**
279 * Generate points randomly accross the volume
280 * @param points [SpacePoint*] Array of 3D points to be filled
281 * @param n [int] Number of points to generat
282 * @param radius [double] Radius of the sphere
283 * @return Number of points generated
284 */
285int simcylinder_generatePoints(SpacePoint * points, int n, double radius, double length) {
286        int i;
287        int testcounter;
288        double x, y, z;
289       
290        // Create points       
291        // To have a uniform density, you want to generate
292        // random points in a box and keep only those that are
293        // within the volume.
294       
295        // Initialize random number generator
296        int seed;   
297        time_t now;
298       
299        time(&now);
300        //seed = 10000;
301       
302        seed = (int)floor(fmod(now,10000));
303        //seed = 10009;
304        srand(seed);   
305        printf("Seed = %i\n", seed);
306       
307        testcounter = 0;
308               
309        memset(points,0,n*sizeof(SpacePoint));
310        for(i=0;i<n;i++) {
311                // Generate in a box centered around zero
312                x = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
313                y = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
314                z = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * length/2.0;
315               
316                // reject those that are not within the volume
317                if( sqrt(x*x+y*y) < radius && fabs(z)<length/2.0) {
318                        points[i].x =  x;
319                        points[i].y =  y;
320                        points[i].z =  z;
321                        testcounter++;
322                } else {
323                        i--;
324                }
325        }
326        //printf("test counter = %d\n", testcounter);
327       
328        // Consistency check
329        if(testcounter != n) {
330                return -1;
331        }
332               
333        return testcounter;
334}
Note: See TracBrowser for help on using the repository browser.