source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2010/PolyRaspberry_v40.ipf @ 762

Last change on this file since 762 was 762, checked in by srkline, 12 years ago

Fixing bugs associated with 2D fitting, making sure that smeared and unsmeared 2D models both fit correctly and report errors to the FitManager?.

Added "_v40" extension to the Raspberry models. These models still need to be added to the Procedure list
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.