source: sasview/src/sas/models/c_extension/c_models/multishell.cpp @ fc18690

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 fc18690 was 79492222, checked in by krzywon, 10 years ago

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

  • Property mode set to 100644
File size: 6.1 KB
Line 
1/**
2        This software was developed by the University of Tennessee as part of the
3        Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
4        project funded by the US National Science Foundation.
5
6        If you use DANSE applications to do scientific research that leads to
7        publication, we ask that you acknowledge the use of the software with the
8        following sentence:
9
10        "This work benefited from DANSE software developed under NSF award DMR-0520547."
11
12        copyright 2008, University of Tennessee
13 */
14
15/**
16 * Scattering model classes
17 * The classes use the IGOR library found in
18 *   sansmodels/src/libigor
19 *
20 */
21
22#include <math.h>
23#include "parameters.hh"
24#include <stdio.h>
25using namespace std;
26#include "multishell.h"
27
28extern "C" {
29#include "libSphere.h"
30}
31
32typedef struct {
33  double scale;
34  double core_radius;
35  double s_thickness;
36  double w_thickness;
37  double core_sld;
38  double shell_sld;
39  double n_pairs;
40  double background;
41
42} MultiShellParameters;
43
44MultiShellModel :: MultiShellModel() {
45  scale      = Parameter(1.0);
46  core_radius     = Parameter(60.0, true);
47  core_radius.set_min(0.0);
48  s_thickness  = Parameter(10.0, true);
49  s_thickness.set_min(0.0);
50  w_thickness   = Parameter(10.0, true);
51  w_thickness.set_min(0.0);
52  core_sld   = Parameter(6.4e-6);
53  shell_sld   = Parameter(4.0e-7);
54  n_pairs   = Parameter(2);
55  background = Parameter(0.0);
56}
57
58/**
59 * Function to evaluate 1D scattering function
60 * The NIST IGOR library is used for the actual calculation.
61 * @param q: q-value
62 * @return: function value
63 */
64double MultiShellModel :: operator()(double q) {
65  double dp[8];
66
67  // Fill parameter array for IGOR library
68  // Add the background after averaging
69  dp[0] = scale();
70  dp[1] = core_radius();
71  dp[2] = s_thickness();
72  dp[3] = w_thickness();
73  dp[4] = core_sld();
74  dp[5] = shell_sld();
75  dp[6] = n_pairs();
76  dp[7] = 0.0;
77
78  // Get the dispersion points for the core radius
79  vector<WeightPoint> weights_core_radius;
80  core_radius.get_weights(weights_core_radius);
81
82  // Get the dispersion points for the s_thickness
83  vector<WeightPoint> weights_s_thickness;
84  s_thickness.get_weights(weights_s_thickness);
85
86  // Get the dispersion points for the w_thickness
87  vector<WeightPoint> weights_w_thickness;
88  w_thickness.get_weights(weights_w_thickness);
89
90  // Perform the computation, with all weight points
91  double sum = 0.0;
92  double norm = 0.0;
93  double vol = 0.0;
94
95  // Loop over radius weight points
96  for(int i=0; i< (int)weights_core_radius.size(); i++) {
97    dp[1] = weights_core_radius[i].value;
98    for(int j=0; j< (int)weights_s_thickness.size(); j++){
99      dp[2] = weights_s_thickness[j].value;
100      for(int k=0; k< (int)weights_w_thickness.size(); k++){
101        dp[3] = weights_w_thickness[k].value;
102        //Un-normalize SphereForm by volume
103        sum += weights_core_radius[i].weight*weights_s_thickness[j].weight
104            *weights_w_thickness[k].weight* MultiShell(dp, q)
105        *pow(weights_core_radius[i].value+dp[6]*weights_s_thickness[j].value+(dp[6]-1)*weights_w_thickness[k].value,3);
106        //Find average volume
107        vol += weights_core_radius[i].weight*weights_s_thickness[j].weight
108            *weights_w_thickness[k].weight
109            *pow(weights_core_radius[i].value+dp[6]*weights_s_thickness[j].value+(dp[6]-1)*weights_w_thickness[k].value,3);
110        norm += weights_core_radius[i].weight*weights_s_thickness[j].weight
111            *weights_w_thickness[k].weight;
112      }
113    }
114  }
115  if (vol != 0.0 && norm != 0.0) {
116    //Re-normalize by avg volume
117    sum = sum/(vol/norm);}
118  return sum/norm + background();
119}
120
121/**
122 * Function to evaluate 2D scattering function
123 * @param q_x: value of Q along x
124 * @param q_y: value of Q along y
125 * @return: function value
126 */
127double MultiShellModel :: operator()(double qx, double qy) {
128  double q = sqrt(qx*qx + qy*qy);
129  return (*this).operator()(q);
130}
131
132/**
133 * Function to evaluate 2D scattering function
134 * @param pars: parameters of the multishell
135 * @param q: q-value
136 * @param phi: angle phi
137 * @return: function value
138 */
139double MultiShellModel :: evaluate_rphi(double q, double phi) {
140  return (*this).operator()(q);
141}
142/**
143 * Function to calculate effective radius
144 * @return: effective radius value
145 */
146double MultiShellModel :: calculate_ER() {
147  MultiShellParameters dp;
148
149  dp.core_radius     = core_radius();
150  dp.s_thickness  = s_thickness();
151  dp.w_thickness  = w_thickness();
152  dp.n_pairs = n_pairs();
153
154  double rad_out = 0.0;
155
156  // Perform the computation, with all weight points
157  double sum = 0.0;
158  double norm = 0.0;
159  if (dp.n_pairs <= 0.0 ){
160    dp.n_pairs = 0.0;
161  }
162
163  // Get the dispersion points for the core radius
164  vector<WeightPoint> weights_core_radius;
165  core_radius.get_weights(weights_core_radius);
166
167  // Get the dispersion points for the s_thickness
168  vector<WeightPoint> weights_s_thickness;
169  s_thickness.get_weights(weights_s_thickness);
170
171  // Get the dispersion points for the w_thickness
172  vector<WeightPoint> weights_w_thickness;
173  w_thickness.get_weights(weights_w_thickness);
174
175  // Loop over major shell weight points
176  for(int i=0; i< (int)weights_s_thickness.size(); i++) {
177    dp.s_thickness = weights_s_thickness[i].value;
178    for(int j=0; j< (int)weights_w_thickness.size(); j++) {
179      dp.w_thickness = weights_w_thickness[j].value;
180      for(int k=0; k< (int)weights_core_radius.size(); k++) {
181        dp.core_radius = weights_core_radius[k].value;
182        sum += weights_s_thickness[i].weight*weights_w_thickness[j].weight
183            * weights_core_radius[k].weight*(dp.core_radius+dp.n_pairs*dp.s_thickness+(dp.n_pairs-1.0)*dp.w_thickness);
184        norm += weights_s_thickness[i].weight*weights_w_thickness[j].weight* weights_core_radius[k].weight;
185      }
186    }
187  }
188  if (norm != 0){
189    //return the averaged value
190    rad_out =  sum/norm;}
191  else{
192    //return normal value
193    rad_out = (dp.core_radius+dp.n_pairs*dp.s_thickness+(dp.n_pairs-1.0)*dp.w_thickness);}
194
195  return rad_out;
196}
197double MultiShellModel :: calculate_VR() {
198  return 1.0;
199}
Note: See TracBrowser for help on using the repository browser.