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