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

#include "modelCalculations.h"
#include "canvas.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <memory.h>
#include <time.h>

/**
 * Initialization function for simulation structure
 */
int canvas_init(CanvasParams *pars) {
        int i;
        
        pars->n_shapes = 0;
        pars->max_shapes = 3; 
        
        // Allocate memory
        pars->shapes = (SpaceObject*) malloc(pars->max_shapes*sizeof(SpaceObject));
        
        if(pars->shapes==NULL) {
                printf("canvas_init: Problem allocating memory for space objects\n");
                return -1;
        }
        
        return 0;
}

int canvas_dealloc(CanvasParams *self) {
        int i;
        printf("dealloc\n");
        /*
        //for(i=0; i<self->n_shapes; i++) {
        for(i=0; i<self->max_shapes; i++) {
                free(self->shapes[i].points);
                free(self->shapes[i].params);
                //free(self->shapes[i].generate);
        }
        */
        free(self->shapes);
        printf("deallocated\n");
        return 0;
}

int canvas_add(CanvasParams *self, int objectCode) {
        int id;
        id = self->n_shapes;
        
        if(objectCode==REALSPACE_SPHERE) {
                initSphereObject( &(self->shapes[id]) );
        }
        
        self->shapes[id].layer = id;
        self->n_shapes++;
        
        return id;
}

int canvas_setParam(CanvasParams *self, int shapeID, int paramID, double value) {
        self->shapes[shapeID].params[paramID]=value;
        self->shapes[shapeID].points_available = 0;
        
        // update volume
        self->shapes[shapeID].volume = (*self->shapes[shapeID].getVolume)(&(self->shapes[shapeID].params));
        
        return 0;
}

double canvas_intensity(CanvasParams *self, double q, double phi) {
        double phase, cos_term, sin_term;
        double qx, qy, vol;
        int i, j, npts;
        double pars[4];
        
        double cos_shape, sin_shape;
        
        double vol_core, vol_shell, vol_sum;
        
        pars[0] = 1.0;
        pars[1] = 40.0;
        
        cos_term = 0.0;
        sin_term = 0.0;
        
        vol_sum = 0.0;
        
        vol = self->shapes[0].volume;
        vol_core = self->shapes[1].volume;
        vol_shell = vol-vol_core;

        npts = 0;
        
        qx = q*cos(phi);
        qy = q*sin(phi);
        
        for(i=0; i<self->n_shapes; i++) {
                if (self->shapes[i].points_available==0) {
                        (*self->shapes[i].generate)( self->shapes[i].points, self->shapes[i].params, self->shapes[i].npts);
                        self->shapes[i].points_available=1;
                }
                npts = 0;
                cos_shape = 0;
                sin_shape = 0;
                for(j=0; j<self->shapes[i].npts; j++) {
                        
                        if(canvas_PointAllowed(self, i, j)==1) {
                                npts++;
                                
                                phase = qx*self->shapes[i].points[j].x + qy*self->shapes[i].points[j].y;
                                //phase = 0.0;
                                
                                cos_shape += cos(phase);
                                sin_shape += sin(phase);
                        }                       
                }

                //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]);


                
                
                // TODO: must find an efficient way to get the exact volume of each SLD regions.
                // For instance, if a volume is completely included in another (and it has a higher layer
                // number, we can simply subtract the two "theoretical" volumes instead of using the 
                // ratio of the points.

                // The following are the three lines we would like to write:
                
                //cos_term += cos_shape * self->shapes[i].volume /npts  * self->shapes[i].params[2];
                //sin_term += sin_shape * self->shapes[i].volume /npts  * self->shapes[i].params[2];
                //vol_sum += self->shapes[i].volume;
                
                // In that case, the volume is approx
                //     vol = (full volume) - sum(volume parts overlapping with other objects)
                // That can easily be done for object completely contained by others.

                cos_term += cos_shape * self->shapes[i].volume /self->shapes[i].npts  * self->shapes[i].params[2];
                sin_term += sin_shape * self->shapes[i].volume /self->shapes[i].npts  * self->shapes[i].params[2];
                vol_sum += self->shapes[i].volume*npts/self->shapes[i].npts;
                
                
                //vol_sum += self->shapes[i].volume;

                // The following is more precise (but implies better investigation of the topology)
                /*
                if(i==0){
                        cos_term += cos_shape * vol_shell/npts  * self->shapes[i].params[2];
                        sin_term += sin_shape * vol_shell/npts  * self->shapes[i].params[2];
                        vol_sum += vol_shell;
                } else {
                        cos_term += cos_shape * vol_core/npts  * self->shapes[i].params[2];
                        sin_term += sin_shape * vol_core/npts  * self->shapes[i].params[2];
                        vol_sum += vol_core;
                }
                */
        }
        
        return 1.0e8*(cos_term*cos_term + sin_term*sin_term)/vol_sum;   
} 

