source: sasview/sansmodels/src/sans/models/c_models/HayterMSA.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 25579e8, checked in by Jae Cho <jhjcho@…>, 16 years ago

Enable structure factors

  • Property mode set to 100644
File size: 4.5 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 "libStructureFactor.h"
31        #include "HayterMSA.h"
32}
33
34HayterMSAStructure :: HayterMSAStructure() {
35        radius      = Parameter(20.75, true);
36        radius.set_min(0.0);
37        charge      = Parameter(19.0, true);
38        charge.set_min(0.0);
39        volfraction = Parameter(0.0192, true);
40        volfraction.set_min(0.0);
41        temperature = Parameter(318.16, true);
42        temperature.set_min(0.0);
43        saltconc   = Parameter(0.0);
44        dielectconst  = Parameter(71.08);
45}
46
47/**
48 * Function to evaluate 1D scattering function
49 * The NIST IGOR library is used for the actual calculation.
50 * @param q: q-value
51 * @return: function value
52 */
53double HayterMSAStructure :: operator()(double q) {
54        double dp[6];
55
56        // Fill parameter array for IGOR library
57        // Add the background after averaging
58        dp[0] = radius();
59        dp[1] = charge();
60        dp[2] = volfraction();
61        dp[3] = temperature();
62        dp[4] = saltconc();
63        dp[5] = dielectconst();
64
65        // Get the dispersion points for the radius
66        vector<WeightPoint> weights_rad;
67        radius.get_weights(weights_rad);
68
69        // Perform the computation, with all weight points
70        double sum = 0.0;
71        double norm = 0.0;
72
73        // Loop over radius weight points
74        for(int i=0; i<weights_rad.size(); i++) {
75                dp[0] = weights_rad[i].value;
76
77                sum += weights_rad[i].weight
78                                * HayterPenfoldMSA(dp, q);
79                norm += weights_rad[i].weight;
80        }
81        return sum/norm ;
82}
83
84/**
85 * Function to evaluate 2D scattering function
86 * @param q_x: value of Q along x
87 * @param q_y: value of Q along y
88 * @return: function value
89 */
90double HayterMSAStructure :: operator()(double qx, double qy) {
91        HayterMSAParameters dp;
92        // Fill parameter array
93        dp.radius      = radius();
94        dp.charge      = charge();
95        dp.volfraction = volfraction();
96        dp.temperature   = temperature();
97        dp.saltconc   = saltconc();
98        dp.dielectconst   = dielectconst();
99
100        // Get the dispersion points for the radius
101        vector<WeightPoint> weights_rad;
102        radius.get_weights(weights_rad);
103
104        // Perform the computation, with all weight points
105        double sum = 0.0;
106        double norm = 0.0;
107
108        // Loop over radius weight points
109        for(int i=0; i<weights_rad.size(); i++) {
110                dp.radius = weights_rad[i].value;
111
112                                        double _ptvalue = weights_rad[i].weight
113                                                * HayterMSA_analytical_2DXY(&dp, qx, qy);
114                                        sum += _ptvalue;
115
116                                        norm += weights_rad[i].weight;
117        }
118        // Averaging in theta needs an extra normalization
119        // factor to account for the sin(theta) term in the
120        // integration (see documentation).
121        return sum/norm;
122}
123
124/**
125 * Function to evaluate 2D scattering function
126 * @param pars: parameters of the cylinder
127 * @param q: q-value
128 * @param phi: angle phi
129 * @return: function value
130 */
131double HayterMSAStructure :: evaluate_rphi(double q, double phi) {
132        double qx = q*cos(phi);
133        double qy = q*sin(phi);
134        return (*this).operator()(qx, qy);
135}
136
137// Testing code
138/*
139int main(void)
140{
141        SquareWellModel c = SquareWellModel();
142
143        printf("I(Qx=%g,Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001));
144        printf("I(Q=%g) = %g\n", 0.001, c(0.001));
145        c.radius.dispersion = new GaussianDispersion();
146        c.radius.dispersion->npts = 100;
147        c.radius.dispersion->width = 5;
148
149        //c.length.dispersion = GaussianDispersion();
150        //c.length.dispersion.npts = 20;
151        //c.length.dispersion.width = 65;
152
153        printf("I(Q=%g) = %g\n", 0.001, c(0.001));
154        printf("I(Q=%g) = %g\n", 0.001, c(0.001));
155        printf("I(Qx=%g, Qy=%g) = %g\n", 0.001, 0.001, c(0.001, 0.001));
156        printf("I(Q=%g,  Phi=%g) = %g\n", 0.00447, .7854, c.evaluate_rphi(sqrt(0.00002), .7854));
157
158
159
160        double i_avg = c(0.01, 0.01);
161        double i_1d = c(sqrt(0.0002));
162
163        printf("\nI(Qx=%g, Qy=%g) = %g\n", 0.01, 0.01, i_avg);
164        printf("I(Q=%g)         = %g\n", sqrt(0.0002), i_1d);
165        printf("ratio %g %g\n", i_avg/i_1d, i_1d/i_avg);
166
167
168        return 0;
169}
170*/
Note: See TracBrowser for help on using the repository browser.