source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/FCC_ParaCrystal_v40.ipf @ 455

Last change on this file since 455 was 455, checked in by srkline, 14 years ago

Two changes: (1) remove threading of Cyl-PolyRad? model until WM can fix the crashing bug related to compiling/threadsafe...
(2) fixed bug in the calculation of the random deviate for MC calculations. The wrong scattering cross section was calculated depending on wavelength. This has been corrected.

File size: 8.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.0
3
4//
5//
6// FCC paracrystal, powder average
7//
8// VERY slow, since the function is so ill-behaved and needs LOTS of quadrature
9// points. Adaptive methods were even slower and troublesom to converge,
10// although in theory they should be a better choice than blindly increasing the number of points.
11//
12// using 150 points for the quadrature
13// 76 points for the smearing is untested
14//
15//
16// Original implementation - Danilo Pozzo
17//              modified and modernized for more efficient integration SRK Nov 2008
18//
19//REFERENCE
20//Hideki Matsuoka etal. Physical Review B, Vol 36 Num 3, p1754 1987   ORIGINAL PAPER
21//Hideki Matsuoka etal. Physical Review B, Vol 41 Num 6, p3854 1990   CORRECTIONS TO PAPER
22//
23////////////////////////////////////////////////////
24
25
26
27Proc PlotFCC_ParaCrystal(num,qmin,qmax)
28        Variable num=100, qmin=0.001, qmax=0.7
29        Prompt num "Enter number of data points for model: "
30        Prompt qmin "Enter minimum q-value (^-1) for model: "
31        Prompt qmax "Enter maximum q-value (^-1) for model: "
32//
33        Make/O/D/n=(num) xwave_FCC_ParaCrystal, ywave_FCC_ParaCrystal
34        xwave_FCC_ParaCrystal =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
35        Make/O/D coef_FCC_ParaCrystal = {1,220,0.06,40,3e-6,6.3e-6,0.0}
36        make/o/t parameters_FCC_ParaCrystal = {"scale","Nearest Neighbor (A)","distortion, g","Sphere Radius (A)","SLD sphere (A-2)","SLD solvent (A-2)", "Background (cm-1)"} 
37        Edit parameters_FCC_ParaCrystal, coef_FCC_ParaCrystal
38       
39        Variable/G root:gNordFCC=150
40       
41        Variable/G root:g_FCC_ParaCrystal
42        g_FCC_ParaCrystal := FCC_ParaCrystal(coef_FCC_ParaCrystal, ywave_FCC_ParaCrystal, xwave_FCC_ParaCrystal)
43        Display ywave_FCC_ParaCrystal vs xwave_FCC_ParaCrystal
44        ModifyGraph marker=29, msize=2, mode=4
45        ModifyGraph grid=1,mirror=2
46        ModifyGraph log=0
47        Label bottom "q (\\S-1\\M) "
48        Label left "I(q) (cm\\S-1\\M)"
49        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
50       
51        AddModelToStrings("FCC_ParaCrystal","coef_FCC_ParaCrystal","FCC_ParaCrystal")
52//
53End
54
55//
56//this macro sets up all the necessary parameters and waves that are
57//needed to calculate the  smeared model function.
58//
59//no input parameters are necessary, it MUST use the experimental q-values
60// from the experimental data read in from an AVE/QSIG data file
61////////////////////////////////////////////////////
62// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
63Proc PlotSmearedFCC_ParaCrystal(str)                                                           
64        String str
65        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
66       
67        // if any of the resolution waves are missing => abort
68        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
69                Abort
70        endif
71       
72        SetDataFolder $("root:"+str)
73       
74        // Setup parameter table for model function
75        Make/O/D smear_coef_FCC_ParaCrystal = {1,220,0.06,40,3e-6,6.3e-6,0.0}
76        make/o/t smear_parameters_FCC_ParaCrystal = {"scale","Nearest Neighbor (A)","distortion, g","Sphere Radius (A)","SLD sphere (A-2)","SLD solvent (A-2)", "Background (cm-1)"}
77        Edit smear_parameters_FCC_ParaCrystal,smear_coef_FCC_ParaCrystal                                        //display parameters in a table
78       
79        // output smeared intensity wave, dimensions are identical to experimental QSIG values
80        // make extra copy of experimental q-values for easy plotting
81        Duplicate/O $(str+"_q") smeared_FCC_ParaCrystal,smeared_qvals
82        SetScale d,0,0,"1/cm",smeared_FCC_ParaCrystal
83       
84        Variable/G gNordFCC = 150               
85        Variable/G gs_FCC_ParaCrystal=0
86        gs_FCC_ParaCrystal := fSmearedFCC_ParaCrystal(smear_coef_FCC_ParaCrystal,smeared_FCC_ParaCrystal,smeared_qvals) //this wrapper fills the STRUCT
87       
88        Display smeared_FCC_ParaCrystal vs smeared_qvals
89        ModifyGraph marker=29,msize=2,mode=4
90        ModifyGraph log=0
91        Label bottom "q (\\S-1\\M)"
92        Label left "I(q) (cm\\S-1\\M)"
93        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
94       
95        SetDataFolder root:
96        AddModelToStrings("SmearedFCC_ParaCrystal","smear_coef_FCC_ParaCrystal","FCC_ParaCrystal")
97End
98
99
100// nothing to change here
101//
102//AAO version, uses XOP if available
103// simply calls the original single point calculation with
104// a wave assignment (this will behave nicely if given point ranges)
105Function FCC_ParaCrystal(cw,yw,xw) : FitFunc
106        Wave cw,yw,xw
107       
108#if exists("FCC_ParaCrystalX")
109        yw = FCC_ParaCrystalX(cw,xw)
110#else
111        yw = fFCC_ParaCrystal(cw,xw)
112#endif
113        return(0)
114End
115
116
117//
118// unsmeared model calculation
119//
120Function fFCC_ParaCrystal(w,x) : FitFunc
121        Wave w
122        Variable x
123       
124//       Input (fitting) variables are not used
125//      you would give them nice names
126        Variable integral,loLim,upLim
127        loLim = 0
128        upLim = 2*Pi
129       
130        Variable/G root:gDumY=0         //root:gDumX=0
131       
132
133        Variable scale,Dnn,gg,Rad,contrast,background,yy,latticeScale,sld,sldSolv
134        scale = w[0]
135        Dnn = w[1] //Nearest neighbor distance A
136        gg = w[2] //Paracrystal distortion factor
137        Rad = w[3] //Sphere radius
138        sld = w[4]
139        sldSolv = w[5]
140        background = w[6]
141               
142        contrast = sld - sldSolv
143       
144        latticeScale = 4*(4/3)*pi*(Rad^3)/((Dnn*(2^0.5))^3)
145
146        NVAR/Z nord=root:gNordFCC
147        if(NVAR_Exists(nord)!=1)
148                nord=20
149        endif
150       
151
152        integral = IntegrateFn_N(Integrand_FCC_Outer,loLim,upLim,w,x,nord)
153       
154       
155        integral *= SphereForm_FCC(Rad,contrast,x)*scale*latticeScale
156//      integral *= scale               //for testing, returns Z(q) only
157
158        integral += background 
159       
160        Return (integral)
161       
162End
163
164// the outer integral is also an integral
165Function Integrand_FCC_Outer(w,x,dum)
166        Wave w
167        Variable x,dum
168               
169        NVAR yy = root:gDumY           
170        yy = dum                                        // save the current dummy yy for use in the inner loop
171        Variable retVal,loLim,upLim
172        //
173        loLim = 0
174        upLim = Pi
175
176        NVAR/Z nord=root:gNordFCC
177        if(NVAR_Exists(nord)!=1)
178                nord=20
179        endif
180       
181        retVal = IntegrateFn_N(Integrand_FCC_Inner,loLim,upLim,w,x,nord)
182       
183        return(retVal)
184End
185
186//returns the value of the integrand of the inner integral
187Function Integrand_FCC_Inner(w,qq,dum)
188        Wave w
189        Variable qq,dum
190       
191        NVAR yy = root:gDumY            //use the yy value from the outer loop
192        Variable xx,retVal
193        xx = dum
194
195        retVal = FCC_Integrand(w,qq,xx,yy)
196       
197        return(retVal)
198End
199
200Function FCC_Integrand(w,qq,xx,yy)
201        Wave w
202        Variable qq,xx,yy
203       
204        Variable retVal,temp1,temp3,aa,Da,Dnn,gg
205        Dnn = w[1] //Nearest neighbor distance A
206        gg = w[2] //Paracrystal distortion factor
207//      aa = Dnn*(2^0.5)                //Danilo's version. As defined in paper, |bi| = a
208        aa = Dnn
209        Da = gg*aa
210       
211        temp1 = qq*qq*Da*Da
212        temp3 = qq*aa
213       
214        retVal = FCCeval(xx,yy,temp1,temp3)
215        retVal /=4*Pi
216       
217        return(retVal)
218end
219
220Function FCCeval(Theta,Phi,temp1,temp3)
221        Variable Theta,Phi,temp1,temp3
222        Variable temp6,temp7,temp8,temp9,temp10
223        Variable result
224       
225        temp6 = sin(Theta)
226        temp7 = sin(Theta)*sin(Phi)+cos(Theta)
227        temp8 = -1*sin(Theta)*cos(Phi)+cos(Theta)
228        temp9 = -1*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)
229        temp10 = exp((-1/8)*temp1*((temp7^2)+(temp8^2)+(temp9^2)))
230        result = ((1-(temp10^2))^3)*temp6/((1-2*temp10*cos(0.5*temp3*(temp7))+(temp10^2))*(1-2*temp10*cos(0.5*temp3*(temp8))+(temp10^2))*(1-2*temp10*cos(0.5*temp3*(temp9))+(temp10^2)))
231       
232        return (result)
233end
234
235Function SphereForm_FCC(radius,delrho,x)                                       
236        Variable radius,delrho,x
237       
238        // variables are:                                                       
239        //[2] radius ()
240        //[3] delrho (-2)
241        //[4] background (cm-1)
242       
243        // calculates scale * f^2/Vol where f=Vol*3*delrho*(sin(qr)-qrcos(qr))/qr^3
244        // and is rescaled to give [=] cm^-1
245       
246        Variable bes,f,vol,f2
247        ////handle q==0 separately
248        If(x==0)
249                f = 4/3*pi*radius^3*delrho*delrho*1e8
250                return(f)
251        Endif
252       
253        bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
254        vol = 4*pi/3*radius^3
255        f = vol*bes*delrho              // [=] 
256        // normalize to single particle volume, convert to 1/cm
257        f2 = f * f / vol * 1.0e8                // [=] 1/cm
258       
259        return (f2)     
260       
261End
262
263
264
265
266///////////////////////////////////////////////////////////////
267// smeared model calculation
268//
269Function SmearedFCC_ParaCrystal(s) : FitFunc
270        Struct ResSmearAAOStruct &s
271
272//      the name of your unsmeared model (AAO) is the first argument
273        Smear_Model_20(FCC_ParaCrystal,s.coefW,s.xW,s.yW,s.resW)
274
275        return(0)
276End
277
278
279/////////////////////////////////////////////////////////////////
280//wrapper to calculate the smeared model as an AAO-Struct
281// fills the struct and calls the ususal function with the STRUCT parameter
282//
283// used only for the dependency, not for fitting
284//
285Function fSmearedFCC_ParaCrystal(coefW,yW,xW)
286        Wave coefW,yW,xW
287       
288        String str = getWavesDataFolder(yW,0)
289        String DF="root:"+str+":"
290       
291        WAVE resW = $(DF+str+"_res")
292       
293        STRUCT ResSmearAAOStruct fs
294        WAVE fs.coefW = coefW   
295        WAVE fs.yW = yW
296        WAVE fs.xW = xW
297        WAVE fs.resW = resW
298       
299        Variable err
300        err = SmearedFCC_ParaCrystal(fs)
301       
302        return (0)
303End
Note: See TracBrowser for help on using the repository browser.