Context Navigation

PlotUtils2D_v40.ipf @ 1001

Last change on this file since 1001 was 1001, checked in by srkline, 7 years ago
A number of simple fixes to make the macros compatible with Igor 7, and in some places back-compatible with Igor 6. There were only a few instances where the IgorVersion? was checked, so it does not caues a huge disruption.
File size: 36.5 KB

Line
1	#pragma rtGlobals=1 // Use modern global access method.
2	#pragma version=4.00
3	#pragma IgorVersion=6.1
4
5	//
6	// TO DO:
7	//
8	// - more intelligent beam stop masking
9	// - constraints
10	// - interactive masking
11	//
12	//
13	// - OCT 2010 - added fDoBinning_QxQy2D(folderStr) that takes the QxQyQz data and bins it to I(Q)
14	// 1D result must still be plotted manually, can be automated later.
15	// - for Scaled Image data (QxQy), the function fDoBinning_Scaled2D (in FFT_Cubes) could be
16	// modified and added to this file. It was meant for an FFT slice, but applies to
17	// 2D model calculations (model_lin) 2D data files as well.
18
19	// Call the procedure BinQxQy_to_1D() to re-bin the QxQy data into I(q). good for testing to check the 2D error propagation
20	//
21	//
22
23
24	//
25	// changed June 2010 to read in the resolution information (now 8! columns)
26	// -- subject to change --
27	//
28	// look for either the old-style 3-column (no resolution information) or the newer 8-column format
29	Proc LoadQxQy()
30
31	SetDataFolder root:
32
33	LoadWave/G/D/W/A
34	String fileName = S_fileName
35	String path = S_Path
36	Variable numCols = V_flag
37
38	String w0,w1,w2,w3,w4,w5,w6,w7
39	String n0,n1,n2,n3,n4,n5,n6,n7
40
41	if(numCols == 8)
42	// put the names of the 8 loaded waves into local names
43	n0 = StringFromList(0, S_waveNames ,";" )
44	n1 = StringFromList(1, S_waveNames ,";" )
45	n2 = StringFromList(2, S_waveNames ,";" )
46	n3 = StringFromList(3, S_waveNames ,";" )
47	n4 = StringFromList(4, S_waveNames ,";" )
48	n5 = StringFromList(5, S_waveNames ,";" )
49	n6 = StringFromList(6, S_waveNames ,";" )
50	n7 = StringFromList(7, S_waveNames ,";" )
51
52	//remove the semicolon AND period from file names
53	w0 = CleanupName((S_fileName + "_qx"),0)
54	w1 = CleanupName((S_fileName + "_qy"),0)
55	w2 = CleanupName((S_fileName + "_i"),0)
56	w3 = CleanupName((S_fileName + "_iErr"),0)
57	w4 = CleanupName((S_fileName + "_qz"),0)
58	w5 = CleanupName((S_fileName + "_sQpl"),0)
59	w6 = CleanupName((S_fileName + "_sQpp"),0)
60	w7 = CleanupName((S_fileName + "_fs"),0)
61
62	String baseStr=w1[0,strlen(w1)-4]
63	if(DataFolderExists("root:"+baseStr))
64	DoAlert 1,"The file "+S_filename+" has already been loaded. Do you want to load the new data file, overwriting the data in memory?"
65	if(V_flag==2) //user selected No, don't load the data
66	SetDataFolder root:
67	KillWaves/Z $n0,$n1,$n2,$n3,$n4,$n5,$n6,$n7 // kill the default waveX that were loaded
68	return //quits the macro
69	endif
70	SetDataFolder $("root:"+baseStr)
71	else
72	NewDataFolder/S $("root:"+baseStr)
73	endif
74
75	//read in the 18 lines of header (18th line starts w/ ASCII... 19th line is blank)
76	Make/O/T/N=18 header
77	Variable refnum,ii
78	string tmpStr=""
79	Open/R refNum as (path+filename)
80	ii=0
81	do
82	tmpStr = ""
83	FReadLine refNum, tmpStr
84	header[ii] = tmpStr
85	ii+=1
86	while(ii < 18)
87	Close refnum
88
89	////overwrite the existing data, if it exists
90	Duplicate/O $("root:"+n0), $w0
91	Duplicate/O $("root:"+n1), $w1
92	Duplicate/O $("root:"+n2), $w2
93	Duplicate/O $("root:"+n3), $w3
94	Duplicate/O $("root:"+n4), $w4
95	Duplicate/O $("root:"+n5), $w5
96	Duplicate/O $("root:"+n6), $w6
97	Duplicate/O $("root:"+n7), $w7
98
99	endif //8-columns
100
101
102	// the very bad case where users don't supply the error data...
103	if(numCols == 3)
104	// put the names of the 3 loaded waves into local names
105	n0 = StringFromList(0, S_waveNames ,";" )
106	n1 = StringFromList(1, S_waveNames ,";" )
107	n2 = StringFromList(2, S_waveNames ,";" )
108
109	//remove the semicolon AND period from file names
110	w0 = CleanupName((S_fileName + "_qx"),0)
111	w1 = CleanupName((S_fileName + "_qy"),0)
112	w2 = CleanupName((S_fileName + "_i"),0)
113	w3 = CleanupName((S_fileName + "_iErr"),0) //make a name for the error wave, to be generated here
114
115	String baseStr=w1[0,strlen(w1)-4]
116	if(DataFolderExists("root:"+baseStr))
117	DoAlert 1,"The file "+S_filename+" has already been loaded. Do you want to load the new data file, overwriting the data in memory?"
118	if(V_flag==2) //user selected No, don't load the data
119	SetDataFolder root:
120	KillWaves/Z $n0,$n1,$n2 // kill the default waveX that were loaded
121	return //quits the macro
122	endif
123	SetDataFolder $("root:"+baseStr)
124	else
125	NewDataFolder/S $("root:"+baseStr)
126	endif
127
128	//read in the 18 lines of header (18th line starts w/ ASCII... 19th line is blank)
129	Make/O/T/N=18 header
130	Variable refnum,ii
131	string tmpStr=""
132	Open/R refNum as (path+filename)
133	ii=0
134	do
135	tmpStr = ""
136	FReadLine refNum, tmpStr
137	header[ii] = tmpStr
138	ii+=1
139	while(ii < 18)
140	Close refnum
141
142	////overwrite the existing data, if it exists
143	Duplicate/O $("root:"+n0), $w0
144	Duplicate/O $("root:"+n1), $w1
145	Duplicate/O $("root:"+n2), $w2
146
147
148
149	// generate my own error wave for I(qx,qy). This is exactly the same procedure that is used in the QxQy_Export function
150	Duplicate/O $("root:"+n0), $w3
151	$w3 = sqrt($w2) //assumes Poisson statistics for each cell (counter)
152	// sw = 0.05*sw // uniform 5% error? tends to favor the low intensity too strongly
153	// get rid of the "bad" errorsby replacing the NaN, Inf, and zero with V_avg
154	// THIS IS EXTREMEMLY IMPORTANT - if this is not done, there are some "bad" values in the
155	// error wave (things that are not numbers) - and this wrecks the smeared model fitting.
156	// It appears to have no effect on the unsmeared model.
157	WaveStats/Q $w3
158	$w3 = numtype($w3[p]) == 0 ? $w3[p] : V_avg
159	$w3 = $w3[p] != 0 ? $w3[p] : V_avg
160
161	endif //3-columns
162
163	// with 4-columns, assume that they are qx,qy,I, err
164	if(numCols == 4)
165	// put the names of the 3 loaded waves into local names
166	n0 = StringFromList(0, S_waveNames ,";" )
167	n1 = StringFromList(1, S_waveNames ,";" )
168	n2 = StringFromList(2, S_waveNames ,";" )
169	n3 = StringFromList(3, S_waveNames ,";" )
170
171	//remove the semicolon AND period from file names
172	w0 = CleanupName((S_fileName + "_qx"),0)
173	w1 = CleanupName((S_fileName + "_qy"),0)
174	w2 = CleanupName((S_fileName + "_i"),0)
175	w3 = CleanupName((S_fileName + "_iErr"),0) //make a name for the error wave, to be generated here
176
177	String baseStr=w1[0,strlen(w1)-4]
178	if(DataFolderExists("root:"+baseStr))
179	DoAlert 1,"The file "+S_filename+" has already been loaded. Do you want to load the new data file, overwriting the data in memory?"
180	if(V_flag==2) //user selected No, don't load the data
181	SetDataFolder root:
182	KillWaves/Z $n0,$n1,$n2,$n3 // kill the default waveX that were loaded
183	return //quits the macro
184	endif
185	SetDataFolder $("root:"+baseStr)
186	else
187	NewDataFolder/S $("root:"+baseStr)
188	endif
189
190	//read in the 18 lines of header (18th line starts w/ ASCII... 19th line is blank)
191	Make/O/T/N=18 header
192	Variable refnum,ii
193	string tmpStr=""
194	Open/R refNum as (path+filename)
195	ii=0
196	do
197	tmpStr = ""
198	FReadLine refNum, tmpStr
199	header[ii] = tmpStr
200	ii+=1
201	while(ii < 18)
202	Close refnum
203
204	////overwrite the existing data, if it exists
205	Duplicate/O $("root:"+n0), $w0
206	Duplicate/O $("root:"+n1), $w1
207	Duplicate/O $("root:"+n2), $w2
208	Duplicate/O $("root:"+n3), $w3
209
210	endif //4-columns
211
212
213	/// do this for all 2D data, whether or not resolution information was read in
214
215	Variable/G gIsLogScale = 0
216
217	Variable num=numpnts($w0)
218	// assume that the Q-grid is "uniform enough" for DISPLAY ONLY
219	// use the 3 original waves for all of the fitting...
220	ConvertQxQy2Mat($w0,$w1,$w2,baseStr+"_mat")
221	Duplicate/O $(baseStr+"_mat"),$(baseStr+"_lin") //keep a linear-scaled version of the data
222
223	PlotQxQy(baseStr) //this sets the data folder back to root:!!
224
225	//clean up
226	SetDataFolder root:
227	KillWaves/Z $n0,$n1,$n2,$n3,$n4,$n5,$n6,$n7
228
229	EndMacro
230
231	//does not seem to need to be flipped at all from the standard QxQy output
232	//
233	// plots the data, as surface0, YellowHot color table
234	//
235	Proc PlotQxQy(str)
236	String str
237	Prompt str,"Pick the data folder containing 2D data",popup,getAList(4)
238
239	PauseUpdate; Silent 1 // Building Gizmo 6 window...
240
241	// Do nothing if the Gizmo XOP is not available.
242	if(exists("NewGizmo")!=4)
243	DoAlert 0, "Gizmo XOP must be installed"
244	return
245	endif
246
247	String data = "root:"+str+":"+str+"_mat"
248
249	NewGizmo/N=$str/T=str /W=(169,44,495,359)
250	ModifyGizmo startRecMacro
251	AppendToGizmo Surface=$data,name=surface0
252
253	// for constant color (data is a blue color, lighter on the bottom)
254	// ModifyGizmo ModifyObject=surface0 property={ surfaceColorType,1}
255	// ModifyGizmo ModifyObject=surface0 property={ srcMode,0}
256	// ModifyGizmo ModifyObject=surface0 property={ frontColor,0.749996,0.811093,1,1}
257	// ModifyGizmo ModifyObject=surface0 property={ backColor,0.250019,0.433326,1,1}
258	//e constant color
259
260	//for color table
261	ModifyGizmo ModifyObject=surface0 property={ srcMode,0}
262	// ModifyGizmo ModifyObject=surface0 property={ surfaceCTab,Blue}
263	ModifyGizmo ModifyObject=surface0 property={ surfaceCTab,ColdWarm}
264	ModifyGizmo ModifyObject=surface0 property={ SurfaceCTABScaling,4}
265	//e color table
266
267	AppendToGizmo Axes=boxAxes,name=axes0
268	ModifyGizmo ModifyObject=axes0,property={0,axisRange,-1,-1,-1,1,-1,-1}
269	ModifyGizmo ModifyObject=axes0,property={1,axisRange,-1,-1,-1,-1,1,-1}
270	ModifyGizmo ModifyObject=axes0,property={2,axisRange,-1,-1,-1,-1,-1,1}
271	ModifyGizmo ModifyObject=axes0,property={3,axisRange,-1,1,-1,-1,1,1}
272	ModifyGizmo ModifyObject=axes0,property={4,axisRange,1,1,-1,1,1,1}
273	ModifyGizmo ModifyObject=axes0,property={5,axisRange,1,-1,-1,1,-1,1}
274	ModifyGizmo ModifyObject=axes0,property={6,axisRange,-1,-1,1,-1,1,1}
275	ModifyGizmo ModifyObject=axes0,property={7,axisRange,1,-1,1,1,1,1}
276	ModifyGizmo ModifyObject=axes0,property={8,axisRange,1,-1,-1,1,1,-1}
277	ModifyGizmo ModifyObject=axes0,property={9,axisRange,-1,1,-1,1,1,-1}
278	ModifyGizmo ModifyObject=axes0,property={10,axisRange,-1,1,1,1,1,1}
279	ModifyGizmo ModifyObject=axes0,property={11,axisRange,-1,-1,1,1,-1,1}
280	ModifyGizmo ModifyObject=axes0,property={-1,axisScalingMode,1}
281	ModifyGizmo ModifyObject=axes0,property={-1,axisColor,0,0,0,1}
282	ModifyGizmo ModifyObject=axes0,property={0,ticks,3}
283	ModifyGizmo ModifyObject=axes0,property={1,ticks,3}
284	ModifyGizmo ModifyObject=axes0,property={2,ticks,3}
285	ModifyGizmo ModifyObject=axes0,property={0,fontScaleFactor,1.5}
286	ModifyGizmo ModifyObject=axes0,property={1,fontScaleFactor,1.5}
287	ModifyGizmo ModifyObject=axes0,property={2,fontScaleFactor,1.5}
288	AppendToGizmo freeAxesCue={0,0,0,1.5},name=freeAxesCue0
289	ModifyGizmo setDisplayList=0, opName=clearColor0, operation=clearColor, data={1,0.917,0.75,1}
290	ModifyGizmo setDisplayList=1, object=surface0
291	ModifyGizmo setDisplayList=2, object=axes0
292	ModifyGizmo setDisplayList=3, object=freeAxesCue0
293	ModifyGizmo SETQUATERNION={0.521287,-0.097088,-0.138769,0.836408}
294	ModifyGizmo autoscaling=1
295	ModifyGizmo currentGroupObject=""
296	ModifyGizmo compile
297
298	ModifyGizmo NamedHook={GizmoContours,WMGizmoContoursNamedHook}
299	ModifyGizmo endRecMacro
300
301	// don't bother with the flat image plot right now
302
303	// Display $yw vs $xw
304	// modifygraph log=0
305	// ModifyGraph mode=3,marker=16,zColor($yw)={$zw,,,YellowHot,0}
306	// ModifyGraph standoff=0
307	// ModifyGraph width={Aspect,1}
308	// ModifyGraph lowTrip=0.001
309
310	End
311
312	// if flag set, data folder will be overwritten (only if it exists)
313	// otherwise, the user will be asked what to do if there is a name conflict
314	//
315	Proc FakeQxQy(minQx,maxQx,minQy,maxQy,numPix,dataFolder,forceOverwrite,doPlot)
316	Variable minQx=-0.1,maxQx=0.1,minQy=-0.1,maxQy=0.1
317	Variable numPix=64
318	String dataFolder="fake2DData"
319	Variable forceOverwrite = 0
320	Variable doPlot = 1
321
322	String baseStr=dataFolder
323
324	if(DataFolderExists("root:"+baseStr) && forceOverwrite)
325	SetDataFolder $("root:"+baseStr)
326	else
327	if(DataFolderExists("root:"+baseStr))
328	DoAlert 1,"Fake data "+baseStr+" has already been created. Do you want to overwrite this fake data set?"
329	if(V_flag==2) //user selected No, don't load the data
330	SetDataFolder root:
331	return //quits the macro
332	endif
333	SetDataFolder $("root:"+baseStr)
334	else
335	NewDataFolder/S $("root:"+baseStr)
336	endif
337	endif
338
339
340	Make/O/D/N=(numPix*numPix) $(baseStr+"_qx"),$(baseStr+"_qy"),$(baseStr+"_i")
341
342	$(baseStr+"_i") = 1
343
344	Make/O/D/N=(numPix) tmpX,tmpY
345
346	variable delX,delY
347	delX = (maxQx - minQx)/numPix
348	delY = (maxQy - minQy)/numPix
349	tmpX = minQx + x*delX
350	tmpY = minQy + x*delY
351
352	// get rid of Q=0 values in the waves (may be a bad actor in model calculations)
353	//
354	tmpX = (tmpX == 0) ? 1e-6 : tmpX
355	tmpY = (tmpY == 0) ? 1e-6 : tmpY
356
357	// X wave varies more rapidly
358	// Y wave is blocks of values
359	Variable ii=0
360	do
361	$(baseStr+"_qy")[iinumPix,(ii+1)numPix] = tmpY[ii]
362	ii+=1
363	while(ii<numPix)
364
365	$(baseStr+"_qx") = tmpX[mod(p,numPix)]
366
367	Variable/G gIsLogScale = 0
368
369	// assume that the Q-grid is "uniform enough" for DISPLAY ONLY
370	// use the 3 original waves for all of the fitting...
371	ConvertQxQy2Mat($(baseStr+"_qx"),$(baseStr+"_qy"),$(baseStr+"_i"),baseStr+"_mat")
372	Duplicate/O $(baseStr+"_mat"),$(baseStr+"_lin") //keep a linear-scaled version of the data
373
374	if(doPlot)
375	PlotQxQy(baseStr) //this sets the data folder back to root:, but may not be executed
376	endif
377	//clean up
378	SetDataFolder root:
379
380	EndMacro
381
382	// this assumes that:
383	// --QxQy data was written out in the format specified by the Igor macros, that is the x varies most rapidly
384	// --the matrix is square!
385	//
386	// probably some other stuff...
387	//
388	Function UpdateQxQy2Mat(Qx,Qy,inten,linMat,mat)
389	Wave Qx,Qy,inten,linMat,mat
390
391	Variable xrows=DimSize(mat, 0 ) //assumes square matrix!!
392
393	String folderStr=GetWavesDataFolder(Qx,1)
394	NVAR/Z gIsLogScale=$(folderStr+"gIsLogScale")
395
396	linMat = inten[q*xrows+p]
397
398	if(gIsLogScale)
399	mat = log(linMat)
400	else
401	mat = linMat
402	endif
403
404	return(0)
405	End
406
407	// this assumes that:
408	// --QxQy data was written out in the format specified by the Igor macros, that is the x varies most rapidly
409	// --the matrix is square!
410	//
411	// probably some other stuff...
412	//
413	Function ConvertQxQy2Mat(Qx,Qy,inten,matStr)
414	Wave Qx,Qy,inten
415	String matStr
416
417	String folderStr=GetWavesDataFolder(Qx,1)
418
419	Variable num=sqrt(numpnts(Qx)) //assumes square matrix, Qx = num x num points long
420	Make/O/D/N=(num,num) $(folderStr + matStr)
421	Wave mat=$matStr
422
423	WaveStats/Q Qx
424	SetScale/I x, V_min, V_max, "", mat
425	WaveStats/Q Qy
426	SetScale/I y, V_min, V_max, "", mat
427
428	Variable xrows=num
429
430	mat = inten[q*xrows+p]
431
432	return(0)
433	End
434
435
436	//str is the full path to the surface to append
437	Proc AppendSurfaceToGizmo(str)
438	String str
439
440	PauseUpdate; Silent 1 // Building Gizmo 6 window...
441
442	AppendToGizmo/Z Surface=$str,name=surface1
443
444	// for a constant color (model is darker on top, lighter on the bottom)
445	//need these two lines, plus a color
446	ModifyGizmo ModifyObject=surface1 property={ surfaceColorType,1}
447	ModifyGizmo ModifyObject=surface1 property={ srcMode,0}
448	//green
449	ModifyGizmo ModifyObject=surface1 property={ frontColor,0.528923,0.882353,0.321584,1}
450	ModifyGizmo ModifyObject=surface1 property={ backColor,0.300221,0.6,1.5259e-05,1}
451	//red
452	// ModifyGizmo ModifyObject=surface1 property={ frontColor,1,0.749996,0.749996,1}
453	// ModifyGizmo ModifyObject=surface1 property={ backColor,1,0.250019,0.250019,1}
454
455	ModifyGizmo ModifyObject=surface1 property={ SurfaceCTABScaling,4}
456
457	// for a color table (maybe not a good choice for vizualization if data uses a color table)
458	// ModifyGizmo ModifyObject=surface1 property={ srcMode,0}
459	// ModifyGizmo ModifyObject=surface1 property={ surfaceCTab,Red}
460
461	// ModifyGizmo setDisplayList=0, object=surface0
462	// ModifyGizmo setDisplayList=1, object=axes0
463	// object 3 is the axisCue
464	ModifyGizmo setDisplayList=4, object=surface1
465	ModifyGizmo SETQUATERNION={0.565517,-0.103105,-0.139134,0.806350}
466	// ModifyGizmo autoscaling=1
467	// ModifyGizmo currentGroupObject=""
468	ModifyGizmo compile
469
470	// ModifyGizmo endRecMacro
471	End
472
473	// would be nice, but I can't get this to work...
474	//
475	//Macro AdjustColorTables()
476	//
477	// ModifyGizmo ModifyObject=surface1 property={ srcMode,0}
478	// ModifyGizmo ModifyObject=surface1 property={ surfaceCTab,Rainbow}
479	// ModifyGizmo setDisplayList=3, object=surface1
480	//// ModifyGizmo SETQUATERNION={0.565517,-0.103105,-0.139134,0.806350}
481	// ModifyGizmo autoscaling=1
482	// ModifyGizmo currentGroupObject=""
483	// ModifyGizmo compile
484	//
485	//end
486
487	Function LogToggle2DButtonProc(ba) : ButtonControl
488	STRUCT WMButtonAction &ba
489
490	switch( ba.eventCode )
491	case 2: // mouse up
492	// click code here
493	String winStr=TopGizmoWindow()
494	if(strlen(winstr)==0)
495	Print "no gizmo window"
496	break
497	endif
498	// the winStr is also the data folder
499	// toggle everything in the data folder
500	ToggleFolderScaling(winStr)
501	break
502	endswitch
503
504	return 0
505	End
506
507	Function Plot2DButtonProc(ba) : ButtonControl
508	STRUCT WMButtonAction &ba
509
510	switch( ba.eventCode )
511	case 2: // mouse up
512	// click code here
513	Execute "PlotQxQy()"
514	break
515	endswitch
516
517	return 0
518	End
519
520	Function Append2DModelButtonProc(ba) : ButtonControl
521	STRUCT WMButtonAction &ba
522
523	switch( ba.eventCode )
524	case 2: // mouse up
525	// click code here
526	//str is the full path to the surface to append
527	String DF=TopGizmoWindow()
528	String str,funcStr,suffix
529	ControlInfo/W=WrapperPanel popup_1
530	funcStr = S_Value
531	suffix = getModelSuffix(funcStr)
532	if(stringmatch(funcStr, "smear*") == 1)
533	suffix = "sm_"+suffix
534	endif
535	str = "root:"+DF+":"+suffix+"_mat"
536
537	Execute "AppendSurfaceToGizmo(\""+str+"\")"
538	// Print str
539	break
540	endswitch
541
542	return 0
543	End
544
545	// toggle the scaling of every matrix in the folder
546	Function ToggleFolderScaling(DF)
547	String DF
548
549	// look for waves DF+"_mat" (the data)
550	// and models prefix + "_mat"
551	SetDataFolder $("root:"+DF)
552	//check the global to see the state of the data
553	NVAR gIsLogScale = gIsLogScale
554
555	String matrixList=WaveList("*_mat",";",""),item
556	Variable ii=0,num=itemsinlist(matrixlist,";"),len
557
558	for(ii=0;ii<num;ii+=1)
559	item = StringFromList(ii,matrixList,";")
560	len = strlen(item)
561	Wave w = $item
562	Wave linW = $(item[0,len-5]+"_lin")
563	if(gisLogScale)
564	//make linear
565	w=linW
566	else
567	//make log scale
568	w=log(linW)
569	w = w[p][q] == -inf ? NaN : w[p][q] //remove the -inf for display
570	endif
571	endfor
572
573	// toggle the global
574	gIsLogScale = !gIsLogScale
575
576	SetDataFolder root:
577	return(0)
578	End
579
580
581	Function Load2DDataButtonProc(ba) : ButtonControl
582	STRUCT WMButtonAction &ba
583
584	switch( ba.eventCode )
585	case 2: // mouse up
586	// click code here
587	Execute "LoadQxQy()"
588	break
589	endswitch
590
591	return 0
592	End
593
594
595	// dispatch the fit
596	Function Do2DFitButtonProc(ba) : ButtonControl
597	STRUCT WMButtonAction &ba
598
599	String folderStr,funcStr,coefStr
600	Variable useCursors,useEps,useConstr
601
602	switch( ba.eventCode )
603	case 2: // mouse up
604	ControlInfo/W=WrapperPanel popup_0
605	folderStr=S_Value
606
607	ControlInfo/W=WrapperPanel popup_1
608	funcStr=S_Value
609
610	ControlInfo/W=WrapperPanel popup_2
611	coefStr=S_Value
612
613	ControlInfo/W=WrapperPanel check_0
614	useCursors=V_Value
615	//
616	// NO CURSORS for 2D waves - force to zero
617	//
618	useCursors = 0
619	//
620	ControlInfo/W=WrapperPanel check_1
621	useEps=V_Value
622	ControlInfo/W=WrapperPanel check_2
623	useConstr=V_Value
624
625	if(!CheckFunctionAndCoef(funcStr,coefStr))
626	DoAlert 0,"The coefficients and function type do not match. Please correct the selections in the popup menus."
627	break
628	endif
629
630	FitWrapper2D(folderStr,funcStr,coefStr,useCursors,useEps,useConstr)
631
632	// DoUpdate (does not work!)
633	//?? why do I need to force an update ??
634	if(!exists("root:"+folderStr+":"+coefStr))
635	Wave w=$coefStr
636	else
637	Wave w=$("root:"+folderStr+":"+coefStr) //smeared coefs in data folder
638	endif
639	w[0] += 1e-6
640	w[0] -= 1e-6
641
642	break
643	endswitch
644
645	return 0
646	End
647
648	// wrapper to do the desired fit, pretty much a clone of FitWrapper, the 1D version
649	//
650	// folderStr is the data folder for the desired data set
651	//
652	// Currently the limitations are:
653	// - I have no error waves for the intensity (fixed 10/2010)
654	// - There is no smeared model (coming soon after 10/2010)
655	// - Cursors can't be used
656	// - the report works OK, but I have little control over the graphics
657	// - the mask is generated here with a default radius of 8 pixels around the beam center
658	//
659	Function FitWrapper2D(folderStr,funcStr,coefStr,useCursors,useEps,useConstr)
660	String folderStr,funcStr,coefStr
661	Variable useCursors,useEps,useConstr
662
663	//These only make sense for the 1D fits, but put them here so keep the look of the dispatching the same
664	Variable useResiduals, useTextBox
665	useResiduals = 0
666	useTextBox = 0
667
668	String suffix=getModelSuffix(funcStr)
669
670	SetDataFolder $("root:"+folderStr)
671	if(!exists(coefStr))
672	// must be unsmeared model, work in the root folder
673	SetDataFolder root:
674	if(!exists(coefStr)) //this should be fine if the coef filter is working, but check anyhow
675	DoAlert 0,"the coefficient and data sets do not match"
676	return 0
677	endif
678	endif
679
680	WAVE cw=$(coefStr)
681	Wave hold=$("Hold_"+suffix)
682	Wave/T lolim=$("LoLim_"+suffix)
683	Wave/T hilim=$("HiLim_"+suffix)
684	Wave eps=$("epsilon_"+suffix)
685
686	// fill a struct instance whether I need one or not
687	// note that the resolution waves may or may not exist, and may or may not be used in the fitting
688	String DF="root:"+folderStr+":"
689
690	WAVE inten=$(DF+folderStr+"_i")
691	WAVE sw=$(DF+folderStr+"_iErr")
692	WAVE qx=$(DF+folderStr+"_qx")
693	WAVE qy=$(DF+folderStr+"_qy")
694	WAVE/Z qz=$(DF+folderStr+"_qz")
695	WAVE/Z sQpl=$(DF+folderStr+"_sQpl")
696	WAVE/Z sQpp=$(DF+folderStr+"_sQpp")
697	WAVE/Z shad=$(DF+folderStr+"_fs")
698
699	//just a dummy - I shouldn't need this
700	Duplicate/O qx resultW
701	resultW=0
702
703	STRUCT ResSmear_2D_AAOStruct s
704	WAVE s.coefW = cw
705	WAVE s.zw = resultW
706	WAVE s.xw[0] = qx
707	WAVE s.xw[1] = qy
708	WAVE/Z s.qz = qz
709	WAVE/Z s.sQpl = sQpl
710	WAVE/Z s.sQpp = sQpp
711	WAVE/Z s.fs = shad
712
713
714	// generate my own mask wave - as a matrix first, then redimension to N*N vector
715	// default mask is two pixels all the way around, (0 is excluded, 1 is included)
716	WAVE DataMat=$(DF+folderStr+"_lin")
717	if(exists(DF+"mask") == 0)
718	Duplicate/O dataMat mask
719	Variable bsRadius=8 //pixels?
720	MakeBSMask(mask,bsRadius)
721	Endif
722
723
724	Duplicate/O inten inten_masked
725	inten_masked = (mask[p][q] == 0) ? NaN : inten[p][q]
726
727
728	//for now, use res will always be 0 for 2D functions
729	Variable useResol=0
730	if(stringmatch(funcStr, "Smear*")) // if it's a smeared function, need a struct
731	useResol=1
732	endif
733
734	// can't use constraints defined as a single text wave for multivariate fits. See the curve fitting help file
735	// and "Contraint Matrix and Vector"
736	//
737	// -- generate a constraint text wave, then the /C flag automatically generates a constraint matrix and vector
738	// -- then use the proper logic to dispatch (can't use /NWOK anymore)
739
740	// if(useConstr)
741	// Print "Constraints not yet implemented"
742	// useConstr = 0
743	// endif
744	// WAVE/Z constr=constr //will be a null reference
745
746	// do not construct constraints for any of the coefficients that are being held
747	// -- this will generate an "unknown error" from the curve fitting
748	Make/O/T/N=0 constr
749	if(useConstr)
750	String constraintExpression
751	Variable i, nPnts=DimSize(lolim, 0),nextRow=0
752	for (i=0; i < nPnts; i += 1)
753	if (strlen(lolim[i]) > 0 && hold[i] == 0)
754	InsertPoints nextRow, 1, constr
755	sprintf constraintExpression, "K%d > %s", i, lolim[i]
756	constr[nextRow] = constraintExpression
757	nextRow += 1
758	endif
759	if (strlen(hilim[i]) > 0 && hold[i] == 0)
760	InsertPoints nextRow, 1, constr
761	sprintf constraintExpression, "K%d < %s", i, hilim[i]
762	constr[nextRow] = constraintExpression
763	nextRow += 1
764	endif
765	endfor
766
767	endif
768
769	if(useCursors)
770	Print "Cursors not yet implemented"
771	useCursors = 0
772	endif
773	///// NO CURSORS for 2D waves
774	//if useCursors, and the data is USANS, need to feed a (reassigned) trimmed matrix to the fit
775	Variable pt1,pt2,newN
776	pt1 = 0
777	pt2 = numpnts(inten)-1
778	// if(useCursors && (dimsize(resW,1) > 4) )
779	// if(pcsr(A) > pcsr(B))
780	// pt1 = pcsr(B)
781	// pt2 = pcsr(A)
782	// else
783	// pt1 = pcsr(A)
784	// pt2 = pcsr(B)
785	// endif
786	// newN = pt2 - pt1 + 1 // +1 includes both cursors in the fit
787	// Make/O/D/N=(newN,newN) $(DF+"crsrResW")
788	// WAVE crsrResW = $(DF+"crsrResW")
789	// crsrResW = resW[p+pt1][q+pt1]
790	// //assign to the struct
791	// WAVE fs.resW = crsrResW
792	// endif
793
794	// create these variables so that FuncFit will set them on exit
795	Variable/G V_FitError=0 //0=no err, 1=error,(2^1+2^0)=3=singular matrix
796	Variable/G V_FitQuitReason=0 //0=ok,1=maxiter,2=user stop,3=no chisq decrease
797
798	// don't use the auto-destination with no flag, it doesn't appear to work correctly
799	// dispatch the fit
800	// FuncFit/H="11110111111"/NTHR=0 Cylinder2D_D :cyl2d_c_txt:coef_Cyl2D_D :cyl2d_c_txt:cyl2d_c_txt_i /X={:cyl2d_c_txt:cyl2d_c_txt_qy,:cyl2d_c_txt:cyl2d_c_txt_qx} /W=:cyl2d_c_txt:sw /I=1 /M=:cyl2d_c_txt:mask /D
801	Variable t1=StopMSTimer(-2)
802	Variable tb = 0 //no textbox
803
804	// /NTHR=1 means just one thread for the fit (since the function itself is threaded)
805	// NTHR = 0 == "Auto" mode, using as many processors as are available (not appropriate here since the function itself is threaded?)
806
807	do
808
809	// now useEps, and useConstr are all handled w/ /NWOK, just like FitWrapper
810	// useCursors needs to have the /C flag in the command for the constraint matrix and vector to be auto-generated
811
812	// if(useResol && useResiduals && useTextBox) //do it all
813	// FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s /R /NWOK
814	// break
815	// endif
816	//
817	// if(useResol && useResiduals) //res + resid
818	// FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s /R /NWOK
819	// break
820	// endif
821	//
822	//
823	// if(useResol && useTextBox) //res + text
824	// FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s /NWOK
825	// break
826	// endif
827
828	if(useResol && useConstr) //res and constraints
829	Print "useRes only"
830	FuncFit/C/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s /NWOK
831	break
832	endif
833
834	if(useResol) //res only
835	Print "useRes only"
836	FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /STRC=s /NWOK
837	break
838	endif
839
840	///// same as above, but all without useResol (no /STRC flag)
841	// if(useResiduals && useTextBox) //resid+ text
842	// FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /R /NWOK
843	// break
844	// endif
845	//
846	// if(useResiduals) //resid
847	// FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /R /NWOK
848	// break
849	// endif
850	//
851	//
852	// if(useTextBox) //text
853	// FuncFit/H=getHStr(hold) /NTHR=0 /TBOX=(tb) $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /NWOK
854	// break
855	// endif
856
857	if(useConstr)
858	FuncFit/C/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /NWOK
859	break
860	Endif
861
862	//just a plain vanilla fit
863
864	FuncFit/H=getHStr(hold) /NTHR=0 $funcStr cw, inten[pt1,pt2] /X={Qx,Qy} /M=mask /W=sw /I=1 /E=eps /C=constr /NWOK
865
866	while(0)
867
868	Print "elapsed 2D fit time = ",(StopMSTimer(-2) - t1)/1e6," seconds = ",(StopMSTimer(-2) - t1)/1e6/60," minutes"
869
870	// append the fit
871	// need to manage duplicate copies
872	// Don't plot the full curve if cursors were used (set fitYw to NaN on entry...)
873	// String traces=TraceNameList("", ";", 1 ) //"" as first parameter == look on the target graph
874	// if(strsearch(traces,"FitYw",0) == -1)
875	// AppendToGraph FitYw vs xw
876	// else
877	// RemoveFromGraph FitYw
878	// AppendToGraph FitYw vs xw
879	// endif
880	// ModifyGraph lsize(FitYw)=2,rgb(FitYw)=(0,0,0)
881
882	// DoUpdate //force update of table and graph with fitted values (why doesn't this work? - the table still does not update)
883
884	// report the results (to the panel?)
885	print "V_chisq = ",V_chisq
886	print cw
887	WAVE w_sigma
888	print w_sigma
889	String resultStr=""
890
891	if(waveexists(W_sigma))
892	//append it to the table, if it's not already there
893	CheckDisplayed/W=WrapperPanel#T0 W_sigma
894	if(V_flag==0)
895	//not there, append it
896	AppendtoTable/W=wrapperPanel#T0 W_sigma
897	else
898	//remove it, and put it back on to make sure it's the right one (do I need to do this?)
899	// -- not really, since any switch of the function menu takes W_Sigma off
900	endif
901	endif
902
903	//now re-write the results
904	sprintf resultStr,"Chi^2 = %g Sqrt(X^2/N) = %g",V_chisq,sqrt(V_chisq/V_Npnts)
905	resultStr = PadString(resultStr,63,0x20)
906	GroupBox grpBox_2 title=resultStr
907	ControlUpdate/W=WrapperPanel grpBox_2
908	sprintf resultStr,"FitErr = %s : FitQuit = %s",W_ErrorMessage(V_FitError),W_QuitMessage(V_FitQuitReason)
909	resultStr = PadString(resultStr,63,0x20)
910	GroupBox grpBox_3 title=resultStr
911	ControlUpdate/W=WrapperPanel grpBox_3
912
913	Variable yesSave=0,yesReport=0
914	ControlInfo/W=WrapperPanel check_4
915	yesReport = V_Value
916	ControlInfo/W=WrapperPanel check_5
917	yesSave = V_Value
918
919
920	if(yesReport)
921	String parStr = getFunctionParams(funcStr)
922	// String parStr=GetWavesDataFolder(cw,1)+ WaveList("param"+suffix, "", "TEXT:1," ) //old way, but doesn't work in 2D folders
923	String topGraph= TopGizmoWindow() //this is the topmost Gizmo (XOP) window
924
925	DoUpdate //force an update of the graph before making a copy of it for the report
926
927	W_GenerateReport(funcStr,folderStr,$parStr,cw,yesSave,V_chisq,W_sigma,V_npnts,V_FitError,V_FitQuitReason,V_startRow,V_endRow,topGraph)
928	endif
929
930	SetDataFolder root:
931	return(0)
932	End
933
934
935	// right now, there are only unsmeared models, but I want experimental points
936	// for the QxQy of the calculation - so get the folder string
937	Function Plot2DFunctionButtonProc(ba) : ButtonControl
938	STRUCT WMButtonAction &ba
939
940	String folderStr,funcStr,coefStr,cmdStr=""
941
942	Variable killWhat=0 //kill nothing as default
943
944	switch( ba.eventCode )
945	case 2: // mouse up
946	ControlInfo/W=WrapperPanel popup_0
947	folderStr=S_Value
948
949	ControlInfo/W=WrapperPanel popup_1
950	funcStr=S_Value
951
952	// check for smeared or smeared function
953	// if(stringmatch(funcStr, "Smear*" )==1)
954	//it's a smeared model
955	// check for the special case of RPA that has an extra parameter
956	// if(strsearch(funcStr, "RPAForm", 0 ,0) == -1)
957	sprintf cmdStr, "Plot%s(\"%s\")",funcStr,folderStr //not RPA
958	// else
959	// sprintf cmdStr, "Plot%s(\"%s\",)",funcStr,folderStr //yes RPA, leave a comma for input
960	// endif
961	// else
962	// // it's not,
963	// sprintf cmdStr, "Plot%s()",funcStr
964	// endif
965
966	//Print cmdStr
967	Execute cmdStr
968
969	//pop the function menu to set the proper coefficients
970	DoWindow/F WrapperPanel
971	STRUCT WMPopupAction pa
972	pa.popStr = funcStr
973	pa.eventcode = 2
974	Function_PopMenuProc(pa)
975
976	KillWhat = 2 //kill just the table, leave the 2d visible for now
977	KillTopGraphAndTable(killWhat) // crude
978
979	break
980	endswitch
981
982	return 0
983	End
984
985	// unused - all of this functionality has been added to the Wrapper Panel
986	Window Plot_2D_Controls()
987	PauseUpdate; Silent 1 // building window...
988	NewPanel /W=(950,44,1250,244)
989	Button button0,pos={172,37},size={70,20},proc=LogToggle2DButtonProc,title="Log/Lin"
990	Button button1,pos={146,9},size={100,20},proc=Plot2DButtonProc,title="Plot 2D Data"
991	Button button2,pos={164,118},size={120,20},proc=Append2DModelButtonProc,title="Append Model"
992	PopupMenu popup_1,pos={9,84},size={132,20},proc=Function_PopMenuProc,title="Function"
993	PopupMenu popup_1,mode=1,popvalue="Cylinder2D",value= #"W_FunctionPopupList()"
994	Button button3,pos={11,9},size={100,20},proc=Load2DDataButtonProc,title="Load 2D Data"
995	Button button4,pos={164,161},size={120,20},proc=Do2DFitButtonProc,title="Do 2D Fit"
996	Button button5,pos={164,85},size={120,20},proc=Plot2DFunctionButtonProc,title="Plot 2D Function"
997	EndMacro
998
999	//
1000	// in Igor 6, Gizmo windows are code 4096
1001	// in Igor 7, Gizmo windows are code 65536
1002	// need to check for version, since code 65536 is not allowed in Igor 6
1003	//
1004	// this function returns only the topmost window from the list
1005	Function/S TopGizmoWindow()
1006	if(IgorVersion() >= 7)
1007	return(StringFromList(0,WinList("*",";","WIN:65536")))
1008	else
1009	return(StringFromList(0,WinList("*",";","WIN:4096")))
1010	endif
1011	end
1012
1013	////
1014	//
1015	// unused proc for testing
1016	Proc Plot2DVsPointNumber(str)
1017	String str
1018	Prompt str,"Pick the data folder containing 2D data",popup,getAList(4)
1019
1020	PauseUpdate; Silent 1 // building window...
1021	String fldrSav0= GetDataFolder(1)
1022	SetDataFolder $("root:"+str)
1023	Display /W=(241,352,810,842) $(str+"_i")
1024	SetDataFolder fldrSav0
1025	ModifyGraph mode=3
1026	ModifyGraph marker=19
1027	ModifyGraph lSize=2
1028	ModifyGraph rgb($(str+"_i"))=(0,0,0)
1029	ModifyGraph msize=1
1030	ModifyGraph grid=1
1031	ModifyGraph log(left)=1
1032	ModifyGraph mirror=2
1033	if(exists("root:"+str+":sw"))
1034	ErrorBars/T=0 $(str+"_i") Y,wave=($("root:"+str+":sw"),$("root:"+str+":sw"))
1035	endif
1036	SetDataFolder fldrSav0
1037	EndMacro
1038
1039	// testing procedure, used when I had x,y swapped and was getting nonsensical results
1040	//Macro CalculateChiSquared(str)
1041	// String str
1042	// Prompt str,"Pick the data folder containing 2D data",popup,getAList(4)
1043	//
1044	//
1045	// String fldrSav0= GetDataFolder(1)
1046	// SetDataFolder $("root:"+str)
1047	//
1048	// Duplicate/O $(str+"_i") chi
1049	// chi = ((zwave_cyl2D_D - $(str+"_i"))/sw )^2
1050	//
1051	// chi = (mask == 1) ? chi : 0
1052	//
1053	// Print sum(chi,-inf,inf)
1054	//
1055	// SetDataFolder fldrSav0
1056	//EndMacro
1057
1058	// NaN or zero == value not used
1059	// any non-zero data == value used for fit (I use 1)
1060	Function MakeBSMask(mask,rad)
1061	Wave mask
1062	Variable rad
1063
1064
1065	// find the center based on the wave scaling
1066	Variable xCtr,yCtr,Qzero=0
1067	xCtr = (Qzero - DimOffset(Mask, 0))/DimDelta(Mask,0)
1068	yCtr = (Qzero - DimOffset(Mask, 1))/DimDelta(Mask,1)
1069
1070	Print xctr,yctr
1071
1072	// Variable center = sqrt(numpnts(mask))/2 -0.5
1073
1074	mask = (sqrt((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1
1075
1076	Variable xDim, yDim
1077	xDim = DimSize(mask,0)
1078	yDim = DimSize(mask,1)
1079	mask[][0] = 0
1080	mask[][1] = 0
1081	mask[][yDim-2] = 0
1082	mask[][yDim-1] = 0
1083	mask[0][] = 0
1084	mask[1][] = 0
1085	mask[xDim-2][] = 0
1086	mask[xDim-1][] = 0
1087
1088	Redimension/N=(xDim*yDim) mask //now 1D
1089
1090
1091	End
1092
1093	// This routine assumes that the 2D data was loaded with the NCNR loader, so that the
1094	// data is in a data folder, and the extensions are known. A more generic form could
1095	// be made too, if needed.
1096	//
1097	// X- need error on I(q) - now calculated both ways, from 2D error propagation, and from deviation from the mean
1098	// X- need to set "proper" number of data points (delta and qMax?)
1099	// X- need to remove points at high Q end
1100	//
1101	// -- like the routines in CircSectAve, start with 500 points, and trim after binning is done.
1102	// you'l end up with < 200 points.
1103	//
1104	// the results are in iBin_qxqy, qBin_qxqy, and eBin_qxqy, in the folder passed
1105	//
1106
1107	Proc BinQxQy_to_1D(folderStr)
1108	String folderStr
1109	Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
1110
1111	fDoBinning_QxQy2D(folderStr)
1112
1113	SetDataFolder $("root:"+folderStr)
1114	Display iBin_qxqy vs qBin_qxqy
1115	ModifyGraph mirror=2,grid=1,log=1
1116	ModifyGraph mode=4,marker=19,msize=2
1117	ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy)
1118	legend
1119
1120	SetDataFolder root:
1121	End
1122
1123
1124
1125	//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
1126	// Wave inten,qx,qy,qz
1127
1128	Function fDoBinning_QxQy2D(folderStr)
1129	String folderStr
1130
1131
1132	SetDataFolder $("root:"+folderStr)
1133
1134	// WAVE inten = $("smeared_sf2D")
1135
1136	WAVE inten = $(folderStr + "_i")
1137	WAVE iErr = $(folderStr + "_iErr")
1138	WAVE qx = $(folderStr + "_qx")
1139	WAVE qy = $(folderStr + "_qy")
1140	WAVE qz = $(folderStr + "_qz")
1141
1142	Variable xDim=numpnts(qx),yDim
1143	Variable ii,jj,delQ
1144	Variable qTot,nq,var,avesq,aveisq
1145	Variable binIndex,val
1146
1147	nq = 500
1148
1149	yDim = XDim
1150	Make/O/D/N=(nq) iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1151	delQ = abs(sqrt(qx[2]^2+qy[2]^2+qz[2]^2) - sqrt(qx[1]^2+qy[1]^2+qz[1]^2)) //use bins of 1 pixel width
1152	qBin_qxqy[] = p* delQ
1153	SetScale/P x,0,delQ,"",qBin_qxqy //allows easy binning
1154
1155	iBin_qxqy = 0
1156	iBin2_qxqy = 0
1157	eBin_qxqy = 0
1158	eBin2D_qxqy = 0
1159	nBin_qxqy = 0 //number of intensities added to each bin
1160
1161	for(ii=0;ii<xDim;ii+=1)
1162	qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1163	binIndex = trunc(x2pnt(qBin_qxqy, qTot))
1164	val = inten[ii]
1165	if (numType(val)==0) //count only the good points, ignore Nan or Inf
1166	iBin_qxqy[binIndex] += val
1167	iBin2_qxqy[binIndex] += val*val
1168	eBin2D_qxqy[binIndex] += iErr[ii]*iErr[ii]
1169	nBin_qxqy[binIndex] += 1
1170	endif
1171	endfor
1172
1173	//calculate errors, just like in CircSectAve.ipf
1174	for(ii=0;ii<nq;ii+=1)
1175	if(nBin_qxqy[ii] == 0)
1176	//no pixels in annuli, data unknown
1177	iBin_qxqy[ii] = 0
1178	eBin_qxqy[ii] = 1
1179	eBin2D_qxqy[ii] = NaN
1180	else
1181	if(nBin_qxqy[ii] <= 1)
1182	//need more than one pixel to determine error
1183	iBin_qxqy[ii] /= nBin_qxqy[ii]
1184	eBin_qxqy[ii] = 1
1185	eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1186	else
1187	//assume that the intensity in each pixel in annuli is normally
1188	// distributed about mean...
1189	iBin_qxqy[ii] /= nBin_qxqy[ii]
1190	avesq = iBin_qxqy[ii]^2
1191	aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1192	var = aveisq-avesq
1193	if(var<=0)
1194	eBin_qxqy[ii] = 1e-6
1195	else
1196	eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1197	endif
1198	// and calculate as it is propagated pixel-by-pixel
1199	eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1200	endif
1201	endif
1202	endfor
1203
1204	eBin2D_qxqy = sqrt(eBin2D_qxqy) // as equation (3) of John's memo
1205
1206	// find the last non-zero point, working backwards
1207	val=nq
1208	do
1209	val -= 1
1210	while(nBin_qxqy[val] == 0)
1211
1212	// print val, nBin_qxqy[val]
1213	DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1214
1215	SetDataFolder root:
1216
1217	return(0)
1218	End

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

Download in other formats:

Original Format