Context Navigation

LakeDesmearing_JB.ipf @ 425

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

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

Download in other formats:

Original Format