source: sasview/src/sas/models/c_extension/c_models/lamellarPS.cpp @ 7bcbeb7

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 7bcbeb7 was 37649847, checked in by richardh, 10 years ago

fix polydisp in lamellarPS; fix S(Q) in that and also lamellarPS_HG

  • Property mode set to 100644
File size: 5.1 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 "parameters.hh"
26#include <stdio.h>
27#include <stdlib.h>
28using namespace std;
29#include "lamellarPS.h"
30
31extern "C" {
32#include "libCylinder.h"
33}
34
35/*LamellarPS_kernel() was moved from libigor to get rid of polydipersity in del(thickness) that we provide from control panel.
36  LamellarPSX  :  calculates the form factor of a lamellar structure - with S(q) effects included
37-------
38------- resolution effects ARE NOT included, but only a CONSTANT default value, not the real q-dependent resolution!!
39
40 */
41static double LamellarPS_kernel(double dp[], double q)
42{
43  double scale,dd,del,sld_bi,sld_sol,contr,NN,Cp,bkg;   //local variables of coefficient wave
44  double inten, qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ;
45  double Pi,Euler,dQDefault,fii;
46  int ii,NNint;
47  Euler = 0.5772156649;   // Euler's constant
48  dQDefault = 0.0;    //[=] 1/A, q-resolution, default value
49  dQ = dQDefault;
50
51  Pi = 4.0*atan(1.0);
52  qval = q;
53
54  scale = dp[0];
55  dd = dp[1];
56  del = dp[2];
57  sld_bi = dp[3];
58  sld_sol = dp[4];
59  NN = trunc(dp[5]);    //be sure that NN is an integer
60  Cp = dp[6];
61  bkg = dp[7];
62
63  contr = sld_bi - sld_sol;
64
65  Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del));
66
67  NNint = (int)NN;    //cast to an integer for the loop
68  ii=0;
69  Sq = 0.0;
70  for(ii=1;ii<=(NNint-1);ii+=1) {
71
72    fii = (double)ii;   //do I really need to do this?
73
74    temp = 0.0;
75    alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler);
76    t1 = 2.0*dQ*dQ*dd*dd*alpha;
77    t2 = 2.0*qval*qval*dd*dd*alpha;
78    t3 = dQ*dQ*dd*dd*ii*ii;
79
80    temp = 1.0-ii/NN;
81    temp *= cos(dd*qval*ii/(1.0+t1));
82    temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
83    temp /= sqrt(1.0+t1);
84
85    Sq += temp;
86  }
87
88  Sq *= 2.0;
89  Sq += 1.0;
90
91  inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
92
93  inten *= 1.0e8;   // 1/A to 1/cm
94
95  return(inten+bkg);
96}
97
98LamellarPSModel :: LamellarPSModel() {
99  scale      = Parameter(1.0);
100  spacing    = Parameter(400.0, true);
101  spacing.set_min(0.0);
102  delta     = Parameter(30.0);
103  delta.set_min(0.0);
104  sld_bi   = Parameter(6.3e-6);
105  sld_sol   = Parameter(1.0e-6);
106  n_plates     = Parameter(20.0);
107  caille = Parameter(0.1);
108  background = Parameter(0.0);
109
110}
111
112/**
113 * Function to evaluate 1D scattering function
114 * The NIST IGOR library is used for the actual calculation.
115 * @param q: q-value
116 * @return: function value
117 */
118double LamellarPSModel :: operator()(double q) {
119  double dp[8];
120
121  // Fill parameter array for IGOR library
122  // Add the background after averaging
123  dp[0] = scale();
124  dp[1] = spacing();
125  dp[2] = delta();
126  dp[3] = sld_bi();
127  dp[4] = sld_sol();
128  dp[5] = n_plates();
129  dp[6] = caille();
130  dp[7] = 0.0;
131
132
133  // Get the dispersion points for spacing and delta (thickness)
134  vector<WeightPoint> weights_spacing;
135  spacing.get_weights(weights_spacing);
136  vector<WeightPoint> weights_delta;
137  delta.get_weights(weights_delta);
138
139  // Perform the computation, with all weight points
140  double sum = 0.0;
141  double norm = 0.0;
142
143  // Loop over short_edgeA weight points
144  for(int i=0; i< (int)weights_spacing.size(); i++) {
145    dp[1] = weights_spacing[i].value;
146    //for(int j=0; j< (int)weights_spacing.size(); j++) {    BUGS fixed March 2015
147    for(int j=0; j< (int)weights_delta.size(); j++) {
148      //dp[2] = weights_delta[i].value;        BUG
149      dp[2] = weights_delta[j].value;
150
151      sum += weights_spacing[i].weight * weights_delta[j].weight * LamellarPS_kernel(dp, q);
152      norm += weights_spacing[i].weight * weights_delta[j].weight;
153    }
154  }
155  return sum/norm + background();
156}
157/**
158 * Function to evaluate 2D scattering function
159 * @param q_x: value of Q along x
160 * @param q_y: value of Q along y
161 * @return: function value
162 */
163double LamellarPSModel :: operator()(double qx, double qy) {
164  double q = sqrt(qx*qx + qy*qy);
165  return (*this).operator()(q);
166}
167
168/**
169 * Function to evaluate 2D scattering function
170 * @param pars: parameters of the cylinder
171 * @param q: q-value
172 * @param phi: angle phi
173 * @return: function value
174 */
175double LamellarPSModel :: evaluate_rphi(double q, double phi) {
176  return (*this).operator()(q);
177}
178/**
179 * Function to calculate effective radius
180 * @return: effective radius value
181 */
182double LamellarPSModel :: calculate_ER() {
183  //NOT implemented yet!!!
184  return 0.0;
185}
186double LamellarPSModel :: calculate_VR() {
187  return 1.0;
188}
Note: See TracBrowser for help on using the repository browser.