source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Sphere_2D_v40.ipf

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

in 2D models, adjusted the point ranges to be integers. non-integer values occasionally causes issues.

changed names of polarization procedures to start with "Pol_" so that they are easier to locate.

added MixedDumbbell_FFT model as another example of a FFT fitting function.
File size: 11.9 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 PlotSphere2D(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_sf2D = {1.,60,1e-6,6.3e-6,0.01}                                           
27        make/o/t parameters_sf2D = {"scale","Radius (A)","SLD sphere (A-2)","SLD solvent","bkgd (cm-1)"}               
28        Edit parameters_sf2D,coef_sf2D                         
29       
30        // generate the triplet representation
31        Duplicate/O $(str+"_qx") xwave_sf2D
32        Duplicate/O $(str+"_qy") ywave_sf2D,zwave_sf2D                 
33               
34        Variable/G g_sf2D=0
35        g_sf2D := Sphere2D(coef_sf2D,zwave_sf2D,xwave_sf2D,ywave_sf2D)  //AAO 2D calculation
36       
37        Display ywave_sf2D vs xwave_sf2D
38        modifygraph log=0
39        ModifyGraph mode=3,marker=16,zColor(ywave_sf2D)={zwave_sf2D,*,*,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_sf2D,ywave_sf2D,zwave_sf2D,"sf2D_mat")
49        Duplicate/O $"sf2D_mat",$"sf2D_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_sf2Dmat=0
54        g_sf2Dmat := UpdateQxQy2Mat(xwave_sf2D,ywave_sf2D,zwave_sf2D,sf2D_lin,sf2D_mat)
55       
56       
57        SetDataFolder root:
58        AddModelToStrings("Sphere2D","coef_sf2D","parameters_sf2D","sf2D")
59End
60
61// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
62Proc PlotSmearedSphere2D(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_sf2D = {1.,60,1e-6,6.3e-6,0.01}                                     
75        make/o/t smear_parameters_sf2D = {"scale","Radius (A)","SLD sphere (A-2)","SLD solvent (A-2)","bkgd (cm-1)"}
76        Edit smear_parameters_sf2D,smear_coef_sf2D                                     
77       
78        Duplicate/O $(str+"_qx") smeared_sf2D   //1d place for the smeared model
79        SetScale d,0,0,"1/cm",smeared_sf2D                                     
80               
81        Variable/G gs_sf2D=0
82        gs_sf2D := fSmearedSphere2D(smear_coef_sf2D,smeared_sf2D)       //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_sf2D,*,*,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_sf2D,"sm_sf2D_mat")
99        Duplicate/O $"sm_sf2D_mat",$"sm_sf2D_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_sf2Dmat=0
104        gs_sf2Dmat := UpdateQxQy2Mat(sm_qx,sm_qy,smeared_sf2D,sm_sf2D_lin,sm_sf2D_mat)
105       
106        SetDataFolder root:
107        AddModelToStrings("SmearedSphere2D","smear_coef_sf2D","smear_parameters_sf2D","sf2D")
108End
109
110//
111//  Fit function that is actually a wrapper to dispatch the calculation to N threads
112//
113// nthreads is 1 or an even number, typically 2
114// it doesn't matter if npt is odd. In this case, fractional point numbers are passed
115// and the wave indexing works just fine - I tested this with test waves of 7 and 8 points
116// and the points "2.5" and "3.5" evaluate correctly as 2 and 3
117//
118Function Sphere2D(cw,zw,xw,yw) : FitFunc
119        Wave cw,zw,xw,yw
120
121#if exists("Sphere_2DX")                        //to hide the function if XOP not installed
122        MultiThread     zw= Sphere_2DX(cw,xw,yw)
123#endif
124
125        return(0)
126End
127
128
129///
130//// keep this section as an example
131//
132
133//      Variable npt=numpnts(yw)
134//      Variable i,nthreads= ThreadProcessorCount
135//      variable mt= ThreadGroupCreate(nthreads)
136//
137////    Variable t1=StopMSTimer(-2)
138//     
139//      for(i=0;i<nthreads;i+=1)
140//      //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
141//              ThreadStart mt,i,Sphere2D_T(cw,zw,xw,yw,(i*npt/nthreads),((i+1)*npt/nthreads-1))
142//      endfor
143//
144//      do
145//              variable tgs= ThreadGroupWait(mt,100)
146//      while( tgs != 0 )
147//
148//      variable dummy= ThreadGroupRelease(mt)
149//     
150////    Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
151//     
152//
153////// end example of threading
154
155
156
157//threaded version of the function
158//ThreadSafe Function Sphere2D_T(cw,zw,xw,yw,p1,p2)
159//      WAVE cw,zw,xw,yw
160//      Variable p1,p2
161//     
162//#if exists("Sphere_2DX")                      //to hide the function if XOP not installed
163//     
164//      zw[p1,p2]= Sphere_2DX(cw,xw,yw)
165//
166//#endif
167//
168//      return 0
169//End
170
171
172//non-threaded version of the function, necessary for the smearing calculation
173// -- the smearing calculation can only calculate (nord) points at a time.
174//
175ThreadSafe Function Sphere2D_noThread(cw,zw,xw,yw)
176        WAVE cw,zw, xw,yw
177       
178#if exists("Sphere_2DX")                        //to hide the function if XOP not installed
179        zw= Sphere_2DX(cw,xw,yw)
180#endif
181
182        return 0
183End
184
185
186
187//// the threaded version must be specifically written, since
188//// FUNCREF can't be passed into  a threaded calc (structures can't be passed either)
189// so in this implementation, the smearing is dispatched as threads to a function that
190// can calculate the function for a range of points in the input qxqyqz. It is important
191// that the worker calls the un-threaded model function (so write one) and that in the (nord x nord)
192// loop, vectors of length (nord) are calculated rather than pointwise, since the model
193// function is AAO.
194// -- makes things rather messy to code individual functions, but I really see no other way
195// given the restrictions of what can be passed to threaded functions.
196//
197//
198// The smearing is handled this way since 1D smearing is 20 x 200 pts = 4000 evaluations
199// and the 2D is (10 x 10) x 16000 pts = 1,600,000 evaluations (if it's done like the 1D, it's 4000x slower)
200//
201//
202// - the threading gives a clean speedup of 2 for N=2, even for this simple calculation
203//  -- 4.8X speedup for N=8 (4 real cores + 4 virtual cores)
204//
205// nord = 5,10,20 allowed
206//
207Function SmearedSphere2D(s)
208        Struct ResSmear_2D_AAOStruct &s
209       
210//// non-threaded, but generic calculation
211//// the last param is nord     
212//      Smear_2DModel_PP(Sphere2D_noThread,s,10)
213
214
215//// the last param is nord
216        SmearedSphere2D_THR(s,10)               
217
218        return(0)
219end
220
221
222Function fSmearedSphere2D(coefW,resultW)
223        Wave coefW,resultW
224       
225        String str = getWavesDataFolder(resultW,0)
226        String DF="root:"+str+":"
227       
228        WAVE qx = $(DF+str+"_qx")
229        WAVE qy = $(DF+str+"_qy")
230        WAVE qz = $(DF+str+"_qz")
231        WAVE sQpl = $(DF+str+"_sQpl")
232        WAVE sQpp = $(DF+str+"_sQpp")
233        WAVE shad = $(DF+str+"_fs")
234       
235        STRUCT ResSmear_2D_AAOStruct s
236        WAVE s.coefW = coefW   
237        WAVE s.zw = resultW     
238        WAVE s.xw[0] = qx
239        WAVE s.xw[1] = qy
240        WAVE s.qz = qz
241        WAVE s.sQpl = sQpl
242        WAVE s.sQpp = sQpp
243        WAVE s.fs = shad
244       
245        Variable err
246        err = SmearedSphere2D(s)
247       
248        return (0)
249End
250
251
252
253
254//
255// this is the threaded version, that dispatches the calculation out to threads
256//
257// must be written specific to each 2D function
258//
259Function SmearedSphere2D_THR(s,nord)
260        Struct ResSmear_2D_AAOStruct &s
261        Variable nord
262       
263        String weightStr,zStr
264       
265// create all of the necessary quadrature waves here - rather than inside a threadsafe function
266        switch(nord)   
267                case 5:         
268                        weightStr="gauss5wt"
269                        zStr="gauss5z"
270                        if (WaveExists($weightStr) == 0)
271                                Make/O/D/N=(nord) $weightStr,$zStr
272                                Make5GaussPoints($weightStr,$zStr)     
273                        endif
274                        break                           
275                case 10:               
276                        weightStr="gauss10wt"
277                        zStr="gauss10z"
278                        if (WaveExists($weightStr) == 0)
279                                Make/O/D/N=(nord) $weightStr,$zStr
280                                Make10GaussPoints($weightStr,$zStr)     
281                        endif
282                        break                           
283                case 20:               
284                        weightStr="gauss20wt"
285                        zStr="gauss20z"
286                        if (WaveExists($weightStr) == 0)
287                                Make/O/D/N=(nord) $weightStr,$zStr
288                                Make20GaussPoints($weightStr,$zStr)     
289                        endif
290                        break
291                default:                                                       
292                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                       
293        endswitch
294       
295        Wave/Z wt = $weightStr
296        Wave/Z xi = $zStr               // wave references to pass
297
298        Variable npt=numpnts(s.xw[0])
299        Variable i,nthreads= ThreadProcessorCount
300        variable mt= ThreadGroupCreate(nthreads)
301
302        Variable t1=StopMSTimer(-2)
303       
304        for(i=0;i<nthreads;i+=1)
305//              Print trunc(i*npt/nthreads),trunc((i+1)*npt/nthreads-1)
306                ThreadStart mt,i,SmearedSphere2D_T(s.coefW,s.xw[0],s.xw[1],s.qz,s.sQpl,s.sQpp,s.fs,s.zw,wt,xi,trunc(i*npt/nthreads),trunc((i+1)*npt/nthreads-1),nord)
307        endfor
308
309        do
310                variable tgs= ThreadGroupWait(mt,100)
311        while( tgs != 0 )
312
313        variable dummy= ThreadGroupRelease(mt)
314       
315// comment out the threading + uncomment this for testing to make sure that the single thread works
316//      nThreads=1
317//      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)
318
319        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
320       
321        return(0)
322end
323
324//
325// - worker function for threads of Sphere2D
326//
327ThreadSafe Function SmearedSphere2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord)
328        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi
329        Variable pt1,pt2,nord
330       
331// now passed in....
332//      Wave wt = $weightStr
333//      Wave xi = $zStr         
334
335        Variable ii,jj,kk,num
336        Variable qx,qy,qz,qval,sx,sy,fs
337        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
338       
339        Variable normFactor,phi,theta,maxSig,numStdDev=3
340       
341/// keep these waves local
342        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW
343       
344        // now just loop over the points as specified
345       
346        answer=0
347       
348        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt
349        Variable qperp_pt,phi_prime,q_prime
350
351        //loop over q-values
352        for(ii=pt1;ii<(pt2+1);ii+=1)
353               
354                qx = qxw[ii]
355                qy = qyw[ii]
356                qz = qzw[ii]
357                qval = sqrt(qx^2+qy^2+qz^2)
358                spl = sxw[ii]
359                spp = syw[ii]
360                fs = fsw[ii]
361
362               
363                normFactor = 2*pi*spl*spp
364               
365                phi = FindPhi(qx,qy)
366               
367                apl = -numStdDev*spl + qval             //parallel = q integration limits
368                bpl = numStdDev*spl + qval
369///             app = -numStdDev*spp + phi              //perpendicular = phi integration limits (WRONG)
370///             bpp = numStdDev*spp + phi
371                app = -numStdDev*spp + 0                //q_perp = 0
372                bpp = numStdDev*spp + 0
373               
374                //make sure the limits are reasonable.
375                if(apl < 0)
376                        apl = 0
377                endif
378                // do I need to specially handle limits when phi ~ 0?
379       
380               
381                sumOut = 0
382                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
383///                     phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2
384                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp
385                       
386                        sumIn=0
387                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
388
389                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2
390                               
391///                             FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi
392
393                                // find QxQy given Qpl and Qperp on the grid
394                                //
395                                q_prime = sqrt(qpl_pt^2+qperp_pt^2)
396                                phi_prime = phi + qperp_pt/q_prime
397                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt)
398                               
399                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not
400                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
401
402                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) )
403///                             res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) )
404                                res_tot[kk] /= normFactor
405//                              res_tot[kk] *= fs
406
407                        endfor
408                       
409                        Sphere2D_noThread(coef,fcnRet,xptw,yptw)                        //fcn passed in is an AAO
410                       
411                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0]
412                        fcnRet *= wt[jj]*wt*res_tot
413                        //
414                        answer += (bpl-apl)/2.0*sum(fcnRet)             //
415                endfor
416
417                answer *= (bpp-app)/2.0
418                zw[ii] = answer
419        endfor
420       
421        return(0)
422end
423
424
Note: See TracBrowser for help on using the repository browser.