source: sasview/sansmodels/prototypes/src/smeared_cylinder.c @ 7df1a50

/**
 * Scattering model for a cylinder
 * @author: Mathieu Doucet / UTK
 */

#include "cylinder.h"
#include "smeared_cylinder.h"
#include <math.h>
#include "libCylinder.h"


/**
 * Function to evaluate 1D scattering function
 * @param pars: parameters of the cylinder
 * @param q: q-value
 * @return: function value
 */
double smeared_cylinder_analytical_1D(SmearCylinderParameters *pars, double q) {
        double dp[5];
        int i_r;
        double r_0, r, step_r, min_r;
        int npts;
        double weight, func, norm;
        
        // Fill paramater array
        dp[0] = pars->scale;
        dp[1] = pars->radius;
        dp[2] = pars->length;
        dp[3] = pars->contrast;
        dp[4] = pars->background;
        
        if(pars->sigma_radius==0) {
                return CylinderForm(dp, q);
        }
        
        // Central value is the current value 
        r_0     = pars->radius;
        
        npts = 100;
        step_r     = 4.0*pars->sigma_radius/npts;
        min_r     = r_0 - 2.0*pars->sigma_radius;
        
        
        norm = 0.0;
        func = 0.0;
        
        for (i_r=0; i_r<100; i_r++) {
                        r = min_r + step_r*i_r;
                        dp[1] = r;
                        
                        // Weigth for that position
                        weight = smeared_cylinder_dist( r, r_0, pars->sigma_radius );
                        norm += weight;
                        
                        // Evaluate I(q) at that r-value
                        func += weight * CylinderForm(dp, q);   
        }
        
        return func/norm;
}
        
        
/**
 * Function to evaluate 2D scattering function
 * @param pars: parameters of the cylinder
 * @param q: q-value
 * @return: function value
 */
double smeared_cylinder_analytical_2D(SmearCylinderParameters *pars, double q, double phi) {
        CylinderParameters cyl_pars;
        int i_theta, i_phi, i_r;
        int n_theta, n_phi, n_r;
        double theta_0, phi_0, r_0;
        int npts;
        double weight_theta, weight_phi, weight_r;
        double min_theta, min_phi, min_r;
        double step_theta, step_phi, step_r;
        double func, norm;
        double n_width = 3.0;
        
        // Fill paramater struct
        cyl_pars.scale      = pars->scale;
        cyl_pars.radius     = pars->radius;
        cyl_pars.length     = pars->length;
        cyl_pars.contrast   = pars->contrast;
        cyl_pars.background = pars->background;
        cyl_pars.cyl_phi    = pars->cyl_phi;
        cyl_pars.cyl_theta  = pars->cyl_theta;
        
        theta_0 = pars->cyl_theta;
        phi_0   = pars->cyl_phi;
        r_0     = pars->radius;
        
        npts = 25;
        step_theta = 2.0*n_width*pars->sigma_theta/npts;
        step_phi   = 2.0*n_width*pars->sigma_phi/npts;
        step_r     = 2.0*n_width*pars->sigma_radius/npts;
        
        if (step_theta>0) {
                n_theta = npts;
        } else {
                n_theta = 1;
        }
        
        if (step_phi>0) {
                n_phi = npts;
        } else {
                n_phi = 1;
        }
        
        if (step_r>0) {
                n_r = npts;
        } else {
                n_r = 1;
        }
        
        
        
        min_theta = theta_0 - n_width*pars->sigma_theta;
        min_phi   = phi_0 - n_width*pars->sigma_phi;
        min_r     = r_0 - n_width*pars->sigma_radius;
        
        func = 0.0;
        norm = 0.0;
        
        for (i_theta=0; i_theta<n_theta; i_theta++) {
                        
                // Weight for that position
                if(pars->sigma_theta>0) {
                        cyl_pars.cyl_theta = min_theta + step_theta*i_theta;
                        weight_theta = smeared_cylinder_dist( cyl_pars.cyl_theta, theta_0, pars->sigma_theta );
                } else {
                        weight_theta = 1.0;
                }
                                                
                for (i_phi=0; i_phi<n_phi; i_phi++) {
                                
                        // Weight for that position
                        if(pars->sigma_phi>0) {
                                cyl_pars.cyl_phi = min_phi + step_phi*i_phi;
                                weight_phi = smeared_cylinder_dist( cyl_pars.cyl_phi, phi_0, pars->sigma_phi );
                        } else {
                                weight_phi = 1.0;
                        }
                
                        for (i_r=0; i_r<n_r; i_r++) {
                                
                                if(pars->sigma_radius>0) {
                                        cyl_pars.radius = min_r + step_r*i_r;
                                        
                                        // Weight for that position
                                        weight_r = smeared_cylinder_dist( cyl_pars.radius, r_0, pars->sigma_radius );
                                } else {
                                        weight_r = 1.0;
                                }
                                        
                                // Evaluate I(q) at that r-value
                                func += weight_theta * weight_r * weight_phi * cylinder_analytical_2D(&cyl_pars, q, phi);       
                                norm += weight_theta * weight_r * weight_phi;
                        }
                }
        }
        

        return func/norm;
}
    
double smeared_cylinder_dist( double x, double mean, double sigma ) {
        double vary;
        double expo;
                
        //return 1.0;
                
        vary = x-mean;
    expo = -vary*vary/(2.0*sigma*sigma);
        return exp(expo);

}