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