source: sasview/sansmodels/src/c_models/sphere.cpp @ 9e9e9a5

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 9e9e9a5 was 5b07138, checked in by Jae Cho <jhjcho@…>, 12 years ago

removed some compile warnings

  • Property mode set to 100644
File size: 7.8 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 "sphere.h"
27
28extern "C" {
29        #include "libSphere.h"
30        #include "libmultifunc/libfunc.h"
31}
32// Convenience parameter structure
33typedef struct {
34    double scale;
35    double radius;
36    double sldSph;
37    double sldSolv;
38    double background;
39    double M0_sld_sph;
40    double M_theta_sph;
41    double M_phi_sph;
42    double M0_sld_solv;
43    double M_theta_solv;
44    double M_phi_solv;
45    double Up_frac_i;
46        double Up_frac_f;
47        double Up_theta;
48} SphereParameters;
49
50SphereModel :: SphereModel() {
51        scale      = Parameter(1.0);
52        radius     = Parameter(20.0, true);
53        radius.set_min(0.0);
54        sldSph   = Parameter(4.0e-6);
55        sldSolv   = Parameter(1.0e-6);
56        background = Parameter(0.0);
57        M0_sld_sph = Parameter(0.0e-6);
58        M_theta_sph = Parameter(0.0);
59        M_phi_sph = Parameter(0.0); 
60        M0_sld_solv = Parameter(0.0e-6);
61        M_theta_solv = Parameter(0.0);
62        M_phi_solv = Parameter(0.0); 
63        Up_frac_i = Parameter(0.5); 
64        Up_frac_f = Parameter(0.5);
65        Up_theta = Parameter(0.0);
66}
67
68/**
69 * Function to evaluate 1D scattering function
70 * The NIST IGOR library is used for the actual calculation.
71 * @param q: q-value
72 * @return: function value
73 */
74double SphereModel :: operator()(double q) {
75        double dp[5];
76
77        // Fill parameter array for IGOR library
78        // Add the background after averaging
79        dp[0] = scale();
80        dp[1] = radius();
81        dp[2] = sldSph();
82        dp[3] = sldSolv();
83        dp[4] = 0.0;
84
85        // Get the dispersion points for the radius
86        vector<WeightPoint> weights_rad;
87        radius.get_weights(weights_rad);
88
89        // Perform the computation, with all weight points
90        double sum = 0.0;
91        double norm = 0.0;
92        double vol = 0.0;
93
94        // Loop over radius weight points
95        for(size_t i=0; i<weights_rad.size(); i++) {
96                dp[1] = weights_rad[i].value;
97
98                //Un-normalize SphereForm by volume
99                sum += weights_rad[i].weight
100                        * SphereForm(dp, q) * pow(weights_rad[i].value,3);
101                //Find average volume
102                vol += weights_rad[i].weight
103                        * pow(weights_rad[i].value,3);
104
105                norm += weights_rad[i].weight;
106        }
107
108        if (vol != 0.0 && norm != 0.0) {
109                //Re-normalize by avg volume
110                sum = sum/(vol/norm);}
111        return sum/norm + background();
112}
113
114/**
115 * Function to evaluate 2D scattering function
116 * @param pars: parameters
117 * @param q: q-value
118 * @param q_x: q_x / q
119 * @param q_y: q_y / q
120 * @return: function value
121 */
122
123static double sphere_analytical_2D_scaled(SphereParameters *pars, double q, double q_x, double q_y) {
124        double dp[5];
125        //convert angle degree to radian
126        dp[0] = 1.0;
127        dp[1] = pars->radius;
128        dp[2] = 0.0;
129        dp[3] = 0.0;
130        dp[4] = 0.0;
131
132    double sldSph = pars->sldSph;
133    double sldSolv = pars->sldSolv;
134        double answer = 0.0;
135        double m_max = pars->M0_sld_sph;
136        double m_max_solv = pars->M0_sld_solv;
137
138        if (m_max < 1.0e-32 && m_max_solv < 1.0e-32){
139                dp[2] = sldSph;
140                dp[3] = sldSolv;
141                answer = SphereForm(dp, q);
142        }
143        else{
144                //double contrast = sldSph - sldSolv;
145                double qx = q_x;
146                double qy = q_y;
147                double s_theta = pars->Up_theta;
148                double m_phi = pars->M_phi_sph;
149                double m_theta = pars->M_theta_sph;
150                double m_phi_solv = pars->M_phi_solv;
151                double m_theta_solv = pars->M_theta_solv;
152                double in_spin = pars->Up_frac_i;
153                double out_spin = pars->Up_frac_f;
154                polar_sld p_sld;
155                polar_sld p_sld_solv;
156                p_sld = cal_msld(1, qx, qy, sldSph, m_max, m_theta, m_phi, 
157                                                        in_spin, out_spin, s_theta);
158                p_sld_solv = cal_msld(1, qx, qy, sldSolv, m_max_solv, m_theta_solv, m_phi_solv, 
159                                                                in_spin, out_spin, s_theta);
160                //up_up
161                if (in_spin > 0.0 && out_spin > 0.0){                   
162                        dp[2] = p_sld.uu;
163                        dp[3] = p_sld_solv.uu;
164                        answer += SphereForm(dp, q);
165                        }
166                //down_down
167                if (in_spin < 1.0 && out_spin < 1.0){
168                        dp[2] = p_sld.dd;
169                        dp[3] = p_sld_solv.dd;
170                        answer += SphereForm(dp, q);
171                        }
172                //up_down
173                if (in_spin > 0.0 && out_spin < 1.0){
174                        dp[2] = p_sld.re_ud;
175                        dp[3] = p_sld_solv.re_ud;
176                        answer += SphereForm(dp, q);   
177                        dp[2] = p_sld.im_ud;
178                        dp[3] = p_sld_solv.im_ud;
179                        answer += SphereForm(dp, q);
180                        }
181                //down_up       
182                if (in_spin < 1.0 && out_spin > 0.0){
183                        dp[2] = p_sld.re_du;
184                        dp[3] = p_sld_solv.re_du;
185                        answer += SphereForm(dp, q);   
186                        dp[2] = p_sld.im_du;
187                        dp[3] = p_sld_solv.im_du;
188                        answer += SphereForm(dp, q);
189                        }
190        }
191       
192        // add in the background
193        answer *= pars->scale;
194        answer += pars->background;
195        return answer;
196}
197
198
199/**
200 * Function to evaluate 2D scattering function
201 * @param pars: parameters
202 * @param q: q-value
203 * @return: function value
204 */
205static double sphere_analytical_2DXY(SphereParameters *pars, double qx, double qy) {
206  double q;
207  q = sqrt(qx*qx+qy*qy);
208  return sphere_analytical_2D_scaled(pars, q, qx/q, qy/q);
209}
210
211
212/**
213 * Function to evaluate 2D scattering function
214 * @param q_x: value of Q along x
215 * @param q_y: value of Q along y
216 * @return: function value
217 */
218double SphereModel :: operator()(double qx, double qy) {
219        SphereParameters dp;
220        dp.scale = scale();
221        dp.radius = radius();
222        dp.sldSph = sldSph();
223        dp.sldSolv = sldSolv();
224        dp.background = 0.0;
225        dp.Up_theta =  Up_theta();
226        dp.M_phi_sph =  M_phi_sph();
227        dp.M_theta_sph =  M_theta_sph();
228        dp.M0_sld_sph =  M0_sld_sph();
229        dp.M_phi_solv =  M_phi_solv();
230        dp.M_theta_solv =  M_theta_solv();
231        dp.M0_sld_solv =  M0_sld_solv();
232        dp.Up_frac_i =  Up_frac_i();
233        dp.Up_frac_f =  Up_frac_f();
234
235        // Get the dispersion points for the radius
236        vector<WeightPoint> weights_rad;
237        radius.get_weights(weights_rad);
238
239        // Perform the computation, with all weight points
240        double sum = 0.0;
241        double norm = 0.0;
242        double vol = 0.0;
243
244        // Loop over radius weight points
245        for(size_t i=0; i<weights_rad.size(); i++) {
246                dp.radius = weights_rad[i].value;
247
248                //Un-normalize SphereForm by volume
249                sum += weights_rad[i].weight
250                        * sphere_analytical_2DXY(&dp, qx, qy) * pow(weights_rad[i].value,3);
251                //Find average volume
252                vol += weights_rad[i].weight
253                        * pow(weights_rad[i].value,3);
254
255                norm += weights_rad[i].weight;
256        }
257
258        if (vol != 0.0 && norm != 0.0) {
259                //Re-normalize by avg volume
260                sum = sum/(vol/norm);}
261        return sum/norm + background();
262}
263
264/**
265 * Function to evaluate 2D scattering function
266 * @param pars: parameters of the sphere
267 * @param q: q-value
268 * @param phi: angle phi
269 * @return: function value
270 */
271double SphereModel :: evaluate_rphi(double q, double phi) {
272        double qx = q*cos(phi);
273        double qy = q*sin(phi);
274        return (*this).operator()(qx, qy);
275}
276
277/**
278 * Function to calculate effective radius
279 * @return: effective radius value
280 */
281double SphereModel :: calculate_ER() {
282        double rad_out = 0.0;
283
284        // Perform the computation, with all weight points
285        double sum = 0.0;
286        double norm = 0.0;
287
288        // Get the dispersion points for the radius
289        vector<WeightPoint> weights_rad;
290        radius.get_weights(weights_rad);
291        // Loop over radius weight points to average the radius value
292        for(size_t i=0; i<weights_rad.size(); i++) {
293                sum += weights_rad[i].weight
294                        * weights_rad[i].value;
295                norm += weights_rad[i].weight;
296        }
297        if (norm != 0){
298                //return the averaged value
299                rad_out =  sum/norm;}
300        else{
301                //return normal value
302                rad_out = radius();}
303
304        return rad_out;
305}
306double SphereModel :: calculate_VR() {
307  return 1.0;
308}
Note: See TracBrowser for help on using the repository browser.