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

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

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

File size: 27.5 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","parameters_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","parameters_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","parameters_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","parameters_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","smear_parameters_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","smear_parameters_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","smear_parameters_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","smear_parameters_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.