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

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

fixed the P*S combinations in Fuzzy Spheres.

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, calculate here, removing the #include
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        Variable I0, L
102        I0 = w[0]
103        L = w[1]
104
105        //set up the coefficient values
106        scale=w[0]
107        rad=w[1]
108        pd=w[2]
109        sig=pd*rad
110        sig_surf = w[3]
111        rho=w[4]
112        rhos=w[5]
113        delrho=rho-rhos
114        I0 = w[6]               //for the Lorentzian
115        L = w[7]
116        bkg=w[8]
117
118       
119        //could use 5 pt quadrature to integrate over the size distribution, since it's a gaussian
120        //currently using 20 pts...
121        Variable va,vb,ii,zi,nord,yy,summ,inten
122        Variable bes,f,vol,f2
123        String weightStr,zStr
124       
125        //select number of gauss points by setting nord=20 or76 points
126        nord = 20
127//      nord = 76
128       
129        weightStr = "gauss"+num2str(nord)+"wt"
130        zStr = "gauss"+num2str(nord)+"z"
131       
132        if (WaveExists($weightStr) == 0) // wave reference is not valid,
133                Make/D/N=(nord) $weightStr,$zStr
134                Wave gauWt = $weightStr
135                Wave gauZ = $zStr               // wave references to pass
136                if(nord==20)
137                        Make20GaussPoints(gauWt,gauZ)
138                else
139                        Make76GaussPoints(gauWt,gauZ)
140                endif   
141        else
142                if(exists(weightStr) > 1)
143                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
144                endif
145                Wave gauWt = $weightStr
146                Wave gauZ = $zStr               // create the wave references
147        endif
148       
149        // end points of integration
150        // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero
151        // +/- 3 sigq catches 99.73% of distrubution
152        // change limits (and spacing of zi) at each evaluation based on R()
153        //integration from va to vb
154        va = -4*sig + rad
155        if (va<0)
156                va=0            //to avoid numerical error when  va<0 (-ve q-value)
157        endif
158        vb = 4*sig +rad
159       
160        summ = 0.0              // initialize integral
161        for(ii=0;ii<nord;ii+=1)
162                // calculate Gauss points on integration interval (r-value for evaluation)
163                zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
164               
165                // calculate sphere scattering
166                //
167                //handle q==0 separately
168                If(xx==0)
169                        f2 = 4/3*pi*zi^3*delrho*delrho*1e8
170                        f2 *= exp(-0.5*sig_surf*sig_surf*xx*xx)
171                        f2 *= exp(-0.5*sig_surf*sig_surf*xx*xx)
172                else
173                        bes = 3*(sin(xx*zi)-xx*zi*cos(xx*zi))/xx^3/zi^3
174                        vol = 4*pi/3*zi^3
175                        f = vol*bes*delrho              // [=] A
176                        f *= exp(-0.5*sig_surf*sig_surf*xx*xx)
177                        // normalize to single particle volume, convert to 1/cm
178                        f2 = f * f / vol * 1.0e8                // [=] 1/cm
179                endif
180       
181                yy = gauWt[ii] *  Gauss_f_distr(sig,rad,zi) * f2
182                yy *= 4*pi/3*zi*zi*zi           //un-normalize by current sphere volume
183               
184                summ += yy              //add to the running total of the quadrature
185        endfor
186// calculate value of integral to return
187        inten = (vb-va)/2.0*summ
188       
189        //re-normalize by polydisperse sphere volume
190        inten /= (4*pi/3*rad^3)*(1+3*pd^2)
191       
192        inten *= scale
193       
194        //Lorentzian term       
195        inten += I0/(1 + (xx*L)^2)
196       
197       
198        inten+=bkg
199       
200        Return(inten)
201End
202
203//wrapper to calculate the smeared model as an AAO-Struct
204// fills the struct and calls the ususal function with the STRUCT parameter
205//
206// used only for the dependency, not for fitting
207//
208Function fSmearedFuzzySpheres(coefW,yW,xW)
209        Wave coefW,yW,xW
210       
211        String str = getWavesDataFolder(yW,0)
212        String DF="root:"+str+":"
213       
214        WAVE resW = $(DF+str+"_res")
215       
216        STRUCT ResSmearAAOStruct fs
217        WAVE fs.coefW = coefW   
218        WAVE fs.yW = yW
219        WAVE fs.xW = xW
220        WAVE fs.resW = resW
221       
222        Variable err
223        err = SmearedFuzzySpheres(fs)
224       
225        return (0)
226End
227
228// this is all there is to the smeared calculation!
229Function SmearedFuzzySpheres(s) :FitFunc
230        Struct ResSmearAAOStruct &s
231
232//      the name of your unsmeared model (AAO) is the first argument
233        Smear_Model_20(FuzzySpheres,s.coefW,s.xW,s.yW,s.resW)
234
235        return(0)
236End
237
238
239
240Function Gauss_f_distr(sig,avg,pt)
241        Variable sig,avg,pt
242       
243        Variable retval
244       
245        retval = (1/ ( sig*sqrt(2*Pi)) )*exp(-(avg-pt)^2/sig^2/2)
246        return(retval)
247End
Note: See TracBrowser for help on using the repository browser.