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

Last change on this file since 679 was 679, checked in by srkline, 12 years ago

added Simulation help file to the list of those that are to be killed before installation.

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 high 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.