source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/BCC_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: 9.7 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.0
3
4//
5//
6// BCC 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// 150 points seems to give reasonable reproduction of the peak heights in the paper.
13// peak locations are correct
14// 76 points of quadrature for the smearing is only a guess, it's not been tested yet.
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 PlotBCC_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_BCC_ParaCrystal, ywave_BCC_ParaCrystal
34        xwave_BCC_ParaCrystal =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
35        Make/O/D coef_BCC_ParaCrystal = {1,220,0.06,40,3e-6,6.3e-6,0.0}
36        make/o/t parameters_BCC_ParaCrystal = {"scale","Nearest Neighbor (A)","distortion, g","Sphere Radius (A)","SLD Sphere (A-2)","SLD Solvent (A-2)", "Background (cm-1)"} 
37        Edit parameters_BCC_ParaCrystal, coef_BCC_ParaCrystal
38       
39        Variable/G root:gNordBCC=150
40       
41        Variable/G root:g_BCC_ParaCrystal
42        g_BCC_ParaCrystal := BCC_ParaCrystal(coef_BCC_ParaCrystal, ywave_BCC_ParaCrystal, xwave_BCC_ParaCrystal)
43        Display ywave_BCC_ParaCrystal vs xwave_BCC_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("BCC_ParaCrystal","coef_BCC_ParaCrystal","parameters_BCC_ParaCrystal","BCC_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 PlotSmearedBCC_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_BCC_ParaCrystal = {1,220,0.06,40,3e-6,6.3e-6,0.0}
76        make/o/t smear_param_BCC_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_BCC_ParaCrystal,smear_coef_BCC_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_BCC_ParaCrystal,smeared_qvals
82        SetScale d,0,0,"1/cm",smeared_BCC_ParaCrystal
83       
84        Variable/G gNordBCC = 150       
85        Variable/G gs_BCC_ParaCrystal=0
86        gs_BCC_ParaCrystal := fSmearedBCC_ParaCrystal(smear_coef_BCC_ParaCrystal,smeared_BCC_ParaCrystal,smeared_qvals) //this wrapper fills the STRUCT
87       
88        Display smeared_BCC_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("SmearedBCC_ParaCrystal","smear_coef_BCC_ParaCrystal","smear_param_BCC_ParaCrystal","BCC_ParaCrystal")
97End
98
99// Threaded version
100// Threaded XOP = 2.4 s
101// non-threaded, non-XOP = 46.8 s
102// = x 19.5 speedup !
103//
104Function BCC_ParaCrystal(cw,yw,xw) : FitFunc
105        Wave cw,yw,xw
106
107/////// NO threading /////////
108//#if exists("BCC_ParaCrystalX")
109//      yw = BCC_ParaCrystalX(cw,xw)
110//#else
111//      yw = fBCC_ParaCrystal(cw,xw)
112//#endif
113
114
115/// THREADING ///////
116
117//      Variable t1=StopMSTimer(-2)
118               
119#if exists("BCC_ParaCrystalX")
120
121        Variable npt=numpnts(yw)
122        Variable i,nthreads= ThreadProcessorCount
123        variable mt= ThreadGroupCreate(nthreads)
124
125        for(i=0;i<nthreads;i+=1)
126        //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
127                ThreadStart mt,i,BCC_ParaCrystal_T(cw,yw,xw,(i*npt/nthreads),((i+1)*npt/nthreads-1))
128        endfor
129
130        do
131                variable tgs= ThreadGroupWait(mt,100)
132        while( tgs != 0 )
133
134        variable dummy= ThreadGroupRelease(mt)
135       
136#else
137                yw = fBCC_ParaCrystal(cw,xw)            //the Igor, non-XOP, non-threaded calculation, messy to make ThreadSafe
138#endif
139
140//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
141
142
143
144        return(0)
145End
146
147
148
149// nothing to change here
150//
151//AAO version, uses XOP if available
152// simply calls the original single point calculation with
153// a wave assignment (this will behave nicely if given point ranges)
154//
155// Threaded Version
156ThreadSafe Function BCC_ParaCrystal_T(cw,yw,xw,p1,p2) : FitFunc
157        Wave cw,yw,xw
158        Variable p1,p2
159
160//      Variable t1=StopMSTimer(-2)
161
162#if exists("BCC_ParaCrystalX")
163        yw[p1,p2] = BCC_ParaCrystalX(cw,xw)
164#else
165        yw[p1,p2] = fBCC_ParaCrystal(cw,xw)
166#endif
167
168//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
169
170        return(0)
171End
172
173
174//
175// unsmeared model calculation
176//
177Function fBCC_ParaCrystal(w,x) : FitFunc
178        Wave w
179        Variable x
180       
181//       Input (fitting) variables are not used
182//      you would give them nice names
183        Variable integral,loLim,upLim
184        loLim = 0
185        upLim = 2*Pi
186       
187        Variable/G root:gDumY=0         //root:gDumX=0
188       
189
190        Variable scale,Dnn,gg,Rad,contrast,background,yy,latticeScale,sld,sldSolv
191        scale = w[0]
192        Dnn = w[1] //Nearest neighbor distance A
193        gg = w[2] //Paracrystal distortion factor
194        Rad = w[3] //Sphere radius
195        sld = w[4]
196        sldSolv = w[5]
197        background = w[6]
198       
199        contrast = sld - sldSolv
200       
201        //Volume fraction calculated from lattice symmetry and sphere radius
202        latticeScale = 2*(4/3)*pi*(Rad^3)/((Dnn/((3/4)^0.5))^3)
203
204        NVAR/Z nord=root:gNordBCC
205        if(NVAR_Exists(nord)!=1)
206                nord=20
207        endif
208       
209
210        integral = IntegrateFn_N(Integrand_BCC_Outer,loLim,upLim,w,x,nord)
211       
212        integral *= SphereForm_BCC(Rad,contrast,x)*scale*latticeScale
213//      integral *= scale                       //testing, returns only Z(q)
214
215        integral += background 
216       
217        Return (integral)
218       
219End
220
221// the outer integral is also an integral
222Function Integrand_BCC_Outer(w,x,dum)
223        Wave w
224        Variable x,dum
225               
226        NVAR yy = root:gDumY           
227        yy = dum                                        // save the current dummy yy for use in the inner loop
228        Variable retVal,loLim,upLim
229        //
230        loLim = 0
231        upLim = Pi
232
233        NVAR/Z nord=root:gNordBCC
234        if(NVAR_Exists(nord)!=1)
235                nord=20
236        endif
237
238        retVal = IntegrateFn_N(Integrand_BCC_Inner,loLim,upLim,w,x,nord)
239       
240        return(retVal)
241End
242
243//returns the value of the integrand of the inner integral
244Function Integrand_BCC_Inner(w,qq,dum)
245        Wave w
246        Variable qq,dum
247       
248        NVAR yy = root:gDumY            //use the yy value from the outer loop
249        Variable xx,retVal
250        xx = dum
251
252        retVal = BCC_Integrand(w,qq,xx,yy)
253       
254        return(retVal)
255End
256
257Function BCC_Integrand(w,qq,xx,yy)
258        Wave w
259        Variable qq,xx,yy
260       
261        Variable retVal,temp1,temp2,temp3,temp4,temp5,aa,Da,Dnn,gg
262        Dnn = w[1] //Nearest neighbor distance A
263        gg = w[2] //Paracrystal distortion factor
264//      aa = Dnn*((4/3)^0.5)            //Danilo's version (the paper states |bi| = a)
265        aa = Dnn
266        Da = gg*aa
267       
268        temp1 = qq*qq*Da*Da
269        temp3 = qq*aa   
270       
271        retVal = BCCeval(xx,yy,temp1,temp3)
272        retVal /=4*Pi
273       
274        return(retVal)
275end
276
277Function BCCeval(Theta,Phi,temp1,temp3)
278        Variable Theta,Phi,temp1,temp3
279        Variable temp6,temp7,temp8,temp9,temp10
280        Variable result
281       
282        temp6 = sin(Theta)
283        temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta)
284        temp8 = -1*sin(Theta)*cos(Phi)-sin(Theta)*sin(Phi)+cos(Theta)
285        temp9 = -1*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)-cos(Theta)
286        temp10 = exp((-1/8)*temp1*((temp7^2)+(temp8^2)+(temp9^2)))
287        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)))
288       
289        return (result)
290end
291
292
293Function SphereForm_BCC(radius,delrho,x)                                       
294        Variable radius,delrho,x
295       
296        // variables are:                                                       
297        //[2] radius ()
298        //[3] delrho (-2)
299        //[4] background (cm-1)
300       
301        // calculates scale * f^2/Vol where f=Vol*3*delrho*(sin(qr)-qrcos(qr))/qr^3
302        // and is rescaled to give [=] cm^-1
303       
304        Variable bes,f,vol,f2
305        //
306        //handle q==0 separately
307        If(x==0)
308                f = 4/3*pi*radius^3*delrho*delrho*1e8
309                return(f)
310        Endif
311       
312        bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
313        vol = 4*pi/3*radius^3
314        f = vol*bes*delrho              // [=] 
315        // normalize to single particle volume, convert to 1/cm
316        f2 = f * f / vol * 1.0e8                // [=] 1/cm
317       
318        return (f2)     
319       
320End
321
322
323
324
325///////////////////////////////////////////////////////////////
326// smeared model calculation
327//
328Function SmearedBCC_ParaCrystal(s) : FitFunc
329        Struct ResSmearAAOStruct &s
330
331//      the name of your unsmeared model (AAO) is the first argument
332        Smear_Model_76(BCC_ParaCrystal,s.coefW,s.xW,s.yW,s.resW)
333
334        return(0)
335End
336
337
338/////////////////////////////////////////////////////////////////
339//wrapper to calculate the smeared model as an AAO-Struct
340// fills the struct and calls the ususal function with the STRUCT parameter
341//
342// used only for the dependency, not for fitting
343//
344Function fSmearedBCC_ParaCrystal(coefW,yW,xW)
345        Wave coefW,yW,xW
346       
347        String str = getWavesDataFolder(yW,0)
348        String DF="root:"+str+":"
349       
350        WAVE resW = $(DF+str+"_res")
351       
352        STRUCT ResSmearAAOStruct fs
353        WAVE fs.coefW = coefW   
354        WAVE fs.yW = yW
355        WAVE fs.xW = xW
356        WAVE fs.resW = resW
357       
358        Variable err
359        err = SmearedBCC_ParaCrystal(fs)
360       
361        return (0)
362End
Note: See TracBrowser for help on using the repository browser.