source: sasview/sansmodels/prototypes/src/canvas.c @ dc65f8c

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

moving a file

  • Property mode set to 100644
File size: 7.2 KB
Line 
1#include "modelCalculations.h"
2#include "canvas.h"
3#include <stdio.h>
4#include <stdlib.h>
5#include <math.h>
6#include <memory.h>
7#include <time.h>
8
9/**
10 * Initialization function for simulation structure
11 */
12int canvas_init(CanvasParams *pars) {
13        int i;
14       
15        pars->n_shapes = 0;
16        pars->max_shapes = 3; 
17       
18        // Allocate memory
19        pars->shapes = (SpaceObject*) malloc(pars->max_shapes*sizeof(SpaceObject));
20       
21        if(pars->shapes==NULL) {
22                printf("canvas_init: Problem allocating memory for space objects\n");
23                return -1;
24        }
25       
26        return 0;
27}
28
29int canvas_dealloc(CanvasParams *self) {
30        int i;
31        printf("dealloc\n");
32        /*
33        //for(i=0; i<self->n_shapes; i++) {
34        for(i=0; i<self->max_shapes; i++) {
35                free(self->shapes[i].points);
36                free(self->shapes[i].params);
37                //free(self->shapes[i].generate);
38        }
39        */
40        free(self->shapes);
41        printf("deallocated\n");
42        return 0;
43}
44
45int canvas_add(CanvasParams *self, int objectCode) {
46        int id;
47        id = self->n_shapes;
48       
49        if(objectCode==REALSPACE_SPHERE) {
50                initSphereObject( &(self->shapes[id]) );
51        }
52       
53        self->shapes[id].layer = id;
54        self->n_shapes++;
55       
56        return id;
57}
58
59int canvas_setParam(CanvasParams *self, int shapeID, int paramID, double value) {
60        self->shapes[shapeID].params[paramID]=value;
61        self->shapes[shapeID].points_available = 0;
62       
63        // update volume
64        self->shapes[shapeID].volume = (*self->shapes[shapeID].getVolume)(&(self->shapes[shapeID].params));
65       
66        return 0;
67}
68
69double canvas_intensity(CanvasParams *self, double q, double phi) {
70        double phase, cos_term, sin_term;
71        double qx, qy, vol;
72        int i, j, npts;
73        double pars[4];
74       
75        double cos_shape, sin_shape;
76       
77        double vol_core, vol_shell, vol_sum;
78       
79        pars[0] = 1.0;
80        pars[1] = 40.0;
81       
82        cos_term = 0.0;
83        sin_term = 0.0;
84       
85        vol_sum = 0.0;
86       
87        vol = self->shapes[0].volume;
88        vol_core = self->shapes[1].volume;
89        vol_shell = vol-vol_core;
90
91        npts = 0;
92       
93        qx = q*cos(phi);
94        qy = q*sin(phi);
95       
96        for(i=0; i<self->n_shapes; i++) {
97                if (self->shapes[i].points_available==0) {
98                        (*self->shapes[i].generate)( self->shapes[i].points, self->shapes[i].params, self->shapes[i].npts);
99                        self->shapes[i].points_available=1;
100                }
101                npts = 0;
102                cos_shape = 0;
103                sin_shape = 0;
104                for(j=0; j<self->shapes[i].npts; j++) {
105                       
106                        if(canvas_PointAllowed(self, i, j)==1) {
107                                npts++;
108                               
109                                phase = qx*self->shapes[i].points[j].x + qy*self->shapes[i].points[j].y;
110                                //phase = 0.0;
111                               
112                                cos_shape += cos(phase);
113                                sin_shape += sin(phase);
114                        }                       
115                }
116
117                //printf("ID %i: cos=%g sin=%g fac=%g\n",i,cos_shape, sin_shape, self->shapes[i].volume /self->shapes[i].npts  * self->shapes[i].params[2]);
118
119
120               
121               
122                // TODO: must find an efficient way to get the exact volume of each SLD regions.
123                // For instance, if a volume is completely included in another (and it has a higher layer
124                // number, we can simply subtract the two "theoretical" volumes instead of using the
125                // ratio of the points.
126
127                // The following are the three lines we would like to write:
128               
129                //cos_term += cos_shape * self->shapes[i].volume /npts  * self->shapes[i].params[2];
130                //sin_term += sin_shape * self->shapes[i].volume /npts  * self->shapes[i].params[2];
131                //vol_sum += self->shapes[i].volume;
132               
133                // In that case, the volume is approx
134                //     vol = (full volume) - sum(volume parts overlapping with other objects)
135                // That can easily be done for object completely contained by others.
136
137                cos_term += cos_shape * self->shapes[i].volume /self->shapes[i].npts  * self->shapes[i].params[2];
138                sin_term += sin_shape * self->shapes[i].volume /self->shapes[i].npts  * self->shapes[i].params[2];
139                vol_sum += self->shapes[i].volume*npts/self->shapes[i].npts;
140               
141               
142                //vol_sum += self->shapes[i].volume;
143
144                // The following is more precise (but implies better investigation of the topology)
145                /*
146                if(i==0){
147                        cos_term += cos_shape * vol_shell/npts  * self->shapes[i].params[2];
148                        sin_term += sin_shape * vol_shell/npts  * self->shapes[i].params[2];
149                        vol_sum += vol_shell;
150                } else {
151                        cos_term += cos_shape * vol_core/npts  * self->shapes[i].params[2];
152                        sin_term += sin_shape * vol_core/npts  * self->shapes[i].params[2];
153                        vol_sum += vol_core;
154                }
155                */
156        }
157       
158        return 1.0e8*(cos_term*cos_term + sin_term*sin_term)/vol_sum;   
159} 
160
161int canvas_PointAllowed(CanvasParams *self, int object_id, int i_pt) {
162        // Check whether space point i_pt is already covered by another object
163        int i;
164       
165        for(i=0; i<self->n_shapes; i++) {
166                // Check the object object_id is underneath
167                if(self->shapes[object_id].layer < self->shapes[i].layer) {
168                        // Is the point overlapping, if so skip it.
169                        if(self->shapes[i].isContained(&self->shapes[object_id].points[i_pt], 
170                                self->shapes[i].params) == 1) {
171                                return 0;
172                        }
173                } 
174        }
175        return 1;
176} 
177
178
179
180// SPHERE object ----------------------------------------------------------------------
181
182
183int initSphereObject(SpaceObject *pars) {
184       
185        pars->x = 0.0;
186        pars->y = 0.0;
187        pars->z = 0.0;
188       
189        pars->o1 = 0.0;
190        pars->o2 = 0.0;
191        pars->o3 = 0.0;
192       
193        /*
194        free(pars->params);
195       
196        printf("allocating parameters\n");
197        pars->params = (double*) malloc(5*sizeof(double));
198        if(pars->params==NULL) {
199                printf("initSphereObject: Problem allocating memory\n");
200                return -1;     
201        }
202        */
203
204        // Scale
205        pars->params[0] = 1.0;
206        // Radius
207        pars->params[1] = 40.0;
208        // Contrast
209        pars->params[2] = 1.0;
210         
211        pars->points_available = 0;
212       
213        pars->npts = 150000;
214        free(pars->points);
215       
216        // Generation function
217        pars->generate = &generateSphereObject;
218        // isContained function
219        pars->isContained = &sphere_IsContained;
220        pars->getVolume = &sphere_getVolume;
221       
222        pars->volume = acos(-1.0)*4.0/3.0*pars->params[1]*pars->params[1]*pars->params[1];
223       
224        pars->layer = 0;
225       
226        return 0;
227       
228}
229
230int sphere_IsContained(SpacePoint *point, double *params) {
231        if( sqrt(point->x*point->x + point->y*point->y + point->z*point->z) < params[1]) {
232          return 1;
233        }
234        return 0;
235}
236
237double sphere_getVolume(double *params) {
238        printf("radius=%g\n",params[1]);
239        return acos(-1.0)*4.0/3.0*params[1]*params[1]*params[1];
240}
241
242
243int generateSphereObject(SpacePoint *points, double *params, int npts) {
244        int i, testcounter;
245        double x, y, z;
246        SpacePoint tmp;
247       
248        printf("radius=%g, npts=%i\n", params[1], npts);
249       
250        testcounter = 0;
251       
252        /*
253        free(points);
254        printf("allocating points\n");
255        points = (SpacePoint*) malloc(npts*sizeof(SpacePoint));
256        if(points==NULL) {
257                printf("generateSphereObject: Problem allocating memory\n");
258                return -1;     
259        }
260        */
261        memset(points,0,npts*sizeof(SpacePoint));
262
263        for(i=0;i<npts;i++) {
264                // Generate in a box centered around zero
265                x = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * params[1];
266                y = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * params[1];
267                z = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * params[1];
268               
269                // reject those that are not within the volume
270//              if( sqrt(x*x+y*y+z*z) < params[1]) {
271                tmp.x = x;
272                tmp.y = y;
273                tmp.z = z;
274                if( sphere_IsContained(&tmp, params) == 1 ) {
275                        points[i].x =  x;
276                        points[i].y =  y;
277                        points[i].z =  z;
278                        testcounter++;
279                } else {
280                        i--;
281                }
282        }
283       
284        // Consistency check
285        if(testcounter != npts) {
286                return -1;
287        }
288       
289        return testcounter;
290       
291}
Note: See TracBrowser for help on using the repository browser.