source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2010/PolyRaspberry.ipf @ 712

Last change on this file since 712 was 712, checked in by ajj, 12 years ago

Adding Raspberry (pickering emulsion type object) scattering model from:

Larson-Smith, Jackson and Pozzo, J. Coll. Int. Sci, 343 (1), 36-41, 2010

File size: 8.4 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3
4// Raspberry particles with polydisperse large sphere
5#include "Raspberry"
6
7Macro PlotPolyRaspberry(num,qmin,qmax)
8        Variable num=500, qmin=1e-5, qmax=0.7
9        Prompt num "Enter number of data points for model: "
10        Prompt qmin "Enter minimum q-value (^-1) for model: "
11        Prompt qmax "Enter maximum q-value (^-1) for model: "
12//
13        Make/O/D/n=(num) xwave_PolyRaspberry, ywave_PolyRaspberry
14        xwave_PolyRaspberry =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
15        Make/O/D coef_PolyRaspberry = {0.05,5000,0.1,-4e-7,0.005,100,0.4,3.5e-6,0,6.3e-6,0.0}                   
16        make/o/t parameters_PolyRaspberry =  {"vf Large","Radius Large (A)","pd Large Sphere","SLD Large sphere (A-2)","vf Small", "Radius Small (A)","surface coverage","SLD Small sphere (A-2)","delta","SLD solvent (A-2)","bkgd (cm-1)"}
17        Edit parameters_PolyRaspberry, coef_PolyRaspberry
18       
19        Variable/G root:g_PolyRaspberry
20        g_PolyRaspberry := PolyRaspberry(coef_PolyRaspberry, ywave_PolyRaspberry, xwave_PolyRaspberry)
21        Display ywave_PolyRaspberry vs xwave_PolyRaspberry
22        ModifyGraph marker=29, msize=2, mode=4
23        ModifyGraph log=1,grid=1,mirror=2
24        Label bottom "q (\\S-1\\M) "
25        Label left "I(q) (cm\\S-1\\M)"
26        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
27       
28        AddModelToStrings("PolyRaspberry","coef_PolyRaspberry","parameters_PolyRaspberry","PolyRaspberry")
29//
30End
31
32
33Macro PlotSmearedPolyRaspberry(str)                                                             
34        String str
35        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
36       
37        // if any of the resolution waves are missing => abort
38        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
39                Abort
40        endif
41       
42        SetDataFolder $("root:"+str)
43       
44        // Setup parameter table for model function
45        Make/O/D smear_coef_PolyRaspberry = {0.05,5000,0.1,-4e-7,0.005,100,0.4,3.5e-6,0,6.3e-6,0.0}
46        make/o/t smear_parameters_PolyRaspberry = {"vf Large","Radius Large (A)","pd Large Sphere","SLD Large sphere (A-2)","vf Small", "Radius Small (A)","surface coverage","SLD Small sphere (A-2)","delta","SLD solvent (A-2)","bkgd (cm-1)"}
47        Edit smear_parameters_PolyRaspberry,smear_coef_PolyRaspberry                                    //display parameters in a table
48       
49        // output smeared intensity wave, dimensions are identical to experimental QSIG values
50        // make extra copy of experimental q-values for easy plotting
51        Duplicate/O $(str+"_q") smeared_PolyRaspberry,smeared_qvals
52        SetScale d,0,0,"1/cm",smeared_PolyRaspberry
53                                       
54        Variable/G gs_PolyRaspberry=0
55        gs_PolyRaspberry := fSmearedPolyRaspberry(smear_coef_PolyRaspberry,smeared_PolyRaspberry,smeared_qvals) //this wrapper fills the STRUCT
56       
57        Display smeared_PolyRaspberry vs smeared_qvals
58        ModifyGraph log=1,marker=29,msize=2,mode=4
59        Label bottom "q (\\S-1\\M)"
60        Label left "I(q) (cm\\S-1\\M)"
61        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
62       
63        SetDataFolder root:
64        AddModelToStrings("SmearedPolyRaspberry","smear_coef_PolyRaspberry","smear_parameters_PolyRaspberry","PolyRaspberry")
65End
66
67
68Function PolyRaspberry(cw,yw,xw) : FitFunc
69        Wave cw,yw,xw
70       
71#if exists("PolyRaspberryX")
72        yw = PolyRaspberryX(cw,xw)
73#else
74        yw = fPolyRaspberry(cw,xw)
75#endif
76        return(0)
77End
78
79
80Function fPolyRaspberry(w,x) : FitFunc
81        Wave w
82        Variable x
83       
84
85        //Set up variables
86        // variables are:                                                       
87        //[0] volume fraction large spheres
88        //[1] radius large sphere ()
89        //[2] polydispersity large sphere
90        //[3] sld large sphere (-2)
91        //[4] volume fraction small spheres
92        //[5] fraction of small spheres at surface
93        //[6] radius small sphere (A)
94        //[7] sld small sphere
95        //[8] small sphere penetration (A)
96        //[9] sld solvent
97        //[10] background (cm-1)
98       
99        Variable vfL,rL,pdL,sldL,vfS,rS,sldS,deltaS,delrhoL,delrhoS,bkg,sldSolv,qval,aSs,fSs   
100        vfL = w[0]
101        rL = w[1]
102        pdL = w[2]
103        sldL = w[3]
104        vfS = w[4]
105        rS = w[5]
106        aSs = w[6]
107        sldS = w[7]
108        deltaS = w[8]
109        sldSolv = w[9]
110        bkg = w[10]
111       
112        delrhoL = abs(sldL - sldSolv)
113        delrhoS = abs(sldS - sldSolv)   
114       
115        qval = x                //rename the input q-value, purely for readability
116       
117        Variable f2     
118        Variable va,vb,ii,zi,nord,yy,summ
119        String weightStr,zStr
120       
121        Variable Np,VL,VS
122       
123        Variable sig = pdL*rL           
124       
125        //select number of gauss points by setting nord=20 or76 points
126//      nord = 20
127        nord = 76
128       
129        weightStr = "gauss"+num2str(nord)+"wt"
130        zStr = "gauss"+num2str(nord)+"z"
131       
132        if (WaveExists($weightStr) == 0) // wave reference is not valid,
133                Make/D/N=(nord) $weightStr,$zStr
134                Wave gauWt = $weightStr
135                Wave gauZ = $zStr               // wave references to pass
136                if(nord==20)
137                        Make20GaussPoints(gauWt,gauZ)
138                else
139                        Make76GaussPoints(gauWt,gauZ)
140                endif   
141        else
142                if(exists(weightStr) > 1)
143                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
144                endif
145                Wave gauWt = $weightStr
146                Wave gauZ = $zStr               // create the wave references
147        endif
148       
149        // end points of integration
150        // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero
151        // +/- 3 sigq catches 99.73% of distrubution
152        // change limits (and spacing of zi) at each evaluation based on R()
153        //integration from va to vb
154        va = -4*sig + rL
155        if (va< 4*rS)
156                va=4*rS         //to avoid numerical error when  va<0 (-ve q-value)
157        endif
158        vb = 4*sig +rL
159
160        //VL = 4*pi/3*rL^3
161        //VS = 4*pi/3*rS^3
162        //Np = vfS*fSs*VL/vfL/VS
163
164       
165        Make/O/N=9 rasp_temp
166        rasp_temp[0] = w[0]
167        rasp_temp[1] = w[1]
168        rasp_temp[2] = delrhoL
169        rasp_temp[3] = w[4]
170        rasp_temp[4] = w[5]
171        rasp_temp[5] = w[6]
172        rasp_temp[6] = delrhoS
173        rasp_temp[7] = w[8]
174        rasp_temp[8] = w[9]
175        rasp_temp[9] = Np
176
177        summ = 0.0              // initialize integral
178        for(ii=0;ii<nord;ii+=1)
179                // calculate Gauss points on integration interval (r-value for evaluation)
180                zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
181                rasp_temp[1] = zi       
182                //calculate scattering
183                yy = gauWt[ii] * RaspGauss_distr(sig,rL,zi) * fRaspberryKernel(rasp_temp,qval)
184                summ += yy              //add to the running total of the quadrature   
185        endfor
186       
187        summ = (vb-va)/2.0*summ
188       
189        //Use average volume of oil droplet, so fraction of particles should be correct...
190        VL = (4*pi/3*rL^3) *(1+3*pdL^2)
191        VS = 4*pi/3*rS^3
192       
193        //Using average volume of oil droplet, should get average Np...
194        Np = aSs*4*(rS/(rL+deltaS))*VL/VS
195        //Np = aSs*4*((rL+deltaS)/rS)^2
196
197        fSs = Np*vfL*VS/vfS/VL
198       
199        f2 = summ
200        f2 += vfS*(1-fSs)*delrhoS^2*VS*fRaspBes(qval,rS)*fRaspBes(qval,rS)
201       
202        f2 *= 1e8               // [=] 1/cm
203       
204        return (f2+bkg) // Scale, then add in the background
205
206End
207
208Function RaspGauss_distr(sig,avg,pt)
209        Variable sig,avg,pt
210       
211        Variable retval
212       
213        retval = (1/ ( sig*sqrt(2*Pi)) )*exp(-(avg-pt)^2/sig^2/2)
214        return(retval)
215End
216
217
218///////////////////////////////////////////////////////////////
219// smeared model calculation
220//
221// you don't need to do anything with this function, as long as
222// your Raspberry works correctly, you get the resolution-smeared
223// version for free.
224//
225// this is all there is to the smeared model calculation!
226Function SmearedPolyRaspberry(s) : FitFunc
227        Struct ResSmearAAOStruct &s
228
229//      the name of your unsmeared model (AAO) is the first argument
230        Smear_Model_20(PolyRaspberry,s.coefW,s.xW,s.yW,s.resW)
231
232        return(0)
233End
234
235///////////////////////////////////////////////////////////////
236
237
238// nothing to change here
239//
240//wrapper to calculate the smeared model as an AAO-Struct
241// fills the struct and calls the ususal function with the STRUCT parameter
242//
243// used only for the dependency, not for fitting
244//
245Function fSmearedPolyRaspberry(coefW,yW,xW)
246        Wave coefW,yW,xW
247       
248        String str = getWavesDataFolder(yW,0)
249        String DF="root:"+str+":"
250       
251        WAVE resW = $(DF+str+"_res")
252       
253        STRUCT ResSmearAAOStruct fs
254        WAVE fs.coefW = coefW   
255        WAVE fs.yW = yW
256        WAVE fs.xW = xW
257        WAVE fs.resW = resW
258       
259        Variable err
260        err = SmearedPolyRaspberry(fs)
261       
262        return (0)
263End
264
265
266// plots the Gauss distribution based on the coefficient values
267// a static calculation, so re-run each time
268//
269Macro PlotRaspDistribution()
270
271        variable pd,avg,zz,maxr,vf
272       
273        if(Exists("coef_pgs")!=1)
274                abort "You need to plot the unsmeared model first to create the coefficient table"
275        Endif
276        pd=coef_PolyRaspberry[2]
277        avg = coef_PolyRaspberry[1]
278        vf = coef_PolyRaspberry[0]
279       
280        Variable vfL,rL,pdL,sldL,vfS,rS,sldS,deltaS,delrhoL,delrhoS,bkg,sldSolv,qval    ,fSs   
281
282        vfL = coef_PolyRaspberry[0]
283        rL = coef_PolyRaspberry[1]
284        pdL = coef_PolyRaspberry[2]
285        sldL = coef_PolyRaspberry[3]
286        vfS = coef_PolyRaspberry[4]
287        rS = coef_PolyRaspberry[5]
288        aSs = coef_PolyRaspberry[6]
289        sldS = coef_PolyRaspberry[7]
290        deltaS = coef_PolyRaspberry[8]
291        sldSolv = coef_PolyRaspberry[9]
292        bkg = coef_PolyRaspberry[10]
293       
294        Make/O/D/N=1000 Rasp_distribution,Rasp_Vf,Rasp_Np,Rasp_VL
295        maxr =  avg*(1+10*pd)
296
297        SetScale/I x, 0, maxr, Rasp_distribution
298        Rasp_distribution = RaspGauss_distr(pd*avg,avg,x)
299        Display Rasp_distribution
300       
301End
Note: See TracBrowser for help on using the repository browser.