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

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

double test_analytical_2D(SimCylinderParameters *pars, double q, double phi) {
        double cyl_x, cyl_y, cyl_z;
        double q_x, q_y, q_z, lenq;
        double theta, alpha, f, vol, sin_val, cos_val;
        double answer;
        
    // Cylinder orientation
    cyl_x = sin(pars->theta) * cos(pars->phi);
    cyl_y = sin(pars->theta) * sin(pars->phi);
    cyl_z = cos(pars->theta);
     
    // q orientation
    
    q_x = cos(phi);
    q_y = sin(phi);
    q_z = 0;
    
    // Length of q vector
    //lenq = sqrt(q_x*q_x + q_y*q_y + q_z*q_z);
    
    // Normalize unit vector q
    //q_x = q_x/lenq;
    //q_y = q_y/lenq;
    //q_z = q_z/lenq;
    
    // Compute the angle btw vector q and the
    // axis of the cylinder
    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
    
    // The following test should always pass
    if (fabs(cos_val)>1.0) {
        printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
        return 0;
    }
    
        alpha = acos( cos_val );
        
        // Call the IGOR library function to get the kernel
        answer = acos(-1.0)/2.0 * CylKernel(q, pars->radius, pars->length/2.0, alpha)/sin(alpha);
        
        //normalize by cylinder volume
        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
    vol = acos(-1.0) * pars->radius * pars->radius * pars->length;
        answer *= vol;
        
        //Scale
        answer *= pars->scale;
        
        //convert to [cm-1]
        //answer *= 1.0e8;
        
        return answer;
}


/// 1D scattering function
double simcylinder_analytical_1D(SimCylinderParameters *pars, double q) {
        /***
         * Things to keep here:
         *      - volume calc
         *  - point generation
         *  
         */
        
        // Check if Rho array is available
        int volume_points;
        int r_points;
        int ptsGenerated;
        
        double r_step;
        double vol;
        double bin_width;
        double twopiq;
        double tmp;
                
        //printf("returning %g %g\n",test_analytical_2D(pars, q, 0.0), pars->length);
        return test_analytical_2D(pars, q, 0.0);
        
        
        // These should be parameters
        vol = 1.0*acos(-1.0)*pars->radius*pars->radius*pars->length;
        
        
        r_points      = pars->calcPars.r_points;
        volume_points = pars->calcPars.volume_points;
        
        if(pars->calcPars.isPointMemAllocated==0) {
                
                // Call modelCalc function here to init_volume
                
                //twopiq        = (2*acos(-1.0)*pars->qmax);
                //tmp           = twopiq*twopiq*twopiq*vol;
                //volume_points = (int) floor(tmp);
                //r_points = (int) floor(10.0*pow(tmp,0.3333));
                //r_points      = (int) floor(pow(tmp,0.3333));
                //printf("v, r = %d, %d\n",volume_points, r_points); 
                
                // TEST
                volume_points = 1000;
                r_points = 1000;
                
                pars->calcPars.volume_points = volume_points;
                pars->calcPars.r_points      = r_points;
                
                // Memory allocation
                pars->calcPars.points = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint));
                if(pars->calcPars.points==NULL) {
                        printf("Problem allocating memory for 1D volume points\n");
                        return -1.0;
                }
                pars->calcPars.isPointMemAllocated=1;
        }

    r_step = sqrt(pars->radius*pars->radius*4.0+pars->length*pars->length)/(double)r_points;
        
        //printf("step=%g  r=%g  l=%g\n", r_step, pars->radius, pars->length);
        if(pars->calcPars.isRhoAvailable==0) {
                // Generate random points accross the volume
                ptsGenerated = simcylinder_generatePoints(pars->calcPars.points, volume_points, pars->radius, pars->length);
                
                // Consistency check
                if(ptsGenerated <= 0) {
                        // Set error code here
                        return 0;
                }
                
                // Allocate memory
                
                pars->calcPars.rho = (double*) malloc(r_points*sizeof(double));
                if(pars->calcPars.rho==NULL){
                        printf("Problem allocating memory for 1D correlation points\n");
                        return 0;
                }
                
                                
                if(modelcalculations_calculatePairCorrelation_1D(pars->calcPars.points, volume_points, pars->calcPars.rho, r_points, r_step)<=0){
                        printf("Error occured!\n");
                        return 0;
                };
                pars->calcPars.isRhoAvailable=1;
        } 
        
        // Calculate I(q,phi) and return that value
        
        return acos(-1.0)/2.0*modelcalculations_calculateIq_1D(pars->calcPars.rho, r_points, r_step, q) * vol;

}

