Context Navigation

LakeDesmearing_JB.ipf @ 561

Last change on this file since 561 was 496, checked in by srkline, 14 years ago

Two changes:
(1) Changed the behavior of DIV file creation. Now a panel is presented that asks for all of the information - run numbers, transmission, XY box, etc. and then does all of the reduction/replacing/saving with one click. Cedric's ticket #195

(2) Some fixes in USANS notebook generation - filtering out junk.

Property eol-style set to native

File size: 52.7 KB

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

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

Download in other formats:

Original Format