source: sasview/src/sas/models/c_extension/c_models/raspberry.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, 9 years ago

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

  • Property mode set to 100644
File size: 6.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#include <math.h>
17#include "parameters.hh"
18#include <stdio.h>
19#include <stdlib.h>
20using namespace std;
21#include "raspberry.h"
22
23// scattering
24// Modified from igor model: JHC - May 04, 2012
25//
26// you should write your function to calculate the intensity
27// for a single q-value (that's the input parameter x)
28// based on the wave (array) of parameters that you send it (w)
29// Ref: J. coll. inter. sci. (2010) vol. 343 (1) pp. 36-41.
30//  model calculation
31//
32static double raspberry_func(double dp[], double q){
33        // variables are:                                                       
34        //[0] volume fraction large spheres
35        //[1] radius large sphere (A)
36        //[2] sld large sphere (A-2)
37        //[3] volume fraction small spheres
38        //[4] radius small sphere (A)
39        //[5] fraction of small spheres at surface
40        //[6] sld small sphere
41        //[7] small sphere penetration (A)
42        //[8] sld solvent
43        //[9] background (cm-1)
44        double vfL, rL, sldL, vfS, rS, sldS, deltaS, delrhoL, delrhoS, bkg, sldSolv, aSs;
45        vfL = dp[0];
46        rL = dp[1];
47        sldL = dp[2];
48        vfS = dp[3];
49        rS = dp[4];
50        aSs = dp[5];
51        sldS = dp[6];
52        deltaS = dp[7];
53        sldSolv = dp[8];
54        bkg = dp[9];
55       
56        delrhoL = fabs(sldL - sldSolv);
57        delrhoS = fabs(sldS - sldSolv); 
58       
59        double VL, VS, Np, f2, fSs;
60        double pi = 4.0*atan(1.0);
61         
62        VL = 4*pi/3*pow(rL,3.0);
63        VS = 4*pi/3*pow(rS,3.0);
64        Np = aSs*4.0*pow(((rL+deltaS)/rS), 2.0);
65        fSs = Np*vfL*VS/vfS/VL;
66       
67        double rasp_temp[8];
68        rasp_temp[0] = dp[0];
69        rasp_temp[1] = dp[1];
70        rasp_temp[2] = delrhoL;
71        rasp_temp[3] = dp[3];
72        rasp_temp[4] = dp[4];
73        rasp_temp[5] = dp[5];
74        rasp_temp[6] = delrhoS;
75        rasp_temp[7] = dp[7];
76
77        f2 = raspberry_kernel(rasp_temp,q);
78        f2+= vfS*(1.0-fSs)*pow(delrhoS, 2)*VS*rasp_bes(q,rS)*rasp_bes(q,rS);
79       
80        // normalize to single particle volume and convert to 1/cm
81        f2 *= 1.0e8;            // [=] 1/cm
82       
83        return (f2+bkg);        // Scale, then add in the background
84}
85
86double raspberry_kernel(double dp[], double q){
87        // variables are:                                                       
88        //[0] volume fraction large spheres
89        //[1] radius large sphere (A)
90        //[2] sld large sphere (A-2)
91        //[3] volume fraction small spheres
92        //[4] radius small sphere (A)
93        //[5] fraction of small spheres at surface
94        //[6] sld small sphere
95        //[7] small sphere penetration (A)
96
97        double vfL, rL, vfS, rS, deltaS;
98        double delrhoL, delrhoS, qval, aSs;
99        vfL = dp[0];
100        rL = dp[1];
101        delrhoL = dp[2];
102        vfS = dp[3];
103        rS = dp[4];
104        aSs = dp[5];
105        delrhoS = dp[6];
106        deltaS = dp[7];
107                       
108        qval = q;               //rename the input q-value, purely for readability
109        double pi = 4.0*atan(1.0);
110               
111        double psiL,psiS,f2;
112        double sfLS,sfSS;
113        double VL,VS,slT,Np,fSs;
114
115        VL = 4.0*pi/3.0*pow(rL,3.0);
116        VS = 4.0*pi/3.0*pow(rS,3.0);
117
118        Np = aSs*4.0*(rS/(rL+deltaS))*VL/VS; 
119
120        fSs = Np*vfL*VS/vfS/VL;
121
122        slT = delrhoL*VL + Np*delrhoS*VS;
123
124        psiL = rasp_bes(qval,rL);
125        psiS = rasp_bes(qval,rS);
126
127        sfLS = psiL*psiS*(sin(qval*(rL+deltaS*rS))/qval/(rL+deltaS*rS));
128        sfSS = psiS*psiS*pow((sin(qval*(rL+deltaS*rS))/qval/(rL+deltaS*rS)),2);
129               
130        f2 = delrhoL*delrhoL*VL*VL*psiL*psiL; 
131        f2 += Np*delrhoS*delrhoS*VS*VS*psiS*psiS; 
132        f2 += Np*(Np-1)*delrhoS*delrhoS*VS*VS*sfSS; 
133        f2 += 2*Np*delrhoL*delrhoS*VL*VS*sfLS;
134        if (f2 != 0.0){
135                f2 = f2/slT/slT;
136        }
137
138        f2 = f2*(vfL*delrhoL*delrhoL*VL + vfS*fSs*Np*delrhoS*delrhoS*VS);
139
140        return f2;
141}
142
143double rasp_bes(double qval, double rad){
144        double retval;
145        retval = 3.0*(sin(qval*rad)-qval*rad*cos(qval*rad))/(qval*qval*qval)/(rad*rad*rad);
146        return retval;
147}
148
149RaspBerryModel :: RaspBerryModel() {
150        volf_Lsph = Parameter(0.05);
151        radius_Lsph = Parameter(5000.0, true);
152        radius_Lsph.set_min(0.0);
153        sld_Lsph = Parameter(-4.0e-7);
154        volf_Ssph = Parameter(0.005);
155        radius_Ssph = Parameter(100.0, true);
156        radius_Ssph.set_min(0.0);
157        surfrac_Ssph = Parameter(0.4);
158        sld_Ssph = Parameter(3.5e-6);
159        delta_Ssph = Parameter(0.0);
160        sld_solv   = Parameter(6.3e-6);
161        background = Parameter(0.0);
162}
163
164/**
165 * Function to evaluate 1D scattering function
166 * The NIST IGOR is used for the actual calculation.
167 * @param q: q-value
168 * @return: function value
169 */
170double RaspBerryModel :: operator()(double q) {
171        double dp[10];
172        // Add the background after averaging
173        dp[0] = volf_Lsph();
174        dp[1] = radius_Lsph();
175        dp[2] = sld_Lsph();
176        dp[3] = volf_Ssph();
177        dp[4] = radius_Ssph();
178        dp[5] = surfrac_Ssph();
179        dp[6] = sld_Ssph();
180        dp[7] = delta_Ssph();
181        dp[8] = sld_solv();
182        dp[9] = 0.0;
183
184        // Get the dispersion points for the radius_Lsph
185        vector<WeightPoint> weights_radius_Lsph;
186        radius_Lsph.get_weights(weights_radius_Lsph);
187
188        // Perform the computation, with all weight points
189        double sum = 0.0;
190        double norm = 0.0;
191        //double norm_vol = 0.0;
192        //double vol = 0.0;
193
194        // Loop over radius_Lsph weight points
195        for(size_t i=0; i<weights_radius_Lsph.size(); i++) {
196                dp[1] = weights_radius_Lsph[i].value;
197                        //Un-normalize by volume
198                        sum += weights_radius_Lsph[i].weight
199                                * raspberry_func(dp, q);// * pow(weights_radius_Lsph[i].value,3.0);
200                        //Find average volume
201                        //vol += weights_radius_Lsph[i].weight
202                        //      * pow(weights_radius_Lsph[i].value,3.0);
203                        norm += weights_radius_Lsph[i].weight;
204                        //norm_vol += weights_radius_Lsph[i].weight;
205                }
206       
207        //if (vol != 0.0 && norm_vol != 0.0) {
208                //Re-normalize by avg volume
209                //sum = sum/(vol/norm_vol);}
210        return sum/norm + background();
211}
212
213/**
214 * Function to evaluate 2D scattering function
215 * @param q_x: value of Q along x
216 * @param q_y: value of Q along y
217 * @return: function value
218 */
219double RaspBerryModel :: operator()(double qx, double qy) {
220        double q = sqrt(qx*qx + qy*qy);
221        return (*this).operator()(q);
222}
223
224/**
225 * Function to evaluate 2D scattering function
226 * @param pars: parameters
227 * @param q: q-value
228 * @param phi: angle phi
229 * @return: function value
230 */
231double RaspBerryModel :: evaluate_rphi(double q, double phi) {
232        return (*this).operator()(q);
233}
234
235/**
236 * Function to calculate effective radius
237 * @return: effective radius value
238 */
239double RaspBerryModel :: calculate_ER() {
240  //NOT implemented yet!!!
241  return 0.0;
242}
243double RaspBerryModel :: calculate_VR() {
244  return 1.0;
245}
Note: See TracBrowser for help on using the repository browser.