1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | #pragma Version=2.20 |
---|
3 | #pragma IgorVersion=6.1 |
---|
4 | |
---|
5 | ////////////////////////////////////////////// |
---|
6 | // Igor conversion: 03 FEB 06 SRK |
---|
7 | // Revision: 03 MAR 06 SRK |
---|
8 | // |
---|
9 | // Program uses Lake's iterative technique, J. A. Lake, Acta Cryst. 23 (1967) 191-4. |
---|
10 | // to DESMEAR Infinite slit smeared USANS DATA. |
---|
11 | // |
---|
12 | // TO SEE DESCRIPTION, CHECK COMPUTER SOFTWARE LOGBOOK, P13,41-46, 50 |
---|
13 | // J. BARKER, August, 2001 |
---|
14 | // Switches from fast to slow convergence near target chi JGB 2/2003 |
---|
15 | // |
---|
16 | // steps are: |
---|
17 | // load |
---|
18 | // mask |
---|
19 | // extrapolate (find the power law of the desmeared set = smeared -1) |
---|
20 | // smooth (optional) |
---|
21 | // desmear based on dQv and chi^2 target |
---|
22 | // save result |
---|
23 | // |
---|
24 | ///////// |
---|
25 | // Waves produced at each step: (Q/I/S with the following extensions) |
---|
26 | // |
---|
27 | // Load: Uses -> nothing (Kills old waves) |
---|
28 | // Produces -> "_exp" and "_exp_orig" |
---|
29 | // |
---|
30 | // Mask: Uses -> "_exp" |
---|
31 | // Produces -> "_msk" |
---|
32 | // |
---|
33 | // Extrapolate: Uses -> nothing |
---|
34 | // Produces -> nothing |
---|
35 | // |
---|
36 | // Smooth: Uses -> "_smth" OR "_msk" OR "_exp" (in that order) |
---|
37 | // Produces -> "_smth" |
---|
38 | // |
---|
39 | // Desmear: Uses -> "_smth" OR "_msk" OR "_exp" (in that order) |
---|
40 | // Produces -> "_dsm" |
---|
41 | // |
---|
42 | //////// |
---|
43 | |
---|
44 | ///***** TO FIX ******* |
---|
45 | // - ? power law + background extrapolation (only useful for X-ray data...) |
---|
46 | // |
---|
47 | // see commented code lines for Igor 4 or Igor 5 |
---|
48 | // - there are only a few options and calls that are not Igor 4 compatible. |
---|
49 | // Igor 4 routines are currently used. |
---|
50 | |
---|
51 | //////////////////////////////////////////////////////////////////////// |
---|
52 | |
---|
53 | // main entry routine |
---|
54 | Proc Desmear() |
---|
55 | |
---|
56 | String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
57 | |
---|
58 | //check for the correct data folder, initialize if necessary |
---|
59 | // |
---|
60 | if(DataFolderExists(USANSFolder+":DSM") == 0) |
---|
61 | Execute "Init_Desmearing()" |
---|
62 | endif |
---|
63 | |
---|
64 | SetDataFolder $(USANSFolder+":DSM") |
---|
65 | //always initialize these global variables |
---|
66 | gStr1 = "" |
---|
67 | gStr2 = "" |
---|
68 | gIterStr = "" |
---|
69 | SetDataFolder root: |
---|
70 | |
---|
71 | DoWindow/F Desmear_Graph |
---|
72 | if(V_flag==0) |
---|
73 | Execute "Desmear_Graph()" |
---|
74 | endif |
---|
75 | End |
---|
76 | |
---|
77 | Proc Init_Desmearing() |
---|
78 | |
---|
79 | String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
80 | |
---|
81 | //set up the folder(s) needed |
---|
82 | NewDataFolder/O $(USANSFolder+":DSM") |
---|
83 | NewDataFolder/O $(USANSFolder+":myGlobals") //in case it wasn't created elsewhere |
---|
84 | |
---|
85 | SetDataFolder $(USANSFolder+":DSM") |
---|
86 | |
---|
87 | String/G gCurFile="" |
---|
88 | |
---|
89 | Variable/G gMaxFastIter = 100 //max number of iter in Fast convergence |
---|
90 | Variable/G gMaxSlowIter = 10000 |
---|
91 | |
---|
92 | Variable/G gNptsExtrap = 10 //points for high q extrapolation |
---|
93 | Variable/G gChi2Target = 1 //chi^2 target |
---|
94 | Variable/G gPowerM = -4 |
---|
95 | Variable/G gDqv = 0.117 //2005 measured slit height - see John |
---|
96 | Variable/G gNq = 1 |
---|
97 | Variable/G gS = 0 // global varaible for Midpnt() |
---|
98 | Variable/G gSmoothFac=0.03 |
---|
99 | |
---|
100 | Variable/G gChi2Final = 0 //chi^2 final |
---|
101 | Variable/G gIterations = 0 //total number of iterations |
---|
102 | |
---|
103 | String/G gStr1 = "" //information strings |
---|
104 | String/G gStr2 = "" |
---|
105 | String/G gIterStr = "" |
---|
106 | Variable/G gChi2Smooth = 0 |
---|
107 | |
---|
108 | Variable/G gFreshMask=1 |
---|
109 | |
---|
110 | SetDataFolder root: |
---|
111 | End |
---|
112 | |
---|
113 | //////////// Lake desmearing routines |
---|
114 | |
---|
115 | // Smears guess Y --> YS using weighting array |
---|
116 | Function DSM_Smear(Y_wave,Ys,NQ,FW) |
---|
117 | Wave Y_wave,Ys |
---|
118 | Variable nq |
---|
119 | Wave FW |
---|
120 | |
---|
121 | Variable ii,jj,summ |
---|
122 | for(ii=0;ii<nq;ii+=1) |
---|
123 | summ=0 |
---|
124 | for(jj=0;jj<nq;jj+=1) |
---|
125 | summ = summ + Y_wave[jj]*FW[ii][jj] |
---|
126 | endfor |
---|
127 | Ys[ii] = summ |
---|
128 | endfor |
---|
129 | |
---|
130 | Return (0) |
---|
131 | End |
---|
132 | |
---|
133 | // CALCULATES CHI^2 |
---|
134 | FUNCTION CHI2(ys_new,ys_exp,w,NQ) |
---|
135 | Wave ys_new,ys_exp,w |
---|
136 | Variable nq |
---|
137 | |
---|
138 | Variable CHI2,summ,ii |
---|
139 | |
---|
140 | SUMm = 0.0 |
---|
141 | for(ii=0;ii<nq;ii+=1) |
---|
142 | // if(numtype(YS_EXP[ii]) == 0 && numtype(YS_NEW[ii]) == 0) |
---|
143 | SUMm=SUMm+W[ii]*(YS_EXP[ii]-YS_NEW[ii])^2 |
---|
144 | // endif |
---|
145 | endfor |
---|
146 | |
---|
147 | CHI2=SUMm/(NQ-1) |
---|
148 | |
---|
149 | RETURN CHI2 |
---|
150 | END |
---|
151 | |
---|
152 | // Routine calculates the weights needed to convert a table |
---|
153 | // representation of a scattering function I(q) into the infinite |
---|
154 | // slit smeared table represented as Is(q) |
---|
155 | // Is(qi) = Sum(j=i-n) Wij*I(qj) |
---|
156 | // Program assumes data is extrapolated from qn to dqv using powerlaw |
---|
157 | // I(q) = Aq^m |
---|
158 | // Required input: |
---|
159 | // N : Number of data points {in Global common block } |
---|
160 | // q(N) : Array of q values {in Global common block } |
---|
161 | // dqv : Limit of slit length {in Global common block } |
---|
162 | // m : powerlaw for extrapolation. |
---|
163 | // |
---|
164 | // Routine output: |
---|
165 | // W(N,N) : Weights array |
---|
166 | // |
---|
167 | // The weights obtained from this routine can be used to calculate |
---|
168 | // Is(q) for any I(q) function. |
---|
169 | // |
---|
170 | // Calculation based upon linear interpolation between points. |
---|
171 | // See p.44, Computer Software Log book. |
---|
172 | // 11/00 JGB |
---|
173 | Function Weights_L(m,FW,Q_exp) //changed input w to FW (FW is two dimensional) |
---|
174 | Variable m |
---|
175 | Wave FW,Q_exp |
---|
176 | |
---|
177 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
178 | |
---|
179 | NVAR dqv = $(USANSFolder+":DSM:gDqv") |
---|
180 | NVAR NN = $(USANSFolder+":DSM:gNq") |
---|
181 | |
---|
182 | // Calculate Remainder fractions and put into separate array. |
---|
183 | Variable lower,ii,ss,jj |
---|
184 | Duplicate/O Q_exp $(USANSFolder+":DSM:R_wave") |
---|
185 | wave r_wave = $(USANSFolder+":DSM:R_wave") |
---|
186 | |
---|
187 | // Make/O/D/N=75 root:DSM:SS_save //debug |
---|
188 | // wave SS_save = root:DSM:SS_save |
---|
189 | |
---|
190 | Print "calculating remainders by integration..." |
---|
191 | for(ii=0;ii<NN-1;ii+=1) |
---|
192 | lower = sqrt(Q_exp[NN-1]^2-Q_exp[ii]^2) |
---|
193 | ss = Qromo(lower,dqv,Q_exp[ii],m) //new parameter list for Qromo call |
---|
194 | // SS_save[ii] = ss |
---|
195 | R_wave[ii] = (Q_exp[NN-1]^(-m)/dqv)*SS |
---|
196 | // Printf "I = %d R_wave[ii] =%g\r",ii,R_wave[ii] |
---|
197 | endfor |
---|
198 | |
---|
199 | lower = 0.0 |
---|
200 | ss = Qromo(lower,dqv,Q_exp[NN-1],m) //new parameter list for Qromo call |
---|
201 | R_wave[NN-1] = (Q_exp[NN-1]^(-m)/dqv)*SS |
---|
202 | // Printf "I = %d R_wave[ii] =%g\r",NN-1,R_wave[NN-1] |
---|
203 | // SS_save[ii] = ss |
---|
204 | |
---|
205 | // Make/O/D/N=(75,75) root:DSM:IG_save //debug |
---|
206 | // wave IG_save=root:DSM:IG_save |
---|
207 | |
---|
208 | // Zero weight matrix... then fill it |
---|
209 | FW = 0 |
---|
210 | Print "calculating full weights...." |
---|
211 | // Calculate weights |
---|
212 | for(ii=0;ii<NN;ii+=1) |
---|
213 | for(jj=ii+1;jj<NN-1;jj+=1) |
---|
214 | FW[ii][jj] = DU(ii,jj)*(1.0+Q_exp[jj]/(Q_exp[jj+1]-Q_exp[jj])) |
---|
215 | FW[ii][jj] -= (1.0/(Q_exp[jj+1]-Q_exp[jj]))*IG(ii,jj) |
---|
216 | FW[ii][jj] -= DU(ii,jj-1)*Q_exp[jj-1]/(Q_exp[jj]-Q_exp[jj-1]) |
---|
217 | FW[ii][jj] += (1.0/(Q_exp[jj]-Q_exp[jj-1]))*IG(ii,jj-1) |
---|
218 | FW[ii][jj] *= (1.0/dqv) |
---|
219 | // Printf "FW[%d][%d] = %g\r",ii,jj,FW[ii][jj] |
---|
220 | endfor |
---|
221 | endfor |
---|
222 | // |
---|
223 | // special case: I=J,I=N |
---|
224 | for(ii=0;ii<NN-1;ii+=1) |
---|
225 | FW[ii][ii] = DU(ii,ii)*(1.0+Q_exp[ii]/(Q_exp[ii+1]-Q_exp[ii])) |
---|
226 | FW[ii][ii] -= (1.0/(Q_exp[ii+1]-Q_exp[ii]))*IG(ii,ii) |
---|
227 | FW[ii][ii] *= (1.0/dqv) |
---|
228 | // Printf "FW[%d][%d] = %g\r",ii,jj,FW[ii][jj] |
---|
229 | // following line corrected for N -> NN-1 since Igor indexes from 0! (Q_exp[NN] DNE!) |
---|
230 | FW[ii][NN-1] = -DU(ii,NN-2)*Q_exp[NN-2]/(Q_exp[NN-1]-Q_exp[NN-2]) |
---|
231 | FW[ii][NN-1] += (1.0/(Q_exp[NN-1]-Q_exp[NN-2]))*IG(ii,NN-2) |
---|
232 | FW[ii][NN-1] *= (1.0/dqv) |
---|
233 | FW[ii][NN-1] += R_wave[ii] |
---|
234 | // Printf "FW[%d][%d] = %g\r",ii,NN-1,FW[ii][NN-1] |
---|
235 | endfor |
---|
236 | // |
---|
237 | // special case: I=J=N |
---|
238 | FW[NN-1][NN-1] = R_wave[NN-1] |
---|
239 | // |
---|
240 | Return (0) |
---|
241 | End |
---|
242 | |
---|
243 | //// |
---|
244 | Function Integrand_dsm(u,qi,m) |
---|
245 | Variable u,qi,m |
---|
246 | |
---|
247 | Variable integrand |
---|
248 | Integrand = (u^2+qi^2)^(m/2) |
---|
249 | return integrand |
---|
250 | end |
---|
251 | |
---|
252 | /// |
---|
253 | Function DU(ii,jj) |
---|
254 | Variable ii,jj |
---|
255 | |
---|
256 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
257 | |
---|
258 | Wave Q_exp = $(USANSFolder+":DSM:Q_exp") |
---|
259 | Variable DU |
---|
260 | |
---|
261 | // Wave DU_save=root:DSM:DU_save |
---|
262 | If (ii == jj) |
---|
263 | DU = sqrt(Q_exp[jj+1]^2-Q_exp[ii]^2) |
---|
264 | Else |
---|
265 | DU = sqrt(Q_exp[jj+1]^2-Q_exp[ii]^2) - sqrt(Q_exp[jj]^2-Q_exp[ii]^2) |
---|
266 | EndIf |
---|
267 | |
---|
268 | // DU_save[ii][jj] = DU |
---|
269 | Return DU |
---|
270 | End |
---|
271 | |
---|
272 | Function IG(ii,jj) |
---|
273 | Variable ii,jj |
---|
274 | |
---|
275 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
276 | |
---|
277 | Wave Q_exp=$(USANSFolder+":DSM:Q_exp") |
---|
278 | Variable IG,UL,UU |
---|
279 | |
---|
280 | // WAVE IG_save = root:DSM:IG_save |
---|
281 | If (ii == jj) |
---|
282 | UL=0.0 |
---|
283 | Else |
---|
284 | UL=sqrt(Q_exp[jj]^2-Q_exp[ii]^2) |
---|
285 | EndIf |
---|
286 | UU=sqrt(Q_exp[jj+1]^2-Q_exp[ii]^2) |
---|
287 | |
---|
288 | //in FORTRAN log = natural log.... |
---|
289 | IG = UU*Q_exp[jj+1]+Q_exp[ii]^2*ln(UU+Q_exp[jj+1]) |
---|
290 | IG -= UL*Q_exp[jj] |
---|
291 | IG -= Q_exp[ii]^2*ln(UL+Q_exp[jj]) |
---|
292 | IG *= 0.5 |
---|
293 | // IG_save[ii][jj] = IG |
---|
294 | |
---|
295 | Return IG |
---|
296 | End |
---|
297 | |
---|
298 | // |
---|
299 | Function FF(ii,jj) |
---|
300 | Variable ii,jj |
---|
301 | |
---|
302 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
303 | |
---|
304 | Variable FF |
---|
305 | NVAR dqv = $(USANSFolder+":DSM:gDqv") |
---|
306 | |
---|
307 | FF = (1.0/dqv)*(0.5+HH(ii,jj)) |
---|
308 | Return FF |
---|
309 | End |
---|
310 | |
---|
311 | // |
---|
312 | Function GG(ii,jj) |
---|
313 | Variable ii,jj |
---|
314 | |
---|
315 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
316 | |
---|
317 | Variable GG |
---|
318 | NVAR dqv = $(USANSFolder+":DSM:gDqv") |
---|
319 | |
---|
320 | GG = (1.0/dqv)*(0.5-HH(ii,jj)) |
---|
321 | Return GG |
---|
322 | End |
---|
323 | // |
---|
324 | Function HH(ii,jj) |
---|
325 | Variable ii,jj |
---|
326 | |
---|
327 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
328 | |
---|
329 | Wave Q_exp=$(USANSFolder+":DSM:Q_exp") |
---|
330 | Variable HH |
---|
331 | |
---|
332 | HH = 0.5*(Q_exp[jj+1]+Q_exp[jj])/(Q_exp[jj+1]-Q_exp[jj]) |
---|
333 | HH -= (1.0/(Q_exp[jj+1]-Q_exp[jj]))*(CC(ii,jj+1)-CC(ii,jj)) |
---|
334 | return HH |
---|
335 | End |
---|
336 | // |
---|
337 | // |
---|
338 | Function CC(ii,jj) |
---|
339 | Variable ii,jj |
---|
340 | |
---|
341 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
342 | |
---|
343 | wave Q_exp = $(USANSFolder+":DSM:Q_exp") |
---|
344 | Variable CC |
---|
345 | |
---|
346 | If (ii == jj) |
---|
347 | CC = 0.0 |
---|
348 | Else |
---|
349 | CC = (Q_exp[jj]-Q_exp[ii])^0.5 |
---|
350 | EndIf |
---|
351 | Return CC |
---|
352 | End |
---|
353 | |
---|
354 | // QROMO is a gerneric NR routine that takes function arguments |
---|
355 | // Call Qromo(Integrand,lower,dqv,Q_exp(N),m,SS,Midpnt) |
---|
356 | // -- here, it is always called with Integrand and Midpnt as the functions |
---|
357 | // -- so rewrite in a simpler way.... |
---|
358 | // |
---|
359 | // SS is the returned value? |
---|
360 | // H_wave, S_wave (original names H,S) |
---|
361 | // |
---|
362 | // modified using c-version in NR (pg 143) |
---|
363 | Function QROMO(A,B,qi,m) |
---|
364 | Variable A,B,qi,m |
---|
365 | |
---|
366 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
367 | |
---|
368 | Variable EPS,JMAX,JMAXP,KM,K |
---|
369 | EPS=1.E-6 |
---|
370 | JMAX=14 |
---|
371 | KM=4 |
---|
372 | K=KM+1 |
---|
373 | |
---|
374 | Make/O/D/N=(JMAX) $(USANSFolder+":DSM:S_wave") |
---|
375 | Make/O/D/N=(JMAX+1) $(USANSFolder+":DSM:H_wave") |
---|
376 | wave S_wave=$(USANSFolder+":DSM:S_wave") |
---|
377 | wave H_wave=$(USANSFolder+":DSM:H_wave") |
---|
378 | S_wave=0 |
---|
379 | H_wave=0 |
---|
380 | |
---|
381 | H_Wave[0] = 1 |
---|
382 | |
---|
383 | variable jj,SS,DSS |
---|
384 | |
---|
385 | for(jj=0;jj<jmax;jj+=1) |
---|
386 | S_wave[jj] = Midpnt(A,B,jj,qi,m) //remove FUNC, always call Integrand from within Midpnt |
---|
387 | IF (jj>=KM) //after 1st 5 points calculated |
---|
388 | POLINT(H_wave,S_wave,jj-KM,KM,0.0,SS,DSS) //ss, dss returned |
---|
389 | IF (ABS(DSS) < EPS*ABS(SS)) |
---|
390 | RETURN ss |
---|
391 | endif |
---|
392 | ENDIF |
---|
393 | // S_wave[jj+1]=S_wave[jj] |
---|
394 | H_wave[jj+1]=H_wave[jj]/9. |
---|
395 | endfor |
---|
396 | |
---|
397 | DoAlert 0,"Too many steps in QROMO" |
---|
398 | return 1 //error if you get here |
---|
399 | END |
---|
400 | |
---|
401 | // |
---|
402 | // see NR pg 109 |
---|
403 | // - Given input arrays xa[1..n] and ya[1..n], and a given value x |
---|
404 | // return a value y and an error estimate dy. if P(x) is the polynomial of |
---|
405 | // degree N-1 such that P(xai) = yai, then the returned value y=P(x) |
---|
406 | // |
---|
407 | // arrays XA[] and YA[] are passed in with an offset of the index |
---|
408 | // of where to start the interpolation |
---|
409 | Function POLINT(XA,YA,offset,N,X,Y,DY) |
---|
410 | Wave XA,YA |
---|
411 | Variable offset,N,X,&Y,&DY |
---|
412 | |
---|
413 | Variable ii,mm,nmax,ns,dif,den,ho,hp,wi,dift |
---|
414 | |
---|
415 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
416 | |
---|
417 | NMAX=10 |
---|
418 | |
---|
419 | Make/O/D/N=(NMAX) $(USANSFolder+":DSM:Ci"),$(USANSFolder+":DSM:Di") |
---|
420 | wave Ci = $(USANSFolder+":DSM:Ci") |
---|
421 | wave Di = $(USANSFolder+":DSM:Di") |
---|
422 | |
---|
423 | NS=1 |
---|
424 | DIF=ABS(X-XA[0]) |
---|
425 | for(ii=0;ii<n;ii+=1) |
---|
426 | DIFT=ABS(X-XA[ii+offset]) |
---|
427 | IF (DIFT < DIF) |
---|
428 | NS=ii |
---|
429 | DIF=DIFT |
---|
430 | ENDIF |
---|
431 | Ci[ii]=YA[ii+offset] |
---|
432 | Di[ii]=YA[ii+offset] |
---|
433 | endfor |
---|
434 | |
---|
435 | Y=YA[NS+offset] |
---|
436 | NS=NS-1 |
---|
437 | for(mm=1;mm<n-1;mm+=1) |
---|
438 | for(ii=0;ii<(n-mm);ii+=1) |
---|
439 | HO=XA[ii+offset]-X |
---|
440 | HP=XA[ii+offset+mm]-X |
---|
441 | Wi=Ci[ii+1]-Di[ii] |
---|
442 | DEN=HO-HP |
---|
443 | IF(DEN == 0.) |
---|
444 | print "den == 0 in POLINT - ERROR!!!" |
---|
445 | endif |
---|
446 | DEN=Wi/DEN |
---|
447 | Di[ii]=HP*DEN |
---|
448 | Ci[ii]=HO*DEN |
---|
449 | endfor //ii |
---|
450 | IF (2*NS < (N-mm) ) |
---|
451 | DY=Ci[NS+1] |
---|
452 | ELSE |
---|
453 | DY=Di[NS] |
---|
454 | NS=NS-1 |
---|
455 | ENDIF |
---|
456 | Y=Y+DY |
---|
457 | endfor //mm |
---|
458 | RETURN (0) //y and dy are returned as pass-by-reference |
---|
459 | END |
---|
460 | |
---|
461 | // |
---|
462 | // FUNC is always Integrand() |
---|
463 | // again, see the c-version, NR pg 142 |
---|
464 | Function MIDPNT(A,B,N,qi,m) |
---|
465 | Variable A,B,N,qi,m |
---|
466 | |
---|
467 | Variable it,tnm,del,ddel,x,summ,jj |
---|
468 | |
---|
469 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
470 | |
---|
471 | NVAR S_ret = $(USANSFolder+":DSM:gS") |
---|
472 | |
---|
473 | IF (N == 0) |
---|
474 | S_ret=(B-A)*Integrand_dsm(0.5*(A+B),qi,m) |
---|
475 | // Print "N==0, S_ret = ",s_ret |
---|
476 | return(S_ret) |
---|
477 | //IT=1 |
---|
478 | ELSE |
---|
479 | //AJJ This code is confusing! |
---|
480 | // it = 1 |
---|
481 | // for(jj=1;jj<n-1;jj+=1) |
---|
482 | // it *= 3 |
---|
483 | // endfor |
---|
484 | //AJJ Equivalent and simpler |
---|
485 | it = 3^(N-1) |
---|
486 | // |
---|
487 | TNM=IT |
---|
488 | DEL=(B-A)/(3.*TNM) |
---|
489 | DDEL=DEL+DEL |
---|
490 | X=A+0.5*DEL |
---|
491 | SUMm=0. |
---|
492 | for(jj=1;jj<=it;jj+=1) |
---|
493 | SUMm=SUMm+Integrand_dsm(X,qi,m) |
---|
494 | X=X+DDEL |
---|
495 | SUMm=SUMm+Integrand_dsm(X,qi,m) |
---|
496 | X=X+DEL |
---|
497 | endfor |
---|
498 | S_ret=(S_ret+(B-A)*SUMm/TNM)/3. |
---|
499 | IT=3*IT |
---|
500 | ENDIF |
---|
501 | return S_ret |
---|
502 | END |
---|
503 | |
---|
504 | //////////// end of Lake desmearing routines |
---|
505 | |
---|
506 | // - Almost everything below here is utility routines for data handling, panel |
---|
507 | // controls, and all the fluff that makes this useable. Note that this is |
---|
508 | // a lot more code than the actual guts of the Lake method! |
---|
509 | // (the guts of the iteration are in the DemsearButtonProc) |
---|
510 | |
---|
511 | |
---|
512 | // add three "fake" resolution columns to the output file so that it |
---|
513 | // looks like a regular SANS data set, but make the sigQ column 100x smaller than Q |
---|
514 | // so that there is effectively no resolution smearing. Set Qbar = Q, and fs = 1 |
---|
515 | // |
---|
516 | // SRK 06 FEB 06 |
---|
517 | // |
---|
518 | Function WriteUSANSDesmeared(fullpath,lo,hi,dialog) |
---|
519 | String fullpath |
---|
520 | Variable lo,hi,dialog //=1 will present dialog for name |
---|
521 | |
---|
522 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
523 | |
---|
524 | String termStr="\r\n" |
---|
525 | String destStr = USANSFolder+":DSM:" |
---|
526 | String formatStr = "%15.6g %15.6g %15.6g %15.6g %15.6g %15.6g"+termStr |
---|
527 | |
---|
528 | Variable refNum,integer,realval |
---|
529 | |
---|
530 | //*****these waves MUST EXIST, or IGOR Pro will crash, with a type 2 error**** |
---|
531 | WAVE Q_dsm =$(destStr + "Q_dsm") |
---|
532 | WAVE I_dsm=$(destStr + "I_dsm") |
---|
533 | WAVE S_dsm=$(destStr + "S_dsm") |
---|
534 | |
---|
535 | //check each wave |
---|
536 | If(!(WaveExists(Q_dsm))) |
---|
537 | Abort "Q_dsm DNExist in WriteUSANSDesmeared()" |
---|
538 | Endif |
---|
539 | If(!(WaveExists(I_dsm))) |
---|
540 | Abort "I_dsm DNExist in WriteUSANSDesmeared()" |
---|
541 | Endif |
---|
542 | If(!(WaveExists(S_dsm))) |
---|
543 | Abort "S_dsm DNExist in WriteUSANSDesmeared()" |
---|
544 | Endif |
---|
545 | |
---|
546 | // 06 FEB 06 SRK |
---|
547 | // make dummy waves to hold the "fake" resolution, and write it as the last 3 columns |
---|
548 | // |
---|
549 | Duplicate/O Q_dsm,res1,res2,res3 |
---|
550 | res3 = 1 // "fake" beamstop shadowing |
---|
551 | res1 /= 100 //make the sigmaQ so small that there is no smearing |
---|
552 | |
---|
553 | if(dialog) |
---|
554 | Open/D refnum as fullpath+".dsm" //won't actually open the file |
---|
555 | If(cmpstr(S_filename,"")==0) |
---|
556 | //user cancel, don't write out a file |
---|
557 | Close/A |
---|
558 | Abort "no data file was written" |
---|
559 | Endif |
---|
560 | fullpath = S_filename |
---|
561 | Endif |
---|
562 | |
---|
563 | //write out partial set? |
---|
564 | Duplicate/O Q_dsm,tq,ti,te |
---|
565 | ti=I_dsm |
---|
566 | te=S_dsm |
---|
567 | if( (lo!=hi) && (lo<hi)) |
---|
568 | redimension/N=(hi-lo+1) tq,ti,te,res1,res2,res3 //lo to hi, inclusive |
---|
569 | tq=Q_dsm[p+lo] |
---|
570 | ti=I_dsm[p+lo] |
---|
571 | te=S_dsm[p+lo] |
---|
572 | endif |
---|
573 | |
---|
574 | //tailor the output given the type of data written out... |
---|
575 | String samStr="",dateStr="",str1,str2 |
---|
576 | |
---|
577 | NVAR m = $(USANSFolder+":DSM:gPowerM") // power law exponent |
---|
578 | NVAR chiFinal = $(USANSFolder+":DSM:gChi2Final") //chi^2 final |
---|
579 | NVAR iter = $(USANSFolder+":DSM:gIterations") //total number of iterations |
---|
580 | |
---|
581 | //get the number of spline passes from the wave note |
---|
582 | String noteStr |
---|
583 | Variable boxPass,SplinePass |
---|
584 | noteStr=note(I_dsm) |
---|
585 | BoxPass = NumberByKey("BOX", noteStr, "=", ";") |
---|
586 | splinePass = NumberByKey("SPLINE", noteStr, "=", ";") |
---|
587 | |
---|
588 | samStr = fullpath |
---|
589 | dateStr="CREATED: "+date()+" at "+time() |
---|
590 | sprintf str1,"Chi^2 = %g PowerLaw m = %4.2f Iterations = %d",chiFinal,m,iter |
---|
591 | sprintf str2,"%d box smooth passes and %d smoothing spline passes",boxPass,splinePass |
---|
592 | |
---|
593 | |
---|
594 | //actually open the file |
---|
595 | Open refNum as fullpath |
---|
596 | |
---|
597 | fprintf refnum,"%s"+termStr,samStr |
---|
598 | fprintf refnum,"%s"+termStr,str1 |
---|
599 | fprintf refnum,"%s"+termStr,str2 |
---|
600 | fprintf refnum,"%s"+termStr,dateStr |
---|
601 | |
---|
602 | wfprintf refnum, formatStr, tq,ti,te,res1,res2,res3 |
---|
603 | |
---|
604 | Close refnum |
---|
605 | |
---|
606 | Killwaves/Z ti,tq,te,res1,res2,res3 |
---|
607 | |
---|
608 | Return(0) |
---|
609 | End |
---|
610 | |
---|
611 | /// procedures to do the extrapolation to high Q |
---|
612 | // very similar to the procedures in the Invariant package |
---|
613 | // |
---|
614 | //create the wave to extrapolate |
---|
615 | // w is the input q-values |
---|
616 | Function DSM_SetExtrWaves(w) |
---|
617 | Wave w |
---|
618 | |
---|
619 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
620 | |
---|
621 | Variable num_extr=25 |
---|
622 | |
---|
623 | SetDataFolder $(USANSFolder+":DSM") |
---|
624 | |
---|
625 | Make/O/D/N=(num_extr) extr_hqq,extr_hqi |
---|
626 | extr_hqi=1 //default values |
---|
627 | |
---|
628 | //set the q-range |
---|
629 | Variable qmax,num |
---|
630 | // qmax=0.03 |
---|
631 | |
---|
632 | num=numpnts(w) |
---|
633 | qmax=6*w[num-1] |
---|
634 | |
---|
635 | extr_hqq = w[num-1] + x * (qmax-w[num-1])/num_extr |
---|
636 | |
---|
637 | SetDataFolder root: |
---|
638 | return(0) |
---|
639 | End |
---|
640 | |
---|
641 | //creates I_ext,Q_ext,S_ext with extended q-values |
---|
642 | // |
---|
643 | // if num_extr == 0 , the input waves are returned as _ext |
---|
644 | // |
---|
645 | // !! uses simple linear extrapolation at low Q - just need something |
---|
646 | // reasonable in the waves to keep from getting smoothing artifacts |
---|
647 | // - uses power law at high q |
---|
648 | Function ExtendToSmooth(qw,iw,sw,nbeg,nend,num_extr) |
---|
649 | Wave qw,iw,sw |
---|
650 | Variable nbeg,nend,num_extr |
---|
651 | |
---|
652 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
653 | |
---|
654 | Variable/G V_FitMaxIters=300 |
---|
655 | Variable/G V_fitOptions=4 //suppress the iteration window |
---|
656 | Variable num_new,qmax,num,qmin |
---|
657 | num=numpnts(qw) |
---|
658 | |
---|
659 | num_new = num + 2*num_extr |
---|
660 | |
---|
661 | // Print "num,num_new",num,num_new |
---|
662 | |
---|
663 | SetDataFolder $(USANSFolder+":DSM") |
---|
664 | Make/O/D/N=(num_new) Q_ext,I_ext,S_ext |
---|
665 | |
---|
666 | if(num_extr == 0) //the extended waves are the input, get out |
---|
667 | Q_ext = qw |
---|
668 | I_ext = iw |
---|
669 | S_ext = sw |
---|
670 | setDatafolder root: |
---|
671 | return (0) |
---|
672 | endif |
---|
673 | |
---|
674 | //make the extensions |
---|
675 | Make/O/D/N=(num_extr) hqq,hqi,lqq,lqi |
---|
676 | hqi=1 //default values |
---|
677 | lqi=0 |
---|
678 | //set the q-range based on the high/low values of the data set |
---|
679 | // qmin=1e-6 |
---|
680 | qmin= qw[0]*0.5 |
---|
681 | qmax = qw[num-1]*2 |
---|
682 | |
---|
683 | hqq = qw[num-1] + x * (qmax-qw[num-1])/num_extr |
---|
684 | lqq = qmin + x * (qw[0]-qmin)/num_extr |
---|
685 | |
---|
686 | //do the fits |
---|
687 | Duplicate/O iw dummy //use this as the destination to suppress fit_ from being appended |
---|
688 | |
---|
689 | //Use simple linear fits line: y = K0+K1*x |
---|
690 | CurveFit/Q line iw[0,(nbeg-1)] /X=qw /W=sw /D=dummy |
---|
691 | Wave W_coef=W_coef |
---|
692 | lqi=W_coef[0]+W_coef[1]*lqq |
---|
693 | // |
---|
694 | // Guinier or Power-law fits |
---|
695 | // |
---|
696 | // Make/O/D G_coef={100,-100} //input |
---|
697 | // FuncFit DSM_Guinier_Fit G_coef iw[0,(nbeg-1)] /X=qw /W=sw /D |
---|
698 | // lqi= DSM_Guinier_Fit(G_coef,lqq) |
---|
699 | |
---|
700 | // Printf "I(q=0) = %g (1/cm)\r",G_coef[0] |
---|
701 | // Printf "Rg = %g (A)\r",sqrt(-3*G_coef[1]) |
---|
702 | |
---|
703 | Make/O/D P_coef={0,1,-4} //input --- (set background to zero and hold fixed) |
---|
704 | CurveFit/Q/N/H="100" Power kwCWave=P_coef iw[(num-1-nend),(num-1)] /X=qw /W=sw /D=dummy |
---|
705 | hqi=P_coef[0]+P_coef[1]*hqq^P_coef[2] |
---|
706 | |
---|
707 | // Printf "Power law exponent = %g\r",P_coef[2] |
---|
708 | // Printf "Pre-exponential = %g\r",P_coef[1] |
---|
709 | |
---|
710 | // concatenate the extensions |
---|
711 | Q_ext[0,(num_extr-1)] = lqq[p] |
---|
712 | Q_ext[num_extr,(num_extr+num-1)] = qw[p-num_extr] |
---|
713 | Q_ext[(num_extr+num),(2*num_extr+num-1)] = hqq[p-num_extr-num] |
---|
714 | I_ext[0,(num_extr-1)] = lqi[p] |
---|
715 | I_ext[num_extr,(num_extr+num-1)] = iw[p-num_extr] |
---|
716 | I_ext[(num_extr+num),(2*num_extr+num-1)] = hqi[p-num_extr-num] |
---|
717 | S_ext[0,(num_extr-1)] = sw[0] |
---|
718 | S_ext[num_extr,(num_extr+num-1)] = sw[p-num_extr] |
---|
719 | S_ext[(num_extr+num),(2*num_extr+num-1)] = sw[num-1] |
---|
720 | |
---|
721 | killwaves/z dummy |
---|
722 | SetDataFolder root: |
---|
723 | return(0) |
---|
724 | End |
---|
725 | |
---|
726 | // pass in the (smeared) experimental data and find the power-law extrapolation |
---|
727 | // 10 points is usually good enough, unless the data is crummy |
---|
728 | // returns the prediction for the exponent |
---|
729 | // |
---|
730 | Function DSM_DoExtrapolate(qw,iw,sw,nend) |
---|
731 | Wave qw,iw,sw |
---|
732 | Variable nend |
---|
733 | |
---|
734 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
735 | |
---|
736 | Setdatafolder $(USANSFolder+":DSM") |
---|
737 | |
---|
738 | // Wave extr_lqi=extr_lqi |
---|
739 | // Wave extr_lqq=extr_lqq |
---|
740 | Wave extr_hqi=extr_hqi |
---|
741 | Wave extr_hqq=extr_hqq |
---|
742 | Variable/G V_FitMaxIters=300 |
---|
743 | Variable/G V_fitOptions=4 //suppress the iteration window |
---|
744 | Variable num=numpnts(iw),retVal |
---|
745 | |
---|
746 | Make/O/D P_coef={0,1,-4} //input |
---|
747 | // Make/O/T Constr={"K2<0","K2 > -20"} |
---|
748 | //(set background to zero and hold fixed) |
---|
749 | CurveFit/H="100" Power kwCWave=P_coef iw[(num-1-nend),(num-1)] /X=qw /W=sw /D |
---|
750 | extr_hqi=P_coef[0]+P_coef[1]*extr_hqq^P_coef[2] |
---|
751 | |
---|
752 | // for the case of data with a background |
---|
753 | // Make/O/D P_coef={0,1,-4} //input |
---|
754 | // //(set background to iw[num-1], let it be a free parameter |
---|
755 | // P_coef[0] = iw[num-1] |
---|
756 | // CurveFit Power kwCWave=P_coef iw[(num-1-nend),(num-1)] /X=qw /W=sw /D |
---|
757 | // extr_hqi=P_coef[0]+P_coef[1]*extr_hqq^P_coef[2] |
---|
758 | |
---|
759 | |
---|
760 | Printf "Smeared Power law exponent = %g\r",P_coef[2] |
---|
761 | Printf "**** For Desmearing, use a Power law exponent of %5.1f\r",P_coef[2]-1 |
---|
762 | |
---|
763 | retVal = P_coef[2]-1 |
---|
764 | SetDataFolder root: |
---|
765 | return(retVal) |
---|
766 | End |
---|
767 | |
---|
768 | Function DSM_Guinier_Fit(w,x) //: FitFunc |
---|
769 | Wave w |
---|
770 | Variable x |
---|
771 | |
---|
772 | //fit data to I(q) = A*exp(B*q^2) |
---|
773 | // (B will be negative) |
---|
774 | //two parameters |
---|
775 | Variable a,b,ans |
---|
776 | a=w[0] |
---|
777 | b=w[1] |
---|
778 | ans = a*exp(b*x*x) |
---|
779 | return(ans) |
---|
780 | End |
---|
781 | |
---|
782 | Proc Desmear_Graph() |
---|
783 | String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
784 | |
---|
785 | PauseUpdate; Silent 1 // building window... |
---|
786 | Display /W=(5,44,408,558) /K=1 |
---|
787 | ModifyGraph cbRGB=(51664,44236,58982) |
---|
788 | DoWindow/C Desmear_Graph |
---|
789 | ControlBar 160 |
---|
790 | // break into tabs |
---|
791 | TabControl DSM_Tab,pos={5,3},size={392,128},proc=DSM_TabProc |
---|
792 | TabControl DSM_Tab,labelBack=(49151,49152,65535),tabLabel(0)="Load" |
---|
793 | TabControl DSM_Tab,tabLabel(1)="Mask",tabLabel(2)="Extrapolate" |
---|
794 | TabControl DSM_Tab,tabLabel(3)="Smooth",tabLabel(4)="Desmear",value= 0 |
---|
795 | |
---|
796 | //always visible - revert and save |
---|
797 | //maybe the wrong place here? |
---|
798 | Button DSMControlA,pos={225,135},size={80,20},proc=DSM_RevertButtonProc,title="Revert" |
---|
799 | Button DSMControlA,help={"Revert the smeared data to its original state and start over"} |
---|
800 | Button DSMControlB,pos={325,135},size={50,20},proc=DSM_SaveButtonProc,title="Save" |
---|
801 | Button DSMControlB,help={"Save the desmeared data set"} |
---|
802 | Button DSMControlC,pos={25,135},size={50,20},proc=DSM_HelpButtonProc,title="Help" |
---|
803 | Button DSMControlC,help={"Show the help file for desmearing"} |
---|
804 | |
---|
805 | // add the controls to each tab ---- all names start with "DSMControl_" |
---|
806 | |
---|
807 | //tab(0) Load - initially visible |
---|
808 | Button DSMControl_0a,pos={23,39},size={80,20},proc=DSM_LoadButtonProc,title="Load Data" |
---|
809 | Button DSMControl_0a,help={"Load slit-smeared USANS data = \".cor\" files"} |
---|
810 | CheckBox DSMControl_0b,pos={26,74},size={80,14},proc=DSM_LoadCheckProc,title="Log Axes?" |
---|
811 | CheckBox DSMControl_0b,help={"Toggle Log/Lin Q display"},value= 1 |
---|
812 | TitleBox DSMControl_0c,pos={120,37},size={104,19},font="Courier",fSize=10 |
---|
813 | TitleBox DSMControl_0c,variable= $(USANSFolder+":DSM:gStr1") |
---|
814 | //second message string not used currently |
---|
815 | // TitleBox DSMControl_0d,pos={120,57},size={104,19},font="Courier",fSize=10 |
---|
816 | // TitleBox DSMControl_0d,variable= root:DSM:gStr2 |
---|
817 | |
---|
818 | //tab(1) Mask |
---|
819 | Button DSMControl_1a,pos={20,35},size={90,20},proc=DSM_MyMaskProc,title="Mask Point" //bMask |
---|
820 | Button DSMControl_1a,help={"Toggles the masking of the selected data point"} |
---|
821 | Button DSMControl_1a,disable=1 |
---|
822 | Button DSMControl_1b,pos={20,65},size={140,20},proc=DSM_MaskGTCursor,title="Mask Q >= Cursor" //bMask |
---|
823 | Button DSMControl_1b,help={"Toggles the masking of all q-values GREATER than the current cursor location"} |
---|
824 | Button DSMControl_1b,disable=1 |
---|
825 | Button DSMControl_1c,pos={20,95},size={140,20},proc=DSM_MaskLTCursor,title="Mask Q <= Cursor" //bMask |
---|
826 | Button DSMControl_1c,help={"Toggles the masking of all q-values LESS than the current cursor location"} |
---|
827 | Button DSMControl_1c,disable=1 |
---|
828 | Button DSMControl_1d,pos={180,35},size={90,20},proc=DSM_ClearMaskProc,title="Clear Mask" //bMask |
---|
829 | Button DSMControl_1d,help={"Clears all mask points"} |
---|
830 | Button DSMControl_1d,disable=1 |
---|
831 | // Button DSMControl_1b,pos={144,66},size={110,20},proc=DSM_MaskDoneButton,title="Done Masking" |
---|
832 | // Button DSMControl_1b,disable=1 |
---|
833 | |
---|
834 | //tab(2) Extrapolate |
---|
835 | Button DSMControl_2a,pos={31,42},size={90,20},proc=DSM_ExtrapolateButtonProc,title="Extrapolate" |
---|
836 | Button DSMControl_2a,help={"Extrapolate the high-q region with a power-law"} |
---|
837 | Button DSMControl_2a,disable=1 |
---|
838 | SetVariable DSMControl_2b,pos={31,70},size={100,15},title="# of points" |
---|
839 | SetVariable DSMControl_2b,help={"Set the number of points for the power-law extrapolation"} |
---|
840 | SetVariable DSMControl_2b,limits={5,100,1},value= $(USANSFolder+":DSM:gNptsExtrap") |
---|
841 | SetVariable DSMControl_2b,disable=1 |
---|
842 | CheckBox DSMControl_2c,pos={157,45},size={105,14},proc=DSM_ExtrapolationCheckProc,title="Show Extrapolation" |
---|
843 | CheckBox DSMControl_2c,help={"Show or hide the high q extrapolation"},value= 1 |
---|
844 | CheckBox DSMControl_2c,disable=1 |
---|
845 | SetVariable DSMControl_2d,pos={31,96},size={150,15},title="Power Law Exponent" |
---|
846 | SetVariable DSMControl_2d,help={"Power Law exponent from the fit = the DESMEARED slope - override as needed"} |
---|
847 | SetVariable DSMControl_2d format="%5.2f" |
---|
848 | SetVariable DSMControl_2d,limits={-inf,inf,0},value= $(USANSFolder+":DSM:gPowerM") |
---|
849 | SetVariable DSMControl_2d,disable=1 |
---|
850 | |
---|
851 | //tab(3) Smooth |
---|
852 | Button DSMControl_3a,pos={34,97},size={70,20},proc=DSM_SmoothButtonProc,title="Smooth" |
---|
853 | Button DSMControl_3a,disable=1 |
---|
854 | //BoxCheck |
---|
855 | CheckBox DSMControl_3b,pos={34,39},size={35,14},title="Box",value= 1 |
---|
856 | CheckBox DSMControl_3b,help={"Use a single pass of 3-point box smoothing"} |
---|
857 | CheckBox DSMControl_3b,disable=1 |
---|
858 | //SSCheck |
---|
859 | CheckBox DSMControl_3c,pos={34,60},size={45,14},title="Spline",value= 0 |
---|
860 | CheckBox DSMControl_3c,help={"Use a single pass of a smoothing spline"} |
---|
861 | CheckBox DSMControl_3c,disable=1 |
---|
862 | //extendCheck |
---|
863 | CheckBox DSMControl_3d,pos={268,60},size={71,14},title="Extend Data" |
---|
864 | CheckBox DSMControl_3d,help={"extends the data at both low q and high q to avoid end effects in smoothing"} |
---|
865 | CheckBox DSMControl_3d,value= 0 |
---|
866 | CheckBox DSMControl_3d,disable=1 |
---|
867 | Button DSMControl_3e,pos={125,97},size={90,20},proc=DSM_SmoothUndoButtonProc,title="Start Over" |
---|
868 | Button DSMControl_3e,help={"Start the smoothing over again without needing to re-mask the data set"} |
---|
869 | Button DSMControl_3e,disable=1 |
---|
870 | SetVariable DSMControl_3f,pos={94,60},size={150,15},title="Smoothing factor" |
---|
871 | SetVariable DSMControl_3f,help={"Smoothing factor for the smoothing spline"} |
---|
872 | SetVariable DSMControl_3f format="%5.4f" |
---|
873 | SetVariable DSMControl_3f,limits={0.01,2,0.01},value= $(USANSFolder+":DSM:gSmoothFac") |
---|
874 | SetVariable DSMControl_3f,disable=1 |
---|
875 | CheckBox DSMControl_3g,pos={268,39},size={90,14},title="Log-scale smoothing?" |
---|
876 | CheckBox DSMControl_3g,help={"Use log-scaled intensity during smoothing (reverts to linear if negative intensity points found)"} |
---|
877 | CheckBox DSMControl_3g,value=0 |
---|
878 | CheckBox DSMControl_3g,disable=1 |
---|
879 | ValDisplay DSMControl_3h pos={235,97},title="Chi^2",size={80,20},value=root:Packages:NIST:USANS:DSM:gChi2Smooth |
---|
880 | ValDisplay DSMControl_3h,help={"This is the Chi^2 value for the smoothed data vs experimental data"} |
---|
881 | ValDisplay DSMControl_3h,disable=1 |
---|
882 | |
---|
883 | //tab(4) Desmear |
---|
884 | Button DSMControl_4a,pos={35,93},size={80,20},proc=DSM_DesmearButtonProc,title="Desmear" |
---|
885 | Button DSMControl_4a,help={"Do the desmearing - the result is in I_dsm"} |
---|
886 | Button DSMControl_4a,disable=1 |
---|
887 | SetVariable DSMControl_4b,pos={35,63},size={120,15},title="Chi^2 target" |
---|
888 | SetVariable DSMControl_4b,help={"Set the targetchi^2 for convergence (recommend chi^2=1)"} |
---|
889 | SetVariable DSMControl_4b,limits={0,inf,0.1},value= $(USANSFolder+":DSM:gChi2Target") |
---|
890 | SetVariable DSMControl_4b,disable=1 |
---|
891 | SetVariable DSMControl_4c,pos={35,35},size={80,15},title="dQv" |
---|
892 | SetVariable DSMControl_4c,help={"Slit height as read in from the data file. 0.117 is the NIST value, override if necessary"} |
---|
893 | SetVariable DSMControl_4c,limits={-inf,inf,0},value= $(USANSFolder+":DSM:gDqv") |
---|
894 | SetVariable DSMControl_4c,disable=1 |
---|
895 | TitleBox DSMControl_4d,pos={160,37},size={104,19},font="Courier",fSize=10 |
---|
896 | TitleBox DSMControl_4d,variable= $(USANSFolder+":DSM:gIterStr") |
---|
897 | TitleBox DSMControl_4d,disable=1 |
---|
898 | |
---|
899 | |
---|
900 | SetDataFolder root: |
---|
901 | EndMacro |
---|
902 | |
---|
903 | // function to control the drawing of buttons in the TabControl on the main panel |
---|
904 | // Naming scheme for the buttons MUST be strictly adhered to... else buttons will |
---|
905 | // appear in odd places... |
---|
906 | // all buttons are named DSMControl_NA where N is the tab number and A is the letter denoting |
---|
907 | // the button's position on that particular tab. |
---|
908 | // in this way, buttons will always be drawn correctly :-) |
---|
909 | // |
---|
910 | Function DSM_TabProc(ctrlName,tab) //: TabControl |
---|
911 | String ctrlName |
---|
912 | Variable tab |
---|
913 | |
---|
914 | String ctrlList = ControlNameList("",";"),item="",nameStr="" |
---|
915 | Variable num = ItemsinList(ctrlList,";"),ii,onTab |
---|
916 | for(ii=0;ii<num;ii+=1) |
---|
917 | //items all start w/"DSMControl_" //11 characters |
---|
918 | item=StringFromList(ii, ctrlList ,";") |
---|
919 | nameStr=item[0,10] |
---|
920 | if(cmpstr(nameStr,"DSMControl_")==0) |
---|
921 | onTab = str2num(item[11]) //12th is a number |
---|
922 | ControlInfo $item |
---|
923 | switch(abs(V_flag)) |
---|
924 | case 1: |
---|
925 | Button $item,disable=(tab!=onTab) |
---|
926 | break |
---|
927 | case 2: |
---|
928 | CheckBox $item,disable=(tab!=onTab) |
---|
929 | break |
---|
930 | case 5: |
---|
931 | SetVariable $item,disable=(tab!=onTab) |
---|
932 | break |
---|
933 | case 10: |
---|
934 | TitleBox $item,disable=(tab!=onTab) |
---|
935 | break |
---|
936 | case 4: |
---|
937 | ValDisplay $item,disable=(tab!=onTab) |
---|
938 | break |
---|
939 | // add more items to the switch if different control types are used |
---|
940 | endswitch |
---|
941 | |
---|
942 | endif |
---|
943 | endfor |
---|
944 | |
---|
945 | if(tab==1) |
---|
946 | DSM_MyMaskProc("") //start maksing if you click on the tab |
---|
947 | else |
---|
948 | DSM_MaskDoneButton("") //masking is done if you click off the tab |
---|
949 | endif |
---|
950 | |
---|
951 | return 0 |
---|
952 | End |
---|
953 | |
---|
954 | Proc AppendSmeared() |
---|
955 | String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
956 | |
---|
957 | SetDataFolder $(USANSFolder+":DSM") |
---|
958 | // if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_exp_orig",0,2) == -1) //Igor 5 |
---|
959 | if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_exp_orig",0) == -1) |
---|
960 | AppendToGraph/W=Desmear_Graph I_exp_orig vs Q_exp_orig |
---|
961 | ModifyGraph mode=4,marker=19 //3 is markers, 4 is markers and lines |
---|
962 | ModifyGraph rgb(I_exp_orig)=(0,0,0) |
---|
963 | ModifyGraph msize=2,grid=1,log=1,mirror=2,standoff=0,tickunit=1 |
---|
964 | ErrorBars/T=0 I_exp_orig Y,wave=(S_exp_orig,S_exp_orig) |
---|
965 | Legend/N=text0/J "\\F'Courier'\\s(I_exp_orig) I_exp_orig" |
---|
966 | Label left "Intensity" |
---|
967 | Label bottom "Q (1/A)" |
---|
968 | endif |
---|
969 | //always update the textbox - kill the old one first |
---|
970 | TextBox/K/N=text1 |
---|
971 | // TextBox/C/N=text1/F=0/A=MT/E=2/X=5.50/Y=0.00 root:DSM:gCurFile //Igor 5 |
---|
972 | TextBox/C/N=text1/F=0/A=MT/E=1/X=5.50/Y=0.00 $(USANSFolder+":DSM:gCurFile") |
---|
973 | End |
---|
974 | |
---|
975 | Proc AppendMask() |
---|
976 | String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
977 | |
---|
978 | // if( strsearch(TraceNameList("Desmear_Graph", "", 1),"MaskData",0,2) == -1) //Igor 5 |
---|
979 | if( strsearch(TraceNameList("Desmear_Graph", "", 1),"MaskData",0) == -1) |
---|
980 | SetDataFolder $(USANSFolder+":DSM:") |
---|
981 | AppendToGraph/W=Desmear_Graph MaskData vs Q_exp_orig |
---|
982 | ModifyGraph mode(MaskData)=3,marker(MaskData)=8,msize(MaskData)=2.5,opaque(MaskData)=1 |
---|
983 | ModifyGraph rgb(MaskData)=(65535,16385,16385) |
---|
984 | setdatafolder root: |
---|
985 | endif |
---|
986 | end |
---|
987 | |
---|
988 | Proc AppendSmoothed() |
---|
989 | String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
990 | |
---|
991 | // if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_smth",0,2) == -1) //Igor 5 |
---|
992 | if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_smth",0) == -1) |
---|
993 | SetDataFolder $(USANSFolder+":DSM:") |
---|
994 | AppendToGraph/W=Desmear_Graph I_smth vs Q_smth |
---|
995 | ModifyGraph/W=Desmear_Graph rgb(I_smth)=(3,52428,1),lsize(I_smth)=2 |
---|
996 | setdatafolder root: |
---|
997 | endif |
---|
998 | end |
---|
999 | |
---|
1000 | Function RemoveSmoothed() |
---|
1001 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1002 | |
---|
1003 | SetDataFolder $(USANSFolder+":DSM:") |
---|
1004 | RemoveFromGraph/W=Desmear_Graph/Z I_smth |
---|
1005 | setdatafolder root: |
---|
1006 | end |
---|
1007 | |
---|
1008 | Function RemoveMask() |
---|
1009 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1010 | |
---|
1011 | SetDataFolder $(USANSFolder+":DSM:") |
---|
1012 | RemoveFromGraph/W=Desmear_Graph/Z MaskData |
---|
1013 | setdatafolder root: |
---|
1014 | end |
---|
1015 | |
---|
1016 | Proc AppendDesmeared() |
---|
1017 | String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1018 | |
---|
1019 | // if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_dsm",0,2) == -1) //Igor 5 |
---|
1020 | if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_dsm",0) == -1) |
---|
1021 | SetDataFolder $(USANSFolder+":DSM:") |
---|
1022 | AppendToGraph/W=Desmear_Graph I_dsm vs Q_dsm |
---|
1023 | ModifyGraph mode(I_dsm)=3,marker(I_dsm)=19 |
---|
1024 | ModifyGraph rgb(I_dsm)=(1,16019,65535),msize(I_dsm)=2 |
---|
1025 | ErrorBars/T=0 I_dsm Y,wave=(S_dsm,S_dsm) |
---|
1026 | setdatafolder root: |
---|
1027 | endif |
---|
1028 | end |
---|
1029 | |
---|
1030 | Function RemoveDesmeared() |
---|
1031 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1032 | |
---|
1033 | SetDataFolder $(USANSFolder+":DSM:") |
---|
1034 | RemoveFromGraph/W=Desmear_Graph/Z I_dsm |
---|
1035 | setdatafolder root: |
---|
1036 | end |
---|
1037 | |
---|
1038 | Function AppendExtrapolation() |
---|
1039 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1040 | |
---|
1041 | // if( strsearch(TraceNameList("Desmear_Graph", "", 1),"extr_hqi",0,2) == -1) //Igor 5 |
---|
1042 | if( strsearch(TraceNameList("Desmear_Graph", "", 1),"extr_hqi",0) == -1) |
---|
1043 | SetDataFolder $(USANSFolder+":DSM:") |
---|
1044 | AppendToGraph/W=Desmear_Graph extr_hqi vs extr_hqq |
---|
1045 | ModifyGraph/W=Desmear_Graph lSize(extr_hqi)=2 |
---|
1046 | setdatafolder root: |
---|
1047 | endif |
---|
1048 | end |
---|
1049 | |
---|
1050 | Function RemoveExtrapolation() |
---|
1051 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1052 | |
---|
1053 | SetDataFolder $(USANSFolder+":DSM:") |
---|
1054 | RemoveFromGraph/W=Desmear_Graph/Z extr_hqi |
---|
1055 | setdatafolder root: |
---|
1056 | end |
---|
1057 | |
---|
1058 | // step (1) - read in the data, and plot it |
---|
1059 | // clear out all of the "old" waves, remove them from the graph first |
---|
1060 | // reads in a fresh copy of the data |
---|
1061 | // |
---|
1062 | // produces Q_exp, I_exp, S_exp waves (and originals "_orig") |
---|
1063 | // add a dummy wave note that can be changed on later steps |
---|
1064 | // |
---|
1065 | Function DSM_LoadButtonProc(ctrlName) : ButtonControl |
---|
1066 | String ctrlName |
---|
1067 | |
---|
1068 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1069 | |
---|
1070 | String qStr,iStr,sStr,sqStr |
---|
1071 | Variable nq,dqv,numBad,val |
---|
1072 | |
---|
1073 | // remove any of the old traces on the graph and delete the waves |
---|
1074 | CleanUpJunk() |
---|
1075 | |
---|
1076 | SetDataFolder root: |
---|
1077 | Execute "A_LoadOneDDataWithName(\"\",0)" |
---|
1078 | SVAR fileStr = root:Packages:NIST:gLastFileName |
---|
1079 | if (cmpstr(fileStr,"") == 0) |
---|
1080 | return(0) //get out if no file selected |
---|
1081 | endif |
---|
1082 | //define the waves that the smoothing will be looking for... |
---|
1083 | SVAR fname = $("root:Packages:NIST:gLastFileName") //this changes as any data is loaded |
---|
1084 | SVAR curFile = $(USANSFolder+":DSM:gCurFile") //keep this for title, save |
---|
1085 | curFile = fname |
---|
1086 | |
---|
1087 | qStr = CleanupName((fName + "_q"),0) //the q-wave |
---|
1088 | iStr = CleanupName((fName + "_i"),0) //the i-wave |
---|
1089 | sStr = CleanupName((fName + "_s"),0) //the s-wave |
---|
1090 | sqStr = CleanupName((fName + "sq"),0) //the sq-wave, which should have -dQv as the elements |
---|
1091 | |
---|
1092 | String DFStr= CleanupName(fname,0) |
---|
1093 | |
---|
1094 | Duplicate/O $("root:"+DFStr+":"+qStr) $(USANSFolder+":DSM:Q_exp ") |
---|
1095 | Duplicate/O $("root:"+DFStr+":"+iStr) $(USANSFolder+":DSM:I_exp") |
---|
1096 | Duplicate/O $("root:"+DFStr+":"+sStr) $(USANSFolder+":DSM:S_exp ") |
---|
1097 | wave Q_exp = $(USANSFolder+":DSM:Q_exp") |
---|
1098 | Wave I_exp = $(USANSFolder+":DSM:I_exp") |
---|
1099 | Wave S_exp = $(USANSFolder+":DSM:S_exp") |
---|
1100 | |
---|
1101 | // remove any negative q-values (and q=0 values!)(and report this) |
---|
1102 | // ? and trim the low q to be >= 3.0e-5 (1/A), below this USANS is not reliable. |
---|
1103 | NumBad = RemoveBadQPoints(Q_exp,I_exp,S_exp,0) |
---|
1104 | SVAR str1 = $(USANSFolder+":DSM:gStr1") |
---|
1105 | sprintf str1,"%d negative q-values were removed",numBad |
---|
1106 | |
---|
1107 | // don't trim off any positive q-values |
---|
1108 | // val = 3e-5 //lowest "good" q-value from USANS |
---|
1109 | // NumBad = RemoveBadQPoints(Q_exp,I_exp,S_exp,val-1e-8) |
---|
1110 | // SVAR str2 = root:DSM:gStr2 |
---|
1111 | // sprintf str2,"%d q-values below q = %g were removed",numBad,val |
---|
1112 | |
---|
1113 | Duplicate/O $(USANSFolder+":DSM:Q_exp") $(USANSFolder+":DSM:Q_exp_orig") |
---|
1114 | Duplicate/O $(USANSFolder+":DSM:I_exp") $(USANSFolder+":DSM:I_exp_orig") |
---|
1115 | Duplicate/O $(USANSFolder+":DSM:S_exp") $(USANSFolder+":DSM:S_exp_orig") |
---|
1116 | wave I_exp_orig = $(USANSFolder+":DSM:I_exp_orig") |
---|
1117 | |
---|
1118 | nq = numpnts($(USANSFolder+":DSM:Q_exp")) |
---|
1119 | |
---|
1120 | dQv = NumVarOrDefault("root:"+DFStr+":USANS_dQv", 0.117 ) |
---|
1121 | // if(WaveExists(sigQ)) //try to read dQv |
---|
1122 | //// dqv = -sigQ[0][0] |
---|
1123 | //// DoAlert 0,"Found dQv value of " + num2str(dqv) |
---|
1124 | // else |
---|
1125 | // dqv = 0.117 |
---|
1126 | // // dqv = 0.037 //old value |
---|
1127 | // DoAlert 0,"Could not find dQv in the data file - using " + num2str(dqv) |
---|
1128 | // endif |
---|
1129 | NVAR gDqv = $(USANSFolder+":DSM:gDqv") //needs to be global for Weights_L() |
---|
1130 | NVAR gNq = $(USANSFolder+":DSM:gNq") |
---|
1131 | //reset the globals |
---|
1132 | gDqv = dqv |
---|
1133 | gNq = nq |
---|
1134 | |
---|
1135 | // append the (blank) wave note to the intensity wave |
---|
1136 | Note I_exp,"BOX=0;SPLINE=0;" |
---|
1137 | Note I_exp_orig,"BOX=0;SPLINE=0;" |
---|
1138 | |
---|
1139 | //draw the graph |
---|
1140 | Execute "AppendSmeared()" |
---|
1141 | |
---|
1142 | SetDataFolder root: |
---|
1143 | End |
---|
1144 | |
---|
1145 | // remove any q-values <= val |
---|
1146 | Function RemoveBadQPoints(qw,iw,sw,val) |
---|
1147 | Wave qw,iw,sw |
---|
1148 | Variable val |
---|
1149 | |
---|
1150 | Variable ii,num,numBad,qval |
---|
1151 | num = numpnts(qw) |
---|
1152 | |
---|
1153 | ii=0 |
---|
1154 | numBad=0 |
---|
1155 | do |
---|
1156 | qval = qw[ii] |
---|
1157 | if(qval <= val) |
---|
1158 | numBad += 1 |
---|
1159 | else //keep the points |
---|
1160 | qw[ii-numBad] = qval |
---|
1161 | iw[ii-numBad] = iw[ii] |
---|
1162 | sw[ii-numBad] = sw[ii] |
---|
1163 | endif |
---|
1164 | ii += 1 |
---|
1165 | while(ii<num) |
---|
1166 | //trim the end of the waves |
---|
1167 | DeletePoints num-numBad, numBad, qw,iw,sw |
---|
1168 | return(numBad) |
---|
1169 | end |
---|
1170 | |
---|
1171 | // if mw = Nan, keep the point, if a numerical value, delete it |
---|
1172 | Function RemoveMaskedPoints(mw,qw,iw,sw) |
---|
1173 | Wave mw,qw,iw,sw |
---|
1174 | |
---|
1175 | Variable ii,num,numBad,mask |
---|
1176 | num = numpnts(qw) |
---|
1177 | |
---|
1178 | ii=0 |
---|
1179 | numBad=0 |
---|
1180 | do |
---|
1181 | mask = mw[ii] |
---|
1182 | if(numtype(mask) != 2) //if not NaN |
---|
1183 | numBad += 1 |
---|
1184 | else //keep the points that are NaN |
---|
1185 | qw[ii-numBad] = qw[ii] |
---|
1186 | iw[ii-numBad] = iw[ii] |
---|
1187 | sw[ii-numBad] = sw[ii] |
---|
1188 | endif |
---|
1189 | ii += 1 |
---|
1190 | while(ii<num) |
---|
1191 | //trim the end of the waves |
---|
1192 | DeletePoints num-numBad, numBad, qw,iw,sw |
---|
1193 | return(numBad) |
---|
1194 | end |
---|
1195 | |
---|
1196 | // produces the _msk waves that have the new number of data points |
---|
1197 | // |
---|
1198 | Function DSM_MaskDoneButton(ctrlName) : ButtonControl |
---|
1199 | String ctrlName |
---|
1200 | |
---|
1201 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1202 | |
---|
1203 | // Variable aExists= strlen(CsrInfo(A)) > 0 //Igor 5 |
---|
1204 | Variable aExists= strlen(CsrWave(A)) > 0 //Igor 4 |
---|
1205 | if(!aExists) |
---|
1206 | return(1) //possibly reverted data, no cursor, no Mask wave |
---|
1207 | endif |
---|
1208 | |
---|
1209 | Duplicate/O $(USANSFolder+":DSM:Q_exp_orig"),$(USANSFolder+":DSM:Q_msk") |
---|
1210 | Duplicate/O $(USANSFolder+":DSM:I_exp_orig"),$(USANSFolder+":DSM:I_msk") |
---|
1211 | Duplicate/O $(USANSFolder+":DSM:S_exp_orig"),$(USANSFolder+":DSM:S_msk") |
---|
1212 | Wave Q_msk=$(USANSFolder+":DSM:Q_msk") |
---|
1213 | Wave I_msk=$(USANSFolder+":DSM:I_msk") |
---|
1214 | Wave S_msk=$(USANSFolder+":DSM:S_msk") |
---|
1215 | |
---|
1216 | //finish up - trim the data sets and reassign the working set |
---|
1217 | Wave MaskData=$(USANSFolder+":DSM:MaskData") |
---|
1218 | |
---|
1219 | RemoveMaskedPoints(MaskData,Q_msk,I_msk,S_msk) |
---|
1220 | |
---|
1221 | //reset the number of points |
---|
1222 | NVAR gNq = $(USANSFolder+":DSM:gNq") |
---|
1223 | gNq = numpnts(Q_msk) |
---|
1224 | |
---|
1225 | Cursor/K A |
---|
1226 | HideInfo |
---|
1227 | |
---|
1228 | return(0) |
---|
1229 | End |
---|
1230 | |
---|
1231 | |
---|
1232 | // not quite the same as revert |
---|
1233 | Function DSM_ClearMaskProc(ctrlName) : ButtonControl |
---|
1234 | String ctrlName |
---|
1235 | |
---|
1236 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1237 | |
---|
1238 | Wave MaskData=$(USANSFolder+":DSM:MaskData") |
---|
1239 | MaskData = NaN |
---|
1240 | |
---|
1241 | return(0) |
---|
1242 | end |
---|
1243 | |
---|
1244 | // when the mask button is pressed, A must be on the graph |
---|
1245 | // Displays MaskData wave on the graph |
---|
1246 | // |
---|
1247 | Function DSM_MyMaskProc(ctrlName) : ButtonControl |
---|
1248 | String ctrlName |
---|
1249 | |
---|
1250 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1251 | |
---|
1252 | Wave data=$(USANSFolder+":DSM:I_exp_orig") |
---|
1253 | |
---|
1254 | // Variable aExists= strlen(CsrInfo(A)) > 0 //Igor 5 |
---|
1255 | Variable aExists= strlen(CsrWave(A)) > 0 //Igor 4 |
---|
1256 | |
---|
1257 | // need to get rid of old smoothed data if data is re-masked |
---|
1258 | Execute "RemoveSmoothed()" |
---|
1259 | SetDataFolder $(USANSFolder+":DSM") |
---|
1260 | Killwaves/Z I_smth,Q_smth,S_smth |
---|
1261 | |
---|
1262 | if(aExists) //mask the selected point |
---|
1263 | // toggle NaN (keep) or Data value (= masked) |
---|
1264 | Wave MaskData |
---|
1265 | MaskData[pcsr(A)] = (numType(MaskData[pcsr(A)])==0) ? NaN : data[pcsr(A)] //if NaN, doesn't plot |
---|
1266 | else |
---|
1267 | Cursor /A=1/H=1/L=1/P/W=Desmear_Graph A I_exp_orig leftx(I_exp_orig) |
---|
1268 | ShowInfo |
---|
1269 | //if the mask wave does not exist, make one |
---|
1270 | if(exists("MaskData") != 1) |
---|
1271 | Duplicate/O Q_exp_orig MaskData |
---|
1272 | MaskData = NaN //use all data |
---|
1273 | endif |
---|
1274 | Execute "AppendMask()" |
---|
1275 | endif |
---|
1276 | |
---|
1277 | SetDataFolder root: |
---|
1278 | |
---|
1279 | return(0) |
---|
1280 | End |
---|
1281 | |
---|
1282 | // when the mask button is pressed, A must be on the graph |
---|
1283 | // Displays MaskData wave on the graph |
---|
1284 | // |
---|
1285 | Function DSM_MaskLTCursor(ctrlName) : ButtonControl |
---|
1286 | String ctrlName |
---|
1287 | |
---|
1288 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1289 | |
---|
1290 | // Variable aExists= strlen(CsrInfo(A)) > 0 //Igor 5 |
---|
1291 | Variable aExists= strlen(CsrWave(A)) > 0 //Igor 4 |
---|
1292 | |
---|
1293 | if(!aExists) |
---|
1294 | return(1) |
---|
1295 | endif |
---|
1296 | // need to get rid of old smoothed data if data is re-masked |
---|
1297 | Execute "RemoveSmoothed()" |
---|
1298 | SetDataFolder $(USANSFolder+":DSM") |
---|
1299 | Killwaves/Z I_smth,Q_smth,S_smth |
---|
1300 | |
---|
1301 | Wave data=I_exp_orig |
---|
1302 | |
---|
1303 | Variable pt,ii |
---|
1304 | pt = pcsr(A) |
---|
1305 | for(ii=pt;ii>=0;ii-=1) |
---|
1306 | // toggle NaN (keep) or Data value (= masked) |
---|
1307 | Wave MaskData |
---|
1308 | MaskData[ii] = (numType(MaskData[ii])==0) ? NaN : data[ii] //if NaN, doesn't plot |
---|
1309 | endfor |
---|
1310 | |
---|
1311 | SetDataFolder root: |
---|
1312 | return(0) |
---|
1313 | End |
---|
1314 | |
---|
1315 | // when the mask button is pressed, A must be on the graph |
---|
1316 | // Displays MaskData wave on the graph |
---|
1317 | // |
---|
1318 | Function DSM_MaskGTCursor(ctrlName) : ButtonControl |
---|
1319 | String ctrlName |
---|
1320 | |
---|
1321 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1322 | |
---|
1323 | |
---|
1324 | // Variable aExists= strlen(CsrInfo(A)) > 0 //Igor 5 |
---|
1325 | Variable aExists= strlen(CsrWave(A)) > 0 //Igor 4 |
---|
1326 | |
---|
1327 | if(!aExists) |
---|
1328 | return(1) |
---|
1329 | endif |
---|
1330 | // need to get rid of old smoothed data if data is re-masked |
---|
1331 | Execute "RemoveSmoothed()" |
---|
1332 | SetDataFolder $(USANSFolder+":DSM") |
---|
1333 | Killwaves/Z I_smth,Q_smth,S_smth |
---|
1334 | |
---|
1335 | Wave data=I_exp_orig |
---|
1336 | |
---|
1337 | Variable pt,ii,endPt |
---|
1338 | endPt=numpnts(MaskData) |
---|
1339 | pt = pcsr(A) |
---|
1340 | for(ii=pt;ii<endPt;ii+=1) |
---|
1341 | // toggle NaN (keep) or Data value (= masked) |
---|
1342 | Wave MaskData |
---|
1343 | MaskData[ii] = (numType(MaskData[ii])==0) ? NaN : data[ii] //if NaN, doesn't plot |
---|
1344 | endfor |
---|
1345 | |
---|
1346 | SetDataFolder root: |
---|
1347 | |
---|
1348 | return(0) |
---|
1349 | End |
---|
1350 | |
---|
1351 | Function CleanUpJunk() |
---|
1352 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1353 | |
---|
1354 | // clean up the old junk on the graph, /Z for no error |
---|
1355 | Execute "RemoveExtrapolation()" |
---|
1356 | Execute "RemoveDesmeared()" |
---|
1357 | Execute "RemoveSmoothed()" |
---|
1358 | Execute "RemoveMask()" |
---|
1359 | |
---|
1360 | //remove the cursor |
---|
1361 | Cursor/K A |
---|
1362 | |
---|
1363 | //always initialize these |
---|
1364 | String/G $(USANSFolder+":DSM:gStr1") = "" |
---|
1365 | String/G $(USANSFolder+":DSM:gStr2") = "" |
---|
1366 | String/G $(USANSFolder+":DSM:gIterStr") = "" |
---|
1367 | |
---|
1368 | // clean up the old waves from smoothing and desmearing steps |
---|
1369 | SetDataFolder $(USANSFolder+":DSM") |
---|
1370 | Killwaves/Z I_smth,I_dsm,I_dsm_sm,Q_smth,Q_dsm,S_smth,S_dsm,Yi_SS,Yq_SS |
---|
1371 | Killwaves/Z Weights,FW,R_wave,S_wave,H_wave,Di,Ci |
---|
1372 | Killwaves/Z Is_old,I_old,err |
---|
1373 | Killwaves/Z Q_ext,I_ext,S_ext,hqq,hqi,lqq,lqi |
---|
1374 | Killwaves/Z MaskData,Q_msk,I_msk,S_msk |
---|
1375 | Killwaves/Z Q_work,I_work,S_work //working waves for desmearing step |
---|
1376 | SetDataFolder root: |
---|
1377 | End |
---|
1378 | |
---|
1379 | // does not alter the data sets - reports a power law |
---|
1380 | // exponent and makes it global so it will automatically |
---|
1381 | // be used during the desmearing |
---|
1382 | // |
---|
1383 | // generates extr_hqi vs extr_hqq that are Appended to the graph |
---|
1384 | Function DSM_ExtrapolateButtonProc(ctrlName) : ButtonControl |
---|
1385 | String ctrlName |
---|
1386 | |
---|
1387 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1388 | |
---|
1389 | NVAR nend = $(USANSFolder+":DSM:gNptsExtrap") |
---|
1390 | NVAR m_pred = $(USANSFolder+":DSM:gPowerM") |
---|
1391 | |
---|
1392 | SetDataFolder $(USANSFolder+":DSM") |
---|
1393 | //use masked data if it exists |
---|
1394 | if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1) |
---|
1395 | wave Qw = $(USANSFolder+":DSM:Q_msk") |
---|
1396 | Wave Iw = $(USANSFolder+":DSM:I_msk") |
---|
1397 | Wave Sw = $(USANSFolder+":DSM:S_msk") |
---|
1398 | else |
---|
1399 | //start from the "_exp" waves |
---|
1400 | if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1) |
---|
1401 | wave Qw = $(USANSFolder+":DSM:Q_exp") |
---|
1402 | Wave Iw = $(USANSFolder+":DSM:I_exp") |
---|
1403 | Wave Sw = $(USANSFolder+":DSM:S_exp") |
---|
1404 | endif |
---|
1405 | endif |
---|
1406 | SetDataFolder root: |
---|
1407 | |
---|
1408 | DSM_SetExtrWaves(Qw) |
---|
1409 | m_pred = DSM_DoExtrapolate(Qw,Iw,Sw,nend) |
---|
1410 | AppendExtrapolation() |
---|
1411 | |
---|
1412 | return(0) |
---|
1413 | End |
---|
1414 | |
---|
1415 | //smooths the data in steps as requested... |
---|
1416 | // |
---|
1417 | // typically do a simple Box smooth first, |
---|
1418 | // then do a smoothing spline, keeping the same number of data points |
---|
1419 | // |
---|
1420 | // chi-squared is reported - so you can see how "bad" the smoothing is |
---|
1421 | // smoothing of the data is done on a log(I) scale if requested |
---|
1422 | // setting doLog variable to 0 will return to linear smoothing |
---|
1423 | // (I see little difference) |
---|
1424 | // |
---|
1425 | // start from the "_smth" waves if they exist |
---|
1426 | // otheriwse start from the working waves |
---|
1427 | // |
---|
1428 | // create Q_smth, I_smth, S_smth |
---|
1429 | // keep track of the smoothing types and passes in the wave note |
---|
1430 | // |
---|
1431 | Function DSM_SmoothButtonProc(ctrlName) : ButtonControl |
---|
1432 | String ctrlName |
---|
1433 | |
---|
1434 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1435 | |
---|
1436 | SetDataFolder $(USANSFolder+":DSM") |
---|
1437 | |
---|
1438 | Variable ii,new_n,pass,nq_ext,offset,doLog=1 |
---|
1439 | String noteStr="" |
---|
1440 | |
---|
1441 | // want log scaling of intensity during smoothing? |
---|
1442 | ControlInfo DSMControl_3g |
---|
1443 | doLog = V_value |
---|
1444 | |
---|
1445 | //// if(exists("I_smth")==1 && exists("Q_smth")==1 && exists("S_smth")==1 && ! freshMask) |
---|
1446 | if(exists("I_smth")==1 && exists("Q_smth")==1 && exists("S_smth")==1) |
---|
1447 | //start from the smoothed waves --- just need the wave references |
---|
1448 | else |
---|
1449 | //start from the "msk", creating smth waves |
---|
1450 | if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1) |
---|
1451 | wave Q_msk,I_msk,S_msk |
---|
1452 | Duplicate/O I_msk,I_smth,Q_smth,S_smth |
---|
1453 | I_smth = I_msk |
---|
1454 | Q_smth = Q_msk |
---|
1455 | S_smth = S_msk |
---|
1456 | else |
---|
1457 | //start from the "_exp" waves |
---|
1458 | if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1) |
---|
1459 | wave Q_exp,I_exp,S_exp |
---|
1460 | Duplicate/O I_exp,I_smth,Q_smth,S_smth |
---|
1461 | I_smth = I_exp |
---|
1462 | Q_smth = Q_exp |
---|
1463 | S_smth = S_exp |
---|
1464 | endif |
---|
1465 | endif |
---|
1466 | endif |
---|
1467 | |
---|
1468 | wave Q_smth,I_smth,S_smth |
---|
1469 | |
---|
1470 | // extend the data to avoid end effects |
---|
1471 | //creates I_ext,Q_ext,S_ext with extended q-values |
---|
1472 | // fit 15 pts at each end, typically add 10 pts to each end |
---|
1473 | // does nothing if num_extr = 0 |
---|
1474 | ControlInfo/W=Desmear_Graph DSMControl_3d //extendCheck |
---|
1475 | if(V_value == 1) |
---|
1476 | ExtendToSmooth(Q_smth,I_smth,S_smth,15,15,10) |
---|
1477 | else |
---|
1478 | ExtendToSmooth(Q_smth,I_smth,S_smth,15,15,0) //don't extend, just use the original data |
---|
1479 | endif |
---|
1480 | |
---|
1481 | //whether extending or not, the working data is "_ext", set by ExtendToSmooth() |
---|
1482 | SetDataFolder $(USANSFolder+":DSM") |
---|
1483 | wave Q_ext,I_ext ,S_ext |
---|
1484 | |
---|
1485 | noteStr=note(I_smth) |
---|
1486 | Note I_ext , noteStr |
---|
1487 | |
---|
1488 | WaveStats/Q I_ext |
---|
1489 | if(V_min<=0) |
---|
1490 | Print "negative itensity found, using linear scale for smoothing" |
---|
1491 | doLog = 0 |
---|
1492 | endif |
---|
1493 | |
---|
1494 | if(doLog) |
---|
1495 | //convert to log scale |
---|
1496 | Duplicate/O I_ext I_log,I_log_err |
---|
1497 | Wave I_log ,I_log_err |
---|
1498 | I_log = log(I_ext) |
---|
1499 | WaveStats/Q I_log |
---|
1500 | offset = 2*(-V_min) |
---|
1501 | I_log += offset |
---|
1502 | I_log_err = S_ext/(2.30*I_ext) |
---|
1503 | I_ext = I_log |
---|
1504 | endif |
---|
1505 | // |
---|
1506 | |
---|
1507 | ControlInfo/W=Desmear_Graph DSMControl_3b //BoxCheck |
---|
1508 | if(V_value == 1) |
---|
1509 | //do a simple Box smooth first - number of points does not change |
---|
1510 | // fills ends with neighboring value (E=3) (n=3 points in smoothing window) |
---|
1511 | Smooth/E=3/B 3, I_ext |
---|
1512 | |
---|
1513 | noteStr=note(I_ext) |
---|
1514 | pass = NumberByKey("BOX", noteStr, "=", ";") |
---|
1515 | noteStr = ReplaceNumberByKey("BOX", noteStr, pass+1, "=", ";") |
---|
1516 | Note/K I_ext |
---|
1517 | Note I_ext , noteStr |
---|
1518 | // Print "box = ",noteStr |
---|
1519 | endif |
---|
1520 | |
---|
1521 | NVAR sParam = gSmoothFac //already in the right DF |
---|
1522 | |
---|
1523 | ControlInfo/W=Desmear_Graph DSMControl_3c //SSCheck |
---|
1524 | if(V_value == 1) |
---|
1525 | // Interpolate2/T=3/N=(nq)/I=2/F=1/SWAV=S_ext/Y=Yi_SS/X=Yq_SS Q_ext, I_ext //Igor 5 |
---|
1526 | // Interpolate/T=3/N=(nq)/I=2/F=1/S=S_ext/Y=Yi_SS/X=Yq_SS I_ext /X=Q_ext //Igor 4 |
---|
1527 | //Igor 4 |
---|
1528 | String str="" |
---|
1529 | nq_ext = numpnts(Q_ext) |
---|
1530 | str = "Interpolate/T=3/N=("+num2str(nq_ext)+")/I=1/F=("+num2str(sParam)+")/Y=Yi_SS/X=Yq_SS I_ext /X=Q_ext" |
---|
1531 | Execute str |
---|
1532 | // Print Str |
---|
1533 | // Interpolate2/T=3/N=(nq_ext)/I=1/F=(sParam)/Y=Yi_SS/X=Yq_SS Q_ext, I_ext |
---|
1534 | // end Igor 4 |
---|
1535 | wave yi_ss = yi_ss // already in the right DF |
---|
1536 | wave yq_ss = yq_ss |
---|
1537 | // reassign the "working" waves with the result of interpolate, which has the same I/Q values |
---|
1538 | I_ext = yi_ss |
---|
1539 | Q_ext = yq_ss |
---|
1540 | |
---|
1541 | noteStr=note(I_ext) |
---|
1542 | pass = NumberByKey("SPLINE", noteStr, "=", ";") |
---|
1543 | noteStr = ReplaceNumberByKey("SPLINE", noteStr, pass+1, "=", ";") |
---|
1544 | Note/K I_ext |
---|
1545 | Note I_ext , noteStr |
---|
1546 | // Print "spline = ",noteStr |
---|
1547 | endif |
---|
1548 | |
---|
1549 | //undo the scaling |
---|
1550 | If(doLog) |
---|
1551 | I_ext -= offset |
---|
1552 | I_ext = 10^I_ext |
---|
1553 | endif |
---|
1554 | |
---|
1555 | // at this point, I_ext has too many points - and we need to just return the |
---|
1556 | // center chunk that is the original q-values |
---|
1557 | // as assign this to the working set for the desmear step |
---|
1558 | // and to the _smth set in case another smoothng pass is desired |
---|
1559 | // Q_smth has not changed, S_smth has not changed |
---|
1560 | I_smth = interp(Q_smth[p], Q_ext, I_ext ) |
---|
1561 | |
---|
1562 | Note/K I_smth |
---|
1563 | Note I_smth , noteStr |
---|
1564 | // Print "end of smoothed, note = ",note(I_smth) |
---|
1565 | |
---|
1566 | Variable nq = numpnts($(USANSFolder+":DSM:Q_smth")) |
---|
1567 | // Print "nq after smoothing = ",nq |
---|
1568 | |
---|
1569 | //reset the global |
---|
1570 | NVAR gNq = gNq |
---|
1571 | gNq = nq |
---|
1572 | //report the chi^2 difference between the smoothed curve and the experimental data |
---|
1573 | NVAR chi2 = gChi2Smooth |
---|
1574 | chi2 = SmoothedChi2(I_smth) |
---|
1575 | |
---|
1576 | Execute "AppendSmoothed()" |
---|
1577 | |
---|
1578 | setDataFolder root: |
---|
1579 | return(0) |
---|
1580 | End |
---|
1581 | |
---|
1582 | // undo the smoothing and start over - useful if you've smoothed |
---|
1583 | // too aggressively and need to back off |
---|
1584 | Function DSM_SmoothUndoButtonProc(ctrlName) : ButtonControl |
---|
1585 | String ctrlName |
---|
1586 | |
---|
1587 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1588 | |
---|
1589 | Execute "RemoveSmoothed()" |
---|
1590 | SetDataFolder $(USANSFolder+":DSM") |
---|
1591 | Killwaves/Z I_smth,Q_smth,S_smth,Q_ext,I_ext,S_ext,Yi_SS,Yq_SS |
---|
1592 | SetDataFolder root: |
---|
1593 | return (0) |
---|
1594 | end |
---|
1595 | |
---|
1596 | //calculate the chi^2 value for the smoothed data |
---|
1597 | Function SmoothedChi2(I_smth) |
---|
1598 | Wave I_smth |
---|
1599 | |
---|
1600 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1601 | |
---|
1602 | //start from the "msk", if they exist |
---|
1603 | if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1) |
---|
1604 | Wave iw = $(USANSFolder+":DSM:I_msk") |
---|
1605 | Wave sw = $(USANSFolder+":DSM:S_msk") |
---|
1606 | else |
---|
1607 | //start from the "_exp" waves |
---|
1608 | if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1) |
---|
1609 | Wave iw = $(USANSFolder+":DSM:I_exp") |
---|
1610 | Wave sw = $(USANSFolder+":DSM:S_exp") |
---|
1611 | endif |
---|
1612 | endif |
---|
1613 | |
---|
1614 | Variable ii,chi2,num |
---|
1615 | chi2=0 |
---|
1616 | num = numpnts(iw) |
---|
1617 | for(ii=0;ii<num;ii+=1) |
---|
1618 | chi2 += (iw[ii] - I_smth[ii])^2 / (sw[ii])^2 |
---|
1619 | endfor |
---|
1620 | Chi2 /= (num-1) |
---|
1621 | return (chi2) |
---|
1622 | end |
---|
1623 | |
---|
1624 | // I_dsm is the desmeared data |
---|
1625 | // |
---|
1626 | // step (7) - desmearing is done, write out the result |
---|
1627 | // |
---|
1628 | Function DSM_SaveButtonProc(ctrlName) : ButtonControl |
---|
1629 | String ctrlName |
---|
1630 | |
---|
1631 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1632 | NVAR useXMLOutput = root:Packages:NIST:gXML_Write |
---|
1633 | |
---|
1634 | String saveStr |
---|
1635 | SVAR curFile = $(USANSFolder+":DSM:gCurFile") |
---|
1636 | saveStr = CleanupName((curFile),0) //the output filename |
---|
1637 | // |
---|
1638 | |
---|
1639 | if (useXMLOutput == 1) |
---|
1640 | WriteXMLUSANSDesmeared(saveStr,0,0,1) |
---|
1641 | else |
---|
1642 | WriteUSANSDesmeared(saveStr,0,0,1) //use the full set (lo=hi=0) and present a dialog |
---|
1643 | endif |
---|
1644 | |
---|
1645 | SetDataFolder root: |
---|
1646 | return(0) |
---|
1647 | End |
---|
1648 | |
---|
1649 | Function DSM_HelpButtonProc(ctrlName) : ButtonControl |
---|
1650 | String ctrlName |
---|
1651 | |
---|
1652 | DisplayHelpTopic/Z/K=1 "Desmearing USANS Data" |
---|
1653 | if(V_flag !=0) |
---|
1654 | DoAlert 0,"The Desmearing USANS Data Help file could not be found" |
---|
1655 | endif |
---|
1656 | return(0) |
---|
1657 | End |
---|
1658 | |
---|
1659 | |
---|
1660 | //toggles the log/lin display of the loaded data set |
---|
1661 | Function DSM_LoadCheckProc(ctrlName,checked) : CheckBoxControl |
---|
1662 | String ctrlName |
---|
1663 | Variable checked |
---|
1664 | |
---|
1665 | ModifyGraph log=(checked) |
---|
1666 | return(0) |
---|
1667 | End |
---|
1668 | |
---|
1669 | |
---|
1670 | Function DSM_ExtrapolationCheckProc(ctrlName,checked) : CheckBoxControl |
---|
1671 | String ctrlName |
---|
1672 | Variable checked |
---|
1673 | |
---|
1674 | if(checked) |
---|
1675 | AppendExtrapolation() |
---|
1676 | else |
---|
1677 | RemoveExtrapolation() |
---|
1678 | endif |
---|
1679 | return(0) |
---|
1680 | End |
---|
1681 | |
---|
1682 | |
---|
1683 | // takes as input the "working waves" |
---|
1684 | // |
---|
1685 | // creates intermediate waves to work on |
---|
1686 | // |
---|
1687 | // output of Q_dsm,I_dsm,S_dsm (and I_dsm_sm, the result of smearing I_dsm) |
---|
1688 | // |
---|
1689 | Function DSM_DesmearButtonProc(ctrlName) : ButtonControl |
---|
1690 | String ctrlName |
---|
1691 | |
---|
1692 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1693 | |
---|
1694 | SetDataFolder $(USANSFolder+":DSM") |
---|
1695 | if(exists("I_smth")==1 && exists("Q_smth")==1 && exists("S_smth")==1) |
---|
1696 | wave Q_smth , I_smth ,S_smth |
---|
1697 | Duplicate/O I_smth,I_work,Q_work,S_work |
---|
1698 | I_work = I_smth |
---|
1699 | Q_work = Q_smth |
---|
1700 | S_work = S_smth |
---|
1701 | else |
---|
1702 | //start from the "msk", creating work waves |
---|
1703 | if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1) |
---|
1704 | wave Q_msk,I_msk,S_msk |
---|
1705 | Duplicate/O I_msk,I_work,Q_work,S_work |
---|
1706 | I_work = I_msk |
---|
1707 | Q_work = Q_msk |
---|
1708 | S_work = S_msk |
---|
1709 | else |
---|
1710 | //start from the "_exp" waves |
---|
1711 | if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1) |
---|
1712 | wave Q_exp,I_exp,S_exp |
---|
1713 | Duplicate/O I_exp,I_work,Q_work,S_work |
---|
1714 | I_work = I_exp |
---|
1715 | Q_work = Q_exp |
---|
1716 | S_work = S_exp |
---|
1717 | endif |
---|
1718 | endif |
---|
1719 | endif |
---|
1720 | //SetDataFolder root: |
---|
1721 | |
---|
1722 | NVAR nq = gNq |
---|
1723 | NVAR m = gPowerM |
---|
1724 | NVAR chi2_target = gChi2Target |
---|
1725 | NVAR maxFastIter = gMaxFastIter |
---|
1726 | NVAR maxSlowIter = gMaxSlowIter |
---|
1727 | |
---|
1728 | // SET WEIGHTING OF EXPERIMENTAL DATA. |
---|
1729 | Duplicate/O Q_work weights |
---|
1730 | Wave weights = weights |
---|
1731 | weights = 1/S_work^2 |
---|
1732 | |
---|
1733 | // calculate weighting array for smearing of data |
---|
1734 | Make/O/D/N=(nq,nq) FW |
---|
1735 | Wave FW |
---|
1736 | Weights_L(m,FW,Q_work) |
---|
1737 | |
---|
1738 | // ^^^^ Iterative desmearing ^^^^* |
---|
1739 | Variable chi2_old,chi2_new,done,iter |
---|
1740 | // FOR 0TH ITERATION, EXPERIMENTAL DATA IS USED FOR Y_OLD, create ys_old for result |
---|
1741 | // y_old = I_old, y_new = I_dsm, I_dsm_sm = ys_new, |
---|
1742 | // duplicate preserves the wave note! |
---|
1743 | Duplicate/O I_work I_old,Is_old,I_dsm,I_dsm_sm |
---|
1744 | Duplicate/O Q_work Q_dsm,S_dsm //sets Q_dsm correctly |
---|
1745 | wave S_dsm,I_old,Is_old,I_dsm,I_dsm_sm |
---|
1746 | I_old = I_work |
---|
1747 | Is_old = 0 |
---|
1748 | I_dsm = 0 |
---|
1749 | I_dsm_sm = 0 |
---|
1750 | |
---|
1751 | // Smear 0TH iter guess Y --> YS |
---|
1752 | DSM_Smear(I_old,Is_old,NQ,FW) |
---|
1753 | chi2_OLD = chi2(Is_old,I_work,weights,NQ) |
---|
1754 | |
---|
1755 | printf "starting chi^2 = %g\r",chi2_old |
---|
1756 | Print "Starting fast convergence..... " |
---|
1757 | |
---|
1758 | done = 0 //false |
---|
1759 | iter = 0 |
---|
1760 | Do // while not done, use Fast convergence |
---|
1761 | I_dsm = I_work * I_old / Is_old // Calculate corrected guess (no need for do-loop) |
---|
1762 | |
---|
1763 | // Smear iter guess I_dsm --> I_dsm_sm and see how well I_dsm_sm matches experimental data |
---|
1764 | DSM_Smear(I_dsm,I_dsm_sm,NQ,FW) |
---|
1765 | chi2_new = chi2(I_dsm_sm,I_work,weights,NQ) |
---|
1766 | |
---|
1767 | // Stop iteration if fit from new iteration has worse fit... |
---|
1768 | If (chi2_new > chi2_old) |
---|
1769 | Done = 1 |
---|
1770 | Endif |
---|
1771 | |
---|
1772 | // Stop iteration if fit is better than target value... |
---|
1773 | If (chi2_new < chi2_target) |
---|
1774 | Done = 1 |
---|
1775 | Else |
---|
1776 | Iter += 1 |
---|
1777 | Printf "Iteration %d, Chi^2(new) = %g\r", Iter,chi2_new |
---|
1778 | I_old = I_dsm |
---|
1779 | Is_old = I_dsm_sm |
---|
1780 | if(iter>maxFastIter) |
---|
1781 | break |
---|
1782 | endif |
---|
1783 | CHI2_OLD = CHI2_NEW |
---|
1784 | EndIf |
---|
1785 | while( !done ) |
---|
1786 | |
---|
1787 | // append I_dsm,I_exp to the graph |
---|
1788 | Execute "AppendDesmeared()" |
---|
1789 | DoUpdate |
---|
1790 | |
---|
1791 | SetDataFolder $(USANSFolder+":DSM") |
---|
1792 | |
---|
1793 | |
---|
1794 | // step (6) - refine the desmearing using slow convergence |
---|
1795 | Print "Starting slow convergence..... " |
---|
1796 | done = 0 // ! reset flag for slow convergence |
---|
1797 | Do |
---|
1798 | I_dsm = I_old + (I_work - Is_old) // Calculate corrected guess |
---|
1799 | |
---|
1800 | // Smear iter guess Y --> YS |
---|
1801 | DSM_Smear(I_dsm,I_dsm_sm,NQ,FW) |
---|
1802 | chi2_new = chi2(I_dsm_sm,I_work,weights,NQ) |
---|
1803 | |
---|
1804 | // Stop iteration if fit from new iteration has worse fit... |
---|
1805 | If (chi2_new > chi2_old) |
---|
1806 | Done = 1 |
---|
1807 | Endif |
---|
1808 | |
---|
1809 | // Stop iteration if fit is better than target value... |
---|
1810 | If (chi2_new < chi2_target) |
---|
1811 | Done = 1 |
---|
1812 | Else |
---|
1813 | Iter += 1 |
---|
1814 | |
---|
1815 | if(mod(iter, 50 ) ==0 ) |
---|
1816 | Printf "Iteration %d, Chi^2(new) = %g\r", Iter,chi2_new |
---|
1817 | DoUpdate |
---|
1818 | endif |
---|
1819 | I_old = I_dsm |
---|
1820 | Is_old = I_dsm_sm |
---|
1821 | |
---|
1822 | CHI2_OLD = CHI2_NEW |
---|
1823 | EndIf |
---|
1824 | if(iter>maxSlowIter) |
---|
1825 | break |
---|
1826 | endif |
---|
1827 | While ( !done ) |
---|
1828 | |
---|
1829 | Printf "Iteration %d, Chi^2(new) = %g\r", Iter,chi2_new |
---|
1830 | |
---|
1831 | // adjust the error |
---|
1832 | SetDataFolder $(USANSFolder+":DSM") |
---|
1833 | Duplicate/O S_work err |
---|
1834 | S_dsm = abs(err/I_work*I_dsm) //proportional error |
---|
1835 | //John simply keeps the same error as was read in from the smeared data set - is this right? |
---|
1836 | // S_dsm = S_Work |
---|
1837 | |
---|
1838 | NVAR gChi = gChi2Final //chi^2 final |
---|
1839 | NVAR gIter = gIterations //total number of iterations |
---|
1840 | gChi = chi2_new |
---|
1841 | gIter = Iter |
---|
1842 | |
---|
1843 | SVAR str = gIterStr |
---|
1844 | sprintf str,"%d iterations required",iter |
---|
1845 | |
---|
1846 | SetDataFolder root: |
---|
1847 | return(0) |
---|
1848 | End |
---|
1849 | |
---|
1850 | Function DSM_RevertButtonProc(ctrlName) : ButtonControl |
---|
1851 | String ctrlName |
---|
1852 | |
---|
1853 | SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder |
---|
1854 | |
---|
1855 | |
---|
1856 | CleanUpJunk() |
---|
1857 | |
---|
1858 | SetDataFolder $(USANSFolder+":DSM") |
---|
1859 | |
---|
1860 | //reset the working waves to the original |
---|
1861 | wave Q_exp_orig,I_exp_orig,S_exp_orig |
---|
1862 | |
---|
1863 | Duplicate/O Q_exp_orig Q_exp |
---|
1864 | Duplicate/O I_exp_orig I_exp |
---|
1865 | Duplicate/O S_exp_orig S_exp |
---|
1866 | // kill the data folder |
---|
1867 | // re-initialize? |
---|
1868 | |
---|
1869 | SetDataFolder root: |
---|
1870 | |
---|
1871 | return(0) |
---|
1872 | End |
---|