Context Navigation

LakeDesmearing_JB.ipf @ 609

Last change on this file since 609 was 570, checked in by srkline, 13 years ago

Change (1):
In preparation for release, updated pragma IgorVersion?=6.1 in all procedures

Change (2):
As a side benefit of requiring 6.1, we can use the MultiThread? keyword to thread any model function we like. The speed benefit is only noticeable on functions that require at least one integration and at least 100 points (resolution smearing is NOT threaded, too many threadSafe issues, too little benefit). I have chosen to use the MultiThread? only on the XOP assignment. In the Igor code there are too many functions that are not explicitly declared threadsafe, making for a mess.

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.1
4
5	//////////////////////////////////////////////
6	// Igor conversion: 03 FEB 06 SRK
7	// Revision: 03 MAR 06 SRK
8	//
9	// Program uses Lake's iterative technique, J. A. Lake, Acta Cryst. 23 (1967) 191-4.
10	// to DESMEAR Infinite slit smeared USANS DATA.
11	//
12	// TO SEE DESCRIPTION, CHECK COMPUTER SOFTWARE LOGBOOK, P13,41-46, 50
13	// J. BARKER, August, 2001
14	// Switches from fast to slow convergence near target chi JGB 2/2003
15	//
16	// steps are:
17	// load
18	// mask
19	// extrapolate (find the power law of the desmeared set = smeared -1)
20	// smooth (optional)
21	// desmear based on dQv and chi^2 target
22	// save result
23	//
24	/////////
25	// Waves produced at each step: (Q/I/S with the following extensions)
26	//
27	// Load: Uses -> nothing (Kills old waves)
28	// Produces -> "_exp" and "_exp_orig"
29	//
30	// Mask: Uses -> "_exp"
31	// Produces -> "_msk"
32	//
33	// Extrapolate: Uses -> nothing
34	// Produces -> nothing
35	//
36	// Smooth: Uses -> "_smth" OR "_msk" OR "_exp" (in that order)
37	// Produces -> "_smth"
38	//
39	// Desmear: Uses -> "_smth" OR "_msk" OR "_exp" (in that order)
40	// Produces -> "_dsm"
41	//
42	////////
43
44	///*** TO FIX *****
45	// - ? power law + background extrapolation (only useful for X-ray data...)
46	//
47	// see commented code lines for Igor 4 or Igor 5
48	// - there are only a few options and calls that are not Igor 4 compatible.
49	// Igor 4 routines are currently used.
50
51	////////////////////////////////////////////////////////////////////////
52
53	// main entry routine
54	Proc Desmear()
55
56	String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
57
58	//check for the correct data folder, initialize if necessary
59	//
60	if(DataFolderExists(USANSFolder+":DSM") == 0)
61	Execute "Init_Desmearing()"
62	endif
63
64	SetDataFolder $(USANSFolder+":DSM")
65	//always initialize these global variables
66	gStr1 = ""
67	gStr2 = ""
68	gIterStr = ""
69	SetDataFolder root:
70
71	DoWindow/F Desmear_Graph
72	if(V_flag==0)
73	Execute "Desmear_Graph()"
74	endif
75	End
76
77	Proc Init_Desmearing()
78
79	String USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
80
81	//set up the folder(s) needed
82	NewDataFolder/O $(USANSFolder+":DSM")
83	NewDataFolder/O $(USANSFolder+":myGlobals") //in case it wasn't created elsewhere
84
85	SetDataFolder $(USANSFolder+":DSM")
86
87	String/G gCurFile=""
88
89	Variable/G gMaxFastIter = 100 //max number of iter in Fast convergence
90	Variable/G gMaxSlowIter = 10000
91
92	Variable/G gNptsExtrap = 10 //points for high q extrapolation
93	Variable/G gChi2Target = 1 //chi^2 target
94	Variable/G gPowerM = -4
95	Variable/G gDqv = 0.117 //2005 measured slit height - see John
96	Variable/G gNq = 1
97	Variable/G gS = 0 // global varaible for Midpnt()
98	Variable/G gSmoothFac=0.03
99
100	Variable/G gChi2Final = 0 //chi^2 final
101	Variable/G gIterations = 0 //total number of iterations
102
103	String/G gStr1 = "" //information strings
104	String/G gStr2 = ""
105	String/G gIterStr = ""
106	Variable/G gChi2Smooth = 0
107
108	Variable/G gFreshMask=1
109
110	SetDataFolder root:
111	End
112
113	//////////// Lake desmearing routines
114
115	// Smears guess Y --> YS using weighting array
116	Function DSM_Smear(Y_wave,Ys,NQ,FW)
117	Wave Y_wave,Ys
118	Variable nq
119	Wave FW
120
121	Variable ii,jj,summ
122	for(ii=0;ii<nq;ii+=1)
123	summ=0
124	for(jj=0;jj<nq;jj+=1)
125	summ = summ + Y_wave[jj]*FW[ii][jj]
126	endfor
127	Ys[ii] = summ
128	endfor
129
130	Return (0)
131	End
132
133	// CALCULATES CHI^2
134	FUNCTION CHI2(ys_new,ys_exp,w,NQ)
135	Wave ys_new,ys_exp,w
136	Variable nq
137
138	Variable CHI2,summ,ii
139
140	SUMm = 0.0
141	for(ii=0;ii<nq;ii+=1)
142	// if(numtype(YS_EXP[ii]) == 0 && numtype(YS_NEW[ii]) == 0)
143	SUMm=SUMm+W[ii]*(YS_EXP[ii]-YS_NEW[ii])^2
144	// endif
145	endfor
146
147	CHI2=SUMm/(NQ-1)
148
149	RETURN CHI2
150	END
151
152	// Routine calculates the weights needed to convert a table
153	// representation of a scattering function I(q) into the infinite
154	// slit smeared table represented as Is(q)
155	// Is(qi) = Sum(j=i-n) Wij*I(qj)
156	// Program assumes data is extrapolated from qn to dqv using powerlaw
157	// I(q) = Aq^m
158	// Required input:
159	// N : Number of data points {in Global common block }
160	// q(N) : Array of q values {in Global common block }
161	// dqv : Limit of slit length {in Global common block }
162	// m : powerlaw for extrapolation.
163	//
164	// Routine output:
165	// W(N,N) : Weights array
166	//
167	// The weights obtained from this routine can be used to calculate
168	// Is(q) for any I(q) function.
169	//
170	// Calculation based upon linear interpolation between points.
171	// See p.44, Computer Software Log book.
172	// 11/00 JGB
173	Function Weights_L(m,FW,Q_exp) //changed input w to FW (FW is two dimensional)
174	Variable m
175	Wave FW,Q_exp
176
177	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
178
179	NVAR dqv = $(USANSFolder+":DSM:gDqv")
180	NVAR NN = $(USANSFolder+":DSM:gNq")
181
182	// Calculate Remainder fractions and put into separate array.
183	Variable lower,ii,ss,jj
184	Duplicate/O Q_exp $(USANSFolder+":DSM:R_wave")
185	wave r_wave = $(USANSFolder+":DSM:R_wave")
186
187	// Make/O/D/N=75 root:DSM:SS_save //debug
188	// wave SS_save = root:DSM:SS_save
189
190	Print "calculating remainders by integration..."
191	for(ii=0;ii<NN-1;ii+=1)
192	lower = sqrt(Q_exp[NN-1]^2-Q_exp[ii]^2)
193	ss = Qromo(lower,dqv,Q_exp[ii],m) //new parameter list for Qromo call
194	// SS_save[ii] = ss
195	R_wave[ii] = (Q_exp[NN-1]^(-m)/dqv)*SS
196	// Printf "I = %d R_wave[ii] =%g\r",ii,R_wave[ii]
197	endfor
198
199	lower = 0.0
200	ss = Qromo(lower,dqv,Q_exp[NN-1],m) //new parameter list for Qromo call
201	R_wave[NN-1] = (Q_exp[NN-1]^(-m)/dqv)*SS
202	// Printf "I = %d R_wave[ii] =%g\r",NN-1,R_wave[NN-1]
203	// SS_save[ii] = ss
204
205	// Make/O/D/N=(75,75) root:DSM:IG_save //debug
206	// wave IG_save=root:DSM:IG_save
207
208	// Zero weight matrix... then fill it
209	FW = 0
210	Print "calculating full weights...."
211	// Calculate weights
212	for(ii=0;ii<NN;ii+=1)
213	for(jj=ii+1;jj<NN-1;jj+=1)
214	FW[ii][jj] = DU(ii,jj)*(1.0+Q_exp[jj]/(Q_exp[jj+1]-Q_exp[jj]))
215	FW[ii][jj] -= (1.0/(Q_exp[jj+1]-Q_exp[jj]))*IG(ii,jj)
216	FW[ii][jj] -= DU(ii,jj-1)*Q_exp[jj-1]/(Q_exp[jj]-Q_exp[jj-1])
217	FW[ii][jj] += (1.0/(Q_exp[jj]-Q_exp[jj-1]))*IG(ii,jj-1)
218	FW[ii][jj] *= (1.0/dqv)
219	// Printf "FW[%d][%d] = %g\r",ii,jj,FW[ii][jj]
220	endfor
221	endfor
222	//
223	// special case: I=J,I=N
224	for(ii=0;ii<NN-1;ii+=1)
225	FW[ii][ii] = DU(ii,ii)*(1.0+Q_exp[ii]/(Q_exp[ii+1]-Q_exp[ii]))
226	FW[ii][ii] -= (1.0/(Q_exp[ii+1]-Q_exp[ii]))*IG(ii,ii)
227	FW[ii][ii] *= (1.0/dqv)
228	// Printf "FW[%d][%d] = %g\r",ii,jj,FW[ii][jj]
229	// following line corrected for N -> NN-1 since Igor indexes from 0! (Q_exp[NN] DNE!)
230	FW[ii][NN-1] = -DU(ii,NN-2)*Q_exp[NN-2]/(Q_exp[NN-1]-Q_exp[NN-2])
231	FW[ii][NN-1] += (1.0/(Q_exp[NN-1]-Q_exp[NN-2]))*IG(ii,NN-2)
232	FW[ii][NN-1] *= (1.0/dqv)
233	FW[ii][NN-1] += R_wave[ii]
234	// Printf "FW[%d][%d] = %g\r",ii,NN-1,FW[ii][NN-1]
235	endfor
236	//
237	// special case: I=J=N
238	FW[NN-1][NN-1] = R_wave[NN-1]
239	//
240	Return (0)
241	End
242
243	////
244	Function Integrand_dsm(u,qi,m)
245	Variable u,qi,m
246
247	Variable integrand
248	Integrand = (u^2+qi^2)^(m/2)
249	return integrand
250	end
251
252	///
253	Function DU(ii,jj)
254	Variable ii,jj
255
256	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
257
258	Wave Q_exp = $(USANSFolder+":DSM:Q_exp")
259	Variable DU
260
261	// Wave DU_save=root:DSM:DU_save
262	If (ii == jj)
263	DU = sqrt(Q_exp[jj+1]^2-Q_exp[ii]^2)
264	Else
265	DU = sqrt(Q_exp[jj+1]^2-Q_exp[ii]^2) - sqrt(Q_exp[jj]^2-Q_exp[ii]^2)
266	EndIf
267
268	// DU_save[ii][jj] = DU
269	Return DU
270	End
271
272	Function IG(ii,jj)
273	Variable ii,jj
274
275	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
276
277	Wave Q_exp=$(USANSFolder+":DSM:Q_exp")
278	Variable IG,UL,UU
279
280	// WAVE IG_save = root:DSM:IG_save
281	If (ii == jj)
282	UL=0.0
283	Else
284	UL=sqrt(Q_exp[jj]^2-Q_exp[ii]^2)
285	EndIf
286	UU=sqrt(Q_exp[jj+1]^2-Q_exp[ii]^2)
287
288	//in FORTRAN log = natural log....
289	IG = 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