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

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

cleanup of TODO items in code, no noteworthy changes

prepare a test release package for the January startup of VSANS, not a general release (since no changes to SANS or USANS)

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