source: sasview/sansmodels/prototypes/src/modelCalculations.c @ d8f3cd8

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

moving a file

  • Property mode set to 100644
File size: 11.3 KB
Line 
1#include "modelCalculations.h"
2#include <stdio.h>
3#include <stdlib.h>
4#include <math.h>
5#include <memory.h>
6#include <time.h>
7
8/**
9 * Initialization function for simulation structure
10 */
11void modelcalculations_init(CalcParameters *pars) {
12        pars->isRhoAvailable         = 0;
13        pars->isPointMemAllocated    = 0;
14        pars->isRhoAvailable_2D      = 0;
15        pars->isPointMemAllocated_2D = 0;
16       
17        pars->volume_points = 0;
18        pars->r_points = 0;
19       
20                pars->timePr_1D = 0;
21        pars->timePr_2D = 0;
22        pars->timeIq_1D = 0;
23        pars->timeIq_2D = 0;
24   
25        pars->errorOccured = 0;
26       
27}
28
29/**
30 * Reset function for simulation structure
31 */
32void modelcalculations_reset(CalcParameters *pars) {
33       modelcalculations_dealloc(pars);
34       modelcalculations_init(pars);
35}
36
37/**
38 * Deallocate memory of simullation structure
39 */
40void modelcalculations_dealloc(CalcParameters *pars) {
41        free(pars->rho);
42        free(pars->points);
43        free(pars->rho_2D);
44        free(pars->points_2D);
45}
46
47/**
48 * Calculate pair correlation for 1D simulation
49 */
50int modelcalculations_calculatePairCorrelation_1D(SpacePoint * points, int volume_points, double * rho, int r_points, double bin_width) {
51        int i,j;
52        double dx,dy,dz,dist;
53        int i_bin;
54        //double bin_width;
55        double delta_t;
56        double average;
57        double closest;
58        time_t start_time;
59        clock_t start;
60        clock_t finish;
61        //struct tm *timeStruct;
62       
63       
64        // Allocate memory
65        /*
66        rho = (double*) malloc(r_points*sizeof(double));
67        if(rho==NULL){
68                printf("Problem allocating memory for 1D correlation points\n");
69                return -1;
70        }
71        */
72       
73        // Clear vector
74        memset(rho,0,r_points*sizeof(double));
75
76        // R bin width
77    //bin_width = 2.0*size/r_points;
78   
79                time(&start_time);
80                start = clock();
81               
82               
83        average = 0;
84        for(i=0;i<volume_points-1;i++) {
85                closest = -1;
86                for(j=i+1;j<volume_points;j++) {
87            dx = (points[i].x-points[j].x);
88            dy = (points[i].y-points[j].y);
89            dz = (points[i].z-points[j].z);
90                        dist = sqrt(dx*dx + dy*dy + dz*dz);
91                       
92                        if(closest<0 || dist<closest) {
93                                closest = dist;
94                        }
95                       
96                        //i_bin = (int)dist/bin_width;
97                        i_bin = (int)floor(dist/bin_width);
98           
99            //rho[i_bin] = rho[i_bin] + 1.0/9000000.0;
100            if(i_bin >= r_points) {
101                printf("problem! %i > %i\n", i_bin, r_points);
102            } else {
103                rho[i_bin] = rho[i_bin] + 1.0;
104            }
105                }
106                average += closest;
107        }
108        average = average/(double)volume_points;
109        printf("average distance %f\n",average);
110       
111                finish = clock();
112        //time(&end_time);
113        //delta_t = difftime(end_time,start_time);
114        delta_t = ((double)(finish-start))/CLOCKS_PER_SEC;
115        printf("------------->PR calc time = %f\n", delta_t);
116    return 1;
117}
118
119/**
120 * Calculate I(q) for 1D simulation
121 */
122double modelcalculations_calculateIq_1D(double * rho, int r_points, double r_step, double q) {
123        int i;
124        double value;
125        //double r_step;
126        //double vol;
127        double sum;
128        double qr;
129        clock_t start;
130        clock_t finish;
131        double delta_t;
132       
133        start = clock();
134       
135    //vol = 4.0*acos(-1.0)/3.0*radius*radius*radius;
136        //r_step = 2.0*radius/((double)(r_points));
137       
138        value = 0.0;
139        sum = 0.0;
140        for(i=1; i<r_points; i++) {
141                qr = q*r_step*(double)i;
142                value = value + rho[i] * sin(qr) / qr;
143                sum = sum + rho[i];
144        }
145       
146       
147        value = value/sum;
148       
149        finish = clock();
150    delta_t = ((double)(finish-start))/CLOCKS_PER_SEC;
151    //printf("------------->IQ calc time = %f\n", delta_t);
152       
153       
154        return value;
155}
156
157/**
158 * Calculate pair correlation function for 2D simulation using
159 * a 3D array to store P(r)
160 */
161int modelcalculations_calculatePairCorrelation_2D_3Darray(SpacePoint * points, int volume_points, float * rho, int r_points, double bin_width) {
162        int i,j;
163        int ix_bin, iy_bin, iz_bin;
164       
165        clock_t start;
166        clock_t finish;
167       
168        // Clear vector
169        memset(rho,0,r_points*r_points*r_points*sizeof(float));
170
171    start = clock();
172        for(i=0;i<volume_points-1;i++) {
173                for(j=i+1;j<volume_points;j++) {
174                       
175            // Add entry to the matrix
176            ix_bin = (int)floor(fabs(points[i].x-points[j].x)/bin_width);
177            iy_bin = (int)floor(fabs(points[i].y-points[j].y)/bin_width);
178            iz_bin = (int)floor(fabs(points[i].z-points[j].z)/bin_width);
179           
180            if(ix_bin < r_points && iy_bin < r_points && iz_bin < r_points) {
181                rho[(ix_bin*r_points+iy_bin)*r_points+iz_bin] += 1.0;
182            } else {
183                printf("Bad point! %i %i %i\n", ix_bin, iy_bin, iz_bin);
184            }
185                       
186                }
187        }
188       
189        finish = clock();
190    printf("-------------> Pair Correlation time = %f\n", ((double)(finish-start))/CLOCKS_PER_SEC);
191    return 0;
192}
193
194/**
195 * Calculate pair correlation function for 2D simulation by storing P(r) in a 2D array.
196 * Allows for rotation of object in space by specifying theta, phi, omage of the beam
197 */
198int modelcalculations_calculatePairCorrelation_2D_vector(SpacePoint * points, int volume_points, float * rho, 
199                        int r_points, double bin_width, double theta_beam, double phi_beam, double omega_beam) {
200        int i,j;
201        int ix_bin, iy_bin;
202        SpacePoint p1, p2;
203       
204        clock_t start;
205        clock_t finish;
206       
207        // Clear vector
208        memset(rho,0,r_points*r_points*sizeof(float));
209        //printf("P(r) with theta=%g phi=%g\n", theta_beam, phi_beam);
210    start = clock();
211        for(i=0;i<volume_points-1;i++) {
212                // Rotate point
213                p1 = modelcalculations_rotate(points[i], theta_beam, phi_beam, omega_beam);
214                //printf("p = %g %g %g -> %g %g %g\n", points[i].x,points[i].y,points[i].z, p1.x,p1.y,p1.z);
215                for(j=i+1;j<volume_points;j++) {
216                        // Rotate point
217                        p2 = modelcalculations_rotate(points[j], theta_beam, phi_beam, omega_beam);
218                       
219                       
220                        // Calculate distance in plane perpendicular to beam (z)
221            ix_bin = (int)floor(fabs(p1.x-p2.x)/bin_width);
222            iy_bin = (int)floor(fabs(p1.y-p2.y)/bin_width);
223            //iz_bin = (int)floor((points[i].z-points[j].z)/bin_width+r_points/2);
224           
225            rho[ix_bin*r_points+iy_bin] += 1.0;
226                       
227                }
228        }
229       
230        finish = clock();
231    printf("-------------> 2D (v) Pair Correlation time = %f\n", ((double)(finish-start))/CLOCKS_PER_SEC);
232    return 1;
233}
234
235/**
236 * Calculate pair correlation function for 2D simulation by storing P(r) in a 2D array
237 */
238int modelcalculations_calculatePairCorrelation_2D(SpacePoint * points, int volume_points, float * rho, int r_points, double bin_width) {
239        int i,j;
240        int ix_bin, iy_bin;
241       
242        clock_t start;
243        clock_t finish;
244       
245        // Clear vector
246        memset(rho,0,r_points*r_points*sizeof(float));
247
248        //return 1;
249       
250    start = clock();
251        for(i=0;i<volume_points-1;i++) {
252                for(j=i+1;j<volume_points;j++) {
253                       
254                        // Calculate distance in plane perpendicular to beam (z)
255            ix_bin = (int)floor(fabs(points[i].x-points[j].x)/bin_width);
256            iy_bin = (int)floor(fabs(points[i].y-points[j].y)/bin_width);
257            //iz_bin = (int)floor(fabs(points[i].z-points[j].z)/bin_width);
258           
259            rho[ix_bin*r_points+iy_bin] += 1.0;
260            //rho[ix_bin*r_points+iy_bin] += fabs(points[i].z-points[j].z);
261                }
262        }
263       
264        finish = clock();
265    printf("-------------> Pair Correlation time = %f\n", ((double)(finish-start))/CLOCKS_PER_SEC);
266   
267    /*
268     * for(i=0;i<r_points;i++) {
269        for(j=0;j<r_points;j++) {
270                printf("Pr(%i, %i) = %g\n", i, j, rho[i*r_points+j]);
271        }       
272       
273    }
274    */
275        return 0;
276}
277
278/**
279 * Calculate I(q) for 2D simulation from 3D array
280 */
281double modelcalculations_calculateIq_2D_3Darray(float * rho, int r_points, double r_step, double q, double phi) {
282        //TODO: make rho an array of ints.
283       
284        int ix,iy,iz;
285       
286        // This should be a parameter
287        double lambda = 1.6;
288        double theta;
289        double value;
290        double sum;
291        // This should also be a parameter, about value of radius
292        double r_max = 1;
293        //double r_step =1;
294       
295        double c1;
296        double c2;
297        double c3;
298        int r2;
299        double f3;
300        double iz_c;
301        clock_t start; 
302        clock_t finish;
303       
304        theta = 2*asin(q*lambda/(4*acos(-1.0)));
305       
306        value = 0.0;
307        sum = 0.0;
308        //c1 = lambda*(cos(theta)-1)*r_step;
309        //c2 = lambda*sin(theta)*cos(phi)*r_step;
310        //c3 = lambda*sin(theta)*sin(phi)*r_step;
311       
312        c1 = 0;
313        c2 = q*cos(phi)*r_step;
314        c3 = q*sin(phi)*r_step;
315       
316        start = clock();
317       
318        // TODO: sum(rho) should be equal to the number of points^2 !
319       
320        for(ix=0; ix<r_points; ix++) {
321                for(iy=0; iy<r_points; iy++) {
322                        r2 = ix*r_points*r_points+iy*r_points;
323                        f3 = c2*(ix-r_points/2) + c3*(iy-r_points/2);
324                        for(iz=0; iz<r_points; iz++) {
325                                iz_c = (double)iz-r_points/2; 
326                                value = value + rho[r2+iz] * cos( c1*iz_c + f3 );
327                                sum = sum + rho[r2+iz];
328                        }
329                }
330        }
331        finish = clock();
332    //printf("-------------> I(Q) time = %f, (%f %f %f)\n", ((double)(finish-start))/CLOCKS_PER_SEC, q, phi,value/sum);
333       
334       
335        value = value /sum;
336        return value;
337}
338
339/**
340 * Pair correlation function for a sphere
341 *      @param r: distance value
342 */
343double pair_corr_sphere(double r) {
344       
345        return r*r*(1.0 - 0.75*r + r*r*r/16.0); 
346}
347
348/**
349 * Calculate I(q) for 2D simulation from 2D array P(r)
350 */
351double modelcalculations_calculateIq_2D(float * rho, int r_points, double r_step, double q, double phi) {
352        //TODO: make rho an array of ints.
353       
354        int ix,iy;
355       
356        // This should be a parameter
357        double lambda = 1.0;
358        double value;
359        double sum;
360        // This should also be a parameter, about value of radius
361        double r_max = 1;
362        //double r_step =1;
363       
364        double c2;
365        double c3;
366        int ibin;
367        clock_t start; 
368        clock_t finish;
369       
370        //theta = 2*asin(q*lambda/(4*acos(-1.0)));
371        value = 0.0;
372        sum = 0.0;
373        c2 = q*cos(phi)*r_step;
374        c3 = q*sin(phi)*r_step;
375        //printf("phi = %g, c2=%g, c3=%g, q=%g, r_step=%g\n",phi, c2, c3,q,r_step);
376       
377        start = clock();
378       
379        for(ix=-r_points+1; ix<r_points; ix++) {
380                for(iy=-r_points+1; iy<r_points; iy++) {
381                       
382                        //value += rho[ix*r_points+iy] * cos( c2*((double)ix+0.5) + c3*((double)iy+0.5) );
383                        //sum += rho[ix*r_points+iy];
384                       
385                        ibin = (int)(floor(sqrt(1.0*ix*ix)))*r_points+(int)(floor(sqrt(1.0*iy*iy)));
386                       
387                        if (ibin<r_points*r_points) {
388                       
389                                //value += rho[ibin] * cos( c2*((double)ix+0.5) + c3*((double)iy+0.5) );
390                       
391                                value += rho[ibin] * cos( c2*((double)ix) + c3*((double)iy) );
392                               
393                                sum += rho[ibin];
394                        } else {
395                                printf("Error computing IQ %i >= %i (%i %i)\n", ibin, r_points*r_points,ix, iy);
396                        };
397                       
398                       
399                        //dx = r_step*((double)ix+0.5);
400                        //dy = r_step*((double)iy+0.5);
401                       
402                        // Only works for sphere of radius = 20
403                        //f3 = pair_corr_sphere(sqrt(dx*dx+dy*dy)/20.0);
404                       
405                        //value += f3 * cos( c2*(dx) + c3*(dy) );
406                        //sum += f3;           
407                }
408        }
409        finish = clock();
410       
411        value = value /sum;
412        return value;
413}
414
415
416/**
417 *  Rotation of a space point
418 */
419SpacePoint modelcalculations_rotate(SpacePoint p, double theta, double phi, double omega) {
420        SpacePoint new_point;
421        double x_1, x_2;
422        double y_1, y_2;
423        double z_1, z_2;
424        // P(r) assumes beam along z-axis. Rotate point accordingly
425       
426       
427        // Omega, around z-axis (doesn't change anything for cylindrical symmetry
428        x_1 = p.x*cos(omega) - p.y*sin(omega);
429        y_1 = p.x*sin(omega) + p.y*cos(omega);
430        z_1 = p.z;     
431       
432        // Theta, around y-axis
433        x_2 = x_1*cos(theta) + z_1*sin(theta);
434        y_2 = y_1;
435        z_2 = -x_1*sin(theta) + z_1*cos(theta);
436       
437        // Phi, around z-axis
438        new_point.x = x_2*cos(phi) - y_2*sin(phi);
439        new_point.y = x_2*sin(phi) + y_2*cos(phi);
440        new_point.z = z_2;
441       
442        return new_point;       
443}
444 
445 
Note: See TracBrowser for help on using the repository browser.