source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/BCC_ParaCrystal_v40.ipf @ 519

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

threading of model

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        return(0)
143End
144
145
146
147// nothing to change here
148//
149//AAO version, uses XOP if available
150// simply calls the original single point calculation with
151// a wave assignment (this will behave nicely if given point ranges)
152//
153// Threaded Version
154ThreadSafe Function BCC_ParaCrystal_T(cw,yw,xw,p1,p2) : FitFunc
155        Wave cw,yw,xw
156        Variable p1,p2
157
158//      Variable t1=StopMSTimer(-2)
159
160#if exists("BCC_ParaCrystalX")
161        yw[p1,p2] = BCC_ParaCrystalX(cw,xw)
162#else
163        yw[p1,p2] = fBCC_ParaCrystal(cw,xw)
164#endif
165
166//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
167
168        return(0)
169End
170
171
172//
173// unsmeared model calculation
174//
175Function fBCC_ParaCrystal(w,x) : FitFunc
176        Wave w
177        Variable x
178       
179//       Input (fitting) variables are not used
180//      you would give them nice names
181        Variable integral,loLim,upLim
182        loLim = 0
183        upLim = 2*Pi
184       
185        Variable/G root:gDumY=0         //root:gDumX=0
186       
187
188        Variable scale,Dnn,gg,Rad,contrast,background,yy,latticeScale,sld,sldSolv
189        scale = w[0]
190        Dnn = w[1] //Nearest neighbor distance A
191        gg = w[2] //Paracrystal distortion factor
192        Rad = w[3] //Sphere radius
193        sld = w[4]
194        sldSolv = w[5]
195        background = w[6]
196       
197        contrast = sld - sldSolv
198       
199        //Volume fraction calculated from lattice symmetry and sphere radius
200        latticeScale = 2*(4/3)*pi*(Rad^3)/((Dnn/((3/4)^0.5))^3)
201
202        NVAR/Z nord=root:gNordBCC
203        if(NVAR_Exists(nord)!=1)
204                nord=20
205        endif
206       
207
208        integral = IntegrateFn_N(Integrand_BCC_Outer,loLim,upLim,w,x,nord)
209       
210        integral *= SphereForm_BCC(Rad,contrast,x)*scale*latticeScale
211//      integral *= scale                       //testing, returns only Z(q)
212
213        integral += background 
214       
215        Return (integral)
216       
217End
218
219// the outer integral is also an integral
220Function Integrand_BCC_Outer(w,x,dum)
221        Wave w
222        Variable x,dum
223               
224        NVAR yy = root:gDumY           
225        yy = dum                                        // save the current dummy yy for use in the inner loop
226        Variable retVal,loLim,upLim
227        //
228        loLim = 0
229        upLim = Pi
230
231        NVAR/Z nord=root:gNordBCC
232        if(NVAR_Exists(nord)!=1)
233                nord=20
234        endif
235
236        retVal = IntegrateFn_N(Integrand_BCC_Inner,loLim,upLim,w,x,nord)
237       
238        return(retVal)
239End
240
241//returns the value of the integrand of the inner integral
242Function Integrand_BCC_Inner(w,qq,dum)
243        Wave w
244        Variable qq,dum
245       
246        NVAR yy = root:gDumY            //use the yy value from the outer loop
247        Variable xx,retVal
248        xx = dum
249
250        retVal = BCC_Integrand(w,qq,xx,yy)
251       
252        return(retVal)
253End
254
255Function BCC_Integrand(w,qq,xx,yy)
256        Wave w
257        Variable qq,xx,yy
258       
259        Variable retVal,temp1,temp2,temp3,temp4,temp5,aa,Da,Dnn,gg
260        Dnn = w[1] //Nearest neighbor distance A
261        gg = w[2] //Paracrystal distortion factor
262//      aa = Dnn*((4/3)^0.5)            //Danilo's version (the paper states |bi| = a)
263        aa = Dnn
264        Da = gg*aa
265       
266        temp1 = qq*qq*Da*Da
267        temp3 = qq*aa   
268       
269        retVal = BCCeval(xx,yy,temp1,temp3)
270        retVal /=4*Pi
271       
272        return(retVal)
273end
274
275Function BCCeval(Theta,Phi,temp1,temp3)
276        Variable Theta,Phi,temp1,temp3
277        Variable temp6,temp7,temp8,temp9,temp10
278        Variable result
279       
280        temp6 = sin(Theta)
281        temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta)
282        temp8 = -1*sin(Theta)*cos(Phi)-sin(Theta)*sin(Phi)+cos(Theta)
283        temp9 = -1*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)-cos(Theta)
284        temp10 = exp((-1/8)*temp1*((temp7^2)+(temp8^2)+(temp9^2)))
285        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)))
286       
287        return (result)
288end
289
290
291Function SphereForm_BCC(radius,delrho,x)                                       
292        Variable radius,delrho,x
293       
294        // variables are:                                                       
295        //[2] radius ()
296        //[3] delrho (-2)
297        //[4] background (cm-1)
298       
299        // calculates scale * f^2/Vol where f=Vol*3*delrho*(sin(qr)-qrcos(qr))/qr^3
300        // and is rescaled to give [=] cm^-1
301       
302        Variable bes,f,vol,f2
303        //
304        //handle q==0 separately
305        If(x==0)
306                f = 4/3*pi*radius^3*delrho*delrho*1e8
307                return(f)
308        Endif
309       
310        bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
311        vol = 4*pi/3*radius^3
312        f = vol*bes*delrho              // [=] 
313        // normalize to single particle volume, convert to 1/cm
314        f2 = f * f / vol * 1.0e8                // [=] 1/cm
315       
316        return (f2)     
317       
318End
319
320
321
322
323///////////////////////////////////////////////////////////////
324// smeared model calculation
325//
326Function SmearedBCC_ParaCrystal(s) : FitFunc
327        Struct ResSmearAAOStruct &s
328
329//      the name of your unsmeared model (AAO) is the first argument
330        Smear_Model_76(BCC_ParaCrystal,s.coefW,s.xW,s.yW,s.resW)
331
332        return(0)
333End
334
335
336/////////////////////////////////////////////////////////////////
337//wrapper to calculate the smeared model as an AAO-Struct
338// fills the struct and calls the ususal function with the STRUCT parameter
339//
340// used only for the dependency, not for fitting
341//
342Function fSmearedBCC_ParaCrystal(coefW,yW,xW)
343        Wave coefW,yW,xW
344       
345        String str = getWavesDataFolder(yW,0)
346        String DF="root:"+str+":"
347       
348        WAVE resW = $(DF+str+"_res")
349       
350        STRUCT ResSmearAAOStruct fs
351        WAVE fs.coefW = coefW   
352        WAVE fs.yW = yW
353        WAVE fs.xW = xW
354        WAVE fs.resW = resW
355       
356        Variable err
357        err = SmearedBCC_ParaCrystal(fs)
358       
359        return (0)
360End
Note: See TracBrowser for help on using the repository browser.