int canvas_PointAllowed(CanvasParams *self, int object_id, int i_pt) {
        // Check whether space point i_pt is already covered by another object
        int i;
        
        for(i=0; i<self->n_shapes; i++) {
                // Check the object object_id is underneath
                if(self->shapes[object_id].layer < self->shapes[i].layer) {
                        // Is the point overlapping, if so skip it.
                        if(self->shapes[i].isContained(&self->shapes[object_id].points[i_pt], 
                                self->shapes[i].params) == 1) {
                                return 0;
                        }
                } 
        }
        return 1;
} 



// SPHERE object ----------------------------------------------------------------------


int initSphereObject(SpaceObject *pars) {
        
        pars->x = 0.0;
        pars->y = 0.0;
        pars->z = 0.0;
        
        pars->o1 = 0.0;
        pars->o2 = 0.0;
        pars->o3 = 0.0;
        
        /*
        free(pars->params);
        
        printf("allocating parameters\n");
        pars->params = (double*) malloc(5*sizeof(double));
        if(pars->params==NULL) {
                printf("initSphereObject: Problem allocating memory\n");
                return -1;      
        }
        */

        // Scale
        pars->params[0] = 1.0;
        // Radius
        pars->params[1] = 40.0;
        // Contrast
        pars->params[2] = 1.0;
         
        pars->points_available = 0;
        
        pars->npts = 150000;
        free(pars->points);
        
        // Generation function
        pars->generate = &generateSphereObject;
        // isContained function
        pars->isContained = &sphere_IsContained;
        pars->getVolume = &sphere_getVolume;
        
        pars->volume = acos(-1.0)*4.0/3.0*pars->params[1]*pars->params[1]*pars->params[1];
        
        pars->layer = 0;
        
        return 0;
        
}

int sphere_IsContained(SpacePoint *point, double *params) {
        if( sqrt(point->x*point->x + point->y*point->y + point->z*point->z) < params[1]) {
          return 1;
        }
        return 0;
}

double sphere_getVolume(double *params) {
        printf("radius=%g\n",params[1]);
        return acos(-1.0)*4.0/3.0*params[1]*params[1]*params[1];
}


int generateSphereObject(SpacePoint *points, double *params, int npts) {
        int i, testcounter;
        double x, y, z;
        SpacePoint tmp;
        
        printf("radius=%g, npts=%i\n", params[1], npts);
        
        testcounter = 0;
        
        /*
        free(points);
        printf("allocating points\n");
        points = (SpacePoint*) malloc(npts*sizeof(SpacePoint));
        if(points==NULL) {
                printf("generateSphereObject: Problem allocating memory\n");
                return -1;      
        }
        */
        memset(points,0,npts*sizeof(SpacePoint));

        for(i=0;i<npts;i++) {
                // Generate in a box centered around zero
                x = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * params[1];
                y = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * params[1];
                z = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * params[1];
                
                // reject those that are not within the volume
//              if( sqrt(x*x+y*y+z*z) < params[1]) {
                tmp.x = x;
                tmp.y = y;
                tmp.z = z;
                if( sphere_IsContained(&tmp, params) == 1 ) {
                        points[i].x =  x;
                        points[i].y =  y;
                        points[i].z =  z;
                        testcounter++;
                } else {
                        i--;
                }
        }
        
        // Consistency check
        if(testcounter != npts) {
                return -1;
        }
        
        return testcounter;
        
}