source: sasview/src/sas/models/c_extension/c_models/linearpearls.cpp @ 3a39c2e

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 3a39c2e was 79492222, checked in by krzywon, 9 years ago

Changed the file and folder names to remove all SANS references.

  • Property mode set to 100644
File size: 6.5 KB
Line 
1
2#include <math.h>
3#include "parameters.hh"
4#include <stdio.h>
5using namespace std;
6#include "linearpearls.h"
7
8extern "C" {
9#include "libmultifunc/libfunc.h"
10}
11
12static double linear_pearls_kernel(double dp[], double q) {
13  //fit parameters
14  double scale = dp[0];
15  double radius = dp[1];
16  double edge_separation = dp[2];
17  double num_pearls = dp[3];
18  double sld_pearl = dp[4];
19  double sld_solv = dp[5];
20  double background = dp[6];
21  //others
22  double psi = 0.0;
23  double n_contrib = 0.0;
24  double form_factor = 0.0;
25  //Pi
26  double pi = 4.0 * atan(1.0);
27  //relative sld
28  double contrast_pearl = sld_pearl - sld_solv;
29  //each volume
30  double pearl_vol = 4.0 /3.0 * pi * pow(radius, 3.0);
31  //total volume
32  double tot_vol = num_pearls * pearl_vol;
33  //mass
34  double m_s = contrast_pearl * pearl_vol;
35  //center to center distance between the neighboring pearls
36  double separation = edge_separation + 2.0 * radius;
37  //integer
38  int num = 0;
39  int n_max = 0;
40  // constraints
41  if (scale<=0 || radius<=0 || edge_separation<0 || num_pearls<=0){
42          return 0.0;
43  }
44  //sine functions of a pearl
45  psi = sin(q * radius);
46  psi -= q * radius * cos(q * radius);
47  psi /= pow((q * radius), 3.0);
48
49  // N pearls contribution
50  n_max = num_pearls - 1;
51  for(num=0; num<=n_max; num++) {
52
53          if (num == 0){
54                  n_contrib = num_pearls;
55                  continue;
56          }
57          n_contrib += (2.0*(num_pearls-double(num))*sinc(q*separation*double(num)));
58  }
59  // form factor for num_pearls
60  form_factor = n_contrib;
61  form_factor *= pow((m_s*psi*3.0), 2.0);
62  form_factor /= (tot_vol * 1.0e-8); // norm by volume and A^-1 to cm^-1
63
64  // scale and background
65  form_factor *= scale;
66  form_factor += background;
67  return (form_factor);
68}
69
70LinearPearlsModel :: LinearPearlsModel() {
71  scale = Parameter(1.0);
72  radius = Parameter(80.0, true);
73  radius.set_min(0.0);
74  edge_separation = Parameter(350.0, true);
75  edge_separation.set_min(0.0);
76  num_pearls = Parameter(3);
77  num_pearls.set_min(0.0);
78  sld_pearl = Parameter(1.0e-06);
79  sld_solv = Parameter(6.3e-06);
80  background = Parameter(0.0);
81}
82
83/**
84 * Function to evaluate 1D LinearPearlsModel function
85 * @param q: q-value
86 * @return: function value
87 */
88double LinearPearlsModel :: operator()(double q) {
89  double dp[7];
90  // Add the background after averaging
91  dp[0] = scale();
92  dp[1] = radius();
93  dp[2] = edge_separation();
94  dp[3] = num_pearls();
95  dp[4] = sld_pearl();
96  dp[5] = sld_solv();
97  dp[6] = 0.0;
98  double pi = 4.0*atan(1.0);
99  // No polydispersion supported in this model.
100  // Get the dispersion points for the radius
101  vector<WeightPoint> weights_radius;
102  radius.get_weights(weights_radius);
103  vector<WeightPoint> weights_edge_separation;
104  edge_separation.get_weights(weights_edge_separation);
105  // Perform the computation, with all weight points
106  double sum = 0.0;
107  double norm = 0.0;
108  double vol = 0.0;
109  double pearl_vol = 0.0;
110  double tot_vol = 0.0;
111  // Loop over core weight points
112  for(size_t i=0; i<weights_radius.size(); i++) {
113    dp[1] = weights_radius[i].value;
114    // Loop over thick_inter0 weight points
115    for(size_t j=0; j<weights_edge_separation.size(); j++) {
116      dp[2] = weights_edge_separation[j].value;
117      pearl_vol = 4.0 /3.0 * pi * pow(dp[1], 3);
118      tot_vol += dp[3] * pearl_vol;
119      //Un-normalize Sphere by volume
120      sum += weights_radius[i].weight * weights_edge_separation[j].weight
121          * linear_pearls_kernel(dp,q) * tot_vol;
122      //Find average volume
123      vol += weights_radius[i].weight * weights_edge_separation[j].weight
124          * tot_vol;
125      norm += weights_radius[i].weight * weights_edge_separation[j].weight;
126    }
127  }
128  if (vol != 0.0 && norm != 0.0) {
129    //Re-normalize by avg volume
130    sum = sum/(vol/norm);}
131  return sum/norm + background();
132}
133
134/**
135 * Function to evaluate 2D LinearPearlsModel function
136 * @param q_x: value of Q along x
137 * @param q_y: value of Q along y
138 * @return: function value
139 */
140double LinearPearlsModel :: operator()(double qx, double qy) {
141  double q = sqrt(qx*qx + qy*qy);
142  return (*this).operator()(q);
143}
144
145/**
146 * Function to evaluate LinearPearlsModel function
147 * @param pars: parameters of the StringOfPearlsModel
148 * @param q: q-value
149 * @param phi: angle phi
150 * @return: function value
151 */
152double LinearPearlsModel :: evaluate_rphi(double q, double phi) {
153  return (*this).operator()(q);
154}
155
156/**
157 * Function to calculate TOTAL radius
158 * Todo: decide whether or not we keep this calculation
159 * @return: effective radius value
160 */
161// No polydispersion supported in this model.
162// Calculate max radius assumming max_radius = effective radius
163// Note that this max radius is not affected by sld of layer, sld of interface, or
164// sld of solvent.
165double LinearPearlsModel :: calculate_ER() {
166        LinearPearlsParameters dp;
167
168  dp.scale = scale();
169  dp.radius = radius();
170  dp.edge_separation = edge_separation();
171  dp.num_pearls = num_pearls();
172
173  double rad_out = 0.0;
174  // Perform the computation, with all weight points
175  double sum = 0.0;
176  double norm = 0.0;
177  double pi = 4.0*atan(1.0);
178  // Get the dispersion points for the radius
179  vector<WeightPoint> weights_radius;
180  radius.get_weights(weights_radius);
181  vector<WeightPoint> weights_edge_separation;
182  edge_separation.get_weights(weights_edge_separation);
183  // Perform the computation, with all weight points
184  double pearl_vol = 0.0;
185  double tot_vol = 0.0;
186  // Loop over core weight points
187  for(size_t i=0; i<weights_radius.size(); i++) {
188    dp.radius = weights_radius[i].value;
189    // Loop over thick_inter0 weight points
190    for(size_t j=0; j<weights_edge_separation.size(); j++) {
191      dp.edge_separation = weights_edge_separation[j].value;
192      pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3);
193      tot_vol = dp.num_pearls * pearl_vol;
194      //Find  volume
195      // This may be a too much approximation
196      //Todo: decided whether or not we keep this calculation
197      sum += weights_radius[i].weight * weights_edge_separation[j].weight
198          * pow(3.0*tot_vol/4.0/pi,0.333333);
199      norm += weights_radius[i].weight * weights_edge_separation[j].weight;
200    }
201  }
202
203  if (norm != 0){
204    //return the averaged value
205    rad_out =  sum/norm;}
206  else{
207    //return normal value
208    pearl_vol = 4.0 /3.0 * pi * pow(dp.radius , 3);
209    tot_vol = dp.num_pearls * pearl_vol;
210
211    rad_out =  pow((3.0*tot_vol/4.0/pi), 0.33333);
212  }
213
214  return rad_out;
215
216}
217double LinearPearlsModel :: calculate_VR() {
218  return 1.0;
219}
Note: See TracBrowser for help on using the repository browser.