source: sans/Dev/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/USANS_SlitSmearing_v40.ipf @ 570

Last change on this file since 570 was 570, checked in by srkline, 13 years ago

Change (1):
In preparation for release, updated pragma IgorVersion?=6.1 in all procedures

Change (2):
As a side benefit of requiring 6.1, we can use the MultiThread? keyword to thread any model function we like. The speed benefit is only noticeable on functions that require at least one integration and at least 100 points (resolution smearing is NOT threaded, too many threadSafe issues, too little benefit). I have chosen to use the MultiThread? only on the XOP assignment. In the Igor code there are too many functions that are not explicitly declared threadsafe, making for a mess.

File size: 13.4 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=3.00
3#pragma IgorVersion=6.1
4
5//Functions for doing USANS Slit smearing by method of weight matrix
6//Routines originally from J Barker fortran code
7//Translated to IGOR by M-H Kim
8//Updated to use IGOR features and integrated into SANS Macros by A J Jackson
9//
10// AJJ
11// July 26 2007 - Modified functions to work with new SANS Analysis Macros
12//     - pass basestr to functions to determine wave names and avoid globals
13//     - pass dQv to functions to  avoid globals
14//     - pass N to CalcR to avoid globals
15
16//////////////
17//
18//// July 2008 SRK
19//
20// For fitting withe cursors, the matrix must be recalculated for the exact range of the data
21// - but the inital dependency of the model (any number of them) was (were) set up to use the full matrix
22// -- so:
23//      $_res is the resolution for the dependencies. It is always the full size of the data set. If cursors are used,
24//                      it is padded with zeroes at the edges to fill.
25//
26// weights_save: is a pristine copy of the resolution matrix for the full data set, as when loaded
27// $_qt, $_qi, $_qs, $_rest: are "trimmed" sets that use the range specified by the cursors. They are created
28//                      at load time and are initially the full data range.
29//
30// there is a wave note attached to $_res that has the current point range for the matrix. This is what the
31//              $_rest matrix is expected to be too, so keep these two matrices in sync.
32//
33//      during fitting, $_rest is used in the structure if cursors are used, since the matrix must be the same dimension (N)
34//                      as the (trimmed) data range. It may be the full data range if the cursors are at the ends of the data set.
35//                      If no cursors are used, then the $_res wave is used in the structure
36//////////
37
38
39
40
41// called only by the main file loader
42//
43Function USANS_CalcWeights(basestr, dQv)
44        String basestr
45        Variable dQv
46
47        Variable/G USANS_N=numpnts($(basestr+"_q"))
48        Variable/G USANS_dQv = dQv
49        String/G root:Packages:NIST:USANS_basestr = basestr             //this is the "current data" for slope and weights
50
51        Make/D/O/N=(USANS_N,USANS_N) $(basestr+"_res")
52        Make/D/O/N=(USANS_N,USANS_N) W1mat
53        Make/D/O/N=(USANS_N,USANS_N) W2mat
54        Make/D/O/N=(USANS_N,USANS_N) Rmat
55        Wave weights = $(basestr+"_res")
56
57        //Variable/G USANS_m = EnterSlope(baseStr)
58        Variable/G USANS_m = -4
59        Variable/G USANS_slope_numpnts = 15
60        String/G trimStr=""                     //null string if NOT using cursors, "t" if yes (and re-calculating)
61       
62        Print "---- Calculating Weighting Matrix for USANS Data ----"
63       
64        EnterSlope(basestr)
65       
66        //Deal with broken Pauseforuser
67        SetDataFolder $("root:"+basestr)
68       
69        print USANS_m
70
71        if (USANS_m == 999)
72                KillWaves/Z  $(basestr+"_res"),W1Mat,W2mat,Rmat
73                return(1)
74        endif
75       
76        Variable tref = startMSTimer
77        print "Calculating W1..."
78        W1mat = (p <= q ) && (q <= USANS_N-2) ? CalcW1(p,q)  : 0
79        print "Calculating W2..."
80        W2mat = (p+1 <= q ) && (q <= USANS_N-1) ?  CalcW2(p,q) : 0
81        print "Calculating Remainders..."
82        Rmat = (q == USANS_N-1) ? CalcR(p) : 0
83//      print "Summing weights..."
84        weights = W1mat + W2mat + Rmat
85        print "Done"
86        Variable ms = stopMSTimer(tref)
87        print "Time elapsed = ", ms/1e6, "s"
88       
89        // put the point range of the matrix in a wave note (this is the full range)
90        String nStr=""
91        sprintf nStr,"P1=%d;P2=%d;",0,USANS_N-1
92        Note weights ,nStr
93       
94        // make a set of "trimmed" data that is currently the full data set
95        // but will be trimmed as needed for use with cursors
96        // this is necessary to handle the base case of cursors that use the full range of the data!
97        Duplicate/O $(baseStr+"_q") $(baseStr+"_qt")
98        Duplicate/O $(baseStr+"_i") $(baseStr+"_it")
99        Duplicate/O $(baseStr+"_s") $(baseStr+"_st")
100        Duplicate/O weights, $(basestr+"_res"+"t")
101
102        // save a copy of the untainted matrix to return to...
103        Duplicate/O weights, weights_save
104       
105        return(0)
106End
107
108// re-calculate the weighting matrix, maybe easier to do in a separate function
109//
110// supposedly I'm in the data folder, not root:, double check to make sure
111//
112Function USANS_RE_CalcWeights(baseStr,pt1,pt2)
113        String baseStr
114        Variable pt1,pt2
115
116        SetDataFolder $("root:"+baseStr)
117        NVAR USANS_dQv = USANS_dQv
118        SVAR trimStr = trimStr
119        trimStr="t"             //yes, we're re-calculating, using trimmed data sets
120       
121        Variable matN
122        matN=pt2-pt1+1          //new size of matrix
123
124        Make/D/O/N=(matN,matN) $(basestr+"_res"+trimStr)                //this is a temporary matrix
125        Make/D/O/N=(matN,matN) W1mat
126        Make/D/O/N=(matN,matN) W2mat
127        Make/D/O/N=(matN,matN) Rmat
128        Wave tmpWeights = $(basestr+"_res"+trimStr)
129        Wave fullWeights = $(basestr+"_res")
130
131        // make a trimmed set of QIS waves to match the size selected by the cursors
132        // == "t" added to the suffix
133        WAVE qval = $(baseStr+"_q")
134        WAVE ival = $(baseStr+"_i")
135        WAVE sval = $(baseStr+"_s")
136        Duplicate/O/R=[pt1,pt2] qval $(baseStr+"_qt")
137        Duplicate/O/R=[pt1,pt2] ival $(baseStr+"_it")
138        Duplicate/O/R=[pt1,pt2] sval $(baseStr+"_st")
139        WAVE qt = $(baseStr+"_qt")                      //these are trimmed based on the cursor points
140        WAVE it = $(baseStr+"_it")
141        WAVE st = $(baseStr+"_st")
142       
143        //Variable/G USANS_m = EnterSlope(baseStr)
144        Variable/G USANS_m = -4
145        Variable/G USANS_slope_numpnts = 15
146        // must reset the global for baseStr to this set, otherwise it will look for the last set loaded
147        String/G root:Packages:NIST:USANS_basestr = basestr             //this is the "current data" for slope and weights
148
149       
150        Print "---- Calculating Weighting Matrix for USANS Data ----"
151       
152        EnterSlope(basestr)
153       
154        //Deal with broken Pauseforuser
155        SetDataFolder $("root:"+basestr)
156       
157        print USANS_m
158
159        if (USANS_m == 999)
160                KillWaves/Z  $(basestr+"_res"),W1Mat,W2mat,Rmat
161                return(1)
162        endif
163       
164        Variable tref = startMSTimer
165        print "Calculating W1..."
166        W1mat = (p <= q ) && (q <= matN-2) ? CalcW1(p,q)  : 0
167        print "Calculating W2..."
168        W2mat = (p+1 <= q ) && (q <= matN-1) ?  CalcW2(p,q) : 0
169        print "Calculating Remainders..."
170        Rmat = (q == matN-1) ? CalcR(p) : 0
171//      print "Summing weights..."
172        tmpWeights = W1mat + W2mat + Rmat
173        print "Done"
174        Variable ms = stopMSTimer(tref)
175        print "Time elapsed = ", ms/1e6, "s"
176       
177        //  - put the smaller matrix into the larger one (padded w/zero)
178        // get the smaller triangular weights into the proper place in the full matrix
179        //
180        // this is necessary since the original dependency (as plotted) was set up using the
181        // full data set in the STRUCT
182        fullWeights = 0
183        fullWeights[pt1,pt2][pt1,pt2] = tmpWeights[p-pt1][q-pt1]
184       
185       
186        // put the point range of the matrix in a wave note
187        String nStr=""
188        sprintf nStr,"P1=%d;P2=%d;",pt1,pt2
189        Note/K fullWeights ,nStr
190       
191        return(0)
192End
193
194// trimStr = "" if the full (original) data set is to be used, == "t" if the set is trimmed to cursors
195Function EnterSlope(baseStr)
196        String baseStr
197       
198//      NVAR USANS_m,USANS_slope_numpnts
199        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
200        NVAR USANS_slope_numpnts = $("root:"+baseStr+":USANS_slope_numpnts")
201        SVAR trimStr = $("root:"+baseStr+":trimStr")
202       
203//      Variable slope=-4
204
205//      Prompt slope "Enter a slope for the file \""+ baseStr + "\""
206//      DoPrompt "Enter Slope", slope
207//              If (V_Flag)
208//                      return (999)                    //return a bogus slope if the user canceled
209//              Endif
210//      print "slope=", slope
211//      return slope
212
213        NewPanel /W=(600,300,1000,700)/N=USANS_Slope as "USANS Slope Extrapolation"
214        SetDrawLayer UserBack
215        Button button_OK,pos={270,360},size={100,20},title="Accept Slope",font="Geneva"
216        Button button_OK,proc=USANS_Slope_ButtonProc
217        Button button_SlopeCalc,pos={270,310},size={100,20},title="Calculate",font="Geneva"
218        Button button_SlopeCalc,proc=USANS_Slope_ButtonProc
219        SetVariable setvar_numpnts, pos={20,310}, size={130,19}, title="# Points",fSize=13
220        SetVariable setvar_numpnts, value= USANS_slope_numpnts
221        SetVariable setvar_Slope,pos={160,310},size={90,19},title="Slope",fSize=13
222        SetVariable setvar_Slope,limits={-inf,inf,0},value= USANS_m
223        SetVariable setvar_slope, proc=USANS_ManualSlope
224       
225        Display/W=(9,6,402,305)/HOST=USANS_Slope $(basestr+"_i"+trimStr) vs $(basestr+"_q"+trimStr)
226        RenameWindow #,SlopePlot
227        ErrorBars/T=0/W=USANS_Slope#SlopePlot $(basestr+"_i"+trimStr), Y wave=($(basestr+"_s"+trimStr),$(basestr+"_s"+trimStr))
228        ModifyGraph log=1
229        ModifyGraph mode=3,msize=3,marker=1,rgb=(0,0,0)
230//      legend
231        SetActiveSubwindow ##
232       
233        //Print TraceNameList("USANS_Slope#SlopePlot",";",1)
234       
235        //USANS_CalculateSlope(basestr,USANS_slope_numpnts)
236       
237        PauseForUser  USANS_Slope       
238       
239       
240//      return slope
241       
242End
243
244Function USANS_CalculateSlope(basestr, nend)
245        String basestr
246        Variable nend
247       
248//      NVAR USANS_m
249        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
250        SVAR trimStr = $("root:"+baseStr+":trimStr")
251       
252        Wave iw = $(basestr+"_i"+trimStr)
253        Wave qw = $(basestr+"_q"+trimStr)
254        Wave sw = $(basestr+"_s"+trimStr)
255
256        Variable num_extr=25
257        // Make extra waves for extrapolation
258        // Taken from DSM_SetExtrWaves
259               
260        Make/O/D/N=(num_extr) extr_hqq,extr_hqi
261        extr_hqi=1              //default values
262
263        //set the q-range
264        Variable qmax,num
265
266        num=numpnts(qw)
267        qmax=6*qw[num-1]
268       
269        // num-1-nend puts the extrapolation over the data set, as well as extending it
270        extr_hqq = qw[num-1-nend] + x * (qmax-qw[num-1-nend])/num_extr
271               
272       
273        // Modifed from DSM_DoExtraploate in LakeDesmearing_JB.ipf
274        // which is part of the USANS Reduction macros
275       
276               
277        //      Wave extr_lqi=extr_lqi
278        //      Wave extr_lqq=extr_lqq
279        //      Wave extr_hqi=extr_hqi
280        //      Wave extr_hqq=extr_hqq
281                Variable/G V_FitMaxIters=300
282                Variable/G V_fitOptions=4               //suppress the iteration window
283                Variable retval
284                num=numpnts(iw)
285       
286                Make/O/D P_coef={0,1,-4}                        //input
287                Make/O/T Constr={"K2<0","K2 > -8"}
288                //(set background to zero and hold fixed)
289                // initial guess
290                P_coef[1] = iw[num-1]/qw[num-1]^P_coef[2]
291               
292                CurveFit/H="100" Power kwCWave=P_coef  iw[(num-1-nend),(num-1)] /X=qw /W=sw /C=constr
293                extr_hqi=P_coef[0]+P_coef[1]*extr_hqq^P_coef[2]
294       
295                AppendToGraph/W=USANS_Slope#SlopePlot extr_hqi vs extr_hqq
296                ModifyGraph/W=USANS_Slope#SlopePlot lsize(extr_hqi)=2
297       
298                Printf "Smeared Power law exponent = %g\r",P_coef[2]
299                Printf "**** For Desmearing, use a Power law exponent of %5.1f\r",P_coef[2]-1
300       
301                retVal = P_coef[2]-1                   
302                return(retVal)
303               
304End
305
306// adjust the slope manually on the graph (not fitted!!) just for visualization
307//
308// can be turned off by removing the action proc in the setVariable
309//
310Function USANS_ManualSlope(ctrlName,varNum,varStr,varName) : SetVariableControl
311        String ctrlName
312        Variable varNum
313        String varStr           
314        String varName
315               
316        SVAR baseStr = root:Packages:NIST:USANS_basestr
317        NVAR nFit = $("root:"+baseStr+":USANS_slope_numpnts")
318        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
319        WAVE/Z extr_hqq = $("root:"+baseStr+":extr_hqq")
320        WAVE/Z extr_hqi = $("root:"+baseStr+":extr_hqi")
321        SVAR trimStr = $("root:"+baseStr+":trimStr")
322       
323        if(waveExists(extr_hqq)==0 || waveExists(extr_hqi)==0)
324                return(0)               // encourage the user to try to fit first...
325        endif
326       
327        Wave iw = $("root:"+basestr+":"+basestr+"_i"+trimStr)
328        Wave qw = $("root:"+basestr+":"+basestr+"_q"+trimStr)
329
330        Variable matchPt,num=numpnts(iw),correctedSlope
331        correctedSlope = USANS_m + 1
332       
333        matchPt = iw[num-1-nFit]/(qw[num-1-nFit]^correctedSlope)
334       
335        extr_hqi=matchPt*extr_hqq^correctedSlope
336       
337        return(0)
338End
339
340Function USANS_Slope_ButtonProc(ctrlName) : ButtonControl
341        String ctrlName
342       
343        SVAR basestr =  root:Packages:NIST:USANS_basestr
344       
345        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
346        NVAR USANS_slope_numpnts = $("root:"+baseStr+":USANS_slope_numpnts")
347       
348        strswitch (ctrlName)
349                case "button_OK":
350                        ControlUpdate/W=USANS_Slope setvar_Slope
351                        DoWindow/K USANS_Slope
352                        break
353                case "button_SlopeCalc":
354                        USANS_m = USANS_CalculateSlope(basestr, USANS_slope_numpnts)
355                        ControlUpdate/W=USANS_Slope setvar_Slope
356                        break
357        endswitch
358         
359
360End
361
362Function CalcW1(i,j)
363        Variable i,j
364       
365        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
366        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
367        SetDataFolder $("root:"+USANS_basestr)
368
369       
370       
371        NVAR dQv = USANS_dQv
372       
373        Variable UU,UL,dqj,rU,rL,wU,wL,dqw
374        Wave Qval = $(USANS_basestr+"_q"+trimStr)
375       
376        UU =sqrt(Qval[j+1]^2-Qval[i]^2)
377        UL = sqrt(Qval[j]^2-Qval[i]^2)
378        dqw = Qval[j+1]-Qval[j]
379        rU = sqrt(UU^2+Qval[i]^2)
380        rL = sqrt(UL^2+Qval[i]^2)
381       
382        wU = (1.0/dQv)*(Qval[j+1]*UU/dqw - 0.5*UU*rU/dqw - 0.5*Qval[i]^2*ln(UU+rU)/dqw )
383        wL = (1.0/dQv)*(Qval[j+1]*UL/dqw - 0.5*UL*rL/dqw - 0.5*Qval[i]^2*ln(UL+rL)/dqw )
384       
385        Return wU-wL
386
387End
388
389Function CalcW2(i,j)
390        Variable i,j
391       
392        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
393        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
394        SetDataFolder $("root:"+USANS_basestr)
395
396
397
398        NVAR dQv = USANS_dQv
399       
400        variable UU,UL,dqw,rU,rL,wU,wL
401       
402        Wave Qval = $(USANS_basestr+"_q"+trimStr)
403
404        UU = sqrt(Qval[j]^2-Qval[i]^2)                 
405        UL = sqrt(Qval[j-1]^2-Qval[i]^2)               
406        dqw = Qval[j]-Qval[j-1]                 
407        rU = sqrt(UU^2+Qval[i]^2)
408        rL = sqrt(UL^2+Qval[i]^2)
409        wU = (1.0/dQv)*( -Qval[j-1]*UU/dqw + 0.5*UU*rU/dqw + 0.5*Qval[i]^2*ln(UU+rU)/dqw )
410        wL = (1.0/dQv)*( -Qval[j-1]*UL/dqw + 0.5*UL*rL/dqw + 0.5*Qval[i]^2*ln(UL+rL)/dqw )
411
412        Return wU-wL
413
414End
415
416Function CalcR(i)
417        Variable i
418
419        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
420        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
421        SetDataFolder $("root:"+USANS_basestr)
422
423
424        NVAR m = USANS_m
425        NVAR N = USANS_N
426        NVAR dQv = USANS_dQv
427       
428        Variable retval
429        Wave Qval = $(USANS_basestr+"_q"+trimStr)
430        Wave Ival = $(USANS_basestr+"_i"+trimStr)
431        Variable/G USANS_intQpt = Qval[i]
432       
433        Variable lower = sqrt(qval[N-1]^2-qval[i]^2)
434        Variable upper = dQv
435
436        if (i == N)
437                lower = 0
438        endif
439       
440        retval = Integrate1D(Remainder,lower,upper)
441       
442        retval *= 1/dQv
443       
444        Return retval
445
446End
447
448Function Remainder(i)
449       
450        Variable i
451       
452        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
453        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
454        SetDataFolder $("root:"+USANS_basestr)
455
456        NVAR m = USANS_m
457        NVAR qi = USANS_intQpt
458        NVAR N = USANS_N
459        WAVE Qval = $(USANS_basestr+"_q"+trimStr)       
460        Variable retVal
461       
462        retVal=Qval[N-1]^(-m)*(i^2+qi^2)^(m/2)
463
464        return retval
465
466End
Note: See TracBrowser for help on using the repository browser.