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

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


/// 1D scattering function
double sphere_fast_analytical_1D(SimSphereFParameters *pars, double q) {
        int j, npts, volume_points;
        double sum, phi, vol;
        
        return simcylinder_fast_analytical_2D(pars, q, 0.0);
        
        //return simcylinder_fast_analytical_1D_average(pars, q);
        
        sum = 0;
        npts = 21;
        volume_points = (int)floor(pars->npoints);
        
        for(j=0; j<npts; j++) {
                phi = acos(-1.0)/npts * j;
                //sum += simcylinder_simple_analytical_2D(pars, q, phi, volume_points);
                sum += simcylinder_fast_analytical_2D(pars, q, phi);
        }

        // Calculate I(q,phi) and return that value
        vol = 4.0/3.0*acos(-1.0)*pars->radius*pars->radius*pars->radius;
        sum = 1.0e8/volume_points/volume_points*vol*sum*acos(-1.0)/2;
        return sum/npts;
}


/// 1D scattering function
double sphere_fast_analytical_2D(SimSphereFParameters *pars, double q, double phi) {
        // Check if Rho array is available
        int volume_points;
        int r_points;
        int ptsGenerated;
        double bin_width;
        double r_step;
        double vol;
        double retval;
        
        int i,j;
        int ix_bin, iy_bin;
        SpacePoint p1, p2;
        
        clock_t start;
        clock_t finish;
        double cos_term;
        double sin_term;
        double qx, qy, qz;
        double phase;
        double cyl_x, cyl_y, cyl_z;
        double q_x, q_y, q_z;
        double cos_val, alpha;
        
        
        volume_points = (int)floor(pars->npoints);
        cos_term = 0.0;
        sin_term = 0.0;
        
        qx = q*cos(phi);
        qy = q*sin(phi);
        qz = 0.0;
                
        // Generate random points accross the volume
        if(pars->calcPars.isPointMemAllocated_2D==0) {
                pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
                if(pars->calcPars.points_2D==NULL) {
                        printf("Problem allocating memory for 2D volume points\n");
                        return -1.0;
                }
                
            ptsGenerated = sphere_fast_generatePoints(pars->calcPars.points_2D, 
                        volume_points, pars->radius, (int)floor(pars->seed));
                
                pars->calcPars.isPointMemAllocated_2D=1;
        }
                  
                        
        for(i=0;i<volume_points-1;i++) {
                
                //p1 = fast_rotate(pars->calcPars.points_2D[i], pars->theta, pars->phi, 0.0);
                p1 = pars->calcPars.points_2D[i];
                phase = qx*p1.x + qy*p1.y;
                cos_term += cos(phase);
                sin_term += sin(phase);
                    
        }
        
        // Calculate I(q,phi) and return that value
        vol = 4.0/3.0*acos(-1.0)*pars->radius*pars->radius*pars->radius;
        
    return 1.0e8/volume_points/volume_points*vol*(cos_term*cos_term + sin_term*sin_term);       
    //return 1.0e8/volume_points/volume_points*vol*(cos_term*cos_term); 
}


 


/**
 * Generate points randomly accross the volume
 * @param points [SpacePoint*] Array of 3D points to be filled
 * @param n [int] Number of points to generat
 * @param radius [double] Radius of the sphere
 * @return Number of points generated
 */
int sphere_fast_generatePoints(SpacePoint * points, int n, double radius, int seed) {
        int i;
        int testcounter;
        double x, y, z;
        
        // Create points        
        // To have a uniform density, you want to generate
        // random points in a box and keep only those that are
        // within the volume.
        
        // Initialize random number generator
        //int seed;    
        time_t now;
        
        time(&now);
        //seed = 10000;
        
        //seed = (int)floor(fmod(now,10000));
        //seed = 10009; 
        srand(seed);    
        //printf("Seed = %i\n", seed);
        
        testcounter = 0;
                
        memset(points,0,n*sizeof(SpacePoint));
        for(i=0;i<n;i++) {
                // Generate in a box centered around zero
                x = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
                y = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
                z = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius;
                
                // reject those that are not within the volume
                if( sqrt(x*x+y*y+z*z) < radius ) {
                        points[i].x =  x;
                        points[i].y =  y;
                        points[i].z =  z;
                        testcounter++;
                } else {
                        i--;
                }
        }
        
        // Consistency check
        if(testcounter != n) {
                return -1;
        }
                
        return testcounter;
}