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

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

adding fuzzy sphere model

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