source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/PeakGauss2D_v40.ipf @ 775

Last change on this file since 775 was 775, checked in by srkline, 12 years ago

Added 2D versions of sphere and GaussPeak? to the package. Probably never need them, but good for testing.

File size: 10.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4//
5// The plotting macro sets up TWO dependencies
6// - one for the triplet calculation
7// - one for a matrix to display, a copy of the triplet
8//
9// For display, there are two copies of the matrix. One matrix is linear, and is a copy of the
10// triplet (which is ALWAYS linear). The other matrix is toggled log/lin for display
11// in the same way the 2D SANS data matrix is handled.
12//
13
14///  REQUIRES XOP for 2D FUNCTIONS
15
16//
17// the calculation is done as for the QxQy data set:
18// three waves XYZ, then converted to a matrix
19//
20Proc PlotPeakGauss2D(str)                                               
21        String str
22        Prompt str,"Pick the data folder containing the 2D data",popup,getAList(4)
23       
24        SetDataFolder $("root:"+str)
25       
26        Make/O/D coef_pkGau2D = {100.0, 0.25,0.005, 1.0}                                               
27        make/o/t parameters_pkGau2D = {"Scale Factor, I0 ", "Peak position (A^-1)", "Std Dev (A^-1)","Incoherent Bgd (cm-1)"}           
28        Edit parameters_pkGau2D,coef_pkGau2D                           
29       
30        // generate the triplet representation
31        Duplicate/O $(str+"_qx") xwave_pkGau2D
32        Duplicate/O $(str+"_qy") ywave_pkGau2D,zwave_pkGau2D                   
33               
34        Variable/G g_pkGau2D=0
35        g_pkGau2D := PeakGauss2D(coef_pkGau2D,zwave_pkGau2D,xwave_pkGau2D,ywave_pkGau2D)        //AAO 2D calculation
36       
37        Display ywave_pkGau2D vs xwave_pkGau2D
38        modifygraph log=0
39        ModifyGraph mode=3,marker=16,zColor(ywave_pkGau2D)={zwave_pkGau2D,*,*,YellowHot,0}
40        ModifyGraph standoff=0
41        ModifyGraph width={Aspect,1}
42        ModifyGraph lowTrip=0.001
43        Label bottom "qx (A\\S-1\\M)"
44        Label left "qy (A\\S-1\\M)"
45        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
46       
47        // generate the matrix representation
48        ConvertQxQy2Mat(xwave_pkGau2D,ywave_pkGau2D,zwave_pkGau2D,"pkGau2D_mat")
49        Duplicate/O $"pkGau2D_mat",$"pkGau2D_lin"               //keep a linear-scaled version of the data
50        // _mat is for display, _lin is the real calculation
51
52        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation
53        Variable/G g_pkGau2Dmat=0
54        g_pkGau2Dmat := UpdateQxQy2Mat(xwave_pkGau2D,ywave_pkGau2D,zwave_pkGau2D,pkGau2D_lin,pkGau2D_mat)
55       
56       
57        SetDataFolder root:
58        AddModelToStrings("PeakGauss2D","coef_pkGau2D","parameters_pkGau2D","pkGau2D")
59End
60
61// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
62Proc PlotSmearedPeakGauss2D(str)                                                               
63        String str
64        Prompt str,"Pick the data folder containing the 2D data",popup,getAList(4)
65       
66        // if any of the resolution waves are missing => abort
67//      if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
68//              Abort
69//      endif
70       
71        SetDataFolder $("root:"+str)
72       
73        // Setup parameter table for model function
74        Make/O/D smear_coef_pkGau2D = {100.0, 0.25,0.005, 1.0}                                 
75        make/o/t smear_parameters_pkGau2D = {"Scale Factor, I0 ", "Peak position (A^-1)", "Std Dev (A^-1)","Incoherent Bgd (cm-1)"}
76        Edit smear_parameters_pkGau2D,smear_coef_pkGau2D                                       
77       
78        Duplicate/O $(str+"_qx") smeared_pkGau2D        //1d place for the smeared model
79        SetScale d,0,0,"1/cm",smeared_pkGau2D                                   
80               
81        Variable/G gs_pkGau2D=0
82        gs_pkGau2D := fSmearedPeakGauss2D(smear_coef_pkGau2D,smeared_pkGau2D)   //this wrapper fills the STRUCT
83
84        Display $(str+"_qy") vs $(str+"_qx")
85        modifygraph log=0
86        ModifyGraph mode=3,marker=16,zColor($(str+"_qy"))={smeared_pkGau2D,*,*,YellowHot,0}
87        ModifyGraph standoff=0
88        ModifyGraph width={Aspect,1}
89        ModifyGraph lowTrip=0.001
90        Label bottom "qx (A\\S-1\\M)"
91        Label left "qy (A\\S-1\\M)"
92        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
93       
94        // generate the matrix representation
95        Duplicate/O $(str+"_qx"), sm_qx
96        Duplicate/O $(str+"_qy"), sm_qy         // I can't use local variables in dependencies, so I need the name (that I can't get)
97       
98        ConvertQxQy2Mat(sm_qx,sm_qy,smeared_pkGau2D,"sm_pkGau2D_mat")
99        Duplicate/O $"sm_pkGau2D_mat",$"sm_pkGau2D_lin"                 //keep a linear-scaled version of the data
100        // _mat is for display, _lin is the real calculation
101
102        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation
103        Variable/G gs_pkGau2Dmat=0
104        gs_pkGau2Dmat := UpdateQxQy2Mat(sm_qx,sm_qy,smeared_pkGau2D,sm_pkGau2D_lin,sm_pkGau2D_mat)
105       
106        SetDataFolder root:
107        AddModelToStrings("SmearedPeakGauss2D","smear_coef_pkGau2D","smear_parameters_pkGau2D","pkGau2D")
108End
109
110
111//threaded version of the function
112ThreadSafe Function PeakGauss2D_T(cw,zw,xw,yw,p1,p2)
113        WAVE cw,zw,xw,yw
114        Variable p1,p2
115       
116#if exists("PeakGauss2DX")                      //to hide the function if XOP not installed
117       
118        zw[p1,p2]= PeakGauss2DX(cw,xw,yw)
119
120#else
121        zw[p1,p2] =  I_PeakGauss2D(cw,xw,yw)
122#endif
123
124        return 0
125End
126
127//
128//  Fit function that is actually a wrapper to dispatch the calculation to N threads
129//
130// not used
131//
132Function PeakGauss2D(cw,zw,xw,yw) : FitFunc
133        Wave cw,zw,xw,yw
134       
135#if exists("PeakGauss2DX")                      //to hide the function if XOP not installed
136        MultiThread zw= PeakGauss2DX(cw,xw,yw)
137#else
138        MultiThread zw = I_PeakGauss2D(cw,xw,yw)
139#endif 
140       
141//      Variable npt=numpnts(yw)
142//      Variable i,nthreads= ThreadProcessorCount
143//      variable mt= ThreadGroupCreate(nthreads)
144//
145////    Variable t1=StopMSTimer(-2)
146//     
147//      for(i=0;i<nthreads;i+=1)
148//      //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
149//              ThreadStart mt,i,PeakGauss2D_T(cw,zw,xw,yw,(i*npt/nthreads),((i+1)*npt/nthreads-1))
150//      endfor
151//
152//      do
153//              variable tgs= ThreadGroupWait(mt,100)
154//      while( tgs != 0 )
155//
156//      variable dummy= ThreadGroupRelease(mt)
157//     
158////    Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
159       
160        return(0)
161End
162
163
164//non-threaded version of the function
165ThreadSafe Function PeakGauss2D_noThread(cw,zw,xw,yw)
166        WAVE cw,zw, xw,yw
167       
168#if exists("PeakGauss2DX")                      //to hide the function if XOP not installed
169        zw= PeakGauss2DX(cw,xw,yw)
170#else
171        zw = I_PeakGauss2D(cw,xw,yw)
172#endif
173
174        return 0
175End
176
177// I think I did this because when I do the quadrature loops I'm calling the AAO with 1-pt waves, so threading
178// would just be a slowdown
179Function SmearedPeakGauss2D(s)
180        Struct ResSmear_2D_AAOStruct &s
181       
182        //no threads
183//      Smear_2DModel_PP(PeakGauss2D_noThread,s,10)
184
185        //or threaded...
186        SmearPeakGauss2D_THR(s,10)
187       
188        return(0)
189end
190
191
192Function fSmearedPeakGauss2D(coefW,resultW)
193        Wave coefW,resultW
194       
195        String str = getWavesDataFolder(resultW,0)
196        String DF="root:"+str+":"
197       
198        WAVE qx = $(DF+str+"_qx")
199        WAVE qy = $(DF+str+"_qy")
200        WAVE qz = $(DF+str+"_qz")
201        WAVE sQpl = $(DF+str+"_sQpl")
202        WAVE sQpp = $(DF+str+"_sQpp")
203        WAVE shad = $(DF+str+"_fs")
204       
205        STRUCT ResSmear_2D_AAOStruct s
206        WAVE s.coefW = coefW   
207        WAVE s.zw = resultW     
208        WAVE s.xw[0] = qx
209        WAVE s.xw[1] = qy
210        WAVE s.qz = qz
211        WAVE s.sQpl = sQpl
212        WAVE s.sQpp = sQpp
213        WAVE s.fs = shad
214       
215        Variable err
216        err = SmearedPeakGauss2D(s)
217       
218        return (0)
219End
220
221ThreadSafe Function I_PeakGauss2D(w,x,y)
222        Wave w
223        Variable x,y
224       
225        Variable retVal,qval
226        qval = sqrt(x^2+y^2)
227               
228        retval = Peak_Gauss_modelX(w,qval)
229//      retval = fPeak_Gauss_model(w,qval)
230       
231        return(retVal)
232End
233
234///////////////
235//
236// this is the threaded version, that dispatches the calculation out to threads
237//
238// must be written specific to each 2D function
239//
240Function SmearPeakGauss2D_THR(s,nord)
241        Struct ResSmear_2D_AAOStruct &s
242        Variable nord
243       
244        String weightStr,zStr
245       
246// create all of the necessary quadrature waves here - rather than inside a threadsafe function
247        switch(nord)   
248                case 5:         
249                        weightStr="gauss5wt"
250                        zStr="gauss5z"
251                        if (WaveExists($weightStr) == 0)
252                                Make/O/D/N=(nord) $weightStr,$zStr
253                                Make5GaussPoints($weightStr,$zStr)     
254                        endif
255                        break                           
256                case 10:               
257                        weightStr="gauss10wt"
258                        zStr="gauss10z"
259                        if (WaveExists($weightStr) == 0)
260                                Make/O/D/N=(nord) $weightStr,$zStr
261                                Make10GaussPoints($weightStr,$zStr)     
262                        endif
263                        break                           
264                case 20:               
265                        weightStr="gauss20wt"
266                        zStr="gauss20z"
267                        if (WaveExists($weightStr) == 0)
268                                Make/O/D/N=(nord) $weightStr,$zStr
269                                Make20GaussPoints($weightStr,$zStr)     
270                        endif
271                        break
272                default:                                                       
273                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                       
274        endswitch
275       
276        Wave/Z wt = $weightStr
277        Wave/Z xi = $zStr               // wave references to pass
278
279        Variable npt=numpnts(s.xw[0])
280        Variable i,nthreads= ThreadProcessorCount
281        variable mt= ThreadGroupCreate(nthreads)
282
283        Variable t1=StopMSTimer(-2)
284       
285        for(i=0;i<nthreads;i+=1)
286        //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
287                ThreadStart mt,i,SmearPeakGauss2D_T(s.coefW,s.xw[0],s.xw[1],s.qz,s.sQpl,s.sQpp,s.fs,s.zw,wt,xi,(i*npt/nthreads),((i+1)*npt/nthreads-1),nord)
288        endfor
289
290        do
291                variable tgs= ThreadGroupWait(mt,100)
292        while( tgs != 0 )
293
294        variable dummy= ThreadGroupRelease(mt)
295       
296// comment out the threading + uncomment this for testing to make sure that the single thread works
297//      nThreads=1
298//      SmearSphere2D_T(s.coefW,s.xw[0],s.xw[1],s.qz,s.sQpl,s.sQpp,s.fs,s.zw,wt,xi,(i*npt/nthreads),((i+1)*npt/nthreads-1),nord)
299
300        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
301       
302        return(0)
303end
304
305//
306// - worker function for threads of PeakGauss2D
307//
308ThreadSafe Function SmearPeakGauss2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord)
309        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi
310        Variable pt1,pt2,nord
311       
312// now passed in....
313//      Wave wt = $weightStr
314//      Wave xi = $zStr         
315
316        Variable ii,jj,kk,num
317        Variable qx,qy,qz,qval,sx,sy,fs
318        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
319       
320        Variable normFactor,phi,theta,maxSig,numStdDev=3
321       
322/// keep these waves local
323        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW
324       
325        // now just loop over the points as specified
326       
327        answer=0
328       
329        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt
330
331        //loop over q-values
332        for(ii=pt1;ii<(pt2+1);ii+=1)
333               
334                qx = qxw[ii]
335                qy = qyw[ii]
336                qz = qzw[ii]
337                qval = sqrt(qx^2+qy^2+qz^2)
338                spl = sxw[ii]
339                spp = syw[ii]
340                fs = fsw[ii]
341               
342                normFactor = 2*pi*spl*spp
343               
344                phi = FindPhi(qx,qy)
345               
346                apl = -numStdDev*spl + qval             //parallel = q integration limits
347                bpl = numStdDev*spl + qval
348                app = -numStdDev*spp + phi              //perpendicular = phi integration limits
349                bpp = numStdDev*spp + phi
350               
351                //make sure the limits are reasonable.
352                if(apl < 0)
353                        apl = 0
354                endif
355                // do I need to specially handle limits when phi ~ 0?
356       
357               
358                sumOut = 0
359                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
360                        phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2
361
362                        sumIn=0
363                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
364
365                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2
366                               
367                                FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi
368                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not
369                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
370                               
371                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) )
372                                res_tot[kk] /= normFactor
373//                              res_tot[kk] *= fs
374
375                        endfor
376                       
377                        PeakGauss2D_noThread(coef,fcnRet,xptw,yptw)                     //fcn passed in is an AAO
378                       
379                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0]
380                        fcnRet *= wt[jj]*wt*res_tot
381                        //
382                        answer += (bpl-apl)/2.0*sum(fcnRet)             //
383                endfor
384
385                answer *= (bpp-app)/2.0
386                zw[ii] = answer
387        endfor
388       
389        return(0)
390end
391
Note: See TracBrowser for help on using the repository browser.