source: sasview/sansmodels/src/sans/models/c_models/pearlnecklace.cpp @ a24f530

#include <math.h>
#include "models.hh"
#include "parameters.hh"
#include <stdio.h>
using namespace std;

extern "C" {
        #include "pearlnecklace.h"
}

PearlNecklaceModel :: PearlNecklaceModel() {
        scale = Parameter(1.0);
        radius = Parameter(80.0, true);
        radius.set_min(0.0);
        edge_separation = Parameter(350.0, true);
        edge_separation.set_min(0.0);
        thick_string = Parameter(2.5, true);
        thick_string.set_min(0.0);
        num_pearls = Parameter(3);
        num_pearls.set_min(0.0);
        sld_pearl = Parameter(1.0e-06);
        sld_string = Parameter(1.0e-06);
        sld_solv = Parameter(6.3e-06);
    background = Parameter(0.0);

}

/**
 * Function to evaluate 1D PearlNecklaceModel function
 * @param q: q-value
 * @return: function value
 */
double PearlNecklaceModel :: operator()(double q) {
        double dp[9];
        // Fill parameter array for IGOR library
        // Add the background after averaging
        dp[0] = scale();
        dp[1] = radius();
        dp[2] = edge_separation();
        dp[3] = thick_string();
        dp[4] = num_pearls();
        dp[5] = sld_pearl();
        dp[6] = sld_string();
        dp[7] = sld_solv();
        dp[8] = 0.0;
        double pi = 4.0*atan(1.0);
        // No polydispersion supported in this model.
        // Get the dispersion points for the radius
        vector<WeightPoint> weights_radius;
        radius.get_weights(weights_radius);
        vector<WeightPoint> weights_edge_separation;
        edge_separation.get_weights(weights_edge_separation);
        // Perform the computation, with all weight points
        double sum = 0.0;
        double norm = 0.0;
        double vol = 0.0;
        double string_vol = 0.0;
        double pearl_vol = 0.0;
        double tot_vol = 0.0;
        // Loop over core weight points
        for(size_t i=0; i<weights_radius.size(); i++) {
                dp[1] = weights_radius[i].value;
                // Loop over thick_inter0 weight points
                for(size_t j=0; j<weights_edge_separation.size(); j++) {
                        dp[2] = weights_edge_separation[j].value;
                        pearl_vol = 4.0 /3.0 * pi * pow(dp[1], 3);
                        string_vol =dp[2] * pi * pow((dp[3] / 2.0), 2);
                        tot_vol = (dp[4] - 1.0) * string_vol;
                        tot_vol += dp[4] * pearl_vol;
                        //Un-normalize Sphere by volume
                        sum += weights_radius[i].weight * weights_edge_separation[j].weight
                                * pearl_necklace_kernel(dp,q) * tot_vol;
                        //Find average volume
                        vol += weights_radius[i].weight * weights_edge_separation[j].weight
                                * tot_vol;
                        norm += weights_radius[i].weight * weights_edge_separation[j].weight;
                }
        }

        if (vol != 0.0 && norm != 0.0) {
                //Re-normalize by avg volume
                sum = sum/(vol/norm);}

        return sum/norm + background();
}

/**
 * Function to evaluate 2D PearlNecklaceModel function
 * @param q_x: value of Q along x
 * @param q_y: value of Q along y
 * @return: function value
 */
double PearlNecklaceModel :: operator()(double qx, double qy) {
        double q = sqrt(qx*qx + qy*qy);
        return (*this).operator()(q);
}

/**
 * Function to evaluate PearlNecklaceModel function
 * @param pars: parameters of the PearlNecklaceModel
 * @param q: q-value
 * @param phi: angle phi
 * @return: function value
 */
double PearlNecklaceModel :: evaluate_rphi(double q, double phi) {
        return (*this).operator()(q);
}

/**
 * Function to calculate TOTAL radius
 * Todo: decide whether or not we keep this calculation
 * @return: effective radius value
 */
// No polydispersion supported in this model.
// Calculate max radius assumming max_radius = effective radius
// Note that this max radius is not affected by sld of layer, sld of interface, or
// sld of solvent.
double PearlNecklaceModel :: calculate_ER() {
        PeralNecklaceParameters dp;

        dp.scale = scale();
        dp.radius = radius();
        dp.edge_separation = edge_separation();
        dp.thick_string = thick_string();
        dp.num_pearls = num_pearls();

        double rad_out = 0.0;
        // Perform the computation, with all weight points
        double sum = 0.0;
        double norm = 0.0;
        double pi = 4.0*atan(1.0);
        // No polydispersion supported in this model.
        // Get the dispersion points for the radius
        vector<WeightPoint> weights_radius;
        radius.get_weights(weights_radius);
        vector<WeightPoint> weights_edge_separation;
        edge_separation.get_weights(weights_edge_separation);
        // Perform the computation, with all weight points
        double string_vol = 0.0;
        double pearl_vol = 0.0;
        double tot_vol = 0.0;
        // Loop over core weight points
        for(size_t i=0; i<weights_radius.size(); i++) {
                dp.radius = weights_radius[i].value;
                // Loop over thick_inter0 weight points
                for(size_t j=0; j<weights_edge_separation.size(); j++) {
                        dp.edge_separation = weights_edge_separation[j].value;
                        pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3);
                        string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2);
                        tot_vol = (dp.num_pearls - 1.0) * string_vol;
                        tot_vol += dp.num_pearls * pearl_vol;
                        //Find  volume
                        // This may be a too much approximation
                        //Todo: decided whether or not we keep this calculation
                        sum += weights_radius[i].weight * weights_edge_separation[j].weight
                                * pow(3.0*tot_vol/4.0/pi,0.333333);
                        norm += weights_radius[i].weight * weights_edge_separation[j].weight;
                }
        }

        if (norm != 0){
                //return the averaged value
                rad_out =  sum/norm;}
        else{
                //return normal value
                pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3);
                string_vol =dp.edge_separation * pi * pow((dp.thick_string / 2.0), 2);
                tot_vol = (dp.num_pearls - 1.0) * string_vol;
                tot_vol += dp.num_pearls * pearl_vol;

                rad_out =  pow((3.0*tot_vol/4.0/pi), 0.33333);
                }

        return rad_out;

}