Context Navigation

LakeDesmearing_JB.ipf @ 1074

Last change on this file since 1074 was 742, checked in by srkline, 13 years ago

fix for ticket #279: calculation of the slope when desmearing.

Now the behavior is for the desmearing panel to use the slope as found when loading the .cor data. For a fresh file, it proceeds as usual. For a previously-loaded file, be sure to NOT re-load the file (<no> from the dialog) so that the slope is not recalculated, and the data is simply plotted. Then the slope information is copied along with the data to the desmearing folder.

The changed behavior is when selecting the "Extrapolate" tab. Now the extrapolation is automatically calculated from the loading parameters (not re-fit) and displayed on the graph. If the extrapolation is good, nothing needs to be done. In fact, since the extrapolation (slope) should be acceptable from the load, the mask, extrapolate, and smooth tabs can all be skipped, proceeding directly to the desmear tab.

Property eol-style set to native

File size: 53.9 KB

Line
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 = 15 //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 = UUQ_exp[jj+1]+Q_exp[ii]^2ln(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
747	/////////for the fit
748	Make/O/D P_coef={0,1,-4} //input
749	Make/O/T Constr={"K2<0","K2 > -8"}
750	//(set background to zero and hold fixed)
751
752	// initial guess
753	P_coef[1] = iw[num-1]/qw[num-1]^P_coef[2]
754
755	CurveFit/H="100" Power kwCWave=P_coef iw[(num-1-nend),(num-1)] /X=qw /W=sw /I=1 /C=constr
756	extr_hqi=P_coef[0]+P_coef[1]*extr_hqq^P_coef[2]
757
758	// for the case of data with a background
759	// Make/O/D P_coef={0,1,-4} //input
760	// //(set background to iw[num-1], let it be a free parameter
761	// P_coef[0] = iw[num-1]
762	// CurveFit Power kwCWave=P_coef iw[(num-1-nend),(num-1)] /X=qw /W=sw /D
763	// extr_hqi=P_coef[0]+P_coef[1]*extr_hqq^P_coef[2]
764
765	// if(checked && if(not already displayed))
766	// AppendToGraph extr_hqi vs extr_hqq
767	// ModifyGraph lsize(extr_hqi)=2
768	// endif
769
770
771	Printf "Smeared Power law exponent = %g\r",P_coef[2]
772	Printf "**** For Desmearing, use a Power law exponent of %5.1f\r",P_coef[2]-1
773
774	retVal = P_coef[2]-1
775	SetDataFolder root:
776	return(retVal)
777	End
778
779	Function DSM_Guinier_Fit(w,x) //: FitFunc
780	Wave w
781	Variable x
782
783	//fit data to I(q) = Aexp(Bq^2)
784	// (B will be negative)
785	//two parameters
786	Variable a,b,ans
787	a=w[0]
788	b=w[1]
789	ans = aexp(bx*x)
790	return(ans)
791	End
792
793	Proc Desmear_Graph()
794	String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
795
796	PauseUpdate; Silent 1 // building window...
797	Display /W=(5,44,408,558) /K=1
798	ModifyGraph cbRGB=(51664,44236,58982)
799	DoWindow/C Desmear_Graph
800	ControlBar 160
801	// break into tabs
802	TabControl DSM_Tab,pos={5,3},size={392,128},proc=DSM_TabProc
803	TabControl DSM_Tab,labelBack=(49151,49152,65535),tabLabel(0)="Load"
804	TabControl DSM_Tab,tabLabel(1)="Mask",tabLabel(2)="Extrapolate"
805	TabControl DSM_Tab,tabLabel(3)="Smooth",tabLabel(4)="Desmear",value= 0
806
807	//always visible - revert and save
808	//maybe the wrong place here?
809	Button DSMControlA,pos={225,135},size={80,20},proc=DSM_RevertButtonProc,title="Revert"
810	Button DSMControlA,help={"Revert the smeared data to its original state and start over"}
811	Button DSMControlB,pos={325,135},size={50,20},proc=DSM_SaveButtonProc,title="Save"
812	Button DSMControlB,help={"Save the desmeared data set"}
813	Button DSMControlC,pos={25,135},size={50,20},proc=DSM_HelpButtonProc,title="Help"
814	Button DSMControlC,help={"Show the help file for desmearing"}
815
816	// add the controls to each tab ---- all names start with "DSMControl_"
817
818	//tab(0) Load - initially visible
819	Button DSMControl_0a,pos={23,39},size={80,20},proc=DSM_LoadButtonProc,title="Load Data"
820	Button DSMControl_0a,help={"Load slit-smeared USANS data = \".cor\" files"}
821	CheckBox DSMControl_0b,pos={26,74},size={80,14},proc=DSM_LoadCheckProc,title="Log Axes?"
822	CheckBox DSMControl_0b,help={"Toggle Log/Lin Q display"},value= 1
823	TitleBox DSMControl_0c,pos={120,37},size={104,19},font="Courier",fSize=10
824	TitleBox DSMControl_0c,variable= $(USANSFolder+":DSM:gStr1")
825	//second message string not used currently
826	// TitleBox DSMControl_0d,pos={120,57},size={104,19},font="Courier",fSize=10
827	// TitleBox DSMControl_0d,variable= root:DSM:gStr2
828
829	//tab(1) Mask
830	Button DSMControl_1a,pos={20,35},size={90,20},proc=DSM_MyMaskProc,title="Mask Point" //bMask
831	Button DSMControl_1a,help={"Toggles the masking of the selected data point"}
832	Button DSMControl_1a,disable=1
833	Button DSMControl_1b,pos={20,65},size={140,20},proc=DSM_MaskGTCursor,title="Mask Q >= Cursor" //bMask
834	Button DSMControl_1b,help={"Toggles the masking of all q-values GREATER than the current cursor location"}
835	Button DSMControl_1b,disable=1
836	Button DSMControl_1c,pos={20,95},size={140,20},proc=DSM_MaskLTCursor,title="Mask Q <= Cursor" //bMask
837	Button DSMControl_1c,help={"Toggles the masking of all q-values LESS than the current cursor location"}
838	Button DSMControl_1c,disable=1
839	Button DSMControl_1d,pos={180,35},size={90,20},proc=DSM_ClearMaskProc,title="Clear Mask" //bMask
840	Button DSMControl_1d,help={"Clears all mask points"}
841	Button DSMControl_1d,disable=1
842	// Button DSMControl_1b,pos={144,66},size={110,20},proc=DSM_MaskDoneButton,title="Done Masking"
843	// Button DSMControl_1b,disable=1
844
845	//tab(2) Extrapolate
846	Button DSMControl_2a,pos={31,42},size={90,20},proc=DSM_ExtrapolateButtonProc,title="Extrapolate"
847	Button DSMControl_2a,help={"Extrapolate the high-q region with a power-law"}
848	Button DSMControl_2a,disable=1
849	SetVariable DSMControl_2b,pos={31,70},size={100,15},title="# of points"
850	SetVariable DSMControl_2b,help={"Set the number of points for the power-law extrapolation"}
851	SetVariable DSMControl_2b,limits={5,100,1},value= $(USANSFolder+":DSM:gNptsExtrap")
852	SetVariable DSMControl_2b,disable=1
853	CheckBox DSMControl_2c,pos={157,45},size={105,14},proc=DSM_ExtrapolationCheckProc,title="Show Extrapolation"
854	CheckBox DSMControl_2c,help={"Show or hide the high q extrapolation"},value= 1
855	CheckBox DSMControl_2c,disable=1
856	SetVariable DSMControl_2d,pos={31,96},size={150,15},title="Power Law Exponent"
857	SetVariable DSMControl_2d,help={"Power Law exponent from the fit = the DESMEARED slope - override as needed"}
858	SetVariable DSMControl_2d format="%5.2f"
859	SetVariable DSMControl_2d,limits={-inf,inf,0},value= $(USANSFolder+":DSM:gPowerM")
860	SetVariable DSMControl_2d,disable=1
861
862	//tab(3) Smooth
863	Button DSMControl_3a,pos={34,97},size={70,20},proc=DSM_SmoothButtonProc,title="Smooth"
864	Button DSMControl_3a,disable=1
865	//BoxCheck
866	CheckBox DSMControl_3b,pos={34,39},size={35,14},title="Box",value= 1
867	CheckBox DSMControl_3b,help={"Use a single pass of 3-point box smoothing"}
868	CheckBox DSMControl_3b,disable=1
869	//SSCheck
870	CheckBox DSMControl_3c,pos={34,60},size={45,14},title="Spline",value= 0
871	CheckBox DSMControl_3c,help={"Use a single pass of a smoothing spline"}
872	CheckBox DSMControl_3c,disable=1
873	//extendCheck
874	CheckBox DSMControl_3d,pos={268,60},size={71,14},title="Extend Data"
875	CheckBox DSMControl_3d,help={"extends the data at both low q and high q to avoid end effects in smoothing"}
876	CheckBox DSMControl_3d,value= 0
877	CheckBox DSMControl_3d,disable=1
878	Button DSMControl_3e,pos={125,97},size={90,20},proc=DSM_SmoothUndoButtonProc,title="Start Over"
879	Button DSMControl_3e,help={"Start the smoothing over again without needing to re-mask the data set"}
880	Button DSMControl_3e,disable=1
881	SetVariable DSMControl_3f,pos={94,60},size={150,15},title="Smoothing factor"
882	SetVariable DSMControl_3f,help={"Smoothing factor for the smoothing spline"}
883	SetVariable DSMControl_3f format="%5.4f"
884	SetVariable DSMControl_3f,limits={0.01,2,0.01},value= $(USANSFolder+":DSM:gSmoothFac")
885	SetVariable DSMControl_3f,disable=1
886	CheckBox DSMControl_3g,pos={268,39},size={90,14},title="Log-scale smoothing?"
887	CheckBox DSMControl_3g,help={"Use log-scaled intensity during smoothing (reverts to linear if negative intensity points found)"}
888	CheckBox DSMControl_3g,value=0
889	CheckBox DSMControl_3g,disable=1
890	ValDisplay DSMControl_3h pos={235,97},title="Chi^2",size={80,20},value=root:Packages:NIST:USANS:DSM:gChi2Smooth
891	ValDisplay DSMControl_3h,help={"This is the Chi^2 value for the smoothed data vs experimental data"}
892	ValDisplay DSMControl_3h,disable=1
893
894	//tab(4) Desmear
895	Button DSMControl_4a,pos={35,93},size={80,20},proc=DSM_DesmearButtonProc,title="Desmear"
896	Button DSMControl_4a,help={"Do the desmearing - the result is in I_dsm"}
897	Button DSMControl_4a,disable=1
898	SetVariable DSMControl_4b,pos={35,63},size={120,15},title="Chi^2 target"
899	SetVariable DSMControl_4b,help={"Set the targetchi^2 for convergence (recommend chi^2=1)"}
900	SetVariable DSMControl_4b,limits={0,inf,0.1},value= $(USANSFolder+":DSM:gChi2Target")
901	SetVariable DSMControl_4b,disable=1
902	SetVariable DSMControl_4c,pos={35,35},size={80,15},title="dQv"
903	SetVariable DSMControl_4c,help={"Slit height as read in from the data file. 0.117 is the NIST value, override if necessary"}
904	SetVariable DSMControl_4c,limits={-inf,inf,0},value= $(USANSFolder+":DSM:gDqv")
905	SetVariable DSMControl_4c,disable=1
906	TitleBox DSMControl_4d,pos={160,37},size={104,19},font="Courier",fSize=10
907	TitleBox DSMControl_4d,variable= $(USANSFolder+":DSM:gIterStr")
908	TitleBox DSMControl_4d,disable=1
909
910
911	SetDataFolder root:
912	EndMacro
913
914	// function to control the drawing of buttons in the TabControl on the main panel
915	// Naming scheme for the buttons MUST be strictly adhered to... else buttons will
916	// appear in odd places...
917	// all buttons are named DSMControl_NA where N is the tab number and A is the letter denoting
918	// the button's position on that particular tab.
919	// in this way, buttons will always be drawn correctly :-)
920	//
921	Function DSM_TabProc(ctrlName,tab) //: TabControl
922	String ctrlName
923	Variable tab
924
925	String ctrlList = ControlNameList("",";"),item="",nameStr=""
926	Variable num = ItemsinList(ctrlList,";"),ii,onTab
927	for(ii=0;ii<num;ii+=1)
928	//items all start w/"DSMControl_" //11 characters
929	item=StringFromList(ii, ctrlList ,";")
930	nameStr=item[0,10]
931	if(cmpstr(nameStr,"DSMControl_")==0)
932	onTab = str2num(item[11]) //12th is a number
933	ControlInfo $item
934	switch(abs(V_flag))
935	case 1:
936	Button $item,disable=(tab!=onTab)
937	break
938	case 2:
939	CheckBox $item,disable=(tab!=onTab)
940	break
941	case 5:
942	SetVariable $item,disable=(tab!=onTab)
943	break
944	case 10:
945	TitleBox $item,disable=(tab!=onTab)
946	break
947	case 4:
948	ValDisplay $item,disable=(tab!=onTab)
949	break
950	// add more items to the switch if different control types are used
951	endswitch
952
953	endif
954	endfor
955
956	if(tab==1)
957	DSM_MyMaskProc("") //start maksing if you click on the tab
958	else
959	DSM_MaskDoneButton("") //masking is done if you click off the tab
960	endif
961
962	if(tab == 2)
963	//calculate the extrapolation when the tab is selected - this re-fits the data, what we want to avoid
964	// DSM_ExtrapolateButtonProc("")
965
966	// OR
967	// use the coefficients from when it was loaded
968	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
969	WAVE P_coef = $(USANSFolder+":DSM:P_coef")
970	wave Qw = $(USANSFolder+":DSM:Q_exp")
971
972	DSM_SetExtrWaves(Qw)
973	Wave extr_hqi=$(USANSFolder+":DSM:extr_hqi")
974	Wave extr_hqq=$(USANSFolder+":DSM:extr_hqq")
975	extr_hqi=P_coef[0]+P_coef[1]*extr_hqq^P_coef[2]
976
977	AppendExtrapolation()
978
979	endif
980
981	return 0
982	End
983
984	Proc AppendSmeared()
985	String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
986
987	SetDataFolder $(USANSFolder+":DSM")
988	// if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_exp_orig",0,2) == -1) //Igor 5
989	if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_exp_orig",0) == -1)
990	AppendToGraph/W=Desmear_Graph I_exp_orig vs Q_exp_orig
991	ModifyGraph mode=4,marker=19 //3 is markers, 4 is markers and lines
992	ModifyGraph rgb(I_exp_orig)=(0,0,0)
993	ModifyGraph msize=2,grid=1,log=1,mirror=2,standoff=0,tickunit=1
994	ErrorBars/T=0 I_exp_orig Y,wave=(S_exp_orig,S_exp_orig)
995	Legend/N=text0/J "\\F'Courier'\\s(I_exp_orig) I_exp_orig"
996	Label left "Intensity"
997	Label bottom "Q (1/A)"
998	endif
999	//always update the textbox - kill the old one first
1000	TextBox/K/N=text1
1001	// TextBox/C/N=text1/F=0/A=MT/E=2/X=5.50/Y=0.00 root:DSM:gCurFile //Igor 5
1002	TextBox/C/N=text1/F=0/A=MT/E=1/X=5.50/Y=0.00 $(USANSFolder+":DSM:gCurFile")
1003	End
1004
1005	Proc AppendMask()
1006	String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1007
1008	// if( strsearch(TraceNameList("Desmear_Graph", "", 1),"MaskData",0,2) == -1) //Igor 5
1009	if( strsearch(TraceNameList("Desmear_Graph", "", 1),"MaskData",0) == -1)
1010	SetDataFolder $(USANSFolder+":DSM:")
1011	AppendToGraph/W=Desmear_Graph MaskData vs Q_exp_orig
1012	ModifyGraph mode(MaskData)=3,marker(MaskData)=8,msize(MaskData)=2.5,opaque(MaskData)=1
1013	ModifyGraph rgb(MaskData)=(65535,16385,16385)
1014	setdatafolder root:
1015	endif
1016	end
1017
1018	Proc AppendSmoothed()
1019	String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1020
1021	// if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_smth",0,2) == -1) //Igor 5
1022	if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_smth",0) == -1)
1023	SetDataFolder $(USANSFolder+":DSM:")
1024	AppendToGraph/W=Desmear_Graph I_smth vs Q_smth
1025	ModifyGraph/W=Desmear_Graph rgb(I_smth)=(3,52428,1),lsize(I_smth)=2
1026	setdatafolder root:
1027	endif
1028	end
1029
1030	Function RemoveSmoothed()
1031	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1032
1033	SetDataFolder $(USANSFolder+":DSM:")
1034	RemoveFromGraph/W=Desmear_Graph/Z I_smth
1035	setdatafolder root:
1036	end
1037
1038	Function RemoveMask()
1039	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1040
1041	SetDataFolder $(USANSFolder+":DSM:")
1042	RemoveFromGraph/W=Desmear_Graph/Z MaskData
1043	setdatafolder root:
1044	end
1045
1046	Proc AppendDesmeared()
1047	String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1048
1049	// if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_dsm",0,2) == -1) //Igor 5
1050	if( strsearch(TraceNameList("Desmear_Graph", "", 1),"I_dsm",0) == -1)
1051	SetDataFolder $(USANSFolder+":DSM:")
1052	AppendToGraph/W=Desmear_Graph I_dsm vs Q_dsm
1053	ModifyGraph mode(I_dsm)=3,marker(I_dsm)=19
1054	ModifyGraph rgb(I_dsm)=(1,16019,65535),msize(I_dsm)=2
1055	ErrorBars/T=0 I_dsm Y,wave=(S_dsm,S_dsm)
1056	setdatafolder root:
1057	endif
1058	end
1059
1060	Function RemoveDesmeared()
1061	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1062
1063	SetDataFolder $(USANSFolder+":DSM:")
1064	RemoveFromGraph/W=Desmear_Graph/Z I_dsm
1065	setdatafolder root:
1066	end
1067
1068	Function AppendExtrapolation()
1069	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1070
1071	// if( strsearch(TraceNameList("Desmear_Graph", "", 1),"extr_hqi",0,2) == -1) //Igor 5
1072	if( strsearch(TraceNameList("Desmear_Graph", "", 1),"extr_hqi",0) == -1)
1073	SetDataFolder $(USANSFolder+":DSM:")
1074	AppendToGraph/W=Desmear_Graph extr_hqi vs extr_hqq
1075	ModifyGraph/W=Desmear_Graph lSize(extr_hqi)=2
1076	setdatafolder root:
1077	endif
1078	end
1079
1080	Function RemoveExtrapolation()
1081	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1082
1083	SetDataFolder $(USANSFolder+":DSM:")
1084	RemoveFromGraph/W=Desmear_Graph/Z extr_hqi
1085	setdatafolder root:
1086	end
1087
1088	// step (1) - read in the data, and plot it
1089	// clear out all of the "old" waves, remove them from the graph first
1090	// reads in a fresh copy of the data
1091	//
1092	// produces Q_exp, I_exp, S_exp waves (and originals "_orig")
1093	// add a dummy wave note that can be changed on later steps
1094	//
1095	Function DSM_LoadButtonProc(ctrlName) : ButtonControl
1096	String ctrlName
1097
1098	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1099
1100	String qStr,iStr,sStr,sqStr
1101	Variable nq,dqv,numBad,val
1102
1103	// remove any of the old traces on the graph and delete the waves
1104	CleanUpJunk()
1105
1106	SetDataFolder root:
1107	Execute "A_LoadOneDDataWithName(\"\",0)"
1108	SVAR fileStr = root:Packages:NIST:gLastFileName
1109	if (cmpstr(fileStr,"") == 0)
1110	return(0) //get out if no file selected
1111	endif
1112	//define the waves that the smoothing will be looking for...
1113	SVAR fname = $("root:Packages:NIST:gLastFileName") //this changes as any data is loaded
1114	SVAR curFile = $(USANSFolder+":DSM:gCurFile") //keep this for title, save
1115	curFile = fname
1116
1117	qStr = CleanupName((fName + "_q"),0) //the q-wave
1118	iStr = CleanupName((fName + "_i"),0) //the i-wave
1119	sStr = CleanupName((fName + "_s"),0) //the s-wave
1120	sqStr = CleanupName((fName + "sq"),0) //the sq-wave, which should have -dQv as the elements
1121
1122	String DFStr= CleanupName(fname,0)
1123
1124	Duplicate/O $("root:"+DFStr+":"+qStr) $(USANSFolder+":DSM:Q_exp")
1125	Duplicate/O $("root:"+DFStr+":"+iStr) $(USANSFolder+":DSM:I_exp")
1126	Duplicate/O $("root:"+DFStr+":"+sStr) $(USANSFolder+":DSM:S_exp")
1127	wave Q_exp = $(USANSFolder+":DSM:Q_exp")
1128	Wave I_exp = $(USANSFolder+":DSM:I_exp")
1129	Wave S_exp = $(USANSFolder+":DSM:S_exp")
1130
1131	// copy over the high q extrapolation information
1132	Duplicate/O $("root:"+DFStr+":P_coef") $(USANSFolder+":DSM:P_coef")
1133	NVAR slope = $("root:"+DFStr+":USANS_m")
1134	NVAR powerM = $(USANSFolder+":DSM:gPowerM")
1135	powerM = slope
1136
1137	// remove any negative q-values (and q=0 values!)(and report this)
1138	// ? and trim the low q to be >= 3.0e-5 (1/A), below this USANS is not reliable.
1139	NumBad = RemoveBadQPoints(Q_exp,I_exp,S_exp,0)
1140	SVAR str1 = $(USANSFolder+":DSM:gStr1")
1141	sprintf str1,"%d negative q-values were removed",numBad
1142
1143	// don't trim off any positive q-values
1144	// val = 3e-5 //lowest "good" q-value from USANS
1145	// NumBad = RemoveBadQPoints(Q_exp,I_exp,S_exp,val-1e-8)
1146	// SVAR str2 = root:DSM:gStr2
1147	// sprintf str2,"%d q-values below q = %g were removed",numBad,val
1148
1149	Duplicate/O $(USANSFolder+":DSM:Q_exp") $(USANSFolder+":DSM:Q_exp_orig")
1150	Duplicate/O $(USANSFolder+":DSM:I_exp") $(USANSFolder+":DSM:I_exp_orig")
1151	Duplicate/O $(USANSFolder+":DSM:S_exp") $(USANSFolder+":DSM:S_exp_orig")
1152	wave I_exp_orig = $(USANSFolder+":DSM:I_exp_orig")
1153
1154	nq = numpnts($(USANSFolder+":DSM:Q_exp"))
1155
1156	dQv = NumVarOrDefault("root:"+DFStr+":USANS_dQv", 0.117 )
1157	// if(WaveExists(sigQ)) //try to read dQv
1158	//// dqv = -sigQ[0][0]
1159	//// DoAlert 0,"Found dQv value of " + num2str(dqv)
1160	// else
1161	// dqv = 0.117
1162	// // dqv = 0.037 //old value
1163	// DoAlert 0,"Could not find dQv in the data file - using " + num2str(dqv)
1164	// endif
1165	NVAR gDqv = $(USANSFolder+":DSM:gDqv") //needs to be global for Weights_L()
1166	NVAR gNq = $(USANSFolder+":DSM:gNq")
1167	//reset the globals
1168	gDqv = dqv
1169	gNq = nq
1170
1171	// append the (blank) wave note to the intensity wave
1172	Note I_exp,"BOX=0;SPLINE=0;"
1173	Note I_exp_orig,"BOX=0;SPLINE=0;"
1174
1175	//draw the graph
1176	Execute "AppendSmeared()"
1177
1178	SetDataFolder root:
1179	End
1180
1181	// remove any q-values <= val
1182	Function RemoveBadQPoints(qw,iw,sw,val)
1183	Wave qw,iw,sw
1184	Variable val
1185
1186	Variable ii,num,numBad,qval
1187	num = numpnts(qw)
1188
1189	ii=0
1190	numBad=0
1191	do
1192	qval = qw[ii]
1193	if(qval <= val)
1194	numBad += 1
1195	else //keep the points
1196	qw[ii-numBad] = qval
1197	iw[ii-numBad] = iw[ii]
1198	sw[ii-numBad] = sw[ii]
1199	endif
1200	ii += 1
1201	while(ii<num)
1202	//trim the end of the waves
1203	DeletePoints num-numBad, numBad, qw,iw,sw
1204	return(numBad)
1205	end
1206
1207	// if mw = Nan, keep the point, if a numerical value, delete it
1208	Function RemoveMaskedPoints(mw,qw,iw,sw)
1209	Wave mw,qw,iw,sw
1210
1211	Variable ii,num,numBad,mask
1212	num = numpnts(qw)
1213
1214	ii=0
1215	numBad=0
1216	do
1217	mask = mw[ii]
1218	if(numtype(mask) != 2) //if not NaN
1219	numBad += 1
1220	else //keep the points that are NaN
1221	qw[ii-numBad] = qw[ii]
1222	iw[ii-numBad] = iw[ii]
1223	sw[ii-numBad] = sw[ii]
1224	endif
1225	ii += 1
1226	while(ii<num)
1227	//trim the end of the waves
1228	DeletePoints num-numBad, numBad, qw,iw,sw
1229	return(numBad)
1230	end
1231
1232	// produces the _msk waves that have the new number of data points
1233	//
1234	Function DSM_MaskDoneButton(ctrlName) : ButtonControl
1235	String ctrlName
1236
1237	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1238
1239	// Variable aExists= strlen(CsrInfo(A)) > 0 //Igor 5
1240	Variable aExists= strlen(CsrWave(A)) > 0 //Igor 4
1241	if(!aExists)
1242	return(1) //possibly reverted data, no cursor, no Mask wave
1243	endif
1244
1245	Duplicate/O $(USANSFolder+":DSM:Q_exp_orig"),$(USANSFolder+":DSM:Q_msk")
1246	Duplicate/O $(USANSFolder+":DSM:I_exp_orig"),$(USANSFolder+":DSM:I_msk")
1247	Duplicate/O $(USANSFolder+":DSM:S_exp_orig"),$(USANSFolder+":DSM:S_msk")
1248	Wave Q_msk=$(USANSFolder+":DSM:Q_msk")
1249	Wave I_msk=$(USANSFolder+":DSM:I_msk")
1250	Wave S_msk=$(USANSFolder+":DSM:S_msk")
1251
1252	//finish up - trim the data sets and reassign the working set
1253	Wave MaskData=$(USANSFolder+":DSM:MaskData")
1254
1255	RemoveMaskedPoints(MaskData,Q_msk,I_msk,S_msk)
1256
1257	//reset the number of points
1258	NVAR gNq = $(USANSFolder+":DSM:gNq")
1259	gNq = numpnts(Q_msk)
1260
1261	Cursor/K A
1262	HideInfo
1263
1264	return(0)
1265	End
1266
1267
1268	// not quite the same as revert
1269	Function DSM_ClearMaskProc(ctrlName) : ButtonControl
1270	String ctrlName
1271
1272	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1273
1274	Wave MaskData=$(USANSFolder+":DSM:MaskData")
1275	MaskData = NaN
1276
1277	return(0)
1278	end
1279
1280	// when the mask button is pressed, A must be on the graph
1281	// Displays MaskData wave on the graph
1282	//
1283	Function DSM_MyMaskProc(ctrlName) : ButtonControl
1284	String ctrlName
1285
1286	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1287
1288	Wave data=$(USANSFolder+":DSM:I_exp_orig")
1289
1290	// Variable aExists= strlen(CsrInfo(A)) > 0 //Igor 5
1291	Variable aExists= strlen(CsrWave(A)) > 0 //Igor 4
1292
1293	// need to get rid of old smoothed data if data is re-masked
1294	Execute "RemoveSmoothed()"
1295	SetDataFolder $(USANSFolder+":DSM")
1296	Killwaves/Z I_smth,Q_smth,S_smth
1297
1298	if(aExists) //mask the selected point
1299	// toggle NaN (keep) or Data value (= masked)
1300	Wave MaskData
1301	MaskData[pcsr(A)] = (numType(MaskData[pcsr(A)])==0) ? NaN : data[pcsr(A)] //if NaN, doesn't plot
1302	else
1303	Cursor /A=1/H=1/L=1/P/W=Desmear_Graph A I_exp_orig leftx(I_exp_orig)
1304	ShowInfo
1305	//if the mask wave does not exist, make one
1306	if(exists("MaskData") != 1)
1307	Duplicate/O Q_exp_orig MaskData
1308	MaskData = NaN //use all data
1309	endif
1310	Execute "AppendMask()"
1311	endif
1312
1313	SetDataFolder root:
1314
1315	return(0)
1316	End
1317
1318	// when the mask button is pressed, A must be on the graph
1319	// Displays MaskData wave on the graph
1320	//
1321	Function DSM_MaskLTCursor(ctrlName) : ButtonControl
1322	String ctrlName
1323
1324	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1325
1326	// Variable aExists= strlen(CsrInfo(A)) > 0 //Igor 5
1327	Variable aExists= strlen(CsrWave(A)) > 0 //Igor 4
1328
1329	if(!aExists)
1330	return(1)
1331	endif
1332	// need to get rid of old smoothed data if data is re-masked
1333	Execute "RemoveSmoothed()"
1334	SetDataFolder $(USANSFolder+":DSM")
1335	Killwaves/Z I_smth,Q_smth,S_smth
1336
1337	Wave data=I_exp_orig
1338
1339	Variable pt,ii
1340	pt = pcsr(A)
1341	for(ii=pt;ii>=0;ii-=1)
1342	// toggle NaN (keep) or Data value (= masked)
1343	Wave MaskData
1344	MaskData[ii] = (numType(MaskData[ii])==0) ? NaN : data[ii] //if NaN, doesn't plot
1345	endfor
1346
1347	SetDataFolder root:
1348	return(0)
1349	End
1350
1351	// when the mask button is pressed, A must be on the graph
1352	// Displays MaskData wave on the graph
1353	//
1354	Function DSM_MaskGTCursor(ctrlName) : ButtonControl
1355	String ctrlName
1356
1357	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1358
1359
1360	// Variable aExists= strlen(CsrInfo(A)) > 0 //Igor 5
1361	Variable aExists= strlen(CsrWave(A)) > 0 //Igor 4
1362
1363	if(!aExists)
1364	return(1)
1365	endif
1366	// need to get rid of old smoothed data if data is re-masked
1367	Execute "RemoveSmoothed()"
1368	SetDataFolder $(USANSFolder+":DSM")
1369	Killwaves/Z I_smth,Q_smth,S_smth
1370
1371	Wave data=I_exp_orig
1372
1373	Variable pt,ii,endPt
1374	endPt=numpnts(MaskData)
1375	pt = pcsr(A)
1376	for(ii=pt;ii<endPt;ii+=1)
1377	// toggle NaN (keep) or Data value (= masked)
1378	Wave MaskData
1379	MaskData[ii] = (numType(MaskData[ii])==0) ? NaN : data[ii] //if NaN, doesn't plot
1380	endfor
1381
1382	SetDataFolder root:
1383
1384	return(0)
1385	End
1386
1387	Function CleanUpJunk()
1388	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1389
1390	// clean up the old junk on the graph, /Z for no error
1391	Execute "RemoveExtrapolation()"
1392	Execute "RemoveDesmeared()"
1393	Execute "RemoveSmoothed()"
1394	Execute "RemoveMask()"
1395
1396	//remove the cursor
1397	Cursor/K A
1398
1399	//always initialize these
1400	String/G $(USANSFolder+":DSM:gStr1") = ""
1401	String/G $(USANSFolder+":DSM:gStr2") = ""
1402	String/G $(USANSFolder+":DSM:gIterStr") = ""
1403
1404	// clean up the old waves from smoothing and desmearing steps
1405	SetDataFolder $(USANSFolder+":DSM")
1406	Killwaves/Z I_smth,I_dsm,I_dsm_sm,Q_smth,Q_dsm,S_smth,S_dsm,Yi_SS,Yq_SS
1407	Killwaves/Z Weights,FW,R_wave,S_wave,H_wave,Di,Ci
1408	Killwaves/Z Is_old,I_old,err
1409	Killwaves/Z Q_ext,I_ext,S_ext,hqq,hqi,lqq,lqi
1410	Killwaves/Z MaskData,Q_msk,I_msk,S_msk
1411	Killwaves/Z Q_work,I_work,S_work //working waves for desmearing step
1412	SetDataFolder root:
1413	End
1414
1415	// does not alter the data sets - reports a power law
1416	// exponent and makes it global so it will automatically
1417	// be used during the desmearing
1418	//
1419	// generates extr_hqi vs extr_hqq that are Appended to the graph
1420	Function DSM_ExtrapolateButtonProc(ctrlName) : ButtonControl
1421	String ctrlName
1422
1423	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1424
1425	NVAR nend = $(USANSFolder+":DSM:gNptsExtrap")
1426	NVAR m_pred = $(USANSFolder+":DSM:gPowerM")
1427
1428	SetDataFolder $(USANSFolder+":DSM")
1429	//use masked data if it exists
1430	if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1)
1431	wave Qw = $(USANSFolder+":DSM:Q_msk")
1432	Wave Iw = $(USANSFolder+":DSM:I_msk")
1433	Wave Sw = $(USANSFolder+":DSM:S_msk")
1434	else
1435	//start from the "_exp" waves
1436	if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1)
1437	wave Qw = $(USANSFolder+":DSM:Q_exp")
1438	Wave Iw = $(USANSFolder+":DSM:I_exp")
1439	Wave Sw = $(USANSFolder+":DSM:S_exp")
1440	endif
1441	endif
1442	SetDataFolder root:
1443
1444	DSM_SetExtrWaves(Qw)
1445	m_pred = DSM_DoExtrapolate(Qw,Iw,Sw,nend)
1446	AppendExtrapolation()
1447
1448	return(0)
1449	End
1450
1451	//smooths the data in steps as requested...
1452	//
1453	// typically do a simple Box smooth first,
1454	// then do a smoothing spline, keeping the same number of data points
1455	//
1456	// chi-squared is reported - so you can see how "bad" the smoothing is
1457	// smoothing of the data is done on a log(I) scale if requested
1458	// setting doLog variable to 0 will return to linear smoothing
1459	// (I see little difference)
1460	//
1461	// start from the "_smth" waves if they exist
1462	// otheriwse start from the working waves
1463	//
1464	// create Q_smth, I_smth, S_smth
1465	// keep track of the smoothing types and passes in the wave note
1466	//
1467	Function DSM_SmoothButtonProc(ctrlName) : ButtonControl
1468	String ctrlName
1469
1470	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1471
1472	SetDataFolder $(USANSFolder+":DSM")
1473
1474	Variable ii,new_n,pass,nq_ext,offset,doLog=1
1475	String noteStr=""
1476
1477	// want log scaling of intensity during smoothing?
1478	ControlInfo DSMControl_3g
1479	doLog = V_value
1480
1481	//// if(exists("I_smth")==1 && exists("Q_smth")==1 && exists("S_smth")==1 && ! freshMask)
1482	if(exists("I_smth")==1 && exists("Q_smth")==1 && exists("S_smth")==1)
1483	//start from the smoothed waves --- just need the wave references
1484	else
1485	//start from the "msk", creating smth waves
1486	if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1)
1487	wave Q_msk,I_msk,S_msk
1488	Duplicate/O I_msk,I_smth,Q_smth,S_smth
1489	I_smth = I_msk
1490	Q_smth = Q_msk
1491	S_smth = S_msk
1492	else
1493	//start from the "_exp" waves
1494	if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1)
1495	wave Q_exp,I_exp,S_exp
1496	Duplicate/O I_exp,I_smth,Q_smth,S_smth
1497	I_smth = I_exp
1498	Q_smth = Q_exp
1499	S_smth = S_exp
1500	endif
1501	endif
1502	endif
1503
1504	wave Q_smth,I_smth,S_smth
1505
1506	// extend the data to avoid end effects
1507	//creates I_ext,Q_ext,S_ext with extended q-values
1508	// fit 15 pts at each end, typically add 10 pts to each end
1509	// does nothing if num_extr = 0
1510	ControlInfo/W=Desmear_Graph DSMControl_3d //extendCheck
1511	if(V_value == 1)
1512	ExtendToSmooth(Q_smth,I_smth,S_smth,15,15,10)
1513	else
1514	ExtendToSmooth(Q_smth,I_smth,S_smth,15,15,0) //don't extend, just use the original data
1515	endif
1516
1517	//whether extending or not, the working data is "_ext", set by ExtendToSmooth()
1518	SetDataFolder $(USANSFolder+":DSM")
1519	wave Q_ext,I_ext ,S_ext
1520
1521	noteStr=note(I_smth)
1522	Note I_ext , noteStr
1523
1524	WaveStats/Q I_ext
1525	if(V_min<=0)
1526	Print "negative itensity found, using linear scale for smoothing"
1527	doLog = 0
1528	endif
1529
1530	if(doLog)
1531	//convert to log scale
1532	Duplicate/O I_ext I_log,I_log_err
1533	Wave I_log ,I_log_err
1534	I_log = log(I_ext)
1535	WaveStats/Q I_log
1536	offset = 2*(-V_min)
1537	I_log += offset
1538	I_log_err = S_ext/(2.30*I_ext)
1539	I_ext = I_log
1540	endif
1541	//
1542
1543	ControlInfo/W=Desmear_Graph DSMControl_3b //BoxCheck
1544	if(V_value == 1)
1545	//do a simple Box smooth first - number of points does not change
1546	// fills ends with neighboring value (E=3) (n=3 points in smoothing window)
1547	Smooth/E=3/B 3, I_ext
1548
1549	noteStr=note(I_ext)
1550	pass = NumberByKey("BOX", noteStr, "=", ";")
1551	noteStr = ReplaceNumberByKey("BOX", noteStr, pass+1, "=", ";")
1552	Note/K I_ext
1553	Note I_ext , noteStr
1554	// Print "box = ",noteStr
1555	endif
1556
1557	NVAR sParam = gSmoothFac //already in the right DF
1558
1559	ControlInfo/W=Desmear_Graph DSMControl_3c //SSCheck
1560	if(V_value == 1)
1561	// Interpolate2/T=3/N=(nq)/I=2/F=1/SWAV=S_ext/Y=Yi_SS/X=Yq_SS Q_ext, I_ext //Igor 5
1562	// 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
1563	//Igor 4
1564	String str=""
1565	nq_ext = numpnts(Q_ext)
1566	str = "Interpolate/T=3/N=("+num2str(nq_ext)+")/I=1/F=("+num2str(sParam)+")/Y=Yi_SS/X=Yq_SS I_ext /X=Q_ext"
1567	Execute str
1568	// Print Str
1569	// Interpolate2/T=3/N=(nq_ext)/I=1/F=(sParam)/Y=Yi_SS/X=Yq_SS Q_ext, I_ext
1570	// end Igor 4
1571	wave yi_ss = yi_ss // already in the right DF
1572	wave yq_ss = yq_ss
1573	// reassign the "working" waves with the result of interpolate, which has the same I/Q values
1574	I_ext = yi_ss
1575	Q_ext = yq_ss
1576
1577	noteStr=note(I_ext)
1578	pass = NumberByKey("SPLINE", noteStr, "=", ";")
1579	noteStr = ReplaceNumberByKey("SPLINE", noteStr, pass+1, "=", ";")
1580	Note/K I_ext
1581	Note I_ext , noteStr
1582	// Print "spline = ",noteStr
1583	endif
1584
1585	//undo the scaling
1586	If(doLog)
1587	I_ext -= offset
1588	I_ext = 10^I_ext
1589	endif
1590
1591	// at this point, I_ext has too many points - and we need to just return the
1592	// center chunk that is the original q-values
1593	// as assign this to the working set for the desmear step
1594	// and to the _smth set in case another smoothng pass is desired
1595	// Q_smth has not changed, S_smth has not changed
1596	I_smth = interp(Q_smth[p], Q_ext, I_ext )
1597
1598	Note/K I_smth
1599	Note I_smth , noteStr
1600	// Print "end of smoothed, note = ",note(I_smth)
1601
1602	Variable nq = numpnts($(USANSFolder+":DSM:Q_smth"))
1603	// Print "nq after smoothing = ",nq
1604
1605	//reset the global
1606	NVAR gNq = gNq
1607	gNq = nq
1608	//report the chi^2 difference between the smoothed curve and the experimental data
1609	NVAR chi2 = gChi2Smooth
1610	chi2 = SmoothedChi2(I_smth)
1611
1612	Execute "AppendSmoothed()"
1613
1614	setDataFolder root:
1615	return(0)
1616	End
1617
1618	// undo the smoothing and start over - useful if you've smoothed
1619	// too aggressively and need to back off
1620	Function DSM_SmoothUndoButtonProc(ctrlName) : ButtonControl
1621	String ctrlName
1622
1623	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1624
1625	Execute "RemoveSmoothed()"
1626	SetDataFolder $(USANSFolder+":DSM")
1627	Killwaves/Z I_smth,Q_smth,S_smth,Q_ext,I_ext,S_ext,Yi_SS,Yq_SS
1628	SetDataFolder root:
1629	return (0)
1630	end
1631
1632	//calculate the chi^2 value for the smoothed data
1633	Function SmoothedChi2(I_smth)
1634	Wave I_smth
1635
1636	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1637
1638	//start from the "msk", if they exist
1639	if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1)
1640	Wave iw = $(USANSFolder+":DSM:I_msk")
1641	Wave sw = $(USANSFolder+":DSM:S_msk")
1642	else
1643	//start from the "_exp" waves
1644	if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1)
1645	Wave iw = $(USANSFolder+":DSM:I_exp")
1646	Wave sw = $(USANSFolder+":DSM:S_exp")
1647	endif
1648	endif
1649
1650	Variable ii,chi2,num
1651	chi2=0
1652	num = numpnts(iw)
1653	for(ii=0;ii<num;ii+=1)
1654	chi2 += (iw[ii] - I_smth[ii])^2 / (sw[ii])^2
1655	endfor
1656	Chi2 /= (num-1)
1657	return (chi2)
1658	end
1659
1660	// I_dsm is the desmeared data
1661	//
1662	// step (7) - desmearing is done, write out the result
1663	//
1664	Function DSM_SaveButtonProc(ctrlName) : ButtonControl
1665	String ctrlName
1666
1667	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1668	NVAR useXMLOutput = root:Packages:NIST:gXML_Write
1669
1670	String saveStr
1671	SVAR curFile = $(USANSFolder+":DSM:gCurFile")
1672	saveStr = CleanupName((curFile),0) //the output filename
1673	//
1674
1675	if (useXMLOutput == 1)
1676	WriteXMLUSANSDesmeared(saveStr,0,0,1)
1677	else
1678	WriteUSANSDesmeared(saveStr,0,0,1) //use the full set (lo=hi=0) and present a dialog
1679	endif
1680
1681	SetDataFolder root:
1682	return(0)
1683	End
1684
1685	Function DSM_HelpButtonProc(ctrlName) : ButtonControl
1686	String ctrlName
1687
1688	DisplayHelpTopic/Z/K=1 "Desmearing USANS Data"
1689	if(V_flag !=0)
1690	DoAlert 0,"The Desmearing USANS Data Help file could not be found"
1691	endif
1692	return(0)
1693	End
1694
1695
1696	//toggles the log/lin display of the loaded data set
1697	Function DSM_LoadCheckProc(ctrlName,checked) : CheckBoxControl
1698	String ctrlName
1699	Variable checked
1700
1701	ModifyGraph log=(checked)
1702	return(0)
1703	End
1704
1705
1706	Function DSM_ExtrapolationCheckProc(ctrlName,checked) : CheckBoxControl
1707	String ctrlName
1708	Variable checked
1709
1710	if(checked)
1711	AppendExtrapolation()
1712	else
1713	RemoveExtrapolation()
1714	endif
1715	return(0)
1716	End
1717
1718
1719	// takes as input the "working waves"
1720	//
1721	// creates intermediate waves to work on
1722	//
1723	// output of Q_dsm,I_dsm,S_dsm (and I_dsm_sm, the result of smearing I_dsm)
1724	//
1725	Function DSM_DesmearButtonProc(ctrlName) : ButtonControl
1726	String ctrlName
1727
1728	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1729
1730	SetDataFolder $(USANSFolder+":DSM")
1731	if(exists("I_smth")==1 && exists("Q_smth")==1 && exists("S_smth")==1)
1732	wave Q_smth , I_smth ,S_smth
1733	Duplicate/O I_smth,I_work,Q_work,S_work
1734	I_work = I_smth
1735	Q_work = Q_smth
1736	S_work = S_smth
1737	else
1738	//start from the "msk", creating work waves
1739	if(exists("I_msk")==1 && exists("Q_msk")==1 && exists("S_msk")==1)
1740	wave Q_msk,I_msk,S_msk
1741	Duplicate/O I_msk,I_work,Q_work,S_work
1742	I_work = I_msk
1743	Q_work = Q_msk
1744	S_work = S_msk
1745	else
1746	//start from the "_exp" waves
1747	if(exists("I_exp")==1 && exists("Q_exp")==1 && exists("S_exp")==1)
1748	wave Q_exp,I_exp,S_exp
1749	Duplicate/O I_exp,I_work,Q_work,S_work
1750	I_work = I_exp
1751	Q_work = Q_exp
1752	S_work = S_exp
1753	endif
1754	endif
1755	endif
1756	//SetDataFolder root:
1757
1758	NVAR nq = gNq
1759	NVAR m = gPowerM
1760	NVAR chi2_target = gChi2Target
1761	NVAR maxFastIter = gMaxFastIter
1762	NVAR maxSlowIter = gMaxSlowIter
1763
1764	// SET WEIGHTING OF EXPERIMENTAL DATA.
1765	Duplicate/O Q_work weights
1766	Wave weights = weights
1767	weights = 1/S_work^2
1768
1769	// calculate weighting array for smearing of data
1770	Make/O/D/N=(nq,nq) FW
1771	Wave FW
1772	Weights_L(m,FW,Q_work)
1773
1774	// ^^^^ Iterative desmearing ^^^^*
1775	Variable chi2_old,chi2_new,done,iter
1776	// FOR 0TH ITERATION, EXPERIMENTAL DATA IS USED FOR Y_OLD, create ys_old for result
1777	// y_old = I_old, y_new = I_dsm, I_dsm_sm = ys_new,
1778	// duplicate preserves the wave note!
1779	Duplicate/O I_work I_old,Is_old,I_dsm,I_dsm_sm
1780	Duplicate/O Q_work Q_dsm,S_dsm //sets Q_dsm correctly
1781	wave S_dsm,I_old,Is_old,I_dsm,I_dsm_sm
1782	I_old = I_work
1783	Is_old = 0
1784	I_dsm = 0
1785	I_dsm_sm = 0
1786
1787	// Smear 0TH iter guess Y --> YS
1788	DSM_Smear(I_old,Is_old,NQ,FW)
1789	chi2_OLD = chi2(Is_old,I_work,weights,NQ)
1790
1791	printf "starting chi^2 = %g\r",chi2_old
1792	Print "Starting fast convergence..... "
1793
1794	done = 0 //false
1795	iter = 0
1796	Do // while not done, use Fast convergence
1797	I_dsm = I_work * I_old / Is_old // Calculate corrected guess (no need for do-loop)
1798
1799	// Smear iter guess I_dsm --> I_dsm_sm and see how well I_dsm_sm matches experimental data
1800	DSM_Smear(I_dsm,I_dsm_sm,NQ,FW)
1801	chi2_new = chi2(I_dsm_sm,I_work,weights,NQ)
1802
1803	// Stop iteration if fit from new iteration has worse fit...
1804	If (chi2_new > chi2_old)
1805	Done = 1
1806	Endif
1807
1808	// Stop iteration if fit is better than target value...
1809	If (chi2_new < chi2_target)
1810	Done = 1
1811	Else
1812	Iter += 1
1813	Printf "Iteration %d, Chi^2(new) = %g\r", Iter,chi2_new
1814	I_old = I_dsm
1815	Is_old = I_dsm_sm
1816	if(iter>maxFastIter)
1817	break
1818	endif
1819	CHI2_OLD = CHI2_NEW
1820	EndIf
1821	while( !done )
1822
1823	// append I_dsm,I_exp to the graph
1824	Execute "AppendDesmeared()"
1825	DoUpdate
1826
1827	SetDataFolder $(USANSFolder+":DSM")
1828
1829
1830	// step (6) - refine the desmearing using slow convergence
1831	Print "Starting slow convergence..... "
1832	done = 0 // ! reset flag for slow convergence
1833	Do
1834	I_dsm = I_old + (I_work - Is_old) // Calculate corrected guess
1835
1836	// Smear iter guess Y --> YS
1837	DSM_Smear(I_dsm,I_dsm_sm,NQ,FW)
1838	chi2_new = chi2(I_dsm_sm,I_work,weights,NQ)
1839
1840	// Stop iteration if fit from new iteration has worse fit...
1841	If (chi2_new > chi2_old)
1842	Done = 1
1843	Endif
1844
1845	// Stop iteration if fit is better than target value...
1846	If (chi2_new < chi2_target)
1847	Done = 1
1848	Else
1849	Iter += 1
1850
1851	if(mod(iter, 50 ) ==0 )
1852	Printf "Iteration %d, Chi^2(new) = %g\r", Iter,chi2_new
1853	DoUpdate
1854	endif
1855	I_old = I_dsm
1856	Is_old = I_dsm_sm
1857
1858	CHI2_OLD = CHI2_NEW
1859	EndIf
1860	if(iter>maxSlowIter)
1861	break
1862	endif
1863	While ( !done )
1864
1865	Printf "Iteration %d, Chi^2(new) = %g\r", Iter,chi2_new
1866
1867	// adjust the error
1868	SetDataFolder $(USANSFolder+":DSM")
1869	Duplicate/O S_work err
1870	S_dsm = abs(err/I_work*I_dsm) //proportional error
1871	//John simply keeps the same error as was read in from the smeared data set - is this right?
1872	// S_dsm = S_Work
1873
1874	NVAR gChi = gChi2Final //chi^2 final
1875	NVAR gIter = gIterations //total number of iterations
1876	gChi = chi2_new
1877	gIter = Iter
1878
1879	SVAR str = gIterStr
1880	sprintf str,"%d iterations required",iter
1881
1882	SetDataFolder root:
1883	return(0)
1884	End
1885
1886	Function DSM_RevertButtonProc(ctrlName) : ButtonControl
1887	String ctrlName
1888
1889	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1890
1891
1892	CleanUpJunk()
1893
1894	SetDataFolder $(USANSFolder+":DSM")
1895
1896	//reset the working waves to the original
1897	wave Q_exp_orig,I_exp_orig,S_exp_orig
1898
1899	Duplicate/O Q_exp_orig Q_exp
1900	Duplicate/O I_exp_orig I_exp
1901	Duplicate/O S_exp_orig S_exp
1902	// kill the data folder
1903	// re-initialize?
1904
1905	SetDataFolder root:
1906
1907	return(0)
1908	End

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format