source: sans/Dev/trunk/NCNR_User_Procedures/SANS/Analysis/Models/USANS_SlitSmearing_v40.ipf @ 347

Last change on this file since 347 was 347, checked in by srkline, 14 years ago

fixed bug in USANS_RecalculateSlope that was looking for the last loaded USANS file when it is supposed to be looking for the selected data set. Global is now reset as needed.

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