source: sans/Dev/trunk/NCNR_User_Procedures/SANS_Analysis/Models/USANS_SlitSmearing_v40.ipf @ 325

Last change on this file since 325 was 325, checked in by ajj, 14 years ago

Adding SANS Analysis Models to new Dev tree

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.