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

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

some of the changes for smearing calculations.

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