Context Navigation

U_CALC.ipf @ 561

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

Added UCALC window to list of windows is package loader

Box Sum (from the marquee) asks for normalized (SAM) or RAW data

SASCALC uses real beamstop size rather than projected sizewhen doing simulation

Fake1DDataFolder() in MultiScatter_MonteCarlo ipf bug fixed where a bogus beamstop diameter (not 1"-4") would cause oddities in the resolution wave, and an (obviously) incorrect smearing of the simulation.

Fixed range 7 defaults in UCALC panel.

File size: 51.5 KB

Line
1	#pragma rtGlobals=1 // Use modern global access method.
2
3
4	// USANS version of SASCALC
5
6	// to simulate the intensity from a USANS experiment for planning
7	// see John's instrument paper:
8	//
9	// J. Appl. Cryst. (2005) 38 1004-1011.
10	//
11	// SRK JUL 2009
12
13
14	// ideas:
15	//
16	// - more presets
17	//
18	// - NEED -
19	// - printable output that makes sense to an instrument scientist at least
20	// for instrument setup
21	//
22	// - make sure that the # points <-> increment relation meshes with ICP
23
24	//
25	// X plot as countrate, not absolute scale
26	// X 3e-5 cutoff
27	// X ? don't plot lowest angle range (but needs to be in the count time)
28	// - need empty beam and empty cell count rate vs. aperture (Cd vs. Gd?)
29	//
30
31
32	// need to add in empty and background corrections to see "reduced" data
33	// or at least compare to what the empty cell would give in the same count time
34	//
35	// -? model the direct beam?? currently the "red" region from -1 to 0.6 is almost entirely
36	// the primary beam, so it's a bit artificial (qmin is really ~ 3e-5)
37	//
38
39	//
40	// Need T_wide, T_rock, I peak for proper absolute scaling, but I don't know if it's really important
41	// to be able to simlulate to this extent
42	//
43
44
45	//#include "MultScatter_MonteCarlo_2D"
46
47
48
49	// Bring the UCALC panel to the front
50	// ALWAYS initializes the folders and variables
51	// then draws the panel if necessary
52	Proc Show_UCALC()
53
54	Init_UCALC()
55	DoWindow/F UCALC
56	if(V_Flag==0)
57	UCALC_Panel()
58	CalcTotalCountTime()
59	ControlUpdate/W=UCALC U_popup0 //force a pop of the function list
60	Endif
61
62	End
63
64
65	//set up the global values for the angles, # points, and count time
66	Proc Init_UCALC()
67
68	//
69	NewDataFolder/O root:Simulation // required to calculate the RandomDeviate
70	NewDataFolder/O root:Packages:NIST:SAS // required to calculate the RandomDeviate
71
72	//
73	NewDataFolder/O root:Packages:NIST:USANS:SIM //for the fake raw data
74	NewDataFolder/O/S root:Packages:NIST:USANS:Globals:U_Sim //for constants, panel, etc
75
76	Variable/G gAngLow1 = -1
77	Variable/G gAngHigh1 = 0.6
78	Variable/G gNumPts1 = 33
79	Variable/G gCtTime1 = 25
80	Variable/G gIncr1 = 0.05
81
82	Variable/G gAngLow2 = 0.7
83	Variable/G gAngHigh2 = 1.9
84	Variable/G gNumPts2 = 13
85	Variable/G gCtTime2 = 100
86	Variable/G gIncr2 = 0.1
87
88	Variable/G gAngLow3 = 2
89	Variable/G gAngHigh3 = 4.8
90	Variable/G gNumPts3 = 15
91	Variable/G gCtTime3 = 300
92	Variable/G gIncr3 = 0.2
93
94	Variable/G gAngLow4 = 5
95	Variable/G gAngHigh4 = 9.5
96	Variable/G gNumPts4 = 10
97	Variable/G gCtTime4 = 600
98	Variable/G gIncr4 = 0.5
99
100	Variable/G gAngLow5 = 10
101	Variable/G gAngHigh5 = 19
102	Variable/G gNumPts5 = 10
103	Variable/G gCtTime5 = 1200
104	Variable/G gIncr5 = 1
105
106	Variable/G gAngLow6 = 20
107	Variable/G gAngHigh6 = 48
108	Variable/G gNumPts6 = 15
109	Variable/G gCtTime6 = 2000
110	Variable/G gIncr6 = 2
111
112	Variable/G gAngLow7 = 50
113	Variable/G gAngHigh7 = 95
114	Variable/G gNumPts7 = 10
115	Variable/G gCtTime7 = 0
116	Variable/G gIncr7 = 3000
117
118	// results, setup values
119	String/G gFuncStr=""
120	String/G gTotTimeStr=""
121	Variable/G gAnalyzerOmega = 7.1e-7 //solid angle of the analyzer, in steradians
122	Variable/G gBeamCurrent=25000 //beam current Ed*I (n/s) for 5/8" diam = 25000 n/s
123	Variable/G gThick=0.1 //sample thickness (cm)
124	Variable/G gSamTrans=0.8
125	Variable/G g_1D_DoABS = 0 //=1 for abs scale, 0 for just counts
126	Variable/G g_1D_PlotCR = 1 //=1 to plot countrate
127	Variable/G g_1D_AddNoise = 1 // add in appropriate noise to simulation
128
129	// a box for the results
130	// ?? maybe not useful to report
131	Variable/G g_1DEstDetCR = 0 // estimated detector count rate
132	Variable/G g_1DTotCts = 0
133	Variable/G g_1DFracScatt= 0 // ??
134	Variable/G g_1DEstTrans = 0 // ? can I calculate this?
135
136	// not on panel yet
137	Variable/G g_Empirical_EMP = 0 // use an emperical model for empty cell subtraction
138	Variable/G g_EmptyLevel = 0.7
139	Variable/G g_BkgLevel = 0.6
140
141	SetDataFolder root:
142
143	End
144
145	//// make the display panel a graph with a control bar just as in SASCALC
146	//// so that the subwindow syntax doesn't break all of the other functionality
147	////
148	//Window UCALC_Panel() : Graph
149	// PauseUpdate; Silent 1 // building window...
150	// Display /W=(55,44,670,850) /K=1
151	// ModifyGraph cbRGB=(36929,50412,31845)
152	// DoWindow/C UCALC
153	// DoWindow/T UCALC,"USANS Simulation"
154	// ControlBar 320
155	//
156	// GroupBox group0,pos={5,0},size={577,159},title="Instrument Setup"
157	// GroupBox group1,pos={5,165},size={240,147},title="Sample Setup"
158	// GroupBox group2,pos={327,165},size={259,147},title="Results"
159	//
160	// PopupMenu popup0,pos={17,18},size={165,20},title="Sample Aperture Diam (in)"
161	// PopupMenu popup0,mode=3,popvalue="0.625",value="0.25;0.50;0.625;0.75;1.0;1.75;2.0;"
162	// PopupMenu popup2,pos={220,18},size={165,20},title="Presets"
163	// PopupMenu popup2,mode=3,popvalue="Long Count",value="Short Count;Medium Count;Long Count;"
164	// PopupMenu popup2,proc=UCALC_PresetPopup
165	//
166	// SetDataFolder root:Packages:NIST:USANS:Globals:U_Sim
167	//
168	// Variable top=44,pt=0,inc=18
169	// SetVariable setvar1a,pos={12,top},size={100,15},title="theta min",value= gAngLow1
170	// SetVariable setvar1b,pos={119,top},size={100,15},title="theta max",value= gAngHigh1
171	// SetVariable setvar1c,pos={227,top},size={100,15},title="increm",value= gIncr1
172	// SetVariable setvar1d,pos={335,top},size={100,15},title="# points",value= gNumPts1
173	// SetVariable setvar1e,pos={443,top},size={100,15},title="count (s)",value= gCtTime1
174	// SetVariable setvar1a,labelBack=(65535,32768,32768)
175	//
176	// pt += inc
177	// SetVariable setvar2a,pos={12,top+pt},size={100,15},title="theta min",value= gAngLow2
178	// SetVariable setvar2b,pos={119,top+pt},size={100,15},title="theta max",value= gAngHigh2
179	// SetVariable setvar2c,pos={227,top+pt},size={100,15},title="increm",value= gIncr2
180	// SetVariable setvar2d,pos={335,top+pt},size={100,15},title="# points",value= gNumPts2
181	// SetVariable setvar2e,pos={443,top+pt},size={100,15},title="count (s)",value= gCtTime2
182	// SetVariable setvar2a labelBack=(65535,65533,32768)
183	//
184	// pt += inc
185	// SetVariable setvar3a,pos={12,top+pt},size={100,15},title="theta min",value= gAngLow3
186	// SetVariable setvar3b,pos={119,top+pt},size={100,15},title="theta max",value= gAngHigh3
187	// SetVariable setvar3c,pos={227,top+pt},size={100,15},title="increm",value= gIncr3
188	// SetVariable setvar3d,pos={335,top+pt},size={100,15},title="# points",value= gNumPts3
189	// SetVariable setvar3e,pos={443,top+pt},size={100,15},title="count (s)",value= gCtTime3
190	// SetVariable setvar3a labelBack=(32769,65535,32768)
191	//
192	// pt += inc
193	// SetVariable setvar4a,pos={12,top+pt},size={100,15},title="theta min",value= gAngLow4
194	// SetVariable setvar4b,pos={119,top+pt},size={100,15},title="theta max",value= gAngHigh4
195	// SetVariable setvar4c,pos={227,top+pt},size={100,15},title="increm",value= gIncr4
196	// SetVariable setvar4d,pos={335,top+pt},size={100,15},title="# points",value= gNumPts4
197	// SetVariable setvar4e,pos={443,top+pt},size={100,15},title="count (s)",value= gCtTime4
198	// SetVariable setvar4a labelBack=(32768,65535,65535)
199	//
200	// pt += inc
201	// SetVariable setvar5a,pos={12,top+pt},size={100,15},title="theta min",value= gAngLow5
202	// SetVariable setvar5b,pos={119,top+pt},size={100,15},title="theta max",value= gAngHigh5
203	// SetVariable setvar5c,pos={227,top+pt},size={100,15},title="increm",value= gIncr5
204	// SetVariable setvar5d,pos={335,top+pt},size={100,15},title="# points",value= gNumPts5
205	// SetVariable setvar5e,pos={443,top+pt},size={100,15},title="count (s)",value= gCtTime5
206	// SetVariable setvar5a labelBack=(32768,54615,65535)
207	//
208	// pt += inc
209	// SetVariable setvar6a,pos={12,top+pt},size={100,15},title="theta min",value= gAngLow6
210	// SetVariable setvar6b,pos={119,top+pt},size={100,15},title="theta max",value= gAngHigh6
211	// SetVariable setvar6c,pos={227,top+pt},size={100,15},title="increm",value= gIncr6
212	// SetVariable setvar6d,pos={335,top+pt},size={100,15},title="# points",value= gNumPts6
213	// SetVariable setvar6e,pos={443,top+pt},size={100,15},title="count (s)",value= gCtTime6
214	// SetVariable setvar6a labelBack=(44253,29492,58982)
215	//
216	//// the action procedures and limits/increments
217	// SetVariable setvar1a proc=ThetaMinSetVarProc //,limits={-2,0,0.1}
218	// SetVariable setvar2a proc=ThetaMinSetVarProc
219	// SetVariable setvar3a proc=ThetaMinSetVarProc
220	// SetVariable setvar4a proc=ThetaMinSetVarProc
221	// SetVariable setvar5a proc=ThetaMinSetVarProc
222	// SetVariable setvar6a proc=ThetaMinSetVarProc
223	//
224	////
225	// SetVariable setvar1b proc=ThetaMaxSetVarProc //,limits={0.4,1,0.1}
226	// SetVariable setvar2b proc=ThetaMaxSetVarProc
227	// SetVariable setvar3b proc=ThetaMaxSetVarProc
228	// SetVariable setvar4b proc=ThetaMaxSetVarProc
229	// SetVariable setvar5b proc=ThetaMaxSetVarProc
230	// SetVariable setvar6b proc=ThetaMaxSetVarProc
231	////
232	// SetVariable setvar1c proc=IncrSetVarProc,limits={0.01,0.1,0.01}
233	// SetVariable setvar2c proc=IncrSetVarProc,limits={0.02,0.2,0.02}
234	// SetVariable setvar3c proc=IncrSetVarProc,limits={0.05,0.4,0.05}
235	// SetVariable setvar4c proc=IncrSetVarProc,limits={0.1,1,0.1}
236	// SetVariable setvar5c proc=IncrSetVarProc,limits={0.5,5,1}
237	// SetVariable setvar6c proc=IncrSetVarProc,limits={1,10,2}
238	////
239	// SetVariable setvar1d proc=NumPtsSetVarProc,limits={2,50,1}
240	// SetVariable setvar2d proc=NumPtsSetVarProc,limits={2,50,1}
241	// SetVariable setvar3d proc=NumPtsSetVarProc,limits={2,50,1}
242	// SetVariable setvar4d proc=NumPtsSetVarProc,limits={2,50,1}
243	// SetVariable setvar5d proc=NumPtsSetVarProc,limits={2,50,1}
244	// SetVariable setvar6d proc=NumPtsSetVarProc,limits={2,50,1}
245	////
246	// SetVariable setvar1e proc=CtTimeSetVarProc,limits={-1,50000,1}
247	// SetVariable setvar2e proc=CtTimeSetVarProc,limits={-1,50000,10}
248	// SetVariable setvar3e proc=CtTimeSetVarProc,limits={-1,50000,10}
249	// SetVariable setvar4e proc=CtTimeSetVarProc,limits={-1,50000,30}
250	// SetVariable setvar5e proc=CtTimeSetVarProc,limits={-1,50000,100}
251	// SetVariable setvar6e proc=CtTimeSetVarProc,limits={-1,50000,100}
252	//
253	// Button button0,pos={255,180},size={60,20},fColor=(65535,65535,0),proc=U_SimPlotButtonProc,title="Plot"
254	// Button button1,pos={260,286},size={50,20},proc=U_SaveButtonProc,title="Save"
255	//
256	////checkbox for "easy" mode
257	// CheckBox check0 title="Simple mode?",pos={400,19},proc=EnterModeCheckProc,value=1
258	// ThetaEditMode(2) //checked on startup
259	//
260	//// instrument setup
261	// SetVariable U_setvar0_1,pos={20,211},size={160,15},title="Thickness (cm)"
262	// SetVariable U_setvar0_1,limits={0,inf,0.1},value= root:Packages:NIST:USANS:Globals:U_Sim:gThick
263	// SetVariable U_setvar0_3,pos={20,235},size={160,15},title="Sample Transmission"
264	// SetVariable U_setvar0_3,limits={0,1,0.01},value= root:Packages:NIST:USANS:Globals:U_Sim:gSamTrans
265	// PopupMenu U_popup0,pos={20,185},size={165,20},proc=Sim_USANS_ModelPopMenuProc,title="Model"
266	// PopupMenu U_popup0,mode=1,value= #"U_FunctionPopupList()"
267	// SetVariable setvar0,pos={20,259},size={120,15},title="Empty Level"
268	// SetVariable setvar0,limits={0,10,0.01},value= root:Packages:NIST:USANS:Globals:U_Sim:g_EmptyLevel
269	// SetVariable setvar0_1,pos={20,284},size={120,15},title="Bkg Level"
270	// SetVariable setvar0_1,limits={0,10,0.01},value= root:Packages:NIST:USANS:Globals:U_Sim:g_BkgLevel
271	//
272	// CheckBox check0_4 title="Show EMP?",pos={160,260},proc=ShowEMPCheckProc,value=0
273	//
274	// CheckBox check0_2,pos={253,239},size={60,14},title="CountRate?",variable= root:Packages:NIST:USANS:Globals:U_Sim:g_1D_PlotCR
275	// CheckBox check0_3,pos={262,264},size={60,14},title="Noise?",variable= root:Packages:NIST:USANS:Globals:U_Sim:g_1D_AddNoise
276	//
277	//// a box for the results
278	// SetVariable totalTime,pos={338,185},size={150,15},title="Count time (h:m)",value= gTotTimeStr
279	//// ValDisplay valdisp0,pos={338,210},size={220,13},title="Total detector counts"
280	//// ValDisplay valdisp0,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:USANS:Globals:U_Sim:g_1DTotCts
281	// ValDisplay valdisp0_2,pos={338,234},size={220,13},title="Fraction of beam scattered"
282	// ValDisplay valdisp0_2,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:USANS:Globals:U_Sim:g_1DFracScatt
283	// ValDisplay valdisp0_3,pos={338,259},size={220,13},title="Estimated transmission"
284	// ValDisplay valdisp0_3,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:USANS:Globals:U_Sim:g_1DEstTrans
285	//
286	//
287	// SetDataFolder root:
288	//
289	//EndMacro
290
291	// ??make the USANS simulation results into a graph and a separate panel.
292	// the control bar at the top makes the whole thing too large, and control bars are limited
293	// to 500 pix wide, which is really tight
294	// left of 1055 is good for Mac, 955 is better for Win
295	//
296	Window UCALC_Panel() : Graph
297	PauseUpdate; Silent 1 // building window...
298	String platform=UpperStr(IgorInfo(2))
299	Variable pos=strsearch(platform,"WINDOWS",0)
300	if(pos >= 0) //windows
301	Display /W=(55,44,855,450) /K=1
302	else //mac
303	Display /W=(55,44,1055,544) /K=1
304	endif
305
306	ModifyGraph cbRGB=(36929,50412,31845)
307	DoWindow/C UCALC
308	DoWindow/T UCALC,"USANS Simulation"
309	ControlBar/L 500
310
311	GroupBox group0,pos={5,1},size={493,177},title="Instrument Setup"
312	GroupBox group1,pos={5,183},size={240,147},title="Sample Setup"
313	GroupBox group2,pos={5,343},size={259,147},title="Results"
314
315	PopupMenu popup0,pos={17,19},size={165,20},title="Sample Aperture Diam (in)"
316	PopupMenu popup0,mode=3,popvalue="0.625",value="0.25;0.50;0.625;0.75;1.0;1.75;2.0;"
317	PopupMenu popup0,proc=UCALC_SampleAperturePopup
318	PopupMenu popup2,pos={220,19},size={165,20},title="Presets"
319	PopupMenu popup2,mode=3,popvalue="Long Count",value="Short Count;Medium Count;Long Count;"
320	PopupMenu popup2,proc=UCALC_PresetPopup
321
322	SetDataFolder root:Packages:NIST:USANS:Globals:U_Sim
323
324	Variable top=46,pt=0,inc=18,left=0//left=533
325	SetVariable setvar1a,pos={left+17,top},size={90,15},title="theta min",value= gAngLow1
326	SetVariable setvar1b,pos={left+113,top},size={89,15},title="theta max",value= gAngHigh1
327	SetVariable setvar1c,pos={left+209,top},size={89,15},title="increm",value= gIncr1
328	SetVariable setvar1d,pos={left+299,top},size={100,15},title="# points",value= gNumPts1
329	SetVariable setvar1e,pos={left+399,top},size={93,15},title="count (s)",value= gCtTime1
330	// SetVariable setvar1a,labelBack=(65535,32768,32768) //old rainbow
331	SetVariable setvar1a,labelBack=(49858,65535,65535)
332
333	pt += inc
334	SetVariable setvar2a,pos={left+17,top+pt},size={90,15},title="theta min",value= gAngLow2
335	SetVariable setvar2b,pos={left+113,top+pt},size={89,15},title="theta max",value= gAngHigh2
336	SetVariable setvar2c,pos={left+209,top+pt},size={89,15},title="increm",value= gIncr2
337	SetVariable setvar2d,pos={left+299,top+pt},size={100,15},title="# points",value= gNumPts2
338	SetVariable setvar2e,pos={left+399,top+pt},size={93,15},title="count (s)",value= gCtTime2
339	// SetVariable setvar2a labelBack=(65535,65533,32768) //old rainbow
340	SetVariable setvar2a labelBack=(21074,8995,21074)
341
342	pt += inc
343	SetVariable setvar3a,pos={left+17,top+pt},size={90,15},title="theta min",value= gAngLow3
344	SetVariable setvar3b,pos={left+113,top+pt},size={89,15},title="theta max",value= gAngHigh3
345	SetVariable setvar3c,pos={left+209,top+pt},size={89,15},title="increm",value= gIncr3
346	SetVariable setvar3d,pos={left+299,top+pt},size={100,15},title="# points",value= gNumPts3
347	SetVariable setvar3e,pos={left+399,top+pt},size={93,15},title="count (s)",value= gCtTime3
348	// SetVariable setvar3a labelBack=(32769,65535,32768) //old rainbow
349	SetVariable setvar3a labelBack=(0,60652,60652)
350
351	pt += inc
352	SetVariable setvar4a,pos={left+17,top+pt},size={90,15},title="theta min",value= gAngLow4
353	SetVariable setvar4b,pos={left+113,top+pt},size={89,15},title="theta max",value= gAngHigh4
354	SetVariable setvar4c,pos={left+209,top+pt},size={89,15},title="increm",value= gIncr4
355	SetVariable setvar4d,pos={left+299,top+pt},size={100,15},title="# points",value= gNumPts4
356	SetVariable setvar4e,pos={left+399,top+pt},size={93,15},title="count (s)",value= gCtTime4
357	// SetVariable setvar4a labelBack=(32768,65535,65535) //old rainbow
358	SetVariable setvar4a labelBack=(0,51400,0)
359
360	pt += inc
361	SetVariable setvar5a,pos={left+17,top+pt},size={90,15},title="theta min",value= gAngLow5
362	SetVariable setvar5b,pos={left+113,top+pt},size={89,15},title="theta max",value= gAngHigh5
363	SetVariable setvar5c,pos={left+209,top+pt},size={89,15},title="increm",value= gIncr5
364	SetVariable setvar5d,pos={left+299,top+pt},size={100,15},title="# points",value= gNumPts5
365	SetVariable setvar5e,pos={left+399,top+pt},size={93,15},title="count (s)",value= gCtTime5
366	// SetVariable setvar5a labelBack=(32768,54615,65535) //old rainbow
367	SetVariable setvar5a labelBack=(59367,49344,0)
368
369	pt += inc
370	SetVariable setvar6a,pos={left+17,top+pt},size={90,15},title="theta min",value= gAngLow6
371	SetVariable setvar6b,pos={left+113,top+pt},size={89,15},title="theta max",value= gAngHigh6
372	SetVariable setvar6c,pos={left+209,top+pt},size={89,15},title="increm",value= gIncr6
373	SetVariable setvar6d,pos={left+299,top+pt},size={100,15},title="# points",value= gNumPts6
374	SetVariable setvar6e,pos={left+399,top+pt},size={93,15},title="count (s)",value= gCtTime6
375	// SetVariable setvar6a labelBack=(44253,29492,58982) //old rainbow
376	SetVariable setvar6a labelBack=(54998,0,0)
377
378	pt += inc
379	SetVariable setvar7a,pos={left+17,top+pt},size={90,15},title="theta min",value= gAngLow7
380	SetVariable setvar7b,pos={left+113,top+pt},size={89,15},title="theta max",value= gAngHigh7
381	SetVariable setvar7c,pos={left+209,top+pt},size={89,15},title="increm",value= gIncr7
382	SetVariable setvar7d,pos={left+299,top+pt},size={100,15},title="# points",value= gNumPts7
383	SetVariable setvar7e,pos={left+399,top+pt},size={93,15},title="count (s)",value= gCtTime7
384	// SetVariable setvar7a labelBack=(44253,29492,58982) //old rainbow
385	SetVariable setvar7a labelBack=(39321,21845,51657)
386
387	// the action procedures and limits/increments
388	SetVariable setvar1a proc=ThetaMinSetVarProc //,limits={-2,0,0.1}
389	SetVariable setvar2a proc=ThetaMinSetVarProc
390	SetVariable setvar3a proc=ThetaMinSetVarProc
391	SetVariable setvar4a proc=ThetaMinSetVarProc
392	SetVariable setvar5a proc=ThetaMinSetVarProc
393	SetVariable setvar6a proc=ThetaMinSetVarProc
394	SetVariable setvar7a proc=ThetaMinSetVarProc
395
396	//
397	SetVariable setvar1b proc=ThetaMaxSetVarProc //,limits={0.4,1,0.1}
398	SetVariable setvar2b proc=ThetaMaxSetVarProc
399	SetVariable setvar3b proc=ThetaMaxSetVarProc
400	SetVariable setvar4b proc=ThetaMaxSetVarProc
401	SetVariable setvar5b proc=ThetaMaxSetVarProc
402	SetVariable setvar6b proc=ThetaMaxSetVarProc
403	SetVariable setvar7b proc=ThetaMaxSetVarProc
404	//
405	SetVariable setvar1c proc=IncrSetVarProc,limits={0.01,0.1,0.01}
406	SetVariable setvar2c proc=IncrSetVarProc,limits={0.02,0.2,0.02}
407	SetVariable setvar3c proc=IncrSetVarProc,limits={0.05,0.4,0.05}
408	SetVariable setvar4c proc=IncrSetVarProc,limits={0.1,1,0.1}
409	SetVariable setvar5c proc=IncrSetVarProc,limits={0.5,5,1}
410	SetVariable setvar6c proc=IncrSetVarProc,limits={1,10,2}
411	SetVariable setvar7c proc=IncrSetVarProc,limits={1,20,2}
412	//
413	SetVariable setvar1d proc=NumPtsSetVarProc,limits={2,50,1}
414	SetVariable setvar2d proc=NumPtsSetVarProc,limits={2,50,1}
415	SetVariable setvar3d proc=NumPtsSetVarProc,limits={2,50,1}
416	SetVariable setvar4d proc=NumPtsSetVarProc,limits={2,50,1}
417	SetVariable setvar5d proc=NumPtsSetVarProc,limits={2,50,1}
418	SetVariable setvar6d proc=NumPtsSetVarProc,limits={2,50,1}
419	SetVariable setvar7d proc=NumPtsSetVarProc,limits={2,50,1}
420	//
421	SetVariable setvar1e proc=CtTimeSetVarProc,limits={-1,50000,1}
422	SetVariable setvar2e proc=CtTimeSetVarProc,limits={-1,50000,10}
423	SetVariable setvar3e proc=CtTimeSetVarProc,limits={-1,50000,10}
424	SetVariable setvar4e proc=CtTimeSetVarProc,limits={-1,50000,30}
425	SetVariable setvar5e proc=CtTimeSetVarProc,limits={-1,50000,100}
426	SetVariable setvar6e proc=CtTimeSetVarProc,limits={-1,50000,100}
427	SetVariable setvar7e proc=CtTimeSetVarProc,limits={-1,50000,100}
428
429	Button button0,pos={left+280,200},size={130,20},fColor=(65535,65535,0),proc=U_SimPlotButtonProc,title="Simulate USANS"
430	CheckBox check0_2,pos={left+280,250},size={60,14},title="CountRate?",variable= root:Packages:NIST:USANS:Globals:U_Sim:g_1D_PlotCR
431	CheckBox check0_3,pos={left+280,270},size={60,14},title="Noise?",variable= root:Packages:NIST:USANS:Globals:U_Sim:g_1D_AddNoise
432	CheckBox check0_4 title="Show EMP?",pos={left+280,290},proc=ShowEMPCheckProc,value=0
433
434
435	//checkbox for "easy" mode
436	CheckBox check0 title="Simple mode?",pos={left+400,19},proc=EnterModeCheckProc,value=1
437	ThetaEditMode(2) //checked on startup
438
439	// sample setup
440	SetVariable U_setvar0_1,pos={left+20,231},size={160,15},title="Thickness (cm)"
441	SetVariable U_setvar0_1,limits={0,inf,0.1},value= root:Packages:NIST:USANS:Globals:U_Sim:gThick
442	SetVariable U_setvar0_3,pos={left+20,255},size={160,15},title="Sample Transmission"
443	SetVariable U_setvar0_3,limits={0,1,0.01},value= root:Packages:NIST:USANS:Globals:U_Sim:gSamTrans
444	PopupMenu U_popup0,pos={left+20,205},size={165,20},proc=Sim_USANS_ModelPopMenuProc,title="Model"
445	PopupMenu U_popup0,mode=1,value= #"U_FunctionPopupList()"
446	SetVariable setvar0,pos={left+20,279},size={120,15},title="Empty Level",disable=2
447	SetVariable setvar0,limits={0,10,0.01},value= root:Packages:NIST:USANS:Globals:U_Sim:g_EmptyLevel
448	SetVariable setvar0_1,pos={left+20,304},size={120,15},title="Bkg Level",disable=2
449	SetVariable setvar0_1,limits={0,10,0.01},value= root:Packages:NIST:USANS:Globals:U_Sim:g_BkgLevel
450
451
452
453	// a box for the results
454	SetVariable totalTime,pos={left+20,370},size={150,15},title="Count time (h:m)",value= gTotTimeStr
455	ValDisplay valdisp0_2,pos={left+20,395},size={220,13},title="Fraction of beam scattered"
456	ValDisplay valdisp0_2,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:USANS:Globals:U_Sim:g_1DFracScatt
457	ValDisplay valdisp0_2,disable=1
458	ValDisplay valdisp0_3,pos={left+20,420},size={220,13},title="Estimated transmission"
459	ValDisplay valdisp0_3,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:USANS:Globals:U_Sim:g_1DEstTrans
460	ValDisplay valdisp0_3,disable=1
461	Button button1,pos={left+20,440},size={150,20},proc=U_SaveButtonProc,title="Save Simulated Data"
462
463
464	SetDataFolder root:
465
466	EndMacro
467
468
469	// changing theta min - hold incr and #, result is new theta max
470	Function ThetaMinSetVarProc(sva) : SetVariableControl
471	STRUCT WMSetVariableAction &sva
472
473	switch( sva.eventCode )
474	case 1: // mouse up
475	case 2: // Enter key
476	case 3: // Live update
477	Variable ThetaMin = sva.dval
478	String ns=CtrlNumber(sva.ctrlName) //control number as a string
479	NVAR incr = $("root:Packages:NIST:USANS:Globals:U_Sim:"+"gIncr"+ns)
480	NVAR NumPts = $("root:Packages:NIST:USANS:Globals:U_Sim:"+"gNumPts"+ns)
481	NVAR thetaMax = $("root:Packages:NIST:USANS:Globals:U_Sim:"+"gAngHigh"+ns)
482
483	thetaMax = ThetaMin + incr*NumPts
484
485	break
486	endswitch
487
488	return 0
489	End
490
491
492	// changing theta max - hold min and incr, result is new # of points
493	// then need to recalculate the total counting time
494	Function ThetaMaxSetVarProc(sva) : SetVariableControl
495	STRUCT WMSetVariableAction &sva
496
497	switch( sva.eventCode )
498	case 1: // mouse up
499	case 2: // Enter key
500	case 3: // Live update
501	Variable ThetaMax = sva.dval
502	String ns=CtrlNumber(sva.ctrlName) //control number as a string
503	NVAR thetaMin = $("root:Packages:NIST:USANS:Globals:U_Sim:"+"gAngLow"+ns)
504	NVAR NumPts = $("root:Packages:NIST:USANS:Globals:U_Sim:"+"gNumPts"+ns)
505	NVAR Incr = $("root:Packages:NIST:USANS:Globals:U_Sim:"+"gIncr"+ns)
506
507	NumPts = trunc( (thetaMax - ThetaMin) / incr )
508
509	CalcTotalCountTime()
510
511	break
512	endswitch
513
514	return 0
515	End
516
517	// changing increment - hold min and #, result is new max
518	Function IncrSetVarProc(sva) : SetVariableControl
519	STRUCT WMSetVariableAction &sva
520
521	switch( sva.eventCode )
522	case 1: // mouse up
523	case 2: // Enter key
524	case 3: // Live update
525	Variable incr = sva.dval
526	String ns=CtrlNumber(sva.ctrlName) //control number as a string
527	NVAR thetaMin = $("root:Packages:NIST:USANS:Globals:U_Sim:"+"gAngLow"+ns)
528	NVAR NumPts = $("root:Packages:NIST:USANS:Globals:U_Sim:"+"gNumPts"+ns)
529	NVAR thetaMax = $("root:Packages:NIST:USANS:Globals:U_Sim:"+"gAngHigh"+ns)
530
531	thetaMax = ThetaMin + incr*NumPts
532
533	break
534	endswitch
535
536	return 0
537	End
538
539	// changing #pts - hold min and incr, result is new max
540	// then recalculate the total count time
541	Function NumPtsSetVarProc(sva) : SetVariableControl
542	STRUCT WMSetVariableAction &sva
543
544	switch( sva.eventCode )
545	case 1: // mouse up
546	case 2: // Enter key
547	case 3: // Live update
548	Variable NumPts = sva.dval
549	String ns=CtrlNumber(sva.ctrlName) //control number as a string
550	NVAR thetaMin = $("root:Packages:NIST:USANS:Globals:U_Sim:"+"gAngLow"+ns)
551	NVAR incr = $("root:Packages:NIST:USANS:Globals:U_Sim:"+"gIncr"+ns)
552	NVAR thetaMax = $("root:Packages:NIST:USANS:Globals:U_Sim:"+"gAngHigh"+ns)
553
554	// old way, like ICP
555	// thetaMax = ThetaMin + incr*NumPts
556
557	// new way - to spread the points out over the specified angle range
558	incr = (thetaMax - thetaMin) / numpts
559
560	CalcTotalCountTime()
561
562	break
563	endswitch
564
565	return 0
566	End
567
568	// changing count time -
569	// then recalculate the total count time
570	Function CtTimeSetVarProc(sva) : SetVariableControl
571	STRUCT WMSetVariableAction &sva
572
573	switch( sva.eventCode )
574	case 1: // mouse up
575	case 2: // Enter key
576	case 3: // Live update
577	CalcTotalCountTime()
578	break
579	endswitch
580
581	return 0
582	End
583
584	// hard-wired for 7 controls
585	// return value is in seconds
586	// global display string is set with hrs:min
587	Function CalcTotalCountTime()
588
589	Variable ii,num,totTime=0
590	String pathStr="root:Packages:NIST:USANS:Globals:U_Sim:"
591	num=7
592
593	for(ii=1;ii<=num;ii+=1)
594	NVAR ctTime = $(pathStr+"gCtTime"+num2str(ii))
595	NVAR numPts = $(pathStr+"gNumPts"+num2str(ii))
596	if(ctTime>0)
597	totTime += ctTime*numPts
598	endif
599	endfor
600	Variable hrs,mins
601	hrs = trunc(totTime/3600)
602	mins = trunc(mod(totTime,3600)/60)
603	// Printf "Counting time (hr:min) = %d:%d\r",hrs,mins
604
605	SVAR str = $(pathStr+"gTotTimeStr")
606	sprintf str,"%d:%d",hrs,mins
607
608	return(totTime)
609	End
610
611	//returns the control number from the name string
612	// all are setvarNa
613	Function/S CtrlNumber(str)
614	String str
615
616	return(str[6])
617	End
618
619	// changes edit mode of the theta min/max boxes and increment for simplified setup
620	// val = 2 = disable
621	// val = 0 = edit enabled
622	Function ThetaEditMode(val)
623	Variable val
624
625	SetVariable setvar1a,win=UCALC,disable=val
626	SetVariable setvar1b,win=UCALC,disable=val
627	SetVariable setvar1c,win=UCALC,disable=val
628
629	SetVariable setvar2a,win=UCALC,disable=val
630	SetVariable setvar2b,win=UCALC,disable=val
631	SetVariable setvar2c,win=UCALC,disable=val
632
633	SetVariable setvar3a,win=UCALC,disable=val
634	SetVariable setvar3b,win=UCALC,disable=val
635	SetVariable setvar3c,win=UCALC,disable=val
636
637	SetVariable setvar4a,win=UCALC,disable=val
638	SetVariable setvar4b,win=UCALC,disable=val
639	SetVariable setvar4c,win=UCALC,disable=val
640
641	SetVariable setvar5a,win=UCALC,disable=val
642	SetVariable setvar5b,win=UCALC,disable=val
643	SetVariable setvar5c,win=UCALC,disable=val
644
645	SetVariable setvar6a,win=UCALC,disable=val
646	SetVariable setvar6b,win=UCALC,disable=val
647	SetVariable setvar6c,win=UCALC,disable=val
648
649	SetVariable setvar7a,win=UCALC,disable=val
650	SetVariable setvar7b,win=UCALC,disable=val
651	SetVariable setvar7c,win=UCALC,disable=val
652	return(0)
653	End
654
655
656
657	Function EnterModeCheckProc(cba) : CheckBoxControl
658	STRUCT WMCheckboxAction &cba
659
660	switch( cba.eventCode )
661	case 2: // mouse up
662	Variable checked = cba.checked
663
664	if(checked)
665	ThetaEditMode(2)
666	else
667	ThetaEditMode(0)
668	endif
669
670	break
671	endswitch
672
673	return 0
674	End
675
676
677	Function ShowEMPCheckProc(cba) : CheckBoxControl
678	STRUCT WMCheckboxAction &cba
679
680	switch( cba.eventCode )
681	case 2: // mouse up
682	Variable checked = cba.checked
683
684	String list,item,popStr,qval,CR
685	Variable OK=1
686
687	if(exists("root:Packages:NIST:USANS:Globals:Q_2p0")==0) //
688	MakeUSANSEmptyWaves()
689	endif
690
691
692	if(checked)
693	// put it on the graph
694	SetDataFolder root:Packages:NIST:USANS:Globals
695	ControlInfo/W=UCALC popup0
696	popStr = S_Value
697	strswitch(popStr) // string switch
698	case "0.25": // execute if case matches expression
699	qval = "Q_0p25"
700	CR = "CR_0p25"
701	break
702	case "0.50":
703	qval = "Q_0p50"
704	CR = "CR_0p50"
705	break
706	case "0.625":
707	qval = "Q_0p625"
708	CR = "CR_0p625"
709	break
710	case "1.75":
711	qval = "Q_1p75"
712	CR = "CR_1p75"
713	break
714	case "2.0":
715	qval = "Q_2p0"
716	CR = "CR_2p0"
717	break
718	default: // optional default expression executed
719	OK=0
720	endswitch
721
722	if(OK)
723	AppendToGraph/W=UCALC $CR vs $qval
724	ModifyGraph marker=19,mode($CR)=4,msize($CR)=3,rgb($CR)=(0,0,0)
725	endif
726
727	else
728	//take it off of the graph
729	SetDataFolder root:Packages:NIST:USANS:Globals
730	list=WaveList("CR*", ";", "WIN:UCALC")
731	item=StringFromList(0, list ,";") //should be one item
732	if(strlen(item) != 0)
733	RemoveFromGraph/W=UCALC $item
734	endif
735	endif
736
737	break
738	endswitch
739
740	SetDataFolder root:
741	return 0
742	End
743
744
745
746	// based on the angle ranges above (with non-zero count times)
747	// plot where the data points would be
748	//
749	Function U_SimPlotButtonProc(ba) : ButtonControl
750	STRUCT WMButtonAction &ba
751
752
753	switch( ba.eventCode )
754	case 2: // mouse up
755	// click code here
756
757	CalcUSANS()
758
759	break
760	endswitch
761
762	return 0
763	End
764
765
766
767	// Fills a fake USANS folder with the concatenated ranges, and converts to Q
768	// root:Packages:NIST:USANS:SIM
769	//
770	// does the work of calculating and smearing the simluated USANS data set
771	Function CalcUSANS()
772
773	Variable num,ii,firstSet=0
774	String pathStr="root:Packages:NIST:USANS:Globals:U_Sim:gCtTime"
775	String fromType="SWAP",toType="SIM"
776	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
777
778	num = 7 //# of angular ranges
779
780	// only try to plot ranges with non-zero count times
781	firstSet=0
782	for(ii=1;ii<=num;ii+=1)
783	NVAR dum = $(pathStr+num2str(ii))
784	if(dum>0)
785	// print "CtTime = ",dum
786	firstSet += 1
787	LoadSimulatedAngleRange(ii,"SWAP") //overwrite what's in the SWAP folder
788
789	// did not do this step
790	//Convert2Countrate("SWAP")
791	//
792
793	if(firstSet==1) //first time, overwrite
794	NewDataWaves("SWAP","SIM")
795	//plus my two waves
796	Duplicate/O $(USANSFolder+":"+fromType+":countingTime"),$(USANSFolder+":"+toType+":countingTime")
797	Duplicate/O $(USANSFolder+":"+fromType+":SetNumber"),$(USANSFolder+":"+toType+":SetNumber")
798
799	else //append to waves in "SIM"
800	AppendDataWaves("SWAP","SIM")
801	// and my two waves
802	ConcatenateData( (USANSFolder+":"+toType+":countingTime"),(USANSFolder+":"+fromType+":countingTime") )
803	ConcatenateData( (USANSFolder+":"+toType+":SetNumber"),(USANSFolder+":"+fromType+":SetNumber") )
804
805	endif
806
807	endif
808	endfor
809
810
811	//sort after all loaded - not by angle, but by Q
812	// Get rid of the negative angles - the smearing integration does not like these!
813	// (may add them back in later, but probably not)
814	UDoAngleSort("SIM")
815
816	ConvertAngle2Qvals("SIM",0)
817
818
819
820	//fill the data with something
821	WAVE qvals = root:Packages:NIST:USANS:SIM:qvals
822	WAVE DetCts = root:Packages:NIST:USANS:SIM:DetCts
823
824
825	//find the Trans Cts for T_Wide
826	// FindTWideCts("EMP")
827
828	//generate a "fake" 1d data folder/set named "Sim_USANS" that can be used for smearing
829	// DOES NOT calculate a matrix, but instead fills 4-5-6 columns with -dQv
830	// so that the trapezoid rule is used
831	//
832
833	FakeUSANSDataFolder(qvals,DetCts,0.117,"Sim_USANS")
834
835	// now calculate the smearing... instead of the counts
836	SVAR funcStr = root:Packages:NIST:USANS:Globals:U_Sim:gFuncStr //set by the popup
837
838	Wave inten = root:Sim_USANS:Sim_USANS_i
839	Wave sigave = root:Sim_USANS:Sim_USANS_s
840	Wave countingTime = root:Sim_USANS:countingTime //counting time per point in the set
841
842	Duplicate/O inten root:Sim_USANS:Smeared_inten // a place for the smeared result for probabilities
843	Wave Smeared_inten = root:Sim_USANS:Smeared_inten
844
845	String coefStr=""
846	Variable sig_sas=0,wavelength = 2.4
847	Variable Imon
848
849
850	NVAR omega = root:Packages:NIST:USANS:Globals:U_Sim:gAnalyzerOmega
851
852	if(exists(funcStr) != 0)
853	FUNCREF SANSModelAAO_proto func=$("fSmeared"+funcStr) //a wrapper for the structure version
854	FUNCREF SANSModelAAO_proto funcUnsmeared=$(funcStr) //unsmeared
855	coefStr = MC_getFunctionCoef(funcStr)
856
857	if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
858	Abort "Function and coefficients do not match. You must plot the unsmeared function before simulation."
859	endif
860
861	// do the smearing calculation
862	func($coefStr,Smeared_inten,qvals)
863
864
865	NVAR thick = root:Packages:NIST:USANS:Globals:U_Sim:gThick
866	NVAR trans = root:Packages:NIST:USANS:Globals:U_Sim:gSamTrans
867	NVAR SimDetCts = root:Packages:NIST:USANS:Globals:U_Sim:g_1DTotCts //summed counts (simulated)
868	NVAR estDetCR = root:Packages:NIST:USANS:Globals:U_Sim:g_1DEstDetCR // estimated detector count rate
869	NVAR fracScat = root:Packages:NIST:USANS:Globals:U_Sim:g_1DFracScatt // fraction of beam captured on detector
870	NVAR estTrans = root:Packages:NIST:USANS:Globals:U_Sim:g_1DEstTrans // estimated transmission of sample
871	// NVAR SimCountTime = root:Packages:NIST:USANS:Globals:U_Sim:gCntTime //counting time used for simulation
872
873	Imon = GetUSANSBeamIntensity() //based on the aperture size, select the beam intensity
874	Print "imon=",imon
875
876	// calculate the scattering cross section simply to be able to estimate the transmission
877	// unfortunately, this calculation is useless in the USANS range. So although it's calculated, the
878	// results are never reported.
879	// -- the main issue is that the integration range is optimized for SANS, and is not useful for USANS
880	// -there are only a few points in the USANS range, and extending it lower caused issues before.
881	CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
882
883	// if(sig_sas > 100)
884	// sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas
885	// endif
886	estTrans = exp(-1thicksig_sas) //thickness and sigma both in units of cm
887	Print "Sig_sas = ",sig_sas
888
889
890	Duplicate/O qvals prob_i
891
892	prob_i = transthickomega*Smeared_inten //probability of a neutron in q-bin(i)
893
894	// Variable P_on = sum(prob_i,-inf,inf)
895	// Print "P_on = ",P_on
896	fracScat = 1-estTrans
897
898	inten = (ImoncountingTime)prob_i
899
900	// do I round to an integer?
901	inten = round(inten)
902
903	// SimDetCts = sum(inten,-inf,inf)
904	// estDetCR = SimDetCts/SimCountTime
905
906
907	NVAR doABS = root:Packages:NIST:USANS:Globals:U_Sim:g_1D_DoABS
908	NVAR plotCR = root:Packages:NIST:USANS:Globals:U_Sim:g_1D_PlotCR
909	NVAR addNoise = root:Packages:NIST:USANS:Globals:U_Sim:g_1D_AddNoise
910
911	sigave = sqrt(inten) // assuming that N is large
912
913	// add in random error in aveint based on the sigave
914	if(addNoise)
915	inten += gnoise(sigave)
916
917	//round to an integer again
918	inten = round(inten)
919	endif
920
921	// convert to absolute scale? Maybe not needed
922	// does nothing yet - need Ipeak, Twide
923	// if(doABS)
924	// Variable kappa = thickomegatransiMonctTime
925	// inten /= kappa
926	// inten /= kappa
927	// endif
928
929	// plot as countrate - maybe easier to visualize, and all of the data overlaps
930	if(plotCR)
931	inten /= countingTime
932	sigave /= countingTime
933	endif
934
935	GraphSIM()
936
937	else
938	//no function plotted, no simulation can be done
939	DoAlert 0,"No function is selected or plotted, so no simulation is done. The default power law function is used."
940
941	inten = U_Power_Law_Model(1e-6,3,0,qvals)
942	// inten = U_SphereForm(1,9000,6e-6,0,qvals)
943
944	GraphSIM()
945
946	endif
947
948
949	end
950
951
952	//sort the data in the "type"folder, based on angle
953	//carry along all associated waves
954	//
955	// ---a duplicate of DoAngleSort(), by modified to
956	// include counting time and setNumber
957	//
958	// also trims the beginning of each data set so that it does not include any negative or zero angles
959	//
960	Function UDoAngleSort(type)
961	String type
962
963	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
964
965	Wave angle = $(USANSFolder+":"+Type+":Angle")
966	Wave detCts = $(USANSFolder+":"+Type+":DetCts")
967	Wave ErrdetCts = $(USANSFolder+":"+Type+":ErrDetCts")
968	Wave MonCts = $(USANSFolder+":"+Type+":MonCts")
969	Wave TransCts = $(USANSFolder+":"+Type+":TransCts")
970	Wave countingTime = $(USANSFolder+":"+Type+":countingTime")
971	Wave SetNumber = $(USANSFolder+":"+Type+":SetNumber")
972
973	Sort Angle DetCts,ErrDetCts,MonCts,TransCts,Angle,countingTime,SetNumber
974
975	Variable ii,num,numBad,ang,val
976	num=numpnts(angle)
977	ii=0
978	numBad=0
979	val = 0 //cutoff value
980	do
981	ang = angle[ii]
982	if(ang <= val)
983	numBad += 1
984	else //keep the points
985	Angle[ii-numBad] = ang
986	DetCts[ii-numBad] = DetCts[ii]
987	ErrDetCts[ii-numBad] = ErrDetCts[ii]
988	MonCts[ii-numBad] = MonCts[ii]
989	TransCts[ii-numBad] = TransCts[ii]
990	countingTime[ii-numBad] = countingTime[ii]
991	SetNumber[ii-numBad] = SetNumber[ii]
992	endif
993	ii += 1
994	while(ii<num)
995	//trim the end of the waves
996	DeletePoints num-numBad, numBad, DetCts,ErrDetCts,MonCts,TransCts,Angle,countingTime,SetNumber
997
998
999	return(0)
1000	End
1001
1002	// a simple default function
1003	//
1004	Function U_SphereForm(scale,radius,delrho,bkg,x)
1005	Variable scale,radius,delrho,bkg
1006	Variable x
1007
1008	// variables are:
1009	//[0] scale
1010	//[1] radius (A)
1011	//[2] delrho (A-2)
1012	//[3] background (cm-1)
1013
1014
1015	// calculates scale * f^2/Vol where f=Vol3delrho*((sin(qr)-qrcos(qr))/qr^3
1016	// and is rescaled to give [=] cm^-1
1017
1018	Variable bes,f,vol,f2
1019	//
1020	//handle q==0 separately
1021	If(x==0)
1022	f = 4/3piradius^3delrhodelrhoscale1e8 + bkg
1023	return(f)
1024	Endif
1025
1026	bes = 3(sin(xradius)-xradiuscos(x*radius))/x^3/radius^3
1027	vol = 4pi/3radius^3
1028	f = volbesdelrho // [=] A
1029	// normalize to single particle volume, convert to 1/cm
1030	f2 = f * f / vol * 1.0e8 // [=] 1/cm
1031
1032	return (scale*f2+bkg) // Scale, then add in the background
1033
1034	End
1035
1036	// better default function
1037	Function U_Power_Law_Model(A,m,bgd,x) : FitFunc
1038	Variable A, m,bgd,x
1039	// Input (fitting) variables are:
1040	//[0] Coefficient
1041	//[1] (-) Power
1042	//[2] incoherent background
1043
1044	// local variables
1045	Variable inten, qval
1046	// x is the q-value for the calculation
1047	qval = x
1048	// do the calculation and return the function value
1049
1050	inten = A*qval^-m + bgd
1051	Return (inten)
1052	End
1053
1054
1055	// mimics LoadBT5File
1056	// creates two other waves to identify the set and the counting time for that set
1057	//
1058	Function LoadSimulatedAngleRange(set,type)
1059	Variable set
1060	String type
1061
1062	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1063
1064	String s=num2str(set)
1065	NVAR angLow = $("root:Packages:NIST:USANS:Globals:U_Sim:gAngLow"+s)
1066	NVAR angHigh = $("root:Packages:NIST:USANS:Globals:U_Sim:gAngHigh"+s)
1067	NVAR numPts = $("root:Packages:NIST:USANS:Globals:U_Sim:gNumPts"+s)
1068	NVAR ctTime = $("root:Packages:NIST:USANS:Globals:U_Sim:gCtTime"+s)
1069	NVAR incr = $("root:Packages:NIST:USANS:Globals:U_Sim:gIncr"+s)
1070
1071	Variable ii, err=0
1072
1073	// generate q-points based on angular range from panel
1074
1075	Make/O/D/N=(numPts) $(USANSFolder+":"+type+":Angle")
1076	Make/O/D/N=(numPts) $(USANSFolder+":"+type+":DetCts")
1077	Make/O/D/N=(numPts) $(USANSFolder+":"+type+":ErrDetCts")
1078	Make/O/D/N=(numPts) $(USANSFolder+":"+type+":MonCts")
1079	Make/O/D/N=(numPts) $(USANSFolder+":"+type+":TransCts")
1080	//
1081	Make/O/D/N=(numPts) $(USANSFolder+":"+type+":CountingTime")
1082	Make/O/D/N=(numPts) $(USANSFolder+":"+type+":SetNumber")
1083
1084	Wave Angle = $(USANSFolder+":"+type+":Angle")
1085	Wave DetCts = $(USANSFolder+":"+type+":DetCts")
1086	Wave ErrDetCts = $(USANSFolder+":"+type+":ErrDetCts")
1087	Wave MonCts = $(USANSFolder+":"+type+":MonCts")
1088	Wave TransCts = $(USANSFolder+":"+type+":TransCts")
1089	Wave countingTime = $(USANSFolder+":"+type+":countingTime")
1090	Wave SetNumber = $(USANSFolder+":"+type+":SetNumber")
1091
1092	countingTime = ctTime
1093	SetNumber = set
1094
1095	for(ii=0;ii<numPts;ii+=1)
1096	Angle[ii] = angLow + ii*incr
1097	DetCts[ii] = set
1098	endfor
1099
1100	//set the wave note for the DetCts
1101	String str=""
1102	str = "FILE:Sim "+s+";"
1103	str += "TIMEPT:"+num2str(ctTime)+";"
1104	str += "PEAKANG:0;" //no value yet
1105	str += "STATUS:;" //no value yet
1106	str += "LABEL:SimSet "+s+";"
1107	str += "PEAKVAL:;" //no value yet
1108	str += "TWIDE:0;" //no value yet
1109	Note DetCts,str
1110
1111	return err // Zero signifies no error.
1112	End
1113
1114	// add SIM data to the graph if it exists and is not already on the graph
1115	//
1116	//
1117	// ** currently, I have changed the graph to not display the angle at the top. I have not yet
1118	// found a way to keep the scaling of the two axes in-sync when the empty cell data is added to the
1119	// graph (versus Q)
1120	//
1121	Function GraphSIM()
1122
1123	// SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1124	// SetDataFolder $(USANSFolder+":SIM")
1125
1126	SetDataFolder root:Sim_USANS
1127	//is it already on the graph?
1128	String list=""
1129	list = Wavelist("Sim_USANS*",";","WIN:UCALC")
1130
1131	if(strlen(list)!=0)
1132	// Print "SIM already on graph"
1133	SetDataFolder root:
1134	return(0)
1135	Endif
1136
1137	//append the data if it exists
1138	If(waveExists($"Sim_USANS_i")==1)
1139	DoWindow/F UCALC
1140
1141	// lines for the no-noise result vs angle
1142	// AppendToGraph/T Smeared_inten vs angle
1143	AppendToGraph Smeared_inten vs Sim_USANS_q
1144	ModifyGraph rgb(Smeared_inten)=(1,12815,52428)
1145	ModifyGraph mode(Smeared_inten)=4,marker(Smeared_inten)=19,msize(Smeared_inten)=2
1146	ModifyGraph tickUnit=1
1147
1148	// colored points for the simulation with noise on top
1149	AppendToGraph Sim_USANS_i vs Sim_USANS_q
1150	ModifyGraph mode(Sim_USANS_i)=3,marker(Sim_USANS_i)=19,msize(Sim_USANS_i)=4
1151	//don't reverse the dbZ21 or the highest angle will have "invisible" light blue error bars
1152	ModifyGraph zColor(Sim_USANS_i)={setNumber,1,7,dBZ21,0} //better for 7 colors
1153	// ModifyGraph zColor(Sim_USANS_i)={setNumber,1,7,Rainbow,0} //force the colors from 1->7
1154	ModifyGraph useMrkStrokeRGB(Sim_USANS_i)=1
1155	ErrorBars/T=0 Sim_USANS_i Y,wave=(Sim_USANS_s,Sim_USANS_s)
1156
1157	ModifyGraph log=1
1158	ModifyGraph mirror(left)=1
1159	// if no top axis, then
1160	// ModifyGraph mirror(bottom)=1
1161	ModifyGraph grid=2
1162	ModifyGraph standoff=0
1163
1164	// to make sure that the scales are the same (but fails on zoom)
1165	NewFreeAxis/O/T top_angle
1166	ModifyFreeAxis/Z top_angle,master= bottom,hook= TransformAngleAxisHook
1167	ModifyGraph lblPos(top_angle)=50,freePos(top_angle)=0
1168	ModifyGraph log(top_angle)=1
1169
1170	SetDrawEnv linefgc= (39321,1,1),dash= 3,linethick= 3
1171	SetDrawEnv xcoord= bottom
1172	DrawLine 3e-05,0.01,3e-05,0.99
1173
1174	Label top_angle "Angle"
1175	Label bottom "Q (1/A)"
1176	Label left "Counts or Count Rate"
1177
1178	Legend
1179	endif
1180
1181	SetDataFolder root:
1182	End
1183
1184	Function TransformAngleAxisHook(s)
1185	STRUCT WMAxisHookStruct &s
1186
1187	s.max= s.max/5.55e-5
1188	s.min= s.min/5.55e-5
1189
1190	return 0
1191	End
1192
1193	//fakes a folder with loaded 1-d usans data, no calculation of the matrix
1194	Function FakeUSANSDataFolder(qval,aveint,dqv,dataFolder)
1195	WAVE qval,aveint
1196	Variable dqv
1197	String dataFolder
1198
1199
1200	String baseStr=dataFolder
1201	if(DataFolderExists("root:"+baseStr))
1202	SetDataFolder $("root:"+baseStr)
1203	else
1204	NewDataFolder/S $("root:"+baseStr)
1205	endif
1206
1207	////overwrite the existing data, if it exists
1208	Duplicate/O qval, $(baseStr+"_q")
1209	Duplicate/O aveint, $(baseStr+"_i")
1210
1211	Duplicate/O qval, $(baseStr+"_s")
1212	Wave sigave = $(baseStr+"_s")
1213
1214	sigave = sqrt(aveint)
1215
1216	// make a resolution matrix for SANS data
1217	Variable np=numpnts(qval)
1218	Make/D/O/N=(np,4) $(baseStr+"_res")
1219	Wave res=$(baseStr+"_res")
1220
1221	res[][0] = -dQv //sigQ
1222	res[][1] = -dQv //qBar
1223	res[][2] = -dQv //fShad
1224	res[][3] = qval[p] //Qvalues
1225
1226	// extra waves of set number and counting time for the simulation
1227	WAVE ctW = root:Packages:NIST:USANS:SIM:countingTime
1228	WAVE setW = root:Packages:NIST:USANS:SIM:setNumber
1229	WAVE ang = root:Packages:NIST:USANS:SIM:Angle
1230	Duplicate/O ctW countingTime
1231	Duplicate/O setW setNumber
1232	Duplicate/O ang Angle
1233
1234
1235	//clean up
1236	SetDataFolder root:
1237
1238	End
1239
1240
1241	Function/S U_FunctionPopupList()
1242	String list,tmp
1243	list = User_FunctionPopupList()
1244
1245	//simplify the display, forcing smeared calculations behind the scenes
1246	tmp = FunctionList("Smear*",";","NPARAMS:1") //smeared dependency calculations
1247	list = RemoveFromList(tmp, list ,";")
1248
1249
1250	if(strlen(list)==0)
1251	list = "No functions plotted"
1252	endif
1253
1254	list = SortList(list)
1255
1256	return(list)
1257	End
1258
1259	Function Sim_USANS_ModelPopMenuProc(pa) : PopupMenuControl
1260	STRUCT WMPopupAction &pa
1261
1262	switch( pa.eventCode )
1263	case 2: // mouse up
1264	Variable popNum = pa.popNum
1265	String popStr = pa.popStr
1266	SVAR gStr = root:Packages:NIST:USANS:Globals:U_Sim:gFuncStr
1267	gStr = popStr
1268
1269	break
1270	endswitch
1271
1272	return 0
1273	End
1274
1275	// if the sample aperture is changed, AND the empty data is displayed, change to the proper data
1276	Function UCALC_SampleAperturePopup(pa) : PopupMenuControl
1277	STRUCT WMPopupAction &pa
1278
1279	switch( pa.eventCode )
1280	case 2: // mouse up
1281	Variable popNum = pa.popNum
1282	String popStr = pa.popStr
1283
1284	ControlInfo/W=UCALC check0_4
1285	if(V_Value==1) //currently checked, need to update
1286
1287	STRUCT WMCheckboxAction cba
1288	cba.checked=0 //"un-check"
1289	cba.eventCode=2
1290	ShowEMPCheckProc(cba)
1291
1292	cba.checked=1 //"re-check"
1293	ShowEMPCheckProc(cba)
1294
1295	endif
1296	break
1297
1298	endswitch
1299
1300	return 0
1301	End
1302
1303
1304	Function UCALC_PresetPopup(pa) : PopupMenuControl
1305	STRUCT WMPopupAction &pa
1306
1307	switch( pa.eventCode )
1308	case 2: // mouse up
1309	Variable popNum = pa.popNum
1310	String popStr = pa.popStr
1311
1312	SetDataFolder root:Packages:NIST:USANS:Globals:U_Sim
1313
1314	NVAR gAngLow1 = gAngLow1
1315	NVAR gAngHigh1 = gAngHigh1
1316	NVAR gNumPts1 = gNumPts1
1317	NVAR gCtTime1 = gCtTime1
1318	NVAR gIncr1 = gIncr1
1319
1320	NVAR gAngLow2 = gAngLow2
1321	NVAR gAngHigh2 = gAngHigh2
1322	NVAR gNumPts2 = gNumPts2
1323	NVAR gCtTime2 = gCtTime2
1324	NVAR gIncr2 = gIncr2
1325
1326	NVAR gAngLow3 = gAngLow3
1327	NVAR gAngHigh3 = gAngHigh3
1328	NVAR gNumPts3 = gNumPts3
1329	NVAR gCtTime3 = gCtTime3
1330	NVAR gIncr3 = gIncr3
1331
1332	NVAR gAngLow4 = gAngLow4
1333	NVAR gAngHigh4 = gAngHigh4
1334	NVAR gNumPts4 = gNumPts4
1335	NVAR gCtTime4 = gCtTime4
1336	NVAR gIncr4 = gIncr4
1337
1338	NVAR gAngLow5 = gAngLow5
1339	NVAR gAngHigh5 = gAngHigh5
1340	NVAR gNumPts5 = gNumPts5
1341	NVAR gCtTime5 = gCtTime5
1342	NVAR gIncr5 = gIncr5
1343
1344	NVAR gAngLow6 = gAngLow6
1345	NVAR gAngHigh6 = gAngHigh6
1346	NVAR gNumPts6 = gNumPts6
1347	NVAR gCtTime6 = gCtTime6
1348	NVAR gIncr6 = gIncr6
1349
1350	NVAR gAngLow7 = gAngLow7
1351	NVAR gAngHigh7 = gAngHigh7
1352	NVAR gNumPts7 = gNumPts7
1353	NVAR gCtTime7 = gCtTime7
1354	NVAR gIncr7 = gIncr7
1355
1356	strswitch(popStr) // string switch
1357	case "Short Count": // execute if case matches expression
1358
1359	gAngLow1 = -1
1360	gAngHigh1 = 0.6
1361	gNumPts1 = 33
1362	gCtTime1 = 10
1363	gIncr1 = 0.05
1364
1365	gAngLow2 = 0.7
1366	gAngHigh2 = 1.9
1367	gNumPts2 = 6
1368	gCtTime2 = 60
1369	gIncr2 = 0.2
1370
1371	gAngLow3 = 2
1372	gAngHigh3 = 4.8
1373	gNumPts3 = 7
1374	gCtTime3 = 120
1375	gIncr3 = 0.4
1376
1377	gAngLow4 = 5
1378	gAngHigh4 = 9.5
1379	gNumPts4 = 5
1380	gCtTime4 = 180
1381	gIncr4 = 0.9
1382
1383	gAngLow5 = 10
1384	gAngHigh5 = 19
1385	gNumPts5 = 5
1386	gCtTime5 = 240
1387	gIncr5 = 1.8
1388
1389	gAngLow6 = 20
1390	gAngHigh6 = 48
1391	gNumPts6 = 15
1392	gCtTime6 = 0
1393	gIncr6 = 2
1394
1395	gAngLow7 = 50
1396	gAngHigh7 = 95
1397	gNumPts7 = 10
1398	gCtTime7 = 0
1399	gIncr7 = 5
1400
1401	break
1402	case "Medium Count":
1403
1404	gAngLow1 = -1
1405	gAngHigh1 = 0.6
1406	gNumPts1 = 33
1407	gCtTime1 = 10
1408	gIncr1 = 0.05
1409
1410	gAngLow2 = 0.7
1411	gAngHigh2 = 1.9
1412	gNumPts2 = 13
1413	gCtTime2 = 60
1414	gIncr2 = 0.1
1415
1416	gAngLow3 = 2
1417	gAngHigh3 = 4.8
1418	gNumPts3 = 15
1419	gCtTime3 = 120
1420	gIncr3 = 0.2
1421
1422	gAngLow4 = 5
1423	gAngHigh4 = 9.5
1424	gNumPts4 = 10
1425	gCtTime4 = 300
1426	gIncr4 = 0.5
1427
1428	gAngLow5 = 10
1429	gAngHigh5 = 19
1430	gNumPts5 = 10
1431	gCtTime5 = 600
1432	gIncr5 = 1
1433
1434	gAngLow6 = 20
1435	gAngHigh6 = 48
1436	gNumPts6 = 15
1437	gCtTime6 = 1200
1438	gIncr6 = 2
1439
1440	gAngLow7 = 50
1441	gAngHigh7 = 95
1442	gNumPts7 = 10
1443	gCtTime7 = 0
1444	gIncr7 = 5
1445
1446	break
1447	case "Long Count":
1448
1449	gAngLow1 = -1
1450	gAngHigh1 = 0.6
1451	gNumPts1 = 33
1452	gCtTime1 = 25
1453	gIncr1 = 0.05
1454
1455	gAngLow2 = 0.7
1456	gAngHigh2 = 1.9
1457	gNumPts2 = 13
1458	gCtTime2 = 100
1459	gIncr2 = 0.1
1460
1461	gAngLow3 = 2
1462	gAngHigh3 = 4.8
1463	gNumPts3 = 15
1464	gCtTime3 = 300
1465	gIncr3 = 0.2
1466
1467	gAngLow4 = 5
1468	gAngHigh4 = 9.5
1469	gNumPts4 = 10
1470	gCtTime4 = 600
1471	gIncr4 = 0.5
1472
1473	gAngLow5 = 10
1474	gAngHigh5 = 19
1475	gNumPts5 = 10
1476	gCtTime5 = 1200
1477	gIncr5 = 1
1478
1479	gAngLow6 = 20
1480	gAngHigh6 = 48
1481	gNumPts6 = 15
1482	gCtTime6 = 2000
1483	gIncr6 = 2
1484
1485	gAngLow7 = 50
1486	gAngHigh7 = 95
1487	gNumPts7 = 10
1488	gCtTime7 = 3000
1489	gIncr7 = 5
1490
1491	break
1492	default:
1493	endswitch
1494
1495	break
1496	endswitch
1497
1498	//update the count time
1499	CalcTotalCountTime()
1500
1501	SetDataFolder root:
1502
1503	return 0
1504	End
1505
1506
1507	// return the beam intensity based on the sample aperture diameter
1508	//
1509	// based on the equation in John's instrument paper
1510	//
1511	Function GetUSANSBeamIntensity()
1512
1513	String popStr
1514	Variable flux,diam,rad
1515
1516	ControlInfo/W=UCALC popup0
1517	popStr = S_Value
1518	diam=str2num(popStr) //in inches
1519	rad = diam/2*25.4 //radius in mm
1520
1521	flux = 662radrad-39.9radradrad+0.70radradrad*rad
1522
1523	return(flux)
1524	End
1525
1526	// based on the angle ranges above (with non-zero count times)
1527	// save fake data points into a fake BT5 data file
1528	//
1529	Function U_SaveButtonProc(ba) : ButtonControl
1530	STRUCT WMButtonAction &ba
1531
1532
1533	switch( ba.eventCode )
1534	case 2: // mouse up
1535	// click code here
1536
1537	Variable num,ii,baseNumber=100,firstSet=0
1538	String pathStr="root:Packages:NIST:USANS:Globals:U_Sim:gCtTime"
1539	SVAR USANSFolder = root:Packages:NIST:USANS:Globals:gUSANSFolder
1540	String baseName="SIMUL"
1541
1542	num = 7 //# of angular ranges
1543
1544	// only try to plot ranges with non-zero count times
1545	firstSet=0
1546	for(ii=1;ii<=num;ii+=1)
1547	NVAR dum = $(pathStr+num2str(ii))
1548	if(dum>0)
1549	firstSet += 1
1550	if(firstSet==1) //first time, ask for base name
1551	Prompt baseName,"Enter a base name for the files"
1552	Prompt baseNumber,"Enter the starting index"
1553	DoPrompt "Enter Information for Save",baseName,baseNumber
1554
1555	if(V_Flag==1) //user canceled
1556	return(1)
1557	endif
1558	endif
1559
1560	SaveFakeUSANS(baseName,baseNumber-1,ii)
1561
1562	endif
1563	endfor
1564
1565	break
1566	endswitch
1567
1568	return 0
1569	End
1570
1571
1572	//
1573	// duplicates the BT5 file format
1574	//
1575	Function SaveFakeUSANS(nameStr,num,set)
1576	String nameStr
1577	Variable num,set
1578
1579	String folder = "root:Packages:NIST:USANS:SIM:"
1580	String pathStr="root:Packages:NIST:USANS:Globals:U_Sim:"
1581	String termStr="\r\n" //VAX uses only <CR> as terminator, but only CRLF seems to FTP correctly to VAX
1582
1583	// WAVE DetCts = $(folder+"DetCts") //these are only dummy values
1584	WAVE DetCts = root:Sim_USANS:Sim_USANS_i
1585	WAVE angle = $(folder+"Angle")
1586	WAVE SetNumber = $(folder+"SetNumber")
1587	WAVE countingTime = $(folder+"countingTime")
1588
1589	NVAR ang1 = $(pathStr+"gAngLow"+num2str(set))
1590	NVAR ang2 = $(pathStr+"gAngHigh"+num2str(set))
1591	NVAR incr = $(pathStr+"gIncr"+num2str(set))
1592
1593	Variable refNum,ii,wavePts,numPts,first,last,ctTime,monCt,transDet
1594	String str,fileStr,dateStr
1595
1596	wavePts = numpnts(angle)
1597	for(ii=0;ii<wavePts;ii+=1)
1598	if(setNumber[ii] == set)
1599	first = ii
1600	break
1601	endif
1602	endfor
1603
1604	for(ii=wavePts-1;ii>=0;ii-=1)
1605	if(setNumber[ii] == set)
1606	last = ii
1607	break
1608	endif
1609	endfor
1610
1611	// Print "First, last = ",first,last
1612
1613	// set up some of the strings needed
1614	fileStr=nameStr+num2str(num+set)+".bt5"
1615	dateStr = date()// + " "+Secs2Time(DateTime,2)
1616	ctTime = countingTime[first]
1617	numPts = last-first+1
1618
1619	MonCt = ctTime*GetUSANSBeamIntensity()
1620
1621	transDet = 999 //bogus and constant value
1622
1623	//actually open the file
1624	Open refNum as fileStr
1625
1626	sprintf str,"'%s' '%s' 'I' %d 1 'TIME' %d 'RAW'",fileStr,dateStr,ctTime,numpts
1627	fprintf refnum,"%s"+termStr,str
1628	fprintf refnum,"%s"+termStr," Filename Date Scan Mon Prf Base #pts Type"
1629
1630	fprintf refnum,"Simulated USANS data"+termStr
1631
1632	fprintf refnum,"%s"+termStr," 0 0 0 0 0 0 0 0.0000 0.00000 0.00000 0.00000 2"
1633	fprintf refnum,"%s"+termStr," Collimation Mosaic Wavelength T-Start Incr. H-field #Det "
1634	fprintf refnum,"%s"+termStr," 1 0.00000 0.00000 0.00000"
1635	fprintf refnum," 2 %9.5f %9.5f %9.5f"+termStr,ang1,incr,ang2
1636	fprintf refnum,"%s"+termStr," 3 10.00000 0.00000 10.00000"
1637	fprintf refnum,"%s"+termStr," 4 10.00000 0.00000 10.00000"
1638	fprintf refnum,"%s"+termStr," 5 0.00000 0.00000 0.00000"
1639	fprintf refnum,"%s"+termStr," 6 0.00000 0.00000 0.00000"
1640	fprintf refnum,"%s"+termStr," Mot: Start Step End"
1641	fprintf refnum,"%s"+termStr," A2 MIN COUNTS "
1642
1643	//loop over the waves, picking out the desired set
1644	//write 2 lines each time
1645	for(ii=first;ii<=last;ii+=1)
1646	sprintf str," %6.3f %6.2f %d",angle[ii],ctTime/60,DetCts[ii]
1647	fprintf refnum,"%s"+termStr,str
1648
1649	sprintf str,"%d,%d,0,%d,0,0,0,0",MonCt,DetCts[ii],transDet
1650	fprintf refnum,"%s"+termStr,str
1651	endfor
1652
1653	Close refnum
1654
1655	return(0)
1656	end

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

Download in other formats:

Original Format