source: sasview/sansmodels/src/sans/models/c_models/oblate.cpp @ 27953d1

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 27953d1 was 975ec8e, checked in by Jae Cho <jhjcho@…>, 15 years ago

working on 2D models. Still need smore corrections and unit tests.

  • Property mode set to 100644
File size: 6.9 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 *      TODO: refactor so that we pull in the old sansmodels.c_extensions
21 */
22
23#include <math.h>
24#include "models.hh"
25#include "parameters.hh"
26#include <stdio.h>
27using namespace std;
28
29extern "C" {
30        #include "libCylinder.h"
31        #include "oblate.h"
32}
33
34OblateModel :: OblateModel() {
35        scale      = Parameter(1.0);
36        major_core     = Parameter(200.0, true);
37        major_core.set_min(0.0);
38        minor_core     = Parameter(20.0, true);
39        minor_core.set_min(0.0);
40        major_shell   = Parameter(250.0, true);
41        major_shell.set_min(0.0);
42        minor_shell    = Parameter(30.0, true);
43        minor_shell.set_min(0.0);
44        contrast   = Parameter(1e-6);
45        sld_solvent = Parameter(6.3e-6);
46        background = Parameter(0.0);
47}
48
49/**
50 * Function to evaluate 1D scattering function
51 * The NIST IGOR library is used for the actual calculation.
52 * @param q: q-value
53 * @return: function value
54 */
55double OblateModel :: operator()(double q) {
56        double dp[8];
57
58        // Fill parameter array for IGOR library
59        // Add the background after averaging
60        dp[0] = scale();
61        dp[1] = major_core();
62        dp[2] = minor_core();
63        dp[3] = major_shell();
64        dp[4] = minor_shell();
65        dp[5] = contrast();
66        dp[6] = sld_solvent();
67        dp[7] = 0.0;
68
69        // Get the dispersion points for the major core
70        vector<WeightPoint> weights_major_core;
71        major_core.get_weights(weights_major_core);
72
73        // Get the dispersion points for the minor core
74        vector<WeightPoint> weights_minor_core;
75        minor_core.get_weights(weights_minor_core);
76
77        // Get the dispersion points for the major shell
78        vector<WeightPoint> weights_major_shell;
79        major_shell.get_weights(weights_major_shell);
80
81        // Get the dispersion points for the minor_shell
82        vector<WeightPoint> weights_minor_shell;
83        minor_shell.get_weights(weights_minor_shell);
84
85
86        // Perform the computation, with all weight points
87        double sum = 0.0;
88        double norm = 0.0;
89
90        // Loop over major core weight points
91        for(int i=0; i<(int)weights_major_core.size(); i++) {
92                dp[1] = weights_major_core[i].value;
93
94                // Loop over minor core weight points
95                for(int j=0; j<(int)weights_minor_core.size(); j++) {
96                        dp[2] = weights_minor_core[j].value;
97
98                        // Loop over major shell weight points
99                        for(int k=0; k<(int)weights_major_shell.size(); k++) {
100                                dp[3] = weights_major_shell[k].value;
101
102                                // Loop over minor shell weight points
103                                for(int l=0; l<(int)weights_minor_shell.size(); l++) {
104                                        dp[4] = weights_minor_shell[l].value;
105
106                                        sum += weights_major_core[i].weight* weights_minor_core[j].weight * weights_major_shell[k].weight
107                                                * weights_minor_shell[l].weight * OblateForm(dp, q);
108                                        norm += weights_major_core[i].weight* weights_minor_core[j].weight * weights_major_shell[k].weight
109                                                        * weights_minor_shell[l].weight;
110                                }
111                        }
112                }
113        }
114        return sum/norm + background();
115}
116
117/**
118 * Function to evaluate 2D scattering function
119 * @param q_x: value of Q along x
120 * @param q_y: value of Q along y
121 * @return: function value
122 */
123/*
124double OblateModel :: operator()(double qx, double qy) {
125        double q = sqrt(qx*qx + qy*qy);
126
127        return (*this).operator()(q);
128}
129*/
130
131/**
132 * Function to evaluate 2D scattering function
133 * @param pars: parameters of the oblate
134 * @param q: q-value
135 * @param phi: angle phi
136 * @return: function value
137 */
138double OblateModel :: evaluate_rphi(double q, double phi) {
139        return (*this).operator()(q);
140}
141
142
143double OblateModel :: operator()(double qx, double qy) {
144        OblateParameters dp;
145        // Fill parameter array
146        dp.scale      = scale();
147        dp.major_core = major_core();
148        dp.minor_core = minor_core();
149        dp.major_shell = major_shell();
150        dp.minor_shell = minor_shell();
151        dp.contrast = contrast();
152        dp.sld_solvent = sld_solvent();
153        dp.background = background();
154        dp.axis_theta = axis_theta();
155        dp.axis_phi = axis_phi();
156
157        // Get the dispersion points for the major core
158        vector<WeightPoint> weights_major_core;
159        major_core.get_weights(weights_major_core);
160
161        // Get the dispersion points for the minor core
162        vector<WeightPoint> weights_minor_core;
163        minor_core.get_weights(weights_minor_core);
164
165        // Get the dispersion points for the major shell
166        vector<WeightPoint> weights_major_shell;
167        major_shell.get_weights(weights_major_shell);
168
169        // Get the dispersion points for the minor shell
170        vector<WeightPoint> weights_minor_shell;
171        minor_shell.get_weights(weights_minor_shell);
172
173
174        // Get angular averaging for theta
175        vector<WeightPoint> weights_theta;
176        axis_theta.get_weights(weights_theta);
177
178        // Get angular averaging for phi
179        vector<WeightPoint> weights_phi;
180        axis_phi.get_weights(weights_phi);
181
182        // Perform the computation, with all weight points
183        double sum = 0.0;
184        double norm = 0.0;
185
186        // Loop over major core weight points
187        for(int i=0; i< (int)weights_major_core.size(); i++) {
188                dp.major_core = weights_major_core[i].value;
189
190                // Loop over minor core weight points
191                for(int j=0; j< (int)weights_minor_core.size(); j++) {
192                        dp.minor_core = weights_minor_core[j].value;
193
194                        // Loop over major shell weight points
195                        for(int k=0; k< (int)weights_major_shell.size(); k++) {
196                                dp.major_shell = weights_major_shell[i].value;
197
198                                // Loop over minor shell weight points
199                                for(int l=0; l< (int)weights_minor_shell.size(); l++) {
200                                        dp.minor_shell = weights_minor_shell[l].value;
201
202                                        // Average over theta distribution
203                                        for(int m=0; m< (int)weights_theta.size(); m++) {
204                                                dp.axis_theta = weights_theta[m].value;
205
206                                                // Average over phi distribution
207                                                for(int n=0; n< (int)weights_phi.size(); n++) {
208                                                        dp.axis_phi = weights_phi[n].value;
209
210                                                        double _ptvalue = weights_major_core[i].weight *weights_minor_core[j].weight
211                                                                * weights_major_shell[k].weight * weights_minor_shell[l].weight
212                                                                * weights_theta[m].weight
213                                                                * weights_phi[n].weight
214                                                                * oblate_analytical_2DXY(&dp, qx, qy);
215                                                        if (weights_theta.size()>1) {
216                                                                _ptvalue *= sin(weights_theta[k].value);
217                                                        }
218                                                        sum += _ptvalue;
219
220                                                        norm += weights_major_core[i].weight *weights_minor_core[j].weight
221                                                                * weights_major_shell[k].weight * weights_minor_shell[l].weight
222                                                                * weights_theta[m].weight * weights_phi[n].weight;
223                                                }
224                                        }
225                                }
226                        }
227                }
228        }
229        // Averaging in theta needs an extra normalization
230        // factor to account for the sin(theta) term in the
231        // integration (see documentation).
232        if (weights_theta.size()>1) norm = norm / asin(1.0);
233        return sum/norm + background();
234}
235
Note: See TracBrowser for help on using the repository browser.