source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2009/FuzzySpheres_v40.ipf @ 656

Last change on this file since 656 was 656, checked in by ajj, 13 years ago

Update to CatVSTable to show sample position, TISANE to fix a couple of startup bugs and FuzzySpheres? to add the lorentzian term.

File size: 7.0 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4
5// plots the form factor of spheres with a Gaussian radius distribution
6// and a fuzzy surface
7//
8// M. Stieger, J. S. Pedersen, P. Lindner, W. Richtering, Langmuir 20 (2004) 7283-7292.
9// M. Stieger, W. Richtering, J. S. Pedersen, P. Lindner, Journal of Chemical Physics 120(13) (2004) 6197-6206.
10//
11// potentially a lorentzian could be added to the low Q, if absolutely necessary
12//
13// SRK JUL 2009
14//
15// Include lorentzian term for *high* Q component of the scattering.
16//
17// AJJ Feb 2010
18
19#include "Lorentz_model_v40"
20
21Proc PlotFuzzySpheres(num,qmin,qmax)
22        Variable num=128,qmin=0.001,qmax=0.7
23        Prompt num "Enter number of data points for model: "
24        Prompt qmin "Enter minimum q-value (A^-1) for model: "
25        Prompt qmax "Enter maximum q-value (A^-1) for model: "
26       
27        Make/O/D/N=(num) xwave_fuzz,ywave_fuzz
28        xwave_fuzz = alog( log(qmin) + x*((log(qmax)-log(qmin))/num) )
29        Make/O/D coef_fuzz = {0.01,60,0.2,10,1e-6,3e-6,1,50,0.001}
30        make/O/T parameters_fuzz = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","interface thickness (A)","SLD sphere (A-2)","SLD solvent (A-2)","Lorentz Scale","Lorentz length","bkg (cm-1 sr-1)"}
31        Edit parameters_fuzz,coef_fuzz
32       
33        Variable/G root:g_fuzz
34        g_fuzz := FuzzySpheres(coef_fuzz,ywave_fuzz,xwave_fuzz)
35        Display ywave_fuzz vs xwave_fuzz
36        ModifyGraph log=1,marker=29,msize=2,mode=4
37        Label bottom "q (A\\S-1\\M)"
38        Label left "Intensity (cm\\S-1\\M)"
39        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
40       
41        AddModelToStrings("FuzzySpheres","coef_fuzz","parameters_fuzz","fuzz")
42End
43
44// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
45Proc PlotSmearedFuzzySpheres(str)                                                               
46        String str
47        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
48       
49        // if any of the resolution waves are missing => abort
50        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
51                Abort
52        endif
53       
54        SetDataFolder $("root:"+str)
55       
56        // Setup parameter table for model function
57        Make/O/D smear_coef_fuzz = {0.01,60,0.2,10,1e-6,3e-6,1,50,0.001}                                       
58        make/o/t smear_parameters_fuzz = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","interface thickness (A)","SLD sphere (A-2)","SLD solvent (A-2)","Lorentz Scale","Lorentz length","bkg (cm-1 sr-1)"}
59        Edit smear_parameters_fuzz,smear_coef_fuzz                                     
60       
61        // output smeared intensity wave, dimensions are identical to experimental QSIG values
62        // make extra copy of experimental q-values for easy plotting
63        Duplicate/O $(str+"_q") smeared_fuzz,smeared_qvals                             
64        SetScale d,0,0,"1/cm",smeared_fuzz                                                     
65                                       
66        Variable/G gs_fuzz=0
67        gs_fuzz := fSmearedFuzzySpheres(smear_coef_fuzz,smeared_fuzz,smeared_qvals)     //this wrapper fills the STRUCT
68       
69        Display smeared_fuzz vs smeared_qvals                                                                   
70        ModifyGraph log=1,marker=29,msize=2,mode=4
71        Label bottom "q (A\\S-1\\M)"
72        Label left "Intensity (cm\\S-1\\M)"
73        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
74       
75        SetDataFolder root:
76        AddModelToStrings("SmearedFuzzySpheres","smear_coef_fuzz","smear_parameters_fuzz","fuzz")
77End
78       
79
80
81
82//AAO version, uses XOP if available
83// simply calls the original single point calculation with
84// a wave assignment (this will behave nicely if given point ranges)
85Function FuzzySpheres(cw,yw,xw) : FitFunc
86        Wave cw,yw,xw
87       
88#if exists("FuzzySpheresX")
89        yw = FuzzySpheresX(cw,xw)
90#else
91        yw = fFuzzySpheres(cw,xw)
92#endif
93        return(0)
94End
95
96Function fFuzzySpheres(w,xx) : FitFunc
97        wave w
98        variable xx
99       
100        Variable scale,rad,pd,sig,rho,rhos,bkg,delrho,sig_surf,lor_sf,lor_len
101       
102        //set up the coefficient values
103        scale=w[0]
104        rad=w[1]
105        pd=w[2]
106        sig=pd*rad
107        sig_surf = w[3]
108        rho=w[4]
109        rhos=w[5]
110        delrho=rho-rhos
111        bkg=w[8]
112
113       
114        //could use 5 pt quadrature to integrate over the size distribution, since it's a gaussian
115        //currently using 20 pts...
116        Variable va,vb,ii,zi,nord,yy,summ,inten
117        Variable bes,f,vol,f2
118        String weightStr,zStr
119       
120        //select number of gauss points by setting nord=20 or76 points
121        nord = 20
122//      nord = 76
123       
124        weightStr = "gauss"+num2str(nord)+"wt"
125        zStr = "gauss"+num2str(nord)+"z"
126       
127        if (WaveExists($weightStr) == 0) // wave reference is not valid,
128                Make/D/N=(nord) $weightStr,$zStr
129                Wave gauWt = $weightStr
130                Wave gauZ = $zStr               // wave references to pass
131                if(nord==20)
132                        Make20GaussPoints(gauWt,gauZ)
133                else
134                        Make76GaussPoints(gauWt,gauZ)
135                endif   
136        else
137                if(exists(weightStr) > 1)
138                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
139                endif
140                Wave gauWt = $weightStr
141                Wave gauZ = $zStr               // create the wave references
142        endif
143       
144        // end points of integration
145        // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero
146        // +/- 3 sigq catches 99.73% of distrubution
147        // change limits (and spacing of zi) at each evaluation based on R()
148        //integration from va to vb
149        va = -4*sig + rad
150        if (va<0)
151                va=0            //to avoid numerical error when  va<0 (-ve q-value)
152        endif
153        vb = 4*sig +rad
154       
155        summ = 0.0              // initialize integral
156        for(ii=0;ii<nord;ii+=1)
157                // calculate Gauss points on integration interval (r-value for evaluation)
158                zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
159               
160                // calculate sphere scattering
161                //
162                //handle q==0 separately
163                If(xx==0)
164                        f2 = 4/3*pi*zi^3*delrho*delrho*1e8
165                        f2 *= exp(-0.5*sig_surf*sig_surf*xx*xx)
166                        f2 *= exp(-0.5*sig_surf*sig_surf*xx*xx)
167                else
168                        bes = 3*(sin(xx*zi)-xx*zi*cos(xx*zi))/xx^3/zi^3
169                        vol = 4*pi/3*zi^3
170                        f = vol*bes*delrho              // [=] A
171                        f *= exp(-0.5*sig_surf*sig_surf*xx*xx)
172                        // normalize to single particle volume, convert to 1/cm
173                        f2 = f * f / vol * 1.0e8                // [=] 1/cm
174                endif
175       
176                yy = gauWt[ii] *  Gauss_f_distr(sig,rad,zi) * f2
177                yy *= 4*pi/3*zi*zi*zi           //un-normalize by current sphere volume
178               
179                summ += yy              //add to the running total of the quadrature
180        endfor
181// calculate value of integral to return
182        inten = (vb-va)/2.0*summ
183       
184        //re-normalize by polydisperse sphere volume
185        inten /= (4*pi/3*rad^3)*(1+3*pd^2)
186       
187        inten *= scale
188       
189        //Lorentzian term
190        Make/O/N=3 tmp_lor
191        tmp_lor[0] = w[6]
192        tmp_lor[1] = w[7]
193        tmp_lor[2] = 0
194       
195        inten+=fLorentz_model(tmp_lor,xx)
196       
197        inten+=bkg
198       
199        Return(inten)
200End
201
202//wrapper to calculate the smeared model as an AAO-Struct
203// fills the struct and calls the ususal function with the STRUCT parameter
204//
205// used only for the dependency, not for fitting
206//
207Function fSmearedFuzzySpheres(coefW,yW,xW)
208        Wave coefW,yW,xW
209       
210        String str = getWavesDataFolder(yW,0)
211        String DF="root:"+str+":"
212       
213        WAVE resW = $(DF+str+"_res")
214       
215        STRUCT ResSmearAAOStruct fs
216        WAVE fs.coefW = coefW   
217        WAVE fs.yW = yW
218        WAVE fs.xW = xW
219        WAVE fs.resW = resW
220       
221        Variable err
222        err = SmearedFuzzySpheres(fs)
223       
224        return (0)
225End
226
227// this is all there is to the smeared calculation!
228Function SmearedFuzzySpheres(s) :FitFunc
229        Struct ResSmearAAOStruct &s
230
231//      the name of your unsmeared model (AAO) is the first argument
232        Smear_Model_20(FuzzySpheres,s.coefW,s.xW,s.yW,s.resW)
233
234        return(0)
235End
236
237
238
239Function Gauss_f_distr(sig,avg,pt)
240        Variable sig,avg,pt
241       
242        Variable retval
243       
244        retval = (1/ ( sig*sqrt(2*Pi)) )*exp(-(avg-pt)^2/sig^2/2)
245        return(retval)
246End
Note: See TracBrowser for help on using the repository browser.