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

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

Several changes:
1) added /I=1 flag in several places (mostly the invariant) so that the error wave would be interpreted as the standard deviation, not 1/s (Jae_Hie pointed this out)
2) put error checking in ProDiv? to warn if the pixel centers are more than 5 pixels from the expected 65,65 for on-center or 105,65 for the offset (usually run at 20 cm offset)
3) commented out the line in WriteQIS that outputs 2D resolution information to QxQy? data. I just don't think it's correct yet, and the 2D resolution smearing is not ready either.

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 /I=1 /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.