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

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

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

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","parameters_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_param_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_param_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","smear_param_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.