source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/Core_and_NShells_v40.ipf @ 633

Last change on this file since 633 was 633, checked in by srkline, 13 years ago

Corrected models to explicitly return proper values for I(q=0). There are some models that just can't be fixed, and these typically return NaN. Some, however, are simply numerically unstable at extreme conditions. Beware.
File size: 22.2 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4////////////////////////////////////////////////
5//
6// this function is for the form factor of a sphere with some
7// number of shells around a central core (currently 1-2-3)
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 done 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////////////////////////////////////////////////
21// Four shell model added and polyCore versions written in separate ipf (polyCore_and_NShells.ipf)
22//
23// Four shell model added at request of H Wacklin to model lipid vesicles with varying deuteration
24//
25// DEC 2008 AJJ
26//////////////////////////////////////////////
27
28//this macro sets up all the necessary parameters and waves that are
29//needed to calculate the model function.
30//
31Proc PlotOneShell(num,qmin,qmax)
32        Variable num=200, qmin=0.001, qmax=0.7
33        Prompt num "Enter number of data points for model: "
34        Prompt qmin "Enter minimum q-value (^-1) for model: "
35        Prompt qmax "Enter maximum q-value (^-1) for model: "
36//
37        Make/O/D/n=(num) xwave_OneShell, ywave_OneShell
38        xwave_OneShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
39        Make/O/D coef_OneShell = {1.,60,6.4e-6,10,1e-6,6.4e-6,0.001}
40        make/o/t parameters_OneShell = {"scale","core radius (A)","Core SLD (A-2)","Shell thickness (A)","Shell SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
41        Edit parameters_OneShell, coef_OneShell
42       
43        Variable/G root:g_OneShell
44        g_OneShell := OneShell(coef_OneShell, ywave_OneShell, xwave_OneShell)
45        Display ywave_OneShell vs xwave_OneShell
46        ModifyGraph marker=29, msize=2, mode=4
47        ModifyGraph log=1,grid=1,mirror=2
48        Label bottom "q (\\S-1\\M) "
49        Label left "I(q) (cm\\S-1\\M)"
50        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
51       
52        AddModelToStrings("OneShell","coef_OneShell","parameters_OneShell","OneShell")
53//
54End
55
56Proc PlotTwoShell(num,qmin,qmax)
57        Variable num=200, qmin=0.001, qmax=0.7
58        Prompt num "Enter number of data points for model: "
59        Prompt qmin "Enter minimum q-value (^-1) for model: "
60        Prompt qmax "Enter maximum q-value (^-1) for model: "
61//
62        Make/O/D/n=(num) xwave_TwoShell, ywave_TwoShell
63        xwave_TwoShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
64        Make/O/D coef_TwoShell = {1.,60,6.4e-6,10,1e-6,10,2e-6,6.4e-6,0.001}
65        make/o/t parameters_TwoShell = {"scale","core radius (A)","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)"}
66        Edit parameters_TwoShell, coef_TwoShell
67       
68        Variable/G root:g_TwoShell
69        g_TwoShell := TwoShell(coef_TwoShell, ywave_TwoShell, xwave_TwoShell)
70        Display ywave_TwoShell vs xwave_TwoShell
71        ModifyGraph marker=29, msize=2, mode=4
72        ModifyGraph log=1,grid=1,mirror=2
73        Label bottom "q (\\S-1\\M) "
74        Label left "I(q) (cm\\S-1\\M)"
75        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
76       
77        AddModelToStrings("TwoShell","coef_TwoShell","parameters_TwoShell","TwoShell")
78//
79End
80
81Proc PlotThreeShell(num,qmin,qmax)
82        Variable num=200, qmin=0.001, qmax=0.7
83        Prompt num "Enter number of data points for model: "
84        Prompt qmin "Enter minimum q-value (^-1) for model: "
85        Prompt qmax "Enter maximum q-value (^-1) for model: "
86//
87        Make/O/D/n=(num) xwave_ThreeShell, ywave_ThreeShell
88        xwave_ThreeShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
89        Make/O/D coef_ThreeShell ={1.,60,6.4e-6,10,1e-6,10,2e-6,10,3e-6,6.4e-6,0.001}
90        make/o/t parameters_ThreeShell = {"scale","core radius (A)","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)"}
91        Edit parameters_ThreeShell, coef_ThreeShell
92       
93        Variable/G root:g_ThreeShell
94        g_ThreeShell := ThreeShell(coef_ThreeShell, ywave_ThreeShell, xwave_ThreeShell)
95        Display ywave_ThreeShell vs xwave_ThreeShell
96        ModifyGraph marker=29, msize=2, mode=4
97        ModifyGraph log=1,grid=1,mirror=2
98        Label bottom "q (\\S-1\\M) "
99        Label left "I(q) (cm\\S-1\\M)"
100        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
101       
102        AddModelToStrings("ThreeShell","coef_ThreeShell","parameters_ThreeShell","ThreeShell")
103//
104End
105
106Proc PlotFourShell(num,qmin,qmax)
107        Variable num=200, qmin=0.001, qmax=0.7
108        Prompt num "Enter number of data points for model: "
109        Prompt qmin "Enter minimum q-value (^-1) for model: "
110        Prompt qmax "Enter maximum q-value (^-1) for model: "
111//
112        Make/O/D/n=(num) xwave_FourShell, ywave_FourShell
113        xwave_FourShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
114        Make/O/D coef_FourShell ={1.,60,6.4e-6,10,1e-6,10,2e-6,10,3e-6,10,4e-6,6.4e-6,0.001}
115        make/o/t parameters_FourShell = {"scale","core radius (A)","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)"}
116        Edit parameters_FourShell, coef_FourShell
117       
118        Variable/G root:g_FourShell
119        g_FourShell := FourShell(coef_FourShell, ywave_FourShell, xwave_FourShell)
120        Display ywave_FourShell vs xwave_FourShell
121        ModifyGraph marker=29, msize=2, mode=4
122        ModifyGraph log=1,grid=1,mirror=2
123        Label bottom "q (\\S-1\\M) "
124        Label left "I(q) (cm\\S-1\\M)"
125        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
126       
127        AddModelToStrings("FourShell","coef_FourShell","parameters_FourShell","FourShell")
128//
129End
130
131
132
133
134//
135//this macro sets up all the necessary parameters and waves that are
136//needed to calculate the  smeared model function.
137//
138//no input parameters are necessary, it MUST use the experimental q-values
139// from the experimental data read in from an AVE/QSIG data file
140////////////////////////////////////////////////////
141// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
142Proc PlotSmearedOneShell(str)                                                           
143        String str
144        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
145       
146        // if any of the resolution waves are missing => abort
147        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
148                Abort
149        endif
150       
151        SetDataFolder $("root:"+str)
152       
153        // Setup parameter table for model function
154        Make/O/D smear_coef_OneShell =  {1.,60,6.4e-6,10,1e-6,6.4e-6,0.001}
155        make/o/t smear_parameters_OneShell =  {"scale","core radius (A)","Core SLD (A-2)","Shell thickness (A)","Shell SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
156        Edit smear_parameters_OneShell,smear_coef_OneShell                                      //display parameters in a table
157       
158        // output smeared intensity wave, dimensions are identical to experimental QSIG values
159        // make extra copy of experimental q-values for easy plotting
160        Duplicate/O $(str+"_q") smeared_OneShell,smeared_qvals
161        SetScale d,0,0,"1/cm",smeared_OneShell
162                                       
163        Variable/G gs_OneShell=0
164        gs_OneShell := fSmearedOneShell(smear_coef_OneShell,smeared_OneShell,smeared_qvals)     //this wrapper fills the STRUCT
165       
166        Display smeared_OneShell vs smeared_qvals
167        ModifyGraph log=1,marker=29,msize=2,mode=4
168        Label bottom "q (\\S-1\\M)"
169        Label left "I(q) (cm\\S-1\\M)"
170        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
171       
172        SetDataFolder root:
173        AddModelToStrings("SmearedOneShell","smear_coef_OneShell","smear_parameters_OneShell","OneShell")
174End
175
176Proc PlotSmearedTwoShell(str)                                                           
177        String str
178        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
179       
180        // if any of the resolution waves are missing => abort
181        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
182                Abort
183        endif
184       
185        SetDataFolder $("root:"+str)
186       
187        // Setup parameter table for model function
188        Make/O/D smear_coef_TwoShell =  {1.,60,6.4e-6,10,1e-6,10,2e-6,6.4e-6,0.001}
189        make/o/t smear_parameters_TwoShell =  {"scale","core radius (A)","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)"}
190        Edit smear_parameters_TwoShell,smear_coef_TwoShell                                      //display parameters in a table
191       
192        // output smeared intensity wave, dimensions are identical to experimental QSIG values
193        // make extra copy of experimental q-values for easy plotting
194        Duplicate/O $(str+"_q") smeared_TwoShell,smeared_qvals
195        SetScale d,0,0,"1/cm",smeared_TwoShell
196                                       
197        Variable/G gs_TwoShell=0
198        gs_TwoShell := fSmearedTwoShell(smear_coef_TwoShell,smeared_TwoShell,smeared_qvals)     //this wrapper fills the STRUCT
199       
200        Display smeared_TwoShell vs smeared_qvals
201        ModifyGraph log=1,marker=29,msize=2,mode=4
202        Label bottom "q (\\S-1\\M)"
203        Label left "I(q) (cm\\S-1\\M)"
204        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
205       
206        SetDataFolder root:
207        AddModelToStrings("SmearedTwoShell","smear_coef_TwoShell","smear_parameters_TwoShell","TwoShell")
208End
209
210
211Proc PlotSmearedThreeShell(str)                                                         
212        String str
213        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
214       
215        // if any of the resolution waves are missing => abort
216        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
217                Abort
218        endif
219       
220        SetDataFolder $("root:"+str)
221       
222        // Setup parameter table for model function
223        Make/O/D smear_coef_ThreeShell =  {1.,60,6.4e-6,10,1e-6,10,2e-6,10,3e-6,6.4e-6,0.001}
224        make/o/t smear_parameters_ThreeShell =  {"scale","core radius (A)","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)"}
225        Edit smear_parameters_ThreeShell,smear_coef_ThreeShell                                  //display parameters in a table
226       
227        // output smeared intensity wave, dimensions are identical to experimental QSIG values
228        // make extra copy of experimental q-values for easy plotting
229        Duplicate/O $(str+"_q") smeared_ThreeShell,smeared_qvals
230        SetScale d,0,0,"1/cm",smeared_ThreeShell
231                                       
232        Variable/G gs_ThreeShell=0
233        gs_ThreeShell := fSmearedThreeShell(smear_coef_ThreeShell,smeared_ThreeShell,smeared_qvals)     //this wrapper fills the STRUCT
234       
235        Display smeared_ThreeShell vs smeared_qvals
236        ModifyGraph log=1,marker=29,msize=2,mode=4
237        Label bottom "q (\\S-1\\M)"
238        Label left "I(q) (cm\\S-1\\M)"
239        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
240       
241        SetDataFolder root:
242        AddModelToStrings("SmearedThreeShell","smear_coef_ThreeShell","smear_parameters_ThreeShell","ThreeShell")
243End
244
245Proc PlotSmearedFourShell(str)                                                         
246        String str
247        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
248       
249        // if any of the resolution waves are missing => abort
250        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
251                Abort
252        endif
253       
254        SetDataFolder $("root:"+str)
255       
256        // Setup parameter table for model function
257        Make/O/D smear_coef_FourShell = {1.,60,6.4e-6,10,1e-6,10,2e-6,10,3e-6,10,4e-6,6.4e-6,0.001}
258        make/o/t smear_parameters_FourShell =  {"scale","core radius (A)","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)"}
259        Edit smear_parameters_FourShell,smear_coef_FourShell                                    //display parameters in a table
260       
261        // output smeared intensity wave, dimensions are identical to experimental QSIG values
262        // make extra copy of experimental q-values for easy plotting
263        Duplicate/O $(str+"_q") smeared_FourShell,smeared_qvals
264        SetScale d,0,0,"1/cm",smeared_FourShell
265                                       
266        Variable/G gs_FourShell=0
267        gs_FourShell := fSmearedFourShell(smear_coef_FourShell,smeared_FourShell,smeared_qvals) //this wrapper fills the STRUCT
268       
269        Display smeared_FourShell vs smeared_qvals
270        ModifyGraph log=1,marker=29,msize=2,mode=4
271        Label bottom "q (\\S-1\\M)"
272        Label left "I(q) (cm\\S-1\\M)"
273        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
274       
275        SetDataFolder root:
276        AddModelToStrings("SmearedFourShell","smear_coef_FourShell","smear_parameters_FourShell","FourShell")
277End
278
279
280// nothing to change here
281//
282//AAO version, uses XOP if available
283// simply calls the original single point calculation with
284// a wave assignment (this will behave nicely if given point ranges)
285Function OneShell(cw,yw,xw) : FitFunc
286        Wave cw,yw,xw
287       
288#if exists("OneShellX")
289        yw = OneShellX(cw,xw)
290#else
291        yw = fOneShell(cw,xw)
292#endif
293        return(0)
294End
295
296Function TwoShell(cw,yw,xw) : FitFunc
297        Wave cw,yw,xw
298       
299#if exists("TwoShellX")
300        yw = TwoShellX(cw,xw)
301#else
302        yw = fTwoShell(cw,xw)
303#endif
304        return(0)
305End
306
307Function ThreeShell(cw,yw,xw) : FitFunc
308        Wave cw,yw,xw
309       
310#if exists("ThreeShellX")
311        yw =ThreeShellX(cw,xw)
312#else
313        yw = fThreeShell(cw,xw)
314#endif
315        return(0)
316End
317
318Function FourShell(cw,yw,xw) : FitFunc
319        Wave cw,yw,xw
320       
321#if exists("FourShellX")
322        yw = FourShellX(cw,xw)
323#else
324        yw = fFourShell(cw,xw)
325#endif
326        return(0)
327End
328
329
330//
331// unsmeared model calculation
332//
333Function fOneShell(w,x) : FitFunc
334        Wave w
335        Variable x
336       
337        // variables are:
338        //[0] scale factor
339        //[1] radius of core []
340        //[2] SLD of the core   [-2]
341        //[3] thickness of the shell    []
342        //[4] SLD of the shell
343        //[5] SLD of the solvent
344        //[6] background        [cm-1]
345       
346        // All inputs are in ANGSTROMS
347        //OUTPUT is normalized by the particle volume, and converted to [cm-1]
348       
349       
350        Variable scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg
351        scale = w[0]
352        rcore = w[1]
353        rhocore = w[2]
354        thick = w[3]
355        rhoshel = w[4]
356        rhosolv = w[5]
357        bkg = w[6]
358       
359        // calculates scale *( f^2 + bkg)
360        Variable bes,f,vol,qr,contr,f2
361       
362        // core first, then add in shell
363        qr=x*rcore
364        contr = rhocore-rhoshel
365       
366        if(qr == 0)
367                bes = 1
368        else
369                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
370        endif
371       
372        vol = 4*pi/3*rcore^3
373        f = vol*bes*contr
374        //now the shell
375        qr=x*(rcore+thick)
376        contr = rhoshel-rhosolv
377       
378        if(qr == 0)
379                bes = 1
380        else
381                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
382        endif
383       
384        vol = 4*pi/3*(rcore+thick)^3
385        f += vol*bes*contr
386       
387        // normalize to particle volume and rescale from [-1] to [cm-1]
388        f2 = f*f/vol*1.0e8
389       
390        //scale if desired
391        f2 *= scale
392        // then add in the background
393        f2 += bkg
394       
395        return (f2)
396End
397
398
399Function fTwoShell(w,x) : FitFunc
400        Wave w
401        Variable x
402       
403        // variables are:
404        //[0] scale factor
405        //[1] radius of core []
406        //[2] SLD of the core   [-2]
407        //[3] thickness of shell 1 []
408        //[4] SLD of shell 1
409        //[5] thickness of shell 2 []
410        //[6] SLD of shell 2
411        //[7] SLD of the solvent
412        //[8] background        [cm-1]
413       
414        // All inputs are in ANGSTROMS
415        //OUTPUT is normalized by the particle volume, and converted to [cm-1]
416       
417       
418        Variable scale,rcore,thick1,thick2,rhocore,rhoshel1,rhoshel2,rhosolv,bkg
419        scale = w[0]
420        rcore = w[1]
421        rhocore = w[2]
422        thick1 = w[3]
423        rhoshel1 = w[4]
424        thick2 = w[5]
425        rhoshel2 = w[6]
426        rhosolv = w[7]
427        bkg = w[8]
428       
429        // calculates scale *( f^2 + bkg)
430        Variable bes,f,vol,qr,contr,f2
431       
432        // core first, then add in shells
433        qr=x*rcore
434        contr = rhocore-rhoshel1
435        if(qr == 0)
436                bes = 1
437        else
438                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
439        endif
440        vol = 4*pi/3*rcore^3
441        f = vol*bes*contr
442        //now the shell (1)
443        qr=x*(rcore+thick1)
444        contr = rhoshel1-rhoshel2
445        if(qr == 0)
446                bes = 1
447        else
448                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
449        endif
450        vol = 4*pi/3*(rcore+thick1)^3
451        f += vol*bes*contr
452        //now the shell (2)
453        qr=x*(rcore+thick1+thick2)
454        contr = rhoshel2-rhosolv
455        if(qr == 0)
456                bes = 1
457        else
458                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
459        endif
460        vol = 4*pi/3*(rcore+thick1+thick2)^3
461        f += vol*bes*contr
462       
463        // normalize to particle volume and rescale from [-1] to [cm-1]
464        f2 = f*f/vol*1.0e8
465       
466        //scale if desired
467        f2 *= scale
468        // then add in the background
469        f2 += bkg
470       
471        return (f2)
472       
473End
474
475Function fThreeShell(w,x) : FitFunc
476        Wave w
477        Variable x
478       
479        // variables are:
480        //[0] scale factor
481        //[1] radius of core []
482        //[2] SLD of the core   [-2]
483        //[3] thickness of shell 1 []
484        //[4] SLD of shell 1
485        //[5] thickness of shell 2 []
486        //[6] SLD of shell 2
487        //[7] SLD of the solvent
488        //[8] background        [cm-1]
489       
490        // All inputs are in ANGSTROMS
491        //OUTPUT is normalized by the particle volume, and converted to [cm-1]
492       
493       
494        Variable scale,rcore,thick1,thick2,thick3,rhoshel1,rhoshel2,rhoshel3
495        Variable rhocore,rhosolv,bkg
496        scale = w[0]
497        rcore = w[1]
498        rhocore = w[2]
499        thick1 = w[3]
500        rhoshel1 = w[4]
501        thick2 = w[5]
502        rhoshel2 = w[6]
503        thick3 = w[7]
504        rhoshel3 = w[8]
505        rhosolv = w[9]
506        bkg = w[10]
507       
508        // calculates scale *( f^2 + bkg)
509        Variable bes,f,vol,qr,contr,f2
510       
511        // core first, then add in shells
512        qr=x*rcore
513        contr = rhocore-rhoshel1
514        if(qr == 0)
515                bes = 1
516        else
517                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
518        endif
519        vol = 4*pi/3*rcore^3
520        f = vol*bes*contr
521        //now the shell (1)
522        qr=x*(rcore+thick1)
523        contr = rhoshel1-rhoshel2
524        if(qr == 0)
525                bes = 1
526        else
527                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
528        endif
529        vol = 4*pi/3*(rcore+thick1)^3
530        f += vol*bes*contr
531        //now the shell (2)
532        qr=x*(rcore+thick1+thick2)
533        contr = rhoshel2-rhoshel3
534        if(qr == 0)
535                bes = 1
536        else
537                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
538        endif
539        vol = 4*pi/3*(rcore+thick1+thick2)^3
540        f += vol*bes*contr
541        //now the shell (3)
542        qr=x*(rcore+thick1+thick2+thick3)
543        contr = rhoshel3-rhosolv
544        if(qr == 0)
545                bes = 1
546        else
547                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
548        endif
549        vol = 4*pi/3*(rcore+thick1+thick2+thick3)^3
550        f += vol*bes*contr
551       
552        // normalize to particle volume and rescale from [-1] to [cm-1]
553        f2 = f*f/vol*1.0e8
554       
555        //scale if desired
556        f2 *= scale
557        // then add in the background
558        f2 += bkg
559       
560        return (f2)
561End
562
563
564Function fFourShell(w,x) : FitFunc
565        Wave w
566        Variable x
567       
568        // variables are:
569        //[0] scale factor
570        //[1] radius of core []
571        //[2] SLD of the core   [-2]
572        //[3] thickness of shell 1 []
573        //[4] SLD of shell 1
574        //[5] thickness of shell 2 []
575        //[6] SLD of shell 2
576        //[7] SLD of the solvent
577        //[8] background        [cm-1]
578       
579        // All inputs are in ANGSTROMS
580        //OUTPUT is normalized by the particle volume, and converted to [cm-1]
581       
582       
583        Variable scale,rcore,thick1,thick2,thick3,thick4
584        Variable rhoshel1,rhoshel2,rhoshel3,rhoshel4
585        Variable rhocore,rhosolv,bkg
586        scale = w[0]
587        rcore = w[1]
588        rhocore = w[2]
589        thick1 = w[3]
590        rhoshel1 = w[4]
591        thick2 = w[5]
592        rhoshel2 = w[6]
593        thick3 = w[7]
594        rhoshel3 = w[8]
595        thick4 = w[9]
596        rhoshel4 = w[10]
597        rhosolv = w[11]
598        bkg = w[12]
599       
600        // calculates scale *( f^2 + bkg)
601        Variable bes,f,vol,qr,contr,f2
602       
603        // core first, then add in shells
604        qr=x*rcore
605        contr = rhocore-rhoshel1
606        if(qr == 0)
607                bes = 1
608        else
609                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
610        endif
611        vol = 4*pi/3*rcore^3
612        f = vol*bes*contr
613        //now the shell (1)
614        qr=x*(rcore+thick1)
615        contr = rhoshel1-rhoshel2
616        if(qr == 0)
617                bes = 1
618        else
619                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
620        endif
621        vol = 4*pi/3*(rcore+thick1)^3
622        f += vol*bes*contr
623        //now the shell (2)
624        qr=x*(rcore+thick1+thick2)
625        contr = rhoshel2-rhoshel3
626        if(qr == 0)
627                bes = 1
628        else
629                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
630        endif
631        vol = 4*pi/3*(rcore+thick1+thick2)^3
632        f += vol*bes*contr
633        //now the shell (3)
634        qr=x*(rcore+thick1+thick2+thick3)
635        contr = rhoshel3-rhoshel4
636        if(qr == 0)
637                bes = 1
638        else
639                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
640        endif
641        vol = 4*pi/3*(rcore+thick1+thick2+thick3)^3
642        f += vol*bes*contr
643        //now the shell (4)
644        qr=x*(rcore+thick1+thick2+thick3+thick4)
645        contr = rhoshel4-rhosolv
646        if(qr == 0)
647                bes = 1
648        else
649                bes = 3*(sin(qr)-qr*cos(qr))/qr^3
650        endif
651        vol = 4*pi/3*(rcore+thick1+thick2+thick3+thick4)^3
652        f += vol*bes*contr
653       
654       
655        // normalize to particle volume and rescale from [-1] to [cm-1]
656        f2 = f*f/vol*1.0e8
657       
658        //scale if desired
659        f2 *= scale
660        // then add in the background
661        f2 += bkg
662       
663        return (f2)
664End
665
666
667
668
669///////////////////////////////////////////////////////////////
670// smeared model calculation
671//
672// you don't need to do anything with this function, as long as
673// your OneShell works correctly, you get the resolution-smeared
674// version for free.
675//
676// this is all there is to the smeared model calculation!
677Function SmearedOneShell(s) : FitFunc
678        Struct ResSmearAAOStruct &s
679
680//      the name of your unsmeared model (AAO) is the first argument
681        Smear_Model_20(OneShell,s.coefW,s.xW,s.yW,s.resW)
682
683        return(0)
684End
685
686Function SmearedTwoShell(s) : FitFunc
687        Struct ResSmearAAOStruct &s
688
689//      the name of your unsmeared model (AAO) is the first argument
690        Smear_Model_20(TwoShell,s.coefW,s.xW,s.yW,s.resW)
691
692        return(0)
693End
694
695Function SmearedThreeShell(s) : FitFunc
696        Struct ResSmearAAOStruct &s
697
698//      the name of your unsmeared model (AAO) is the first argument
699        Smear_Model_20(ThreeShell,s.coefW,s.xW,s.yW,s.resW)
700
701        return(0)
702End
703
704Function SmearedFourShell(s) : FitFunc
705        Struct ResSmearAAOStruct &s
706
707//      the name of your unsmeared model (AAO) is the first argument
708        Smear_Model_20(FourShell,s.coefW,s.xW,s.yW,s.resW)
709
710        return(0)
711End
712
713
714
715///////////////////////////////////////////////////////////////
716
717
718// nothing to change here
719//
720//wrapper to calculate the smeared model as an AAO-Struct
721// fills the struct and calls the ususal function with the STRUCT parameter
722//
723// used only for the dependency, not for fitting
724//
725Function fSmearedOneShell(coefW,yW,xW)
726        Wave coefW,yW,xW
727       
728        String str = getWavesDataFolder(yW,0)
729        String DF="root:"+str+":"
730       
731        WAVE resW = $(DF+str+"_res")
732       
733        STRUCT ResSmearAAOStruct fs
734        WAVE fs.coefW = coefW   
735        WAVE fs.yW = yW
736        WAVE fs.xW = xW
737        WAVE fs.resW = resW
738       
739        Variable err
740        err = SmearedOneShell(fs)
741       
742        return (0)
743End
744
745//wrapper to calculate the smeared model as an AAO-Struct
746// fills the struct and calls the ususal function with the STRUCT parameter
747//
748// used only for the dependency, not for fitting
749//
750Function fSmearedTwoShell(coefW,yW,xW)
751        Wave coefW,yW,xW
752       
753        String str = getWavesDataFolder(yW,0)
754        String DF="root:"+str+":"
755       
756        WAVE resW = $(DF+str+"_res")
757       
758        STRUCT ResSmearAAOStruct fs
759        WAVE fs.coefW = coefW   
760        WAVE fs.yW = yW
761        WAVE fs.xW = xW
762        WAVE fs.resW = resW
763       
764        Variable err
765        err = SmearedTwoShell(fs)
766       
767        return (0)
768End
769
770//wrapper to calculate the smeared model as an AAO-Struct
771// fills the struct and calls the ususal function with the STRUCT parameter
772//
773// used only for the dependency, not for fitting
774//
775Function fSmearedThreeShell(coefW,yW,xW)
776        Wave coefW,yW,xW
777       
778        String str = getWavesDataFolder(yW,0)
779        String DF="root:"+str+":"
780       
781        WAVE resW = $(DF+str+"_res")
782       
783        STRUCT ResSmearAAOStruct fs
784        WAVE fs.coefW = coefW   
785        WAVE fs.yW = yW
786        WAVE fs.xW = xW
787        WAVE fs.resW = resW
788       
789        Variable err
790        err = SmearedThreeShell(fs)
791       
792        return (0)
793End
794
795//wrapper to calculate the smeared model as an AAO-Struct
796// fills the struct and calls the ususal function with the STRUCT parameter
797//
798// used only for the dependency, not for fitting
799//
800Function fSmearedFourShell(coefW,yW,xW)
801        Wave coefW,yW,xW
802       
803        String str = getWavesDataFolder(yW,0)
804        String DF="root:"+str+":"
805       
806        WAVE resW = $(DF+str+"_res")
807       
808        STRUCT ResSmearAAOStruct fs
809        WAVE fs.coefW = coefW   
810        WAVE fs.yW = yW
811        WAVE fs.xW = xW
812        WAVE fs.resW = resW
813       
814        Variable err
815        err = SmearedFourShell(fs)
816       
817        return (0)
818End
Note: See TracBrowser for help on using the repository browser.