source: sasview/sansmodels/src/sans/models/c_models/parallelepiped.cpp @ 6c5239c

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 6c5239c was 8e36cdd, checked in by Jae Cho <jhjcho@…>, 15 years ago

fixed and tested ppmodel 2d and cleaned it up.

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