/// 1D scattering function
double simcylinder_analytical_2D(SimCylinderParameters *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;
        
        // These should be parameters
        //r_points      = pars->calcPars.r_points;
        //volume_points = pars->calcPars.volume_points;
        volume_points = 50000;
        r_points = 100;
        if(pars->calcPars.isPointMemAllocated_2D==0) {
                //volume_points = 1000;
                //r_points = 100;
                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;
                }
                pars->calcPars.isPointMemAllocated_2D=1;
        }
                  
    r_step = sqrt(pars->radius*pars->radius*4.0+pars->length*pars->length)/(double)r_points;
        if(pars->calcPars.isRhoAvailable_2D==0) {
            // Initialize random number generator 
                int seed;    
                seed = 10000;
                srand(seed);
                
                // Generate random points accross the volume
                ptsGenerated = simcylinder_generatePoints(pars->calcPars.points_2D, volume_points, pars->radius, pars->length);
                
                // Calculate correlation function
                pars->calcPars.rho_2D = (float*) malloc(r_points*r_points*sizeof(float));
                if(pars->calcPars.rho_2D==NULL){
                        printf("Problem allocating memory for 2D correlations points\n");
                        return -1.0;
                }
                //if(modelcalculations_calculatePairCorrelation_2D(pars->calcPars.points_2D, volume_points, pars->calcPars.rho_2D, r_points, r_step)==NULL){
                if(modelcalculations_calculatePairCorrelation_2D_vector(pars->calcPars.points_2D, 
                        volume_points, pars->calcPars.rho_2D, r_points, r_step,
                        pars->theta, pars->phi,0.0)==0){
                        return 0;
                };
                pars->calcPars.isRhoAvailable_2D=1;
        } 
        // Calculate I(q,phi) and return that value
        vol = acos(-1.0)*pars->radius*pars->radius*pars->length;
        
        //printf("in ana_2D %f %f\n",q, phi);
        retval = modelcalculations_calculateIq_2D(pars->calcPars.rho_2D, r_points, r_step, q, phi)*vol; 
        //printf("I=%g %f %f\n",retval, q, pars->theta);
        return acos(-1.0)/2.0*retval;
}

/// 1D scattering function
double simcylinder_analytical_2D_3Darray(SimCylinderParameters *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;
        
        // These should be parameters
        volume_points = 5000;
        r_points = 100;
        
        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;
                }
                pars->calcPars.isPointMemAllocated_2D=1;
        }
                  
    r_step = sqrt(pars->radius*pars->radius*4.0+pars->length*pars->length)/(double)r_points;
        
        if(pars->calcPars.isRhoAvailable_2D==0) {
            // Initialize random number generator 
                int seed;    
                seed = 10000;
                srand(seed);
                 
                // Generate random points accross the volume
                ptsGenerated = simcylinder_generatePoints(pars->calcPars.points_2D, volume_points, pars->radius, pars->length);
                
                // Calculate correlation function
                // TODO: Check memory leaks
                printf("A\n");
                pars->calcPars.rho_2D = (float*) malloc(r_points*r_points*r_points*sizeof(float));
                if(pars->calcPars.rho_2D==NULL){
                        printf("Problem allocating memory for 2D correlations points\n");
                        return -1.0;
                }
                printf("B\n");
                if(modelcalculations_calculatePairCorrelation_2D_3Darray(pars->calcPars.points_2D, 
                        volume_points, pars->calcPars.rho_2D, r_points, r_step)<0){
                        return 0;
                };
                printf("C\n");
                pars->calcPars.isRhoAvailable_2D=1;
        } 
        // Calculate I(q,phi) and return that value
        vol = acos(-1.0)*pars->radius*pars->radius*pars->length;
        
        printf("in ana_2D %f %f\n",q, phi);
        return modelcalculations_calculateIq_2D_3Darray(pars->calcPars.rho_2D, r_points, r_step, q, phi)*vol;
        
        
}


/**
 * 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 simcylinder_generatePoints(SpacePoint * points, int n, double radius, double length) {
        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) * length/2.0;
                
                // reject those that are not within the volume
                if( sqrt(x*x+y*y) < radius && fabs(z)<length/2.0) {
                        points[i].x =  x;
                        points[i].y =  y;
                        points[i].z =  z;
                        testcounter++;
                } else {
                        i--;
                }
        }
        //printf("test counter = %d\n", testcounter);
        
        // Consistency check
        if(testcounter != n) {
                return -1;
        }
                
        return testcounter;
}