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

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

Added model documentation for FuzzySpheres? and CSParallelpiped

Added axis labels to data loaded through XML loader

Writing of XML USANS data is now properly handling cursors

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 (A)","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 (A)","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.