source: sasview/src/sas/models/c_extension/c_models/flexcyl_ellipX.cpp @ 79492222

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

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

  • Property mode set to 100644
File size: 6.4 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 "flexcyl_ellipX.h"
27
28extern "C" {
29#include "libCylinder.h"
30#include "libStructureFactor.h"
31}
32
33typedef struct {
34  double scale;
35  double length;
36  double kuhn_length;
37  double radius;
38  double axis_ratio;
39  double sldCyl;
40  double sldSolv;
41  double background;
42} FlexCylEXParameters;
43
44FlexCylEllipXModel :: FlexCylEllipXModel() {
45  scale      = Parameter(1.0);
46  length     = Parameter(1000.0, true);
47  length.set_min(0.0);
48  kuhn_length = Parameter(100.0, true);
49  kuhn_length.set_min(0.0);
50  radius  = Parameter(20.0, true);
51  radius.set_min(0.0);
52  axis_ratio = Parameter(1.5);
53  axis_ratio.set_min(0.0);
54  sldCyl   = Parameter(1.0e-6);
55  sldSolv   = Parameter(6.3e-6);
56  background = Parameter(0.0001);
57}
58
59/**
60 * Function to evaluate 1D scattering function
61 * The NIST IGOR library is used for the actual calculation.
62 * @param q: q-value
63 * @return: function value
64 */
65double FlexCylEllipXModel :: operator()(double q) {
66  double dp[8];
67
68  // Fill parameter array for IGOR library
69  // Add the background after averaging
70  dp[0] = scale();
71  dp[1] = length();
72  dp[2] = kuhn_length();
73  dp[3] = radius();
74  dp[4] = axis_ratio();
75  dp[5] = sldCyl();
76  dp[6] = sldSolv();
77  dp[7] = 0.0;
78
79  // Get the dispersion points for the length
80  vector<WeightPoint> weights_len;
81  length.get_weights(weights_len);
82
83  // Get the dispersion points for the kuhn_length
84  vector<WeightPoint> weights_kuhn;
85  kuhn_length.get_weights(weights_kuhn);
86
87  // Get the dispersion points for the radius
88  vector<WeightPoint> weights_rad;
89  radius.get_weights(weights_rad);
90
91  // Get the dispersion points for the axis_ratio
92  vector<WeightPoint> weights_ratio;
93  axis_ratio.get_weights(weights_ratio);
94
95  // Perform the computation, with all weight points
96  double sum = 0.0;
97  double norm = 0.0;
98  double vol = 0.0;
99
100  // Loop over length weight points
101  for(int i=0; i< (int)weights_len.size(); i++) {
102    dp[1] = weights_len[i].value;
103
104    // Loop over kuhn_length weight points
105    for(int j=0; j< (int)weights_kuhn.size(); j++) {
106      dp[2] = weights_kuhn[j].value;
107
108      // Loop over radius weight points
109      for(int k=0; k< (int)weights_rad.size(); k++) {
110        dp[3] = weights_rad[k].value;
111        // Loop over axis_ratio weight points
112        for(int l=0; l< (int)weights_ratio.size(); l++) {
113          dp[4] = weights_ratio[l].value;
114
115          //Un-normalize by volume
116          sum += weights_len[i].weight * weights_kuhn[j].weight*weights_rad[k].weight
117              * weights_ratio[l].weight * FlexCyl_Ellip(dp, q)
118          * (pow(weights_rad[k].value,2.0) * weights_ratio[l].value * weights_len[i].value);
119          //Find weighted volume
120          vol += weights_rad[k].weight * weights_kuhn[j].weight
121              * weights_len[i].weight * weights_ratio[l].weight
122              *pow(weights_rad[k].value,2.0)* weights_ratio[l].weight*weights_len[i].value;
123          norm += weights_len[i].weight * weights_kuhn[j].weight
124              *weights_rad[k].weight* weights_ratio[l].weight;
125        }
126      }
127    }
128  }
129  if (vol != 0.0 && norm != 0.0) {
130    //Re-normalize by avg volume
131    sum = sum/(vol/norm);}
132
133  return sum/norm + background();
134}
135
136/**
137 * Function to evaluate 2D scattering function
138 * @param q_x: value of Q along x
139 * @param q_y: value of Q along y
140 * @return: function value
141 */
142double FlexCylEllipXModel :: operator()(double qx, double qy) {
143  double q = sqrt(qx*qx + qy*qy);
144  return (*this).operator()(q);
145}
146
147/**
148 * Function to evaluate 2D scattering function
149 * @param pars: parameters of the triaxial ellipsoid
150 * @param q: q-value
151 * @param phi: angle phi
152 * @return: function value
153 */
154double FlexCylEllipXModel :: evaluate_rphi(double q, double phi) {
155  //double qx = q*cos(phi);
156  //double qy = q*sin(phi);
157  return (*this).operator()(q);
158}
159/**
160 * Function to calculate effective radius
161 * @return: effective radius value
162 */
163double FlexCylEllipXModel :: calculate_ER() {
164  FlexCylEXParameters dp;
165
166  dp.radius  = radius();
167  dp.length     = length();
168  dp.axis_ratio  = axis_ratio();
169
170  double rad_out = 0.0;
171  double suf_rad = sqrt(dp.radius*dp.radius*dp.axis_ratio );
172  // Perform the computation, with all weight points
173  double sum = 0.0;
174  double norm = 0.0;
175
176  // Get the dispersion points for the total length
177  vector<WeightPoint> weights_length;
178  length.get_weights(weights_length);
179
180  // Get the dispersion points for minor radius
181  vector<WeightPoint> weights_radius ;
182  radius.get_weights(weights_radius);
183
184  // Get the dispersion points for axis ratio = major_radius/minor_radius
185  vector<WeightPoint> weights_ratio ;
186  axis_ratio.get_weights(weights_ratio);
187
188  // Loop over major shell weight points
189  for(int i=0; i< (int)weights_length.size(); i++) {
190    dp.length = weights_length[i].value;
191    for(int k=0; k< (int)weights_radius.size(); k++) {
192      dp.radius = weights_radius[k].value;
193      // Loop over axis_ratio weight points
194      for(int l=0; l< (int)weights_ratio.size(); l++) {
195        dp.axis_ratio = weights_ratio[l].value;
196        suf_rad = sqrt(dp.radius  * dp.radius  * dp.axis_ratio);
197        //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
198        sum +=weights_length[i].weight * weights_radius[k].weight
199            * weights_ratio[l].weight *DiamCyl(dp.length,suf_rad)/2.0;
200        norm += weights_length[i].weight* weights_radius[k].weight* weights_ratio[l].weight;
201      }
202    }
203  }
204  if (norm != 0){
205    //return the averaged value
206    rad_out =  sum/norm;}
207  else{
208    //return normal value
209    //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
210    rad_out = DiamCyl(dp.length,suf_rad)/2.0;}
211
212  return rad_out;
213}
214double FlexCylEllipXModel :: calculate_VR() {
215  return 1.0;
216}
Note: See TracBrowser for help on using the repository browser.