source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v4.00/USANS_SlitSmearing_v40.ipf @ 309

Last change on this file since 309 was 309, checked in by srkline, 15 years ago

Changes to how USANS resolution matrix is handled in fitting, both regular and global.
(1) - resolution matirx is re-calculated if the data is masked (with cursors)
(2) - if the matrix is appropriate, as stored in the wave note, it is not recalculated
(3) - global fit does the same, being sure to recalculate the matrix before the Global fit has started
(4) - the global fit report now has pararameter names rather than "coef_n". will need to be tweaked in future to accomodate different fit functions, and to improve the formatting.

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