source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/PolyCore_and_NShells_v40.ipf @ 444

Last change on this file since 444 was 444, checked in by srkline, 14 years ago

Changes to the analysis package to add a few more model functions. Documentation and XOPs are to follow later.

General n-point gaussian quadrature has been added to GaussUtils? by including a Gauss-Laguere point generator from Numerical Recipes. See the Paracrystal models for an example, since they needed more than the 76 point quadrature.

File size: 27.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.0
3
4////////////////////////////////////////////////
5//
6// this function is for the form factor of a  with some
7// number of shells around a central core (currently 1-2-3-4)
8//
9// monodisperse and polydisperse (and smeared) versions are included
10// - for the polydisperse models, only polydispersity of the core is taken
11// into account, and dPolyOne numerically. for a Schulz distribution, this
12// should be possible to do analytically, whith a great savings in computation
13// time.
14//
15// It may also be useful to think of scenarios where the layers as well are
16// polydisperse - to break up the very regular spacing of the layers, which
17// is not a very natural structure.
18//
19// 03 MAR 04 SRK
20// 07 AUG 08 AJJ - redone for new version of software
21// 08 AUG 08 AJJ - bug fixed three and four shell models.
22////////////////////////////////////////////////
23
24#include "Core_and_NShells_v40"
25
26//this macro sets up all the necessary parameters and waves that are
27//needed to calculate the model function.
28//
29Proc PlotPolyOneShell(num,qmin,qmax)
30        Variable num=200, qmin=0.001, qmax=0.7
31        Prompt num "Enter number of data points for model: "
32        Prompt qmin "Enter minimum q-value (^-1) for model: "
33        Prompt qmax "Enter maximum q-value (^-1) for model: "
34//
35        Make/O/D/n=(num) xwave_PolyOneShell, ywave_PolyOneShell
36        xwave_PolyOneShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
37        Make/O/D coef_PolyOneShell = {1.,60,0.1,6.4e-6,10,1e-6,6.4e-6,0.001}
38        make/o/t parameters_PolyOneShell = {"scale","core radius (A)","Core Polydispersity(0,1)","Core SLD (A-2)","Shell thickness (A)","Shell SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
39        Edit parameters_PolyOneShell, coef_PolyOneShell
40       
41        Variable/G root:g_PolyOneShell
42        g_PolyOneShell := PolyOneShell(coef_PolyOneShell, ywave_PolyOneShell, xwave_PolyOneShell)
43        Display ywave_PolyOneShell vs xwave_PolyOneShell
44        ModifyGraph marker=29, msize=2, mode=4
45        ModifyGraph log=1,grid=1,mirror=2
46        Label bottom "q (\\S-1\\M) "
47        Label left "I(q) (cm\\S-1\\M)"
48        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
49       
50        AddModelToStrings("PolyOneShell","coef_PolyOneShell","PolyOneShell")
51//
52End
53
54Proc PlotPolyTwoShell(num,qmin,qmax)
55        Variable num=200, qmin=0.001, qmax=0.7
56        Prompt num "Enter number of data points for model: "
57        Prompt qmin "Enter minimum q-value (^-1) for model: "
58        Prompt qmax "Enter maximum q-value (^-1) for model: "
59//
60        Make/O/D/n=(num) xwave_PolyTwoShell, ywave_PolyTwoShell
61        xwave_PolyTwoShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
62        Make/O/D coef_PolyTwoShell = {1.,60,0.1,6.4e-6,10,1e-6,10,2e-6,6.4e-6,0.001}
63        make/o/t parameters_PolyTwoShell = {"scale","core radius (A)","Core Polydispersity(0,1)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
64        Edit parameters_PolyTwoShell, coef_PolyTwoShell
65       
66        Variable/G root:g_PolyTwoShell
67        g_PolyTwoShell := PolyTwoShell(coef_PolyTwoShell, ywave_PolyTwoShell, xwave_PolyTwoShell)
68        Display ywave_PolyTwoShell vs xwave_PolyTwoShell
69        ModifyGraph marker=29, msize=2, mode=4
70        ModifyGraph log=1,grid=1,mirror=2
71        Label bottom "q (\\S-1\\M) "
72        Label left "I(q) (cm\\S-1\\M)"
73        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
74       
75        AddModelToStrings("PolyTwoShell","coef_PolyTwoShell","PolyTwoShell")
76//
77End
78
79Proc PlotPolyThreeShell(num,qmin,qmax)
80        Variable num=200, qmin=0.001, qmax=0.7
81        Prompt num "Enter number of data points for model: "
82        Prompt qmin "Enter minimum q-value (^-1) for model: "
83        Prompt qmax "Enter maximum q-value (^-1) for model: "
84//
85        Make/O/D/n=(num) xwave_PolyThreeShell, ywave_PolyThreeShell
86        xwave_PolyThreeShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
87        Make/O/D coef_PolyThreeShell ={1.,60,0.1,6.4e-6,10,1e-6,10,2e-6,10,3e-6,6.4e-6,0.001}
88        make/o/t parameters_PolyThreeShell = {"scale","core radius (A)","Core Polydispersity(0,1)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Shell 3 thickness","Shell 3 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
89        Edit parameters_PolyThreeShell, coef_PolyThreeShell
90       
91        Variable/G root:g_PolyThreeShell
92        g_PolyThreeShell := PolyThreeShell(coef_PolyThreeShell, ywave_PolyThreeShell, xwave_PolyThreeShell)
93        Display ywave_PolyThreeShell vs xwave_PolyThreeShell
94        ModifyGraph marker=29, msize=2, mode=4
95        ModifyGraph log=1,grid=1,mirror=2
96        Label bottom "q (\\S-1\\M) "
97        Label left "I(q) (cm\\S-1\\M)"
98        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
99       
100        AddModelToStrings("PolyThreeShell","coef_PolyThreeShell","PolyThreeShell")
101//
102End
103
104Proc PlotPolyFourShell(num,qmin,qmax)
105        Variable num=200, qmin=0.001, qmax=0.7
106        Prompt num "Enter number of data points for model: "
107        Prompt qmin "Enter minimum q-value (^-1) for model: "
108        Prompt qmax "Enter maximum q-value (^-1) for model: "
109//
110        Make/O/D/n=(num) xwave_PolyFourShell, ywave_PolyFourShell
111        xwave_PolyFourShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
112        Make/O/D coef_PolyFourShell ={1.,60,0.1,6.4e-6,10,1e-6,10,2e-6,5,3e-6,10,4e-6,6.4e-6,0.001}
113        make/o/t parameters_PolyFourShell = {"scale","core radius (A)","Core Polydispersity(0,1)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Shell 3 thickness","Shell 3 SLD (A-2)","Shell 4 thickness","Shell 4 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
114        Edit parameters_PolyFourShell, coef_PolyFourShell
115       
116        Variable/G root:g_PolyFourShell
117        g_PolyFourShell := PolyFourShell(coef_PolyFourShell, ywave_PolyFourShell, xwave_PolyFourShell)
118        Display ywave_PolyFourShell vs xwave_PolyFourShell
119        ModifyGraph marker=29, msize=2, mode=4
120        ModifyGraph log=1,grid=1,mirror=2
121        Label bottom "q (\\S-1\\M) "
122        Label left "I(q) (cm\\S-1\\M)"
123        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
124       
125        AddModelToStrings("PolyFourShell","coef_PolyFourShell","PolyFourShell")
126//
127End
128
129
130
131
132//
133//this macro sets up all the necessary parameters and waves that are
134//needed to calculate the  smeared model function.
135//
136//no input parameters are necessary, it MUST use the experimental q-values
137// from the experimental data read in from an AVE/QSIG data file
138////////////////////////////////////////////////////
139// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
140Proc PlotSmearedPolyOneShell(str)                                                               
141        String str
142        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
143       
144        // if any of the resolution waves are missing => abort
145        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
146                Abort
147        endif
148       
149        SetDataFolder $("root:"+str)
150       
151        // Setup parameter table for model function
152        Make/O/D smear_coef_PolyOneShell =  {1.,60,0.1,6.4e-6,10,1e-6,6.4e-6,0.001}
153        make/o/t smear_parameters_PolyOneShell =  {"scale","core radius (A)","Core Polydispersity(0,1)","Core SLD (A-2)","Shell thickness (A)","Shell SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
154        Edit smear_parameters_PolyOneShell,smear_coef_PolyOneShell                                      //display parameters in a table
155       
156        // output smeared intensity wave, dimensions are identical to experimental QSIG values
157        // make extra copy of experimental q-values for easy plotting
158        Duplicate/O $(str+"_q") smeared_PolyOneShell,smeared_qvals
159        SetScale d,0,0,"1/cm",smeared_PolyOneShell
160                                       
161        Variable/G gs_PolyOneShell=0
162        gs_PolyOneShell := fSmearedPolyOneShell(smear_coef_PolyOneShell,smeared_PolyOneShell,smeared_qvals)     //this wrapper fills the STRUCT
163       
164        Display smeared_PolyOneShell vs smeared_qvals
165        ModifyGraph log=1,marker=29,msize=2,mode=4
166        Label bottom "q (\\S-1\\M)"
167        Label left "I(q) (cm\\S-1\\M)"
168        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
169       
170        SetDataFolder root:
171        AddModelToStrings("SmearedPolyOneShell","smear_coef_PolyOneShell","PolyOneShell")
172End
173
174Proc PlotSmearedPolyTwoShell(str)                                                               
175        String str
176        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
177       
178        // if any of the resolution waves are missing => abort
179        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
180                Abort
181        endif
182       
183        SetDataFolder $("root:"+str)
184       
185        // Setup parameter table for model function
186        Make/O/D smear_coef_PolyTwoShell =  {1.,60,0.1,6.4e-6,10,1e-6,10,2e-6,6.4e-6,0.001}
187        make/o/t smear_parameters_PolyTwoShell =  {"scale","core radius (A)","Core Polydispersity(0,1)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
188        Edit smear_parameters_PolyTwoShell,smear_coef_PolyTwoShell                                      //display parameters in a table
189       
190        // output smeared intensity wave, dimensions are identical to experimental QSIG values
191        // make extra copy of experimental q-values for easy plotting
192        Duplicate/O $(str+"_q") smeared_PolyTwoShell,smeared_qvals
193        SetScale d,0,0,"1/cm",smeared_PolyTwoShell
194                                       
195        Variable/G gs_PolyTwoShell=0
196        gs_PolyTwoShell := fSmearedPolyTwoShell(smear_coef_PolyTwoShell,smeared_PolyTwoShell,smeared_qvals)     //this wrapper fills the STRUCT
197       
198        Display smeared_PolyTwoShell vs smeared_qvals
199        ModifyGraph log=1,marker=29,msize=2,mode=4
200        Label bottom "q (\\S-1\\M)"
201        Label left "I(q) (cm\\S-1\\M)"
202        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
203       
204        SetDataFolder root:
205        AddModelToStrings("SmearedPolyTwoShell","smear_coef_PolyTwoShell","PolyTwoShell")
206End
207
208
209Proc PlotSmearedPolyThreeShell(str)                                                             
210        String str
211        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
212       
213        // if any of the resolution waves are missing => abort
214        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
215                Abort
216        endif
217       
218        SetDataFolder $("root:"+str)
219       
220        // Setup parameter table for model function
221        Make/O/D smear_coef_PolyThreeShell =  {1.,60,0.1,6.4e-6,10,1e-6,10,2e-6,5,3e-6,6.4e-6,0.001}
222        make/o/t smear_parameters_PolyThreeShell =  {"scale","core radius (A)","Core Polydispersity(0,1)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Shell 3 thickness","Shell 3 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
223        Edit smear_parameters_PolyThreeShell,smear_coef_PolyThreeShell                                  //display parameters in a table
224       
225        // output smeared intensity wave, dimensions are identical to experimental QSIG values
226        // make extra copy of experimental q-values for easy plotting
227        Duplicate/O $(str+"_q") smeared_PolyThreeShell,smeared_qvals
228        SetScale d,0,0,"1/cm",smeared_PolyThreeShell
229                                       
230        Variable/G gs_PolyThreeShell=0
231        gs_PolyThreeShell := fSmearedPolyThreeShell(smear_coef_PolyThreeShell,smeared_PolyThreeShell,smeared_qvals)     //this wrapper fills the STRUCT
232       
233        Display smeared_PolyThreeShell vs smeared_qvals
234        ModifyGraph log=1,marker=29,msize=2,mode=4
235        Label bottom "q (\\S-1\\M)"
236        Label left "I(q) (cm\\S-1\\M)"
237        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
238       
239        SetDataFolder root:
240        AddModelToStrings("SmearedPolyThreeShell","smear_coef_PolyThreeShell","PolyThreeShell")
241End
242
243Proc PlotSmearedPolyFourShell(str)                                                             
244        String str
245        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
246       
247        // if any of the resolution waves are missing => abort
248        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
249                Abort
250        endif
251       
252        SetDataFolder $("root:"+str)
253       
254        // Setup parameter table for model function
255        Make/O/D smear_coef_PolyFourShell = {1.,60,0.1,6.4e-6,10,1e-6,10,2e-6,5,3e-6,10,4e-6,6.4e-6,0.001}
256        make/o/t smear_parameters_PolyFourShell =  {"scale","core radius (A)","Core Polydispersity(0,1)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Shell 3 thickness","Shell 3 SLD (A-2)","Shell 4 thickness","Shell 4 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
257        Edit smear_parameters_PolyFourShell,smear_coef_PolyFourShell                                    //display parameters in a table
258       
259        // output smeared intensity wave, dimensions are identical to experimental QSIG values
260        // make extra copy of experimental q-values for easy plotting
261        Duplicate/O $(str+"_q") smeared_PolyFourShell,smeared_qvals
262        SetScale d,0,0,"1/cm",smeared_PolyFourShell
263                                       
264        Variable/G gs_PolyFourShell=0
265        gs_PolyFourShell := fSmearedPolyFourShell(smear_coef_PolyFourShell,smeared_PolyFourShell,smeared_qvals) //this wrapper fills the STRUCT
266       
267        Display smeared_PolyFourShell vs smeared_qvals
268        ModifyGraph log=1,marker=29,msize=2,mode=4
269        Label bottom "q (\\S-1\\M)"
270        Label left "I(q) (cm\\S-1\\M)"
271        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
272       
273        SetDataFolder root:
274        AddModelToStrings("SmearedPolyFourShell","smear_coef_PolyFourShell","PolyFourShell")
275End
276
277
278// nothing to change here
279//
280//AAO version, uses XOP if available
281// simply calls the original single point calculation with
282// a wave assignment (this will behave nicely if given point ranges)
283Function PolyOneShell(cw,yw,xw) : FitFunc
284        Wave cw,yw,xw
285       
286#if exists("PolyOneShellX")
287        yw = PolyOneShellX(cw,xw)
288#else
289        yw = fPolyOneShell(cw,xw)
290#endif
291        return(0)
292End
293
294Function PolyTwoShell(cw,yw,xw) : FitFunc
295        Wave cw,yw,xw
296       
297#if exists("PolyTwoShellX")
298        yw = PolyTwoShellX(cw,xw)
299#else
300        yw = fPolyTwoShell(cw,xw)
301#endif
302        return(0)
303End
304
305Function PolyThreeShell(cw,yw,xw) : FitFunc
306        Wave cw,yw,xw
307       
308#if exists("PolyThreeShellX")
309        yw =PolyThreeShellX(cw,xw)
310#else
311        yw = fPolyThreeShell(cw,xw)
312#endif
313        return(0)
314End
315
316Function PolyFourShell(cw,yw,xw) : FitFunc
317        Wave cw,yw,xw
318       
319#if exists("PolyFourShellX")
320        yw = PolyFourShellX(cw,xw)
321#else
322        yw = fPolyFourShell(cw,xw)
323#endif
324        return(0)
325End
326
327
328//
329// unsmeared model calculation
330//
331Function fPolyOneShell(w,x) : FitFunc
332        Wave w
333        Variable x
334
335        Variable scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg,pd,zz
336        scale = w[0]
337        rcore = w[1]
338        pd = w[2]
339        rhocore = w[3]
340        thick = w[4]
341        rhoshel = w[5]
342        rhosolv = w[6]
343        bkg = w[7]
344       
345        zz = (1/pd)^2-1         //polydispersity of the core only
346//
347// local variables
348        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq
349        Variable answer,zp1,zp2,zp3,vpoly
350        String weightStr,zStr
351
352        //select number of gauss points by setting nord=20 or76 points
353        NVAR/Z gNord=gNord
354        if(! NVAR_Exists(gNord) )
355                nord=76         //use 76 pts as default
356        else
357                if( (gNord == 20) || (gNord ==76) )
358                        nord = gNord    // should only allow 20 or 76 points
359                else
360                        abort "global value gNord in SchulzSpheres must be either 20 or 76"
361                endif
362        endif
363       
364        weightStr = "gauss"+num2str(nord)+"wt"
365        zStr = "gauss"+num2str(nord)+"z"
366       
367        if (WaveExists($weightStr) == 0) // wave reference is not valid,
368                Make/D/N=(nord) $weightStr,$zStr
369                Wave gauWt = $weightStr
370                Wave gauZ = $zStr               // wave references to pass
371                if(nord==20)
372                        Make20GaussPoints(gauWt,gauZ)
373                else
374                        Make76GaussPoints(gauWt,gauZ)
375                endif   
376        else
377                if(exists(weightStr) > 1)
378                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
379                endif
380                Wave gauWt = $weightStr
381                Wave gauZ = $zStr               // create the wave references
382        endif
383       
384// set up the integration end points and weights
385// limits are technically 0-inf, but wisely choose non-zero region of distribution
386        Variable range=8        //multiples of the std. dev. from the mean
387        va = rcore*(1-range*pd)
388        if (va<0)
389                va=0            //otherwise numerical error when pd >= 0.3, making a<0
390        endif
391        If(pd>0.3)
392                range = range + (pd-0.3)*18             //stretch upper range to account for skewed tail
393        Endif
394        vb = rcore*(1+range*pd) // is this far enough past avg radius?
395
396//temp set scale=1 and bkg=0 for quadrature calc
397        Make/O/D/N=7 temp_1sf
398        temp_1sf[0] = 1
399        temp_1sf[1] = w[1]              //the core radius will be changed in the loop
400        temp_1sf[2] = w[3]
401        temp_1sf[3] = w[4]
402        temp_1sf[4] = w[5]
403        temp_1sf[5] = w[6]
404        temp_1sf[6] = 0
405       
406// evaluate at Gauss points
407        summ = 0.0              // initialize integral
408        for(ii=0;ii<nord;ii+=1)
409                zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
410                temp_1sf[1] = zi               
411                yyy = gauWt[ii] * Schulz_Point_Nsf(zi,rcore,zz) *fOneShell(temp_1sf,x)
412                //un-normalize by volume
413                yyy *= 4*pi/3*(zi+thick)^3
414                summ += yyy
415        endfor
416// calculate value of integral to return
417   answer = (vb-va)/2.0*summ
418       
419   //re-normalize by the average volume
420   zp1 = zz + 1.
421   zp2 = zz + 2.
422   zp3 = zz + 3.
423   vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(rcore+thick)^3
424        answer /= vpoly
425//scale
426        answer *= scale
427// add in the background
428        answer += bkg
429
430        Return (answer)
431End
432
433
434Function fPolyTwoShell(w,x) : FitFunc
435        Wave w
436        Variable x
437
438        Variable scale,rcore,rhocore,rhosolv,bkg,pd,zz
439        Variable thick1,thick2
440        Variable rhoshel1,rhoshel2
441        scale = w[0]
442        rcore = w[1]
443        pd = w[2]
444        rhocore = w[3]
445        thick1 = w[4]
446        rhoshel1 = w[5]
447        thick2 = w[6]
448        rhoshel2 = w[7]
449        rhosolv = w[8]
450        bkg = w[9]
451       
452        zz = (1/pd)^2-1         //polydispersity of the core only
453//
454// local variables
455        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq
456        Variable answer,zp1,zp2,zp3,vpoly
457        String weightStr,zStr
458
459        //select number of gauss points by setting nord=20 or76 points
460        NVAR/Z gNord=gNord
461        if(! NVAR_Exists(gNord) )
462                nord=76         //use 76 pts as default
463        else
464                if( (gNord == 20) || (gNord ==76) )
465                        nord = gNord    // should only allow 20 or 76 points
466                else
467                        abort "global value gNord in SchulzSpheres must be either 20 or 76"
468                endif
469        endif
470       
471        weightStr = "gauss"+num2str(nord)+"wt"
472        zStr = "gauss"+num2str(nord)+"z"
473       
474        if (WaveExists($weightStr) == 0) // wave reference is not valid,
475                Make/D/N=(nord) $weightStr,$zStr
476                Wave gauWt = $weightStr
477                Wave gauZ = $zStr               // wave references to pass
478                if(nord==20)
479                        Make20GaussPoints(gauWt,gauZ)
480                else
481                        Make76GaussPoints(gauWt,gauZ)
482                endif   
483        else
484                if(exists(weightStr) > 1)
485                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
486                endif
487                Wave gauWt = $weightStr
488                Wave gauZ = $zStr               // create the wave references
489        endif
490       
491// set up the integration end points and weights
492// limits are technically 0-inf, but wisely choose non-zero region of distribution
493        Variable range=8        //multiples of the std. dev. from the mean
494        va = rcore*(1-range*pd)
495        if (va<0)
496                va=0            //otherwise numerical error when pd >= 0.3, making a<0
497        endif
498        If(pd>0.3)
499                range = range + (pd-0.3)*18             //stretch upper range to account for skewed tail
500        Endif
501        vb = rcore*(1+range*pd) // is this far enough past avg radius?
502
503//temp set scale=1 and bkg=0 for quadrature calc
504        Make/O/D/N=9 temp_2sf
505        temp_2sf[0] = 1
506        temp_2sf[1] = w[1]              //the core radius will be changed in the loop
507        temp_2sf[2] = w[3]
508        temp_2sf[3] = w[4]
509        temp_2sf[4] = w[5]
510        temp_2sf[5] = w[6]
511        temp_2sf[6] = w[7]
512        temp_2sf[7] = w[8]
513        temp_2sf[8] = 0
514       
515// evaluate at Gauss points
516        summ = 0.0              // initialize integral
517        for(ii=0;ii<nord;ii+=1)
518                zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
519                temp_2sf[1] = zi               
520                yyy = gauWt[ii] * Schulz_Point_Nsf(zi,rcore,zz) * fTwoShell(temp_2sf,x)
521                //un-normalize by volume
522                yyy *= 4*pi/3*(zi+thick1+thick2)^3
523                summ += yyy
524        endfor
525// calculate value of integral to return
526   answer = (vb-va)/2.0*summ
527       
528   //re-normalize by the average volume
529   zp1 = zz + 1.
530   zp2 = zz + 2.
531   zp3 = zz + 3.
532   vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(rcore+thick1+thick2)^3
533        answer /= vpoly
534//scale
535        answer *= scale
536// add in the background
537        answer += bkg
538
539        Return (answer)
540End
541
542Function fPolyThreeShell(w,x) : FitFunc
543        Wave w
544        Variable x
545
546        Variable scale,rcore,rhocore,rhosolv,bkg,pd,zz
547        Variable thick1,thick2,thick3
548        Variable rhoshel1,rhoshel2,rhoshel3
549        scale = w[0]
550        rcore = w[1]
551        pd = w[2]
552        rhocore = w[3]
553        thick1 = w[4]
554        rhoshel1 = w[5]
555        thick2 = w[6]
556        rhoshel2 = w[7]
557        thick3 = w[8]
558        rhoshel3 = w[9]
559        rhosolv = w[10]
560        bkg = w[11]
561       
562        zz = (1/pd)^2-1         //polydispersity of the core only
563//
564// local variables
565        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq
566        Variable answer,zp1,zp2,zp3,vpoly
567        String weightStr,zStr
568
569        //select number of gauss points by setting nord=20 or76 points
570        NVAR/Z gNord=gNord
571        if(! NVAR_Exists(gNord) )
572                nord=76         //use 76 pts as default
573        else
574                if( (gNord == 20) || (gNord ==76) )
575                        nord = gNord    // should only allow 20 or 76 points
576                else
577                        abort "global value gNord in SchulzSpheres must be either 20 or 76"
578                endif
579        endif
580       
581        weightStr = "gauss"+num2str(nord)+"wt"
582        zStr = "gauss"+num2str(nord)+"z"
583       
584        if (WaveExists($weightStr) == 0) // wave reference is not valid,
585                Make/D/N=(nord) $weightStr,$zStr
586                Wave gauWt = $weightStr
587                Wave gauZ = $zStr               // wave references to pass
588                if(nord==20)
589                        Make20GaussPoints(gauWt,gauZ)
590                else
591                        Make76GaussPoints(gauWt,gauZ)
592                endif   
593        else
594                if(exists(weightStr) > 1)
595                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
596                endif
597                Wave gauWt = $weightStr
598                Wave gauZ = $zStr               // create the wave references
599        endif
600       
601// set up the integration end points and weights
602// limits are technically 0-inf, but wisely choose non-zero region of distribution
603        Variable range=8        //multiples of the std. dev. from the mean
604        va = rcore*(1-range*pd)
605        if (va<0)
606                va=0            //otherwise numerical error when pd >= 0.3, making a<0
607        endif
608        If(pd>0.3)
609                range = range + (pd-0.3)*18             //stretch upper range to account for skewed tail
610        Endif
611        vb = rcore*(1+range*pd) // is this far enough past avg radius?
612
613//temp set scale=1 and bkg=0 for quadrature calc
614        Make/O/D/N=11 temp_3sf
615        temp_3sf[0] = 1
616        temp_3sf[1] = w[1]              //the core radius will be changed in the loop
617        temp_3sf[2] = w[3]
618        temp_3sf[3] = w[4]
619        temp_3sf[4] = w[5]
620        temp_3sf[5] = w[6]
621        temp_3sf[6] = w[7]
622        temp_3sf[7] = w[8]
623        temp_3sf[8] = w[9]
624        temp_3sf[9] = w[10]
625        temp_3sf[10] = 0
626       
627// evaluate at Gauss points
628        summ = 0.0              // initialize integral
629        for(ii=0;ii<nord;ii+=1)
630                zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
631                temp_3sf[1] = zi               
632                yyy = gauWt[ii] * Schulz_Point_Nsf(zi,rcore,zz) * fThreeShell(temp_3sf,x)
633                //un-normalize by volume
634                yyy *= 4*pi/3*(zi+thick1+thick2+thick3)^3
635                summ += yyy
636        endfor
637// calculate value of integral to return
638   answer = (vb-va)/2.0*summ
639       
640   //re-normalize by the average volume
641   zp1 = zz + 1.
642   zp2 = zz + 2.
643   zp3 = zz + 3.
644   vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(rcore+thick1+thick2+thick3)^3
645        answer /= vpoly
646//scale
647        answer *= scale
648// add in the background
649        answer += bkg
650
651        Return (answer)
652End
653
654
655Function fPolyFourShell(w,x) : FitFunc
656        Wave w
657        Variable x
658
659        Variable scale,rcore,rhocore,rhosolv,bkg,pd,zz
660        Variable thick1,thick2,thick3,thick4
661        Variable rhoshel1,rhoshel2,rhoshel3,rhoshel4
662        scale = w[0]
663        rcore = w[1]
664        pd = w[2]
665        rhocore = w[3]
666        thick1 = w[4]
667        rhoshel1 = w[5]
668        thick2 = w[6]
669        rhoshel2 = w[7]
670        thick3 = w[8]
671        rhoshel3 = w[9]
672        thick4 = w[10]
673        rhoshel4 = w[11]
674        rhosolv = w[12]
675        bkg = w[13]
676       
677        zz = (1/pd)^2-1         //polydispersity of the core only
678//
679// local variables
680        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq
681        Variable answer,zp1,zp2,zp3,vpoly
682        String weightStr,zStr
683
684        //select number of gauss points by setting nord=20 or76 points
685        NVAR/Z gNord=gNord
686        if(! NVAR_Exists(gNord) )
687                nord=76         //use 76 pts as default
688        else
689                if( (gNord == 20) || (gNord ==76) )
690                        nord = gNord    // should only allow 20 or 76 points
691                else
692                        abort "global value gNord in SchulzSpheres must be either 20 or 76"
693                endif
694        endif
695       
696        weightStr = "gauss"+num2str(nord)+"wt"
697        zStr = "gauss"+num2str(nord)+"z"
698       
699        if (WaveExists($weightStr) == 0) // wave reference is not valid,
700                Make/D/N=(nord) $weightStr,$zStr
701                Wave gauWt = $weightStr
702                Wave gauZ = $zStr               // wave references to pass
703                if(nord==20)
704                        Make20GaussPoints(gauWt,gauZ)
705                else
706                        Make76GaussPoints(gauWt,gauZ)
707                endif   
708        else
709                if(exists(weightStr) > 1)
710                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
711                endif
712                Wave gauWt = $weightStr
713                Wave gauZ = $zStr               // create the wave references
714        endif
715       
716// set up the integration end points and weights
717// limits are technically 0-inf, but wisely choose non-zero region of distribution
718        Variable range=8        //multiples of the std. dev. from the mean
719        va = rcore*(1-range*pd)
720        if (va<0)
721                va=0            //otherwise numerical error when pd >= 0.3, making a<0
722        endif
723        If(pd>0.3)
724                range = range + (pd-0.3)*18             //stretch upper range to account for skewed tail
725        Endif
726        vb = rcore*(1+range*pd) // is this far enough past avg radius?
727
728//temp set scale=1 and bkg=0 for quadrature calc
729        Make/O/D/N=13 temp_4sf
730        temp_4sf[0] = 1
731        temp_4sf[1] = w[1]              //the core radius will be changed in the loop
732        temp_4sf[2] = w[3]
733        temp_4sf[3] = w[4]
734        temp_4sf[4] = w[5]
735        temp_4sf[5] = w[6]
736        temp_4sf[6] = w[7]
737        temp_4sf[7] = w[8]
738        temp_4sf[8] = w[9]
739        temp_4sf[9] = w[10]
740        temp_4sf[10] = w[11]
741        temp_4sf[11] = w[12]
742        temp_4sf[12] = 0
743       
744// evaluate at Gauss points
745        summ = 0.0              // initialize integral
746        for(ii=0;ii<nord;ii+=1)
747                zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
748                temp_4sf[1] = zi               
749                yyy = gauWt[ii] * Schulz_Point_Nsf(zi,rcore,zz) * fFourShell(temp_4sf,x)
750                //un-normalize by volume
751                yyy *= 4*pi/3*(zi+thick1+thick2+thick3+thick4)^3
752                summ += yyy
753        endfor
754// calculate value of integral to return
755   answer = (vb-va)/2.0*summ
756       
757   //re-normalize by the average volume
758   zp1 = zz + 1.
759   zp2 = zz + 2.
760   zp3 = zz + 3.
761   vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(rcore+thick1+thick2+thick3+thick4)^3
762        answer /= vpoly
763//scale
764        answer *= scale
765// add in the background
766        answer += bkg
767
768        Return (answer)
769End
770
771
772
773
774///////////////////////////////////////////////////////////////
775// smeared model calculation
776//
777// you don't need to do anything with this function, as long as
778// your PolyOneShell works correctly, you get the resolution-smeared
779// version for free.
780//
781// this is all there is to the smeared model calculation!
782Function SmearedPolyOneShell(s) : FitFunc
783        Struct ResSmearAAOStruct &s
784
785//      the name of your unsmeared model (AAO) is the first argument
786        Smear_Model_20(PolyOneShell,s.coefW,s.xW,s.yW,s.resW)
787
788        return(0)
789End
790
791Function SmearedPolyTwoShell(s) : FitFunc
792        Struct ResSmearAAOStruct &s
793
794//      the name of your unsmeared model (AAO) is the first argument
795        Smear_Model_20(PolyTwoShell,s.coefW,s.xW,s.yW,s.resW)
796
797        return(0)
798End
799
800Function SmearedPolyThreeShell(s) : FitFunc
801        Struct ResSmearAAOStruct &s
802
803//      the name of your unsmeared model (AAO) is the first argument
804        Smear_Model_20(PolyThreeShell,s.coefW,s.xW,s.yW,s.resW)
805
806        return(0)
807End
808
809Function SmearedPolyFourShell(s) : FitFunc
810        Struct ResSmearAAOStruct &s
811
812//      the name of your unsmeared model (AAO) is the first argument
813        Smear_Model_20(PolyFourShell,s.coefW,s.xW,s.yW,s.resW)
814
815        return(0)
816End
817
818
819
820///////////////////////////////////////////////////////////////
821
822
823// nothing to change here
824//
825//wrapper to calculate the smeared model as an AAO-Struct
826// fills the struct and calls the ususal function with the STRUCT parameter
827//
828// used only for the dependency, not for fitting
829//
830Function fSmearedPolyOneShell(coefW,yW,xW)
831        Wave coefW,yW,xW
832       
833        String str = getWavesDataFolder(yW,0)
834        String DF="root:"+str+":"
835       
836        WAVE resW = $(DF+str+"_res")
837       
838        STRUCT ResSmearAAOStruct fs
839        WAVE fs.coefW = coefW   
840        WAVE fs.yW = yW
841        WAVE fs.xW = xW
842        WAVE fs.resW = resW
843       
844        Variable err
845        err = SmearedPolyOneShell(fs)
846       
847        return (0)
848End
849
850Function fSmearedPolyTwoShell(coefW,yW,xW)
851        Wave coefW,yW,xW
852       
853        String str = getWavesDataFolder(yW,0)
854        String DF="root:"+str+":"
855       
856        WAVE resW = $(DF+str+"_res")
857       
858        STRUCT ResSmearAAOStruct fs
859        WAVE fs.coefW = coefW   
860        WAVE fs.yW = yW
861        WAVE fs.xW = xW
862        WAVE fs.resW = resW
863       
864        Variable err
865        err = SmearedPolyTwoShell(fs)
866       
867        return (0)
868End
869
870Function fSmearedPolyThreeShell(coefW,yW,xW)
871        Wave coefW,yW,xW
872       
873        String str = getWavesDataFolder(yW,0)
874        String DF="root:"+str+":"
875       
876        WAVE resW = $(DF+str+"_res")
877       
878        STRUCT ResSmearAAOStruct fs
879        WAVE fs.coefW = coefW   
880        WAVE fs.yW = yW
881        WAVE fs.xW = xW
882        WAVE fs.resW = resW
883       
884        Variable err
885        err = SmearedPolyThreeShell(fs)
886       
887        return (0)
888End
889
890Function fSmearedPolyFourShell(coefW,yW,xW)
891        Wave coefW,yW,xW
892       
893        String str = getWavesDataFolder(yW,0)
894        String DF="root:"+str+":"
895       
896        WAVE resW = $(DF+str+"_res")
897       
898        STRUCT ResSmearAAOStruct fs
899        WAVE fs.coefW = coefW   
900        WAVE fs.yW = yW
901        WAVE fs.xW = xW
902        WAVE fs.resW = resW
903       
904        Variable err
905        err = SmearedPolyFourShell(fs)
906       
907        return (0)
908End
909
910//calculate the normalized distribution
911Static Function Schulz_Point_Nsf(x,avg,zz)
912        Variable x,avg,zz
913       
914        Variable dr
915       
916        dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1))
917        return (exp(dr))
918End
Note: See TracBrowser for help on using the repository browser.