source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Ellipsoid2D_v40.ipf @ 771

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

Adding new model functions to the picker list

Adding resolution-smeared functions to the 2D analysis models

Constraints are now functional during a 2D fit

vesib.cor example data is now 6-column

File size: 12.4 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4//
5// !!! FOR THE ELLIPSOID, THE ANGLE THETA IS DEFINED FROM ????
6//
7// The plotting macro sets up TWO dependencies
8// - one for the triplet calculation
9// - one for a matrix to display, a copy of the triplet
10//
11// For display, there are two copies of the matrix. One matrix is linear, and is a copy of the
12// triplet (which is ALWAYS linear). The other matrix is toggled log/lin for display
13// in the same way the 2D SANS data matrix is handled.
14//
15
16///  REQUIRES DANSE XOP for 2D FUNCTIONS
17
18//
19// the calculation is done as for the QxQy data set:
20// three waves XYZ, then converted to a matrix
21//
22Proc PlotEllipsoid2D(str)                                               
23        String str
24        Prompt str,"Pick the data folder containing the 2D data",popup,getAList(4)
25
26        if (!exists("Ellipsoid_2DX"))
27                Abort "You must have the SANSAnalysis XOP installed to use 2D models"
28        endif
29       
30        SetDataFolder $("root:"+str)
31       
32        // Setup parameter table for model function
33//      make/O/T/N=12 parameters_Ellip2D
34//      Make/O/D/N=12 coef_Ellip2D
35        make/O/T/N=12 parameters_Ellip2D
36        Make/O/D/N=12 coef_Ellip2D
37       
38        coef_Ellip2D[0] = 1.0
39        coef_Ellip2D[1] = 20.0
40        coef_Ellip2D[2] = 60.0
41        coef_Ellip2D[3] = 1.0e-6
42        coef_Ellip2D[4] = 6.3e-6
43        coef_Ellip2D[5] = 0.0
44        coef_Ellip2D[6] = 1.57
45        coef_Ellip2D[7] = 0.0
46        coef_Ellip2D[8] = 0.0
47        coef_Ellip2D[9] = 0.0
48        coef_Ellip2D[10] = 0.0
49        coef_Ellip2D[11] = 0.0
50        // hard-wire the number of integration points
51//      coef_Ellip2D[12] = 10
52       
53        parameters_Ellip2D[0] = "Scale"
54        parameters_Ellip2D[1] = "Radius_a (rotation axis)"
55        parameters_Ellip2D[2] = "Radius_b"
56        parameters_Ellip2D[3] = "SLD cylinder (A^-2)"
57        parameters_Ellip2D[4] = "SLD solvent"
58        parameters_Ellip2D[5] = "Background"
59        parameters_Ellip2D[6] = "Axis Theta"
60        parameters_Ellip2D[7] = "Axis Phi"
61        parameters_Ellip2D[8] = "Sigma of polydisp in R_a [Angstrom]"
62        parameters_Ellip2D[9] = "Sigma of polydisp in R_b [Angstrom]"
63        parameters_Ellip2D[10] = "Sigma of polydisp in Theta [rad]"
64        parameters_Ellip2D[11] = "Sigma of polydisp in Phi [rad]"
65       
66//      parameters_Ellip2D[12] = "Num of polydisp points"
67
68       
69        Edit parameters_Ellip2D,coef_Ellip2D                                   
70       
71        // generate the triplet representation
72        Duplicate/O $(str+"_qx") xwave_Ellip2D
73        Duplicate/O $(str+"_qy") ywave_Ellip2D,zwave_Ellip2D                   
74               
75        Variable/G g_Ellip2D=0
76        g_Ellip2D := Ellipsoid2D(coef_Ellip2D,zwave_Ellip2D,xwave_Ellip2D,ywave_Ellip2D)        //AAO 2D calculation
77       
78        Display ywave_Ellip2D vs xwave_Ellip2D
79        modifygraph log=0
80        ModifyGraph mode=3,marker=16,zColor(ywave_Ellip2D)={zwave_Ellip2D,*,*,YellowHot,0}
81        ModifyGraph standoff=0
82        ModifyGraph width={Aspect,1}
83        ModifyGraph lowTrip=0.001
84        Label bottom "qx (A\\S-1\\M)"
85        Label left "qy (A\\S-1\\M)"
86        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
87       
88        // generate the matrix representation
89        ConvertQxQy2Mat(xwave_Ellip2D,ywave_Ellip2D,zwave_Ellip2D,"Ellip2D_mat")
90        Duplicate/O $"Ellip2D_mat",$"Ellip2D_lin"               //keep a linear-scaled version of the data
91        // _mat is for display, _lin is the real calculation
92
93        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation
94        Variable/G g_Ellip2Dmat=0
95        g_Ellip2Dmat := UpdateQxQy2Mat(xwave_Ellip2D,ywave_Ellip2D,zwave_Ellip2D,Ellip2D_lin,Ellip2D_mat)
96       
97       
98        SetDataFolder root:
99        AddModelToStrings("Ellipsoid2D","coef_Ellip2D","parameters_Ellip2D","Ellip2D")
100End
101
102
103// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
104Proc PlotSmearedEllipsoid2D(str)                                               
105        String str
106        Prompt str,"Pick the data folder containing the 2D data",popup,getAList(4)
107
108        if (!exists("Ellipsoid_2DX"))
109                Abort "You must have the SANSAnalysis XOP installed to use 2D models"
110        endif
111       
112        SetDataFolder $("root:"+str)
113       
114        // Setup parameter table for model function
115        make/O/T/N=12 smear_parameters_Ellip2D
116        Make/O/D/N=12 smear_coef_Ellip2D
117       
118        smear_coef_Ellip2D[0] = 1.0
119        smear_coef_Ellip2D[1] = 20.0
120        smear_coef_Ellip2D[2] = 60.0
121        smear_coef_Ellip2D[3] = 1.0e-6
122        smear_coef_Ellip2D[4] = 6.3e-6
123        smear_coef_Ellip2D[5] = 0.0
124        smear_coef_Ellip2D[6] = 1.57
125        smear_coef_Ellip2D[7] = 0.0
126        smear_coef_Ellip2D[8] = 0.0
127        smear_coef_Ellip2D[9] = 0.0
128        smear_coef_Ellip2D[10] = 0.0
129        smear_coef_Ellip2D[11] = 0.0
130        // hard-wire the number of integration points
131//      smear_coef_Ellip2D[12] = 10
132       
133        smear_parameters_Ellip2D[0] = "Scale"
134        smear_parameters_Ellip2D[1] = "Radius_a (rotation axis)"
135        smear_parameters_Ellip2D[2] = "Radius_b"
136        smear_parameters_Ellip2D[3] = "SLD cylinder (A^-2)"
137        smear_parameters_Ellip2D[4] = "SLD solvent"
138        smear_parameters_Ellip2D[5] = "Background"
139        smear_parameters_Ellip2D[6] = "Axis Theta"
140        smear_parameters_Ellip2D[7] = "Axis Phi"
141        smear_parameters_Ellip2D[8] = "Sigma of polydisp in R_a [Angstrom]"
142        smear_parameters_Ellip2D[9] = "Sigma of polydisp in R_b [Angstrom]"
143        smear_parameters_Ellip2D[10] = "Sigma of polydisp in Theta [rad]"
144        smear_parameters_Ellip2D[11] = "Sigma of polydisp in Phi [rad]"
145       
146//      smear_parameters_Ellip2D[12] = "Num of polydisp points"
147
148       
149        Edit smear_parameters_Ellip2D,smear_coef_Ellip2D                                       
150       
151        // generate the triplet representation
152        Duplicate/O $(str+"_qx") smeared_Ellip2D
153        SetScale d,0,0,"1/cm",smeared_Ellip2D                                   
154               
155        Variable/G gs_Ellip2D=0
156        gs_Ellip2D := fSmearedEllipsoid2D(smear_coef_Ellip2D,smeared_Ellip2D)           //wrapper to fill the STRUCT
157       
158        Display $(str+"_qy") vs $(str+"_qx")
159        modifygraph log=0
160        ModifyGraph mode=3,marker=16,zColor($(str+"_qy"))={smeared_Ellip2D,*,*,YellowHot,0}
161        ModifyGraph standoff=0
162        ModifyGraph width={Aspect,1}
163        ModifyGraph lowTrip=0.001
164        Label bottom "qx (A\\S-1\\M)"
165        Label left "qy (A\\S-1\\M)"
166        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
167       
168        // generate the matrix representation
169        Duplicate/O $(str+"_qx"), sm_qx
170        Duplicate/O $(str+"_qy"), sm_qy         // I can't use local variables in dependencies, so I need the name (that I can't get)
171       
172        // generate the matrix representation
173        ConvertQxQy2Mat(sm_qx,sm_qy,smeared_Ellip2D,"sm_Ellip2D_mat")
174        Duplicate/O $"sm_Ellip2D_mat",$"sm_Ellip2D_lin"                 //keep a linear-scaled version of the data
175        // _mat is for display, _lin is the real calculation
176
177        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation
178        Variable/G gs_Ellip2Dmat=0
179        gs_Ellip2Dmat := UpdateQxQy2Mat(sm_qx,sm_qy,smeared_Ellip2D,sm_Ellip2D_lin,sm_Ellip2D_mat)
180       
181        SetDataFolder root:
182        AddModelToStrings("SmearedEllipsoid2D","smear_coef_Ellip2D","smear_parameters_Ellip2D","Ellip2D")
183End
184
185
186
187//
188//  Fit function that is actually a wrapper to dispatch the calculation to N threads
189//
190// nthreads is 1 or an even number, typically 2
191// it doesn't matter if npt is odd. In this case, fractional point numbers are passed
192// and the wave indexing works just fine - I tested this with test waves of 7 and 8 points
193// and the points "2.5" and "3.5" evaluate correctly as 2 and 3
194//
195Function Ellipsoid2D(cw,zw,xw,yw) : FitFunc
196        Wave cw,zw,xw,yw
197       
198#if exists("Ellipsoid_2DX")                     //to hide the function if XOP not installed
199
200        Make/O/D/N=13 Ellip2D_tmp
201        Ellip2D_tmp[0,11] = cw
202        Ellip2D_tmp[12] = 25
203        Ellip2D_tmp[5] = 0              //pass in a zero background and add it in later
204       
205        MultiThread zw= Ellipsoid_2DX(Ellip2D_tmp,xw,yw) + cw[5]
206       
207#endif
208       
209        return(0)
210End
211
212
213/////////////////////smeared functions //////////////////////
214
215Function SmearedEllipsoid2D(s)
216        Struct ResSmear_2D_AAOStruct &s
217       
218//// non-threaded, but generic calculation
219//// the last param is nord     
220//      Smear_2DModel_PP(Ellipsoid2D_noThread,s,10)
221
222
223//// the last param is nord
224        SmearedEllipsoid2D_THR(s,10)           
225
226        return(0)
227end
228
229// for the plot dependency only
230Function fSmearedEllipsoid2D(coefW,resultW)
231        Wave coefW,resultW
232       
233        String str = getWavesDataFolder(resultW,0)
234        String DF="root:"+str+":"
235       
236        WAVE qx = $(DF+str+"_qx")
237        WAVE qy = $(DF+str+"_qy")
238        WAVE qz = $(DF+str+"_qz")
239        WAVE sQpl = $(DF+str+"_sQpl")
240        WAVE sQpp = $(DF+str+"_sQpp")
241        WAVE shad = $(DF+str+"_fs")
242       
243        STRUCT ResSmear_2D_AAOStruct s
244        WAVE s.coefW = coefW   
245        WAVE s.zw = resultW     
246        WAVE s.xw[0] = qx
247        WAVE s.xw[1] = qy
248        WAVE s.qz = qz
249        WAVE s.sQpl = sQpl
250        WAVE s.sQpp = sQpp
251        WAVE s.fs = shad
252       
253        Variable err
254        err = SmearedEllipsoid2D(s)
255       
256        return (0)
257End
258
259//
260// NON-THREADED IMPLEMENTATION
261// -- same as threaded, but no MultiThread KW
262//
263ThreadSafe Function Ellipsoid2D_noThread(cw,zw,xw,yw)
264        Wave cw,zw,xw,yw
265       
266        //      Variable t1=StopMSTimer(-2)
267
268#if exists("Ellipsoid_2DX")                     //to hide the function if XOP not installed
269
270        Make/O/D/N=13 Ellip2D_tmp
271        Ellip2D_tmp[0,11] = cw
272        Ellip2D_tmp[12] = 25
273        Ellip2D_tmp[5] = 0              //pass in a zero background and add it in later
274       
275        zw= Ellipsoid_2DX(Ellip2D_tmp,xw,yw) + cw[5]
276       
277#endif
278
279//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
280       
281        return(0)
282End
283
284//
285// this is the threaded version, that dispatches the calculation out to threads
286//
287// must be written specific to each 2D function
288//
289Function SmearedEllipsoid2D_THR(s,nord)
290        Struct ResSmear_2D_AAOStruct &s
291        Variable nord
292       
293        String weightStr,zStr
294       
295// create all of the necessary quadrature waves here - rather than inside a threadsafe function
296        switch(nord)   
297                case 5:         
298                        weightStr="gauss5wt"
299                        zStr="gauss5z"
300                        if (WaveExists($weightStr) == 0)
301                                Make/O/D/N=(nord) $weightStr,$zStr
302                                Make5GaussPoints($weightStr,$zStr)     
303                        endif
304                        break                           
305                case 10:               
306                        weightStr="gauss10wt"
307                        zStr="gauss10z"
308                        if (WaveExists($weightStr) == 0)
309                                Make/O/D/N=(nord) $weightStr,$zStr
310                                Make10GaussPoints($weightStr,$zStr)     
311                        endif
312                        break                           
313                case 20:               
314                        weightStr="gauss20wt"
315                        zStr="gauss20z"
316                        if (WaveExists($weightStr) == 0)
317                                Make/O/D/N=(nord) $weightStr,$zStr
318                                Make20GaussPoints($weightStr,$zStr)     
319                        endif
320                        break
321                default:                                                       
322                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                       
323        endswitch
324       
325        Wave/Z wt = $weightStr
326        Wave/Z xi = $zStr               // wave references to pass
327
328        Variable npt=numpnts(s.xw[0])
329        Variable i,nthreads= ThreadProcessorCount
330        variable mt= ThreadGroupCreate(nthreads)
331
332        Variable t1=StopMSTimer(-2)
333       
334        for(i=0;i<nthreads;i+=1)
335//              Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
336                ThreadStart mt,i,SmearedEllipsoid2D_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)
337        endfor
338
339        do
340                variable tgs= ThreadGroupWait(mt,100)
341        while( tgs != 0 )
342
343        variable dummy= ThreadGroupRelease(mt)
344       
345// comment out the threading + uncomment this for testing to make sure that the single thread works
346//      nThreads=1
347//      SmearedEllipsoid2D_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)
348
349        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
350       
351        return(0)
352end
353
354//
355// - worker function for threads of Sphere2D
356//
357ThreadSafe Function SmearedEllipsoid2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord)
358        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi
359        Variable pt1,pt2,nord
360       
361// now passed in....
362//      Wave wt = $weightStr
363//      Wave xi = $zStr         
364
365        Variable ii,jj,kk,num
366        Variable qx,qy,qz,qval,sx,sy,fs
367        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
368       
369        Variable normFactor,phi,theta,maxSig,numStdDev=3
370       
371/// keep these waves local
372        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW
373       
374        // now just loop over the points as specified
375       
376        answer=0
377       
378        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt
379
380        //loop over q-values
381        for(ii=pt1;ii<(pt2+1);ii+=1)
382               
383                qx = qxw[ii]
384                qy = qyw[ii]
385                qz = qzw[ii]
386                qval = sqrt(qx^2+qy^2+qz^2)
387                spl = sxw[ii]
388                spp = syw[ii]
389                fs = fsw[ii]
390               
391                normFactor = 2*pi*spl*spp
392               
393                phi = FindPhi(qx,qy)
394               
395                apl = -numStdDev*spl + qval             //parallel = q integration limits
396                bpl = numStdDev*spl + qval
397                app = -numStdDev*spp + phi              //perpendicular = phi integration limits
398                bpp = numStdDev*spp + phi
399               
400                //make sure the limits are reasonable.
401                if(apl < 0)
402                        apl = 0
403                endif
404                // do I need to specially handle limits when phi ~ 0?
405       
406               
407                sumOut = 0
408                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
409                        phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2
410
411                        sumIn=0
412                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
413
414                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2
415                               
416                                FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi
417                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not
418                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
419                               
420                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) )
421                                res_tot[kk] /= normFactor
422//                              res_tot[kk] *= fs
423
424                        endfor
425                       
426                        Ellipsoid2D_noThread(coef,fcnRet,xptw,yptw)                     //fcn passed in is an AAO
427                       
428                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0]
429                        fcnRet *= wt[jj]*wt*res_tot
430                        //
431                        answer += (bpl-apl)/2.0*sum(fcnRet)             //
432                endfor
433
434                answer *= (bpp-app)/2.0
435                zw[ii] = answer
436        endfor
437       
438        return(0)
439end
440
441
Note: See TracBrowser for help on using the repository browser.