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