source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_USANS_SlitSmearing_v40.ipf @ 1078

Last change on this file since 1078 was 1078, checked in by srkline, 5 years ago

more utilities and bug fixes to handle:
(1) generation of DIV files
(2) generation and loading of slit-smeared VSANS data sets

Loading of VSANS data sets and generating the resolution smearing matrix required adding a conditional compile statement to the PlotUtilsMacro? file, to distinguish VSANS data - crudely done now by checking to see that VSANS is defined (which will NOT be the case for a "normal" Analysis experiment. This will need to be revisited in the future.

File size: 14.1 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 V_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        V_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) ? V_CalcW1(p,q)  : 0
79        print "Calculating W2..."
80        W2mat = (p+1 <= q ) && (q <= USANS_N-1) ?  V_CalcW2(p,q) : 0
81        print "Calculating Remainders..."
82        Rmat = (q == USANS_N-1) ? V_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 V_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        WAVE dQvW = $(baseStr+"_dQv")
137        Duplicate/O/R=[pt1,pt2] qval $(baseStr+"_qt")
138        Duplicate/O/R=[pt1,pt2] ival $(baseStr+"_it")
139        Duplicate/O/R=[pt1,pt2] sval $(baseStr+"_st")
140        Duplicate/O/R=[pt1,pt2] dQvW $(baseStr+"_dQvt")
141        WAVE qt = $(baseStr+"_qt")                      //these are trimmed based on the cursor points
142        WAVE it = $(baseStr+"_it")
143        WAVE st = $(baseStr+"_st")
144        WAVE dQvtW = $(baseStr+"_dQvt")
145       
146        //Variable/G USANS_m = EnterSlope(baseStr)
147        Variable/G USANS_m = -4
148        Variable/G USANS_slope_numpnts = 15
149        // must reset the global for baseStr to this set, otherwise it will look for the last set loaded
150        String/G root:Packages:NIST:USANS_basestr = basestr             //this is the "current data" for slope and weights
151
152       
153        Print "---- Calculating Weighting Matrix for USANS Data ----"
154       
155        V_EnterSlope(basestr)
156       
157        //Deal with broken Pauseforuser
158        SetDataFolder $("root:"+basestr)
159       
160        print USANS_m
161
162        if (USANS_m == 999)
163                KillWaves/Z  $(basestr+"_res"),W1Mat,W2mat,Rmat
164                return(1)
165        endif
166       
167        Variable tref = startMSTimer
168        print "Calculating W1..."
169        W1mat = (p <= q ) && (q <= matN-2) ? V_CalcW1(p,q)  : 0
170        print "Calculating W2..."
171        W2mat = (p+1 <= q ) && (q <= matN-1) ?  V_CalcW2(p,q) : 0
172        print "Calculating Remainders..."
173        Rmat = (q == matN-1) ? V_CalcR(p) : 0
174//      print "Summing weights..."
175        tmpWeights = W1mat + W2mat + Rmat
176        print "Done"
177        Variable ms = stopMSTimer(tref)
178        print "Time elapsed = ", ms/1e6, "s"
179       
180        //  - put the smaller matrix into the larger one (padded w/zero)
181        // get the smaller triangular weights into the proper place in the full matrix
182        //
183        // this is necessary since the original dependency (as plotted) was set up using the
184        // full data set in the STRUCT
185        fullWeights = 0
186        fullWeights[pt1,pt2][pt1,pt2] = tmpWeights[p-pt1][q-pt1]
187       
188       
189        // put the point range of the matrix in a wave note
190        String nStr=""
191        sprintf nStr,"P1=%d;P2=%d;",pt1,pt2
192        Note/K fullWeights ,nStr
193       
194        return(0)
195End
196
197// trimStr = "" if the full (original) data set is to be used, == "t" if the set is trimmed to cursors
198Function V_EnterSlope(baseStr)
199        String baseStr
200       
201//      NVAR USANS_m,USANS_slope_numpnts
202        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
203        NVAR USANS_slope_numpnts = $("root:"+baseStr+":USANS_slope_numpnts")
204        SVAR trimStr = $("root:"+baseStr+":trimStr")
205       
206//      Variable slope=-4
207
208//      Prompt slope "Enter a slope for the file \""+ baseStr + "\""
209//      DoPrompt "Enter Slope", slope
210//              If (V_Flag)
211//                      return (999)                    //return a bogus slope if the user canceled
212//              Endif
213//      print "slope=", slope
214//      return slope
215
216        NewPanel /W=(600,300,1000,700)/N=USANS_Slope as "USANS Slope Extrapolation"
217        SetDrawLayer UserBack
218        Button button_OK,pos={270,360},size={100,20},title="Accept Slope",font="Geneva"
219        Button button_OK,proc=V_USANS_Slope_ButtonProc
220        Button button_SlopeCalc,pos={270,310},size={100,20},title="Calculate",font="Geneva"
221        Button button_SlopeCalc,proc=V_USANS_Slope_ButtonProc
222        SetVariable setvar_numpnts, pos={20,310}, size={130,19}, title="# Points",fSize=13
223        SetVariable setvar_numpnts, value= USANS_slope_numpnts
224        SetVariable setvar_Slope,pos={160,310},size={90,19},title="Slope",fSize=13
225        SetVariable setvar_Slope,limits={-inf,inf,0},value= USANS_m
226        SetVariable setvar_slope, proc=V_USANS_ManualSlope
227       
228        Display/W=(9,6,402,305)/HOST=USANS_Slope $(basestr+"_i"+trimStr) vs $(basestr+"_q"+trimStr)
229        RenameWindow #,SlopePlot
230        ErrorBars/T=0/W=USANS_Slope#SlopePlot $(basestr+"_i"+trimStr), Y wave=($(basestr+"_s"+trimStr),$(basestr+"_s"+trimStr))
231        ModifyGraph log=1
232        ModifyGraph mode=3,msize=3,marker=1,rgb=(0,0,0)
233//      legend
234        SetActiveSubwindow ##
235       
236        //Print TraceNameList("USANS_Slope#SlopePlot",";",1)
237       
238        //USANS_CalculateSlope(basestr,USANS_slope_numpnts)
239       
240        PauseForUser  USANS_Slope       
241       
242       
243//      return slope
244       
245End
246
247Function V_USANS_CalculateSlope(basestr, nend)
248        String basestr
249        Variable nend
250       
251//      NVAR USANS_m
252        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
253        SVAR trimStr = $("root:"+baseStr+":trimStr")
254       
255        Wave iw = $(basestr+"_i"+trimStr)
256        Wave qw = $(basestr+"_q"+trimStr)
257        Wave sw = $(basestr+"_s"+trimStr)
258
259        Variable num_extr=25
260        // Make extra waves for extrapolation
261        // Taken from DSM_SetExtrWaves
262               
263        Make/O/D/N=(num_extr) extr_hqq,extr_hqi
264        extr_hqi=1              //default values
265
266        //set the q-range
267        Variable qmax,num
268
269        num=numpnts(qw)
270        qmax=6*qw[num-1]
271       
272        // num-1-nend puts the extrapolation over the data set, as well as extending it
273        extr_hqq = qw[num-1-nend] + x * (qmax-qw[num-1-nend])/num_extr
274               
275       
276        // Modifed from DSM_DoExtraploate in LakeDesmearing_JB.ipf
277        // which is part of the USANS Reduction macros
278       
279               
280        //      Wave extr_lqi=extr_lqi
281        //      Wave extr_lqq=extr_lqq
282        //      Wave extr_hqi=extr_hqi
283        //      Wave extr_hqq=extr_hqq
284                Variable/G V_FitMaxIters=300
285                Variable/G V_fitOptions=4               //suppress the iteration window
286                Variable retval
287                num=numpnts(iw)
288       
289                Make/O/D P_coef={0,1,-4}                        //input
290                Make/O/T Constr={"K2<0","K2 > -8"}
291                //(set background to zero and hold fixed)
292                // initial guess
293                P_coef[1] = iw[num-1]/qw[num-1]^P_coef[2]
294               
295                CurveFit/H="100" Power kwCWave=P_coef  iw[(num-1-nend),(num-1)] /X=qw /I=1 /W=sw /C=constr
296                extr_hqi=P_coef[0]+P_coef[1]*extr_hqq^P_coef[2]
297       
298                AppendToGraph/W=USANS_Slope#SlopePlot extr_hqi vs extr_hqq
299                ModifyGraph/W=USANS_Slope#SlopePlot lsize(extr_hqi)=2
300       
301                Printf "Smeared Power law exponent = %g\r",P_coef[2]
302                Printf "**** For Desmearing, use a Power law exponent of %5.1f\r",P_coef[2]-1
303       
304                retVal = P_coef[2]-1                   
305                return(retVal)
306               
307End
308
309// adjust the slope manually on the graph (not fitted!!) just for visualization
310//
311// can be turned off by removing the action proc in the setVariable
312//
313Function V_USANS_ManualSlope(ctrlName,varNum,varStr,varName) : SetVariableControl
314        String ctrlName
315        Variable varNum
316        String varStr           
317        String varName
318               
319        SVAR baseStr = root:Packages:NIST:USANS_basestr
320        NVAR nFit = $("root:"+baseStr+":USANS_slope_numpnts")
321        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
322        WAVE/Z extr_hqq = $("root:"+baseStr+":extr_hqq")
323        WAVE/Z extr_hqi = $("root:"+baseStr+":extr_hqi")
324        SVAR trimStr = $("root:"+baseStr+":trimStr")
325       
326        if(waveExists(extr_hqq)==0 || waveExists(extr_hqi)==0)
327                return(0)               // encourage the user to try to fit first...
328        endif
329       
330        Wave iw = $("root:"+basestr+":"+basestr+"_i"+trimStr)
331        Wave qw = $("root:"+basestr+":"+basestr+"_q"+trimStr)
332
333        Variable matchPt,num=numpnts(iw),correctedSlope
334        correctedSlope = USANS_m + 1
335       
336        matchPt = iw[num-1-nFit]/(qw[num-1-nFit]^correctedSlope)
337       
338        extr_hqi=matchPt*extr_hqq^correctedSlope
339       
340        return(0)
341End
342
343Function V_USANS_Slope_ButtonProc(ctrlName) : ButtonControl
344        String ctrlName
345       
346        SVAR basestr =  root:Packages:NIST:USANS_basestr
347       
348        NVAR USANS_m = $("root:"+baseStr+":USANS_m")
349        NVAR USANS_slope_numpnts = $("root:"+baseStr+":USANS_slope_numpnts")
350       
351        strswitch (ctrlName)
352                case "button_OK":
353                        ControlUpdate/W=USANS_Slope setvar_Slope
354                        DoWindow/K USANS_Slope
355                        break
356                case "button_SlopeCalc":
357                        USANS_m = V_USANS_CalculateSlope(basestr, USANS_slope_numpnts)
358                        ControlUpdate/W=USANS_Slope setvar_Slope
359                        break
360        endswitch
361         
362
363End
364
365
366//TODO:
367// -- does this work for different values of dQv?
368// - if so, which index fo I need to use? i or j?
369//
370Function V_CalcW1(i,j)
371        Variable i,j
372       
373        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
374        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
375        SetDataFolder $("root:"+USANS_basestr)
376
377       
378// replaced with a wave
379//      NVAR dQv = USANS_dQv
380       
381        Variable UU,UL,dqj,rU,rL,wU,wL,dqw
382        Wave Qval = $(USANS_basestr+"_q"+trimStr)
383        Wave dQv = $(USANS_basestr+"_dQv"+trimStr)
384       
385        UU =sqrt(Qval[j+1]^2-Qval[i]^2)
386        UL = sqrt(Qval[j]^2-Qval[i]^2)
387        dqw = Qval[j+1]-Qval[j]
388        rU = sqrt(UU^2+Qval[i]^2)
389        rL = sqrt(UL^2+Qval[i]^2)
390       
391        wU = (1.0/dQv[i])*(Qval[j+1]*UU/dqw - 0.5*UU*rU/dqw - 0.5*Qval[i]^2*ln(UU+rU)/dqw )
392        wL = (1.0/dQv[i])*(Qval[j+1]*UL/dqw - 0.5*UL*rL/dqw - 0.5*Qval[i]^2*ln(UL+rL)/dqw )
393       
394        Return wU-wL
395
396End
397
398//TODO:
399// -- does this work for different values of dQv?
400// - if so, which index fo I need to use? i or j?
401//
402Function V_CalcW2(i,j)
403        Variable i,j
404       
405        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
406        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
407        SetDataFolder $("root:"+USANS_basestr)
408
409
410// replaced with a wave
411//      NVAR dQv = USANS_dQv
412       
413        variable UU,UL,dqw,rU,rL,wU,wL
414       
415        Wave Qval = $(USANS_basestr+"_q"+trimStr)
416        Wave dQv = $(USANS_basestr+"_dQv"+trimStr)
417
418        UU = sqrt(Qval[j]^2-Qval[i]^2)                 
419        UL = sqrt(Qval[j-1]^2-Qval[i]^2)               
420        dqw = Qval[j]-Qval[j-1]                 
421        rU = sqrt(UU^2+Qval[i]^2)
422        rL = sqrt(UL^2+Qval[i]^2)
423        wU = (1.0/dQv[i])*( -Qval[j-1]*UU/dqw + 0.5*UU*rU/dqw + 0.5*Qval[i]^2*ln(UU+rU)/dqw )
424        wL = (1.0/dQv[i])*( -Qval[j-1]*UL/dqw + 0.5*UL*rL/dqw + 0.5*Qval[i]^2*ln(UL+rL)/dqw )
425
426        Return wU-wL
427
428End
429
430//TODO:
431// -- does this work for different values of dQv?
432//
433Function V_CalcR(i)
434        Variable i
435
436        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
437        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
438        SetDataFolder $("root:"+USANS_basestr)
439
440
441        NVAR m = USANS_m
442        NVAR N = USANS_N
443// replaced with a wave
444//      NVAR dQv = USANS_dQv
445       
446        Variable retval
447        Wave Qval = $(USANS_basestr+"_q"+trimStr)
448        Wave Ival = $(USANS_basestr+"_i"+trimStr)
449        Variable/G USANS_intQpt = Qval[i]
450        Wave dQv = $(USANS_basestr+"_dQv"+trimStr)
451
452       
453        Variable lower = sqrt(qval[N-1]^2-qval[i]^2)
454        Variable upper
455        upper = dQv[i]
456
457        if (i == N)
458                lower = 0
459        endif
460       
461        retval = Integrate1D(V_Remainder,lower,upper)
462       
463        retval *= 1/dQv[i]
464       
465        Return retval
466
467End
468
469Function V_Remainder(i)
470       
471        Variable i
472       
473        SVAR USANS_basestr=root:Packages:NIST:USANS_basestr
474        SVAR trimStr = $("root:"+USANS_baseStr+":trimStr")
475        SetDataFolder $("root:"+USANS_basestr)
476
477        NVAR m = USANS_m
478        NVAR qi = USANS_intQpt
479        NVAR N = USANS_N
480        WAVE Qval = $(USANS_basestr+"_q"+trimStr)       
481        Variable retVal
482       
483        retVal=Qval[N-1]^(-m)*(i^2+qi^2)^(m/2)
484
485        return retval
486
487End
488
Note: See TracBrowser for help on using the repository browser.