Context Navigation

SASCALC.ipf @ 564

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

A variety of changes...

fixed model strings for Sum_Model
/Z flag for RTHelp, Transmission help, FIT help, MRED help
added help and done buttons to SASCALC
added comments in MonteCarlo?
restored the 2-line buttons in USANS panels, now supported in Igor 6.1 (COR_Graph and Main_USANS)
added simple text notebook of UCALC configuration, and a save of PNG of the whole panel
added help and done buttons to UCALC

added links in SANS_Reduction_Help to the simluation help file (not added to SVN yet)

File size: 65.0 KB

Line
1	#pragma rtGlobals=1 // Use modern global access method.
2	#pragma version=1.0
3	#pragma IgorVersion = 6.0
4	#pragma ModuleName=SASCALC
5
6	// SASCALC.ipf
7	//
8	// 04 OCT 2006 SRK
9	// 30 OCT 2006 SRK - corrected beamDiameter size in attenuator calculation (Bh vs Bm)
10	// 11 DEC 2006 SRK - added 2.5 cm A1 option for NG3@7guide config
11	// 09 MAR 2007 SRK - now appends text of each frozen configuration for printing
12	// - colorized frozen traces so that they aren't all red (unfrozen is black)
13	// 19 MAR 2007 SRK - corrections added for projected BS diameter at anode plane
14	// 11 APR 2007 SRK - default aperture offset of 5 cm added to match VAX implementation
15	// nn AUG 2007 SRK - added defulat sample aperture size
16	// added option of lenses, approximated beamDiam=sourceDiam, BSdiam=1"
17	// Lens flux, trans is not corrected for lens/prism transmission
18	// Lenses can still be inserted in incorrect cases, and are not automatically taken out
19	// 27 JAN 2009 SRK - Changed behavior of Lens checkbox. Now, it SETS parameters as needed for proper
20	// configuration. 17.2 can be typed in for lens/prism on NG3. Invalid conditions
21	// will automatically uncheck the box
22	//
23	// calculate what q-values you get based on the instruments settings
24	// that you would typically input to SASCALC
25	// calculation is for (80A radius spheres) * (beam stop shadow factor)
26	// or a Debye function (Rg=50A) * (beam stop shadow factor)
27	// - NOT true intensity, not counts, just a display
28	//
29	// To Do:
30	//
31	// Optional:
32	// - freeze configurations with a user defined tag
33	// - different model functions (+ change simple parameters)
34	// - resolution smeared models
35	// - "simulation" of data and error bars given a model and a total number of detector counts
36	// - streamline code (globals needed in panel vs. wave needed for calculation)
37	//
38	// - there is a lot of re-calculation of things (a consequence of the fake-OO) that could be streamlined
39	//
40	// Done:
41	// - include resolution effects (includes BS effect, smeared model)
42	// - (data = 1) then multiply by a typical form factor
43	// - masked two pixels around the edge, as default
44	// - conversion from # guides to SSD from sascalc
45	// - show multiple configurations at once
46	// - interactive graphics on panel
47	// - full capabilities of SASCALC
48	// - correct beamstop diameter
49	// - get the sample/huber position properly incorporated (ssd and sdd)
50	// - get SDD value to update when switching NG7->NG3 and illegal value
51	// - disallow 6 guides at NG3 (post a warning)
52	//
53	//
54
55	Proc SASCALC()
56	DoWindow/F SASCALC
57	if(V_flag==0)
58	S_initialize_space()
59	initNG3() //start life as NG3
60	Sascalc_Panel()
61	ReCalculateInten(1) //will use defaults
62	Endif
63
64	// now a checkbox as needed
65	// DoWindow/F MC_SASCALC
66	// if(V_flag==0)
67	// MC_SASCALC()
68	// endif
69	End
70
71	Proc S_initialize_space()
72	NewDataFolder/O root:Packages
73	NewDataFolder/O root:Packages:NIST
74	NewDataFolder/O root:Packages:NIST:SAS
75
76	Make/O/D/N=23 root:Packages:NIST:SAS:integersRead
77	Make/O/D/N=52 root:Packages:NIST:SAS:realsRead
78	Make/O/T/N=11 root:Packages:NIST:SAS:textRead
79	// data
80	Make/O/D/N=(128,128) root:Packages:NIST:SAS:data,root:Packages:NIST:SAS:linear_data
81	Make/O/D/N=2 root:Packages:NIST:SAS:aveint,root:Packages:NIST:SAS:qval,root:Packages:NIST:SAS:sigave
82	root:Packages:NIST:SAS:data = 1
83	root:Packages:NIST:SAS:linear_data = 1
84	// fill w/default values
85	S_fillDefaultHeader(root:Packages:NIST:SAS:integersRead,root:Packages:NIST:SAS:realsRead,root:Packages:NIST:SAS:textRead)
86
87	// other variables
88	// -(hard coded right now - look for NVAR declarations)
89	Variable/G root:Packages:NIST:SAS:gBinWidth=1
90	Variable/G root:Packages:NIST:SAS:gisLogScale=0
91	String/G root:Packages:NIST:SAS:FileList = "SASCALC"
92
93	// for the panel
94	Variable/G root:Packages:NIST:SAS:gInst=3 //or 7 for NG7
95	Variable/G root:Packages:NIST:SAS:gNg=0
96	Variable/G root:Packages:NIST:SAS:gTable=2 //2=chamber, 1=table
97	Variable/G root:Packages:NIST:SAS:gDetDist=1000 //sample chamber to detector in cm
98	Variable/G root:Packages:NIST:SAS:gSSD=1632 //!!SSD in cm fo 0 guides (derived from Ng)
99	Variable/G root:Packages:NIST:SAS:gOffset=0
100	Variable/G root:Packages:NIST:SAS:gSamAp=1.27 //samAp diameter in cm
101	Variable/G root:Packages:NIST:SAS:gLambda=6
102	Variable/G root:Packages:NIST:SAS:gDeltaLambda=0.125 //default value
103	String/G root:Packages:NIST:SAS:gSourceApString = "1.43 cm;2.54 cm;3.81 cm;"
104	String/G root:Packages:NIST:SAS:gDeltaLambdaStr = "0.109;0.125;0.236;" //ng3 defaults
105	String/G root:Packages:NIST:SAS:gApPopStr = "1/16\";1/8\";3/16\";1/4\";5/16\";3/8\";7/16\";1/2\";9/16\";5/8\";11/16\";3/4\";other;"
106	Variable/G root:Packages:NIST:SAS:gSamApOther = 10 //non-standard aperture diameter, in mm
107	Variable/G root:Packages:NIST:SAS:gUsingLenses = 0 //0=no lenses, 1=lenses(or prisms)
108	Variable/G root:Packages:NIST:SAS:gModelOffsetFactor = 1
109
110	// for the MC simulation
111	Variable/G root:Packages:NIST:SAS:doSimulation =0 // == 1 if 1D simulated data, 0 if other from the checkbox
112	Variable/G root:Packages:NIST:SAS:gRanDateTime=datetime
113	Variable/G root:Packages:NIST:SAS:gImon = 10000
114	Variable/G root:Packages:NIST:SAS:gThick = 0.1
115	Variable/G root:Packages:NIST:SAS:gSig_incoh = 0.1
116	String/G root:Packages:NIST:SAS:gFuncStr = ""
117	Variable/G root:Packages:NIST:SAS:gR2 = 2.54/2
118	Variable/G root:Packages:NIST:SAS:gSamTrans=0.8 //for 1D, default value
119	Variable/G root:Packages:NIST:SAS:gCntTime = 300
120	Variable/G root:Packages:NIST:SAS:gDoMonteCarlo = 0
121	Variable/G root:Packages:NIST:SAS:gUse_MC_XOP = 1 //set to zero to use Igor code
122	Variable/G root:Packages:NIST:SAS:gBeamStopIn = 1 //set to zero for beamstop out (transmission)
123	Variable/G root:Packages:NIST:SAS:gRawCounts = 0
124	Variable/G root:Packages:NIST:SAS:gSaveIndex = 100
125	String/G root:Packages:NIST:SAS:gSavePrefix = "SIMUL"
126	Make/O/D/N=10 root:Packages:NIST:SAS:results = 0
127	Make/O/T/N=10 root:Packages:NIST:SAS:results_desc = {"total X-section (1/cm)","SAS X-section (1/cm)","number that scatter","number that reach detector","avg # times scattered","fraction single coherent","fraction multiple coherent","fraction multiple scattered","fraction transmitted","detector counts w/beamstop"}
128
129	Variable/G root:Packages:NIST:SAS:g_1DTotCts = 0 //summed counts (simulated)
130	Variable/G root:Packages:NIST:SAS:g_1DEstDetCR = 0 // estimated detector count rate
131	Variable/G root:Packages:NIST:SAS:g_1DFracScatt = 0 // fraction of beam captured on detector
132	Variable/G root:Packages:NIST:SAS:g_1DEstTrans = 0 // estimated transmission of sample
133	Variable/G root:Packages:NIST:SAS:g_1D_DoABS = 1
134	Variable/G root:Packages:NIST:SAS:g_1D_AddNoise = 1
135
136	//tick labels for SDD slider
137	//userTicks={tvWave,tlblWave }
138	Make/O/D/N=5 root:Packages:NIST:SAS:tickSDDNG3,root:Packages:NIST:SAS:tickSDDNG7
139	Make/O/T/N=5 root:Packages:NIST:SAS:lblSDDNG3,root:Packages:NIST:SAS:lblSDDNG7
140	root:Packages:NIST:SAS:tickSDDNG3 = {133,400,700,1000,1317}
141	root:Packages:NIST:SAS:lblSDDNG3 = {"133","400","700","1000","1317"}
142	root:Packages:NIST:SAS:tickSDDNG7 = {100,450,800,1150,1530}
143	root:Packages:NIST:SAS:lblSDDNG7 = {"100","450","800","1150","1530"}
144
145	//for the fake dependency
146	Variable/G root:Packages:NIST:SAS:gTouched=0
147	Variable/G root:Packages:NIST:SAS:gCalculate=0
148	//for plotting
149	Variable/G root:Packages:NIST:SAS:gFreezeCount=1 //start the count at 1 to keep Jeff happy
150	Variable/G root:Packages:NIST:SAS:gDoTraceOffset=0 // (1==Yes, offset 2^n), 0==turn off the offset
151
152	End
153
154	Function initNG3()
155
156	SetDataFolder root:Packages:NIST:SAS
157
158	Variable/G instrument = 3
159	Variable/G s12 = 54.8
160	Variable/G d_det = 0.5
161	Variable/G a_pixel = 0.5
162	Variable/G del_r = 0.5
163	Variable/G det_width = 64.0
164	Variable/G phi_0 = 2.95e13
165	Variable/G lambda_t = 5.50
166	Variable/G l2r_lower = 132.3
167	Variable/G l2r_upper = 1317
168	Variable/G lambda_lower = 2.5
169	Variable/G lambda_upper = 20.0
170	Variable/G d_upper = 25.0
171	Variable/G bs_factor = 1.05
172	Variable/G t1 = 0.63
173	Variable/G t2 = 1.0
174	Variable/G t3 = 0.75
175	Variable/G l_gap = 100.0
176	Variable/G guide_width = 6.0
177	Variable/G idmax = 100.0
178	Variable/G b = 0.023
179	Variable/G c = 0.023
180
181	//fwhm values (new variables) (+3, 0, -3, calibrated 2009)
182	Variable/G fwhm_narrow = 0.109
183	Variable/G fwhm_mid = 0.125
184	Variable/G fwhm_wide = 0.236
185
186	//source apertures (cm)
187	Variable/G a1_0_0 = 1.43
188	Variable/G a1_0_1 = 2.54
189	Variable/G a1_0_2 = 3.81
190	Variable/G a1_7_0 = 2.5 // after the polarizer
191	Variable/G a1_7_1 = 5.0
192	Variable/G a1_def = 5.00
193
194	//default configuration values
195	// ng = 0
196	// a1 = 3.81
197	// pos_table = 2
198	// l2r = 1310.0
199	// a2 = 1.27
200	// det_off = 0.0
201	// lambda = 6.0
202	// lambda_width = 0.15
203	Variable/G l2diff = 5.0
204	//
205	SetDataFolder root:
206	end
207
208	Function initNG7()
209
210	SetDataFolder root:Packages:NIST:SAS
211
212	Variable/G instrument = 7
213	Variable/G s12 = 54.8
214	Variable/G d_det = 0.5
215	Variable/G a_pixel = 0.5
216	Variable/G del_r = 0.5
217	Variable/G det_width = 64.0
218	Variable/G phi_0 = 2.3e13
219	Variable/G lambda_t = 5.50
220	Variable/G l2r_lower = 100
221	Variable/G l2r_upper = 1531
222	Variable/G lambda_lower = 4.0
223	Variable/G lambda_upper = 20.0
224	Variable/G d_upper = 25.0
225	Variable/G bs_factor = 1.05
226	Variable/G t1 = 0.63
227	Variable/G t2 = 0.7
228	Variable/G t3 = 0.75
229	Variable/G l_gap = 188.0
230	Variable/G guide_width = 5.0
231	Variable/G idmax = 100.0
232	Variable/G b = 0.028
233	Variable/G c = 0.028
234
235	//fwhm values (new variables)
236	Variable/G fwhm_narrow = 0.09
237	Variable/G fwhm_mid = 0.11
238	Variable/G fwhm_wide = 0.22
239
240	//source apertures (cm)
241	Variable/G a1_0_0 = 1.43
242	Variable/G a1_0_1 = 2.22
243	Variable/G a1_0_2 = 3.81
244	Variable/G a1_7_0 = 5.0 // don't apply to NG7
245	Variable/G a1_7_1 = 5.0
246	Variable/G a1_def = 5.00
247
248	//default configuration values
249	// ng = 0
250	// a1 = 2.22
251	// pos_table = 2
252	// l2r = 1530.0
253	// a2 = 1.27
254	// det_off = 0.0
255	// lambda = 6.0
256	// lambda_width = 0.11
257	Variable/G l2diff = 5.0
258	//
259	SetDataFolder root:
260	end
261
262	Function S_fillDefaultHeader(iW,rW,tW)
263	Wave iW,rW
264	Wave/T tW
265
266	// text wave
267	// don't need anything
268
269	// integer wave
270	// don't need anything
271
272	// real wave
273	rw = 0
274	rw[16] = 64 // beamcenter X (pixels)
275	rw[17] = 64 // beamcenter Y
276
277	rw[10] = 5.08 //detector resolution (5mm) and calibration constants (linearity)
278	rw[11] = 10000
279	rw[12] = 0
280	rw[13] = 5.08
281	rw[14] = 10000
282	rw[15] = 0
283
284	rw[20] = 65 // det size in cm
285	rw[18] = 6 // SDD in meters (=L2)
286
287	rw[26] = 6 //lambda in Angstroms
288	rw[4] = 1 //transmission
289
290	rw[21] = 76.2 //BS diameter in mm
291	rw[23] = 50 //A1 diameter in mm
292	rw[24] = 12.7 //A2 diameter in mm
293	rw[25] = 7.02 //L1 distance in meters (derived from number of guides)
294	rw[27] = 0.11 //DL/L wavelength spread
295
296	return(0)
297	End
298
299
300	Window SASCALC_Panel()
301
302	PauseUpdate; Silent 1 // building window...
303
304	// if I make the graph a subwindow in a panel, it breaks the "append 1d" from the wrapper
305	// NewPanel/W=(5,44,463,570)/K=1 as "SASCALC"
306	// DoWindow/C SASCALC
307	// ModifyPanel cbRGB=(49151,53155,65535)
308	//
309	//
310	// String fldrSav0= GetDataFolder(1)
311	// SetDataFolder root:Packages:NIST:SAS:
312	//
313	// Display/HOST=#/W=(5,200,463,570) aveint vs qval
314	// ModifyGraph mode=3
315	// ModifyGraph marker=19
316	// ModifyGraph rgb=(0,0,0)
317	// Modifygraph log=1
318	// Modifygraph grid=1
319	// Modifygraph mirror=2
320	// ModifyGraph msize(aveint)=2
321	// ErrorBars/T=0 aveint Y,wave=(sigave,sigave)
322	// Label bottom, "Q (1/A)"
323	// Label left, "Relative Intensity"
324	// legend
325	//
326	// RenameWindow #,G_aveint
327	// SetActiveSubwindow ##
328
329	//// end panel commands
330
331	/// draw as a graph
332	String fldrSav0= GetDataFolder(1)
333	SetDataFolder root:Packages:NIST:SAS:
334
335	Display/W=(5,44,463,570)/K=1 aveint vs qval as "SASCALC"
336	DoWindow/C SASCALC
337	ModifyGraph cbRGB=(49151,53155,65535)
338	ModifyGraph mode=3
339	ModifyGraph marker=19
340	ModifyGraph rgb=(0,0,0)
341	Modifygraph log=1
342	Modifygraph grid=1
343	Modifygraph mirror=2
344	ModifyGraph msize(aveint)=2
345	ErrorBars/T=0 aveint Y,wave=(sigave,sigave)
346	Label bottom, "Q (1/A)"
347	Label left, "Relative Intensity"
348	legend
349
350	ControlBar/T 200
351
352	SetDataFolder fldrSav0
353
354	Slider SC_Slider,pos={11,46},size={150,45},proc=GuideSliderProc,live=0
355	Slider SC_Slider,limits={0,8,1},variable= root:Packages:NIST:SAS:gNg,vert= 0//,thumbColor= (1,16019,65535)
356	Slider SC_Slider_1,pos={234,45},size={150,45},proc=DetDistSliderProc,live=0
357	Slider SC_Slider_1,tkLblRot= 90,userTicks={root:Packages:NIST:SAS:tickSDDNG3,root:Packages:NIST:SAS:lblSDDNG3 }
358	Slider SC_Slider_1,limits={133,1317,1},variable= root:Packages:NIST:SAS:gDetDist,vert= 0//,thumbColor= (1,16019,65535)
359	Slider SC_Slider_2,pos={394,21},size={47,65},proc=OffsetSliderProc,live=0,ticks=4
360	Slider SC_Slider_2,limits={0,25,1},variable= root:Packages:NIST:SAS:gOffset//,thumbColor= (1,16019,65535)
361	CheckBox checkNG3,pos={20,19},size={36,14},proc=SelectInstrumentCheckProc,title="NG3"
362	CheckBox checkNG3,value=1,mode=1
363	CheckBox checkNG7,pos={66,19},size={36,14},proc=SelectInstrumentCheckProc,title="NG7"
364	CheckBox checkNG7,value=0,mode=1
365	CheckBox checkChamber,pos={172,48},size={57,14},proc=TableCheckProc,title="Chamber"
366	CheckBox checkChamber,value=1,mode=1
367	CheckBox checkHuber,pos={172,27},size={44,14},proc=TableCheckProc,title="Huber"
368	CheckBox checkHuber,value=0,mode=1
369	PopupMenu popup0,pos={6,94},size={76,20},proc=SourceAperturePopMenuProc
370	PopupMenu popup0,mode=1,popvalue="3.81 cm",value= root:Packages:NIST:SAS:gSourceApString
371	PopupMenu popup0_1,pos={172,72},size={49,20},proc=SampleAperturePopMenuProc
372	PopupMenu popup0_1,mode=8,popvalue="1/2\"",value= root:Packages:NIST:SAS:gApPopStr
373	SetVariable setvar0,pos={301,107},size={130,15},title="Det Dist (cm)",proc=SDDSetVarProc
374	SetVariable setvar0,limits={133,1317,1},value=root:Packages:NIST:SAS:gDetDist
375	SetVariable setvar0_1,pos={321,129},size={110,15},title="Offset (cm)",proc=OffsetSetVarProc
376	SetVariable setvar0_1,limits={0,25,1},value= root:Packages:NIST:SAS:gOffset
377	SetVariable setvar0_2,pos={6,130},size={90,15},title="Lambda",proc=LambdaSetVarProc
378	SetVariable setvar0_2,limits={4,20,0.1},value= root:Packages:NIST:SAS:gLambda
379	PopupMenu popup0_2,pos={108,127},size={55,20},proc=DeltaLambdaPopMenuProc
380	PopupMenu popup0_2,mode=2,value= root:Packages:NIST:SAS:gDeltaLambdaStr
381	Button FreezeButton title="Freeze",size={60,20},pos={180,166}
382	Button FreezeButton proc=FreezeButtonProc
383	Button ClearButton title="Clear",size={60,20},pos={250,166}
384	Button ClearButton proc=S_ClearButtonProc
385	GroupBox group0,pos={6,1},size={108,36},title="Instrument"
386	SetDataFolder fldrSav0
387
388	SetVariable setvar0_3,pos={140,94},size={110,15},title="Diam (mm)",disable=1
389	SetVariable setvar0_3,limits={0,100,0.1},value= root:Packages:NIST:SAS:gSamApOther,proc=SampleApOtherSetVarProc
390
391	CheckBox checkLens,pos={6,155},size={44,14},proc=LensCheckProc,title="Lenses?"
392	CheckBox checkLens,value=root:Packages:NIST:SAS:gUsingLenses
393
394	CheckBox checkSim,pos={6,175},size={44,14},proc=SimCheckProc,title="Simulation?"
395	CheckBox checkSim,value=0
396
397	// help, done buttons
398	Button SC_helpButton,pos={340,166},size={25,20},proc=showSASCALCHelp,title="?"
399	Button SC_helpButton,help={"Show help file for simulation of SANS Data"}
400	Button SC_DoneButton,pos={380,166},size={50,20},proc=SASCALCDoneButton,title="Done"
401	Button SC_DoneButton,help={"This button will close the panel"}
402
403	// set up a fake dependency to trigger recalculation
404	//root:Packages:NIST:SAS:gCalculate := ReCalculateInten(root:Packages:NIST:SAS:gTouched)
405	EndMacro
406
407	//clean up
408	Proc SASCALCDoneButton(ctrlName): ButtonControl
409	String ctrlName
410	DoWindow/K SASCALC
411	DoWindow/K Trial_Configuration
412	DoWindow/K Saved_Configurations
413	DoWindow/K MC_SASCALC
414	DoWindow/K Sim_1D_Panel
415	end
416
417	//
418	Proc showSASCALCHelp(ctrlName): ButtonControl
419	String ctrlName
420	DisplayHelpTopic/K=1/Z "SASCALC"
421	if(V_flag !=0)
422	DoAlert 0,"The SANS Simulation Help file could not be found"
423	endif
424	end
425
426	// based on the instrument selected...
427	//set the SDD range
428	// set the source aperture popup (based on NGx and number of guides)
429	// set the wavelength spread popup
430	//
431	Function UpdateControls()
432	//poll the controls on the panel, and change needed values
433	Variable isNG3,Ng,mode
434	ControlInfo/W=SASCALC checkNG3
435	isNG3=V_value
436	ControlInfo/W=SASCALC SC_slider
437	Ng=V_value
438	SVAR A1str= root:Packages:NIST:SAS:gSourceApString// = "1.43 cm;2.54 cm;3.81 cm;"
439	SVAR dlStr = root:Packages:NIST:SAS:gDeltaLambdaStr
440	if(isNG3)
441	switch(ng)
442	case 0:
443	ControlInfo/W=SASCALC popup0
444	mode=V_value
445	A1str="1.43 cm;2.54 cm;3.81 cm;"
446	break
447	case 6:
448	A1str = "! 6 Guides invalid;"
449	mode=1
450	break
451	case 7: // switched order in 2009 to keep 5 cm as default, 2.5 cm for polarizer
452	A1Str = "5.00 cm;2.50 cm;"
453	mode = 1
454	break
455	default:
456	A1str = "5 cm;"
457	mode=1
458	endswitch
459	//wavelength spread
460	dlStr = "0.109;0.125;0.236;" //updated calibration 2009
461	//detector limits
462	SetVariable setvar0,limits={133,1317,1}
463	NVAR detDist=root:Packages:NIST:SAS:gDetDist
464	if(detDist < 133 )
465	detDist = 133
466	elseif (detDist > 1317 )
467	detDist = 1317
468	endif
469	Slider SC_Slider_1,limits={133,1317,1},userTicks={root:Packages:NIST:SAS:tickSDDNG3,root:Packages:NIST:SAS:lblSDDNG3 }
470	Slider SC_Slider_1,variable=root:Packages:NIST:SAS:gDetDist //forces update
471	else //ng7
472	switch(ng)
473	case 0:
474	ControlInfo/W=SASCALC popup0
475	mode=V_value
476	A1str="1.43 cm;2.22 cm;3.81 cm;"
477	break
478	default:
479	A1str = "5.08 cm;"
480	mode=1
481	endswitch
482
483	dlStr = "0.09;0.11;0.22;"
484	Slider SC_Slider_1,limits={100,1531,1},userTicks={root:Packages:NIST:SAS:tickSDDNG7,root:Packages:NIST:SAS:lblSDDNG7 }
485	SetVariable setvar0,limits={100,1531,1}
486	Slider SC_Slider_1,variable=root:Packages:NIST:SAS:gDetDist //forces update
487	endif
488	ControlUpdate popup0
489	PopupMenu popup0,mode=mode //source Ap
490	ControlInfo/W=SASCALC popup0
491	SourceAperturePopMenuProc("",0,S_Value) //send popNum==0 so recalculation won't be done
492
493	ControlUpdate popup0_2 // delta lambda pop
494	ControlInfo/W=SASCALC popup0_2
495	DeltaLambdaPopMenuProc("",0,S_Value) //sets the global and the wave
496	End
497
498	//changing the number of guides changes the SSD
499	// the source aperture popup may need to be updated
500	//
501	Function GuideSliderProc(ctrlName,sliderValue,event)
502	String ctrlName
503	Variable sliderValue
504	Variable event // bit field: bit 0: value set, 1: mouse down, 2: mouse up, 3: mouse moved
505
506	Variable recalc=0
507
508	if(event %& 0x1) // bit 0, value set
509	if(cmpstr(ctrlName,"") != 0) //here by direct action, so do LensCheck and recalculate
510	recalc=1
511	LensCheckProc("",2) //make sure lenses are deselected
512	endif
513	sourceToSampleDist() //updates the SSD global and wave
514	//change the sourceAp popup, SDD range, etc
515	UpdateControls()
516	ReCalculateInten(recalc)
517	endif
518	return 0
519	End
520
521	// changing the detector position changes the SDD
522	//
523	Function DetDistSliderProc(ctrlName,sliderValue,event)
524	String ctrlName
525	Variable sliderValue
526	Variable event // bit field: bit 0: value set, 1: mouse down, 2: mouse up, 3: mouse moved
527
528	Variable recalc=0
529	if(event %& 0x1) // bit 0, value set
530	if(cmpstr(ctrlName,"") != 0)
531	recalc=1
532	LensCheckProc("",2) //make sure lenses are only selected for valid configurations
533	endif
534	sampleToDetectorDist() //changes the SDD and wave (DetDist is the global)
535	ReCalculateInten(recalc)
536	endif
537
538	return 0
539	End
540
541	// change the offset
542	// - changes the beamcenter (x,y) position too
543	Function OffsetSliderProc(ctrlName,sliderValue,event)
544	String ctrlName
545	Variable sliderValue
546	Variable event // bit field: bit 0: value set, 1: mouse down, 2: mouse up, 3: mouse moved
547
548	if(event %& 0x1) // bit 0, value set
549	detectorOffset()
550	ReCalculateInten(1)
551	endif
552
553	return 0
554	End
555
556	// changing the instrument -
557	// re-initialize the global values
558	// update the controls to appropriate limits/choices
559	//
560	Function SelectInstrumentCheckProc(ctrlName,checked) : CheckBoxControl
561	String ctrlName
562	Variable checked
563
564	if(cmpstr(ctrlName,"checkNG3")==0)
565	checkBox checkNG3, value=1
566	checkBox checkNG7, value=0
567	initNG3()
568	else
569	checkBox checkNG3, value=0
570	checkBox checkNG7, value=1
571	initNG7()
572	endif
573	LensCheckProc("",2) //check if lenses are still valid (they won't be)
574	UpdateControls()
575	ReCalculateInten(1)
576	End
577
578	//changing the table position
579	// changes the SSD and SDD
580	Function TableCheckProc(ctrlName,checked) : CheckBoxControl
581	String ctrlName
582	Variable checked
583
584	NVAR table=root:Packages:NIST:SAS:gTable
585	if(cmpstr(ctrlName,"checkHuber")==0)
586	checkBox checkHuber, value=1
587	checkBox checkChamber, value=0
588	table=1 //in Huber position
589	else
590	checkBox checkHuber, value=0
591	checkBox checkChamber, value=1
592	table = 2 //in Sample chamber
593	endif
594	sampleToDetectorDist()
595	sourceToSampleDist() //update
596	ReCalculateInten(1)
597	End
598
599
600	//lenses (or prisms) in/out changes resolution
601	//
602	// passing in a checked == 2 will do a check of the configuration without
603	// affecting the box state
604	//
605	// if an invalid configuration is detected, the box is unchecked
606	// and the lens flag is set to zero (=out).
607	//
608	// When necessary controls are "popped", a ctrlName="" is passed to signal the control
609	// to NOT recalculate the intensity, and to NOT (recursively) call LensCheckProc again
610	//
611	// currently, the 17.2 A for lens/prism @ ng3 must be typed in
612	//
613	Function LensCheckProc(ctrlName,checked) : CheckBoxControl
614	String ctrlName
615	Variable checked
616
617	// don't let the box get checked if the conditions are wrong
618	// lambda != 8.09,8.4,17.2
619	// Ng != 0
620	NVAR lens = root:Packages:NIST:SAS:gUsingLenses
621	NVAR ng = root:Packages:NIST:SAS:gNg
622	NVAR lam = root:Packages:NIST:SAS:gLambda
623	NVAR dist = root:Packages:NIST:SAS:gDetDist
624	NVAR instrument = root:Packages:NIST:SAS:instrument
625	Wave rw=root:Packages:NIST:SAS:realsRead
626
627	// directly uncheck the box, just set the flag and get out
628	if(checked == 0)
629	lens = 0
630	CheckBox checkLens,value=0
631	rw[28]=0 //flag for lenses out
632	ReCalculateInten(1)
633	return(0)
634	endif
635
636	// check the box, enforce the proper conditions
637	if(checked == 1)
638	lens = checked
639	if(instrument == 3)
640	dist = 1317
641	DetDistSliderProc("",1317,1)
642
643	lam = 8.4
644	LambdaSetVarProc("",8.4,"8.4","")
645
646	ng=0
647	GuideSliderProc("",0,1) //this updates the controls to the new # of guides
648
649	PopupMenu popup0,mode=1,popvalue="1.43 cm" //first item in source aperture menu
650
651	PopupMenu popup0_2,mode=2 //deltaLambda
652	ControlInfo popup0_2
653	DeltaLambdaPopMenuProc("",0,S_value) //zero as 2nd param skips recalculation
654	else
655	dist = 1531
656	DetDistSliderProc("",1531,1)
657
658	lam = 8.09
659	LambdaSetVarProc("",8.09,"8.09","")
660
661	ng=0
662	GuideSliderProc("",0,1)
663	PopupMenu popup0,mode=1,popvalue="1.43 cm" //first item
664
665	PopupMenu popup0_2,mode=2 //deltaLambda
666	ControlInfo popup0_2
667	DeltaLambdaPopMenuProc("",0,S_value) //zero as 2nd param skips recalculation
668	endif
669	rw[28]=1 //flag for lenses in (not the true number, but OK)
670	ReCalculateInten(1)
671	endif
672
673	// this is my internal check to see if conditions are still valid
674	// I'll uncheck as needed
675	if(checked == 2)
676
677	// source aperture must be 1.43 cm diameter
678	// number of guides must be zero
679	Variable a1 = sourceApertureDiam()
680	if(a1 != 1.43 \|\| Ng !=0)
681	lens = 0
682	CheckBox checkLens,value=0
683	rw[28]=0 //flag for lenses out
684	return(0)
685	endif
686
687	// instrument specific distance requirements
688	if(instrument == 3 && dist != 1317)
689	lens = 0
690	CheckBox checkLens,value=0
691	rw[28]=0 //flag for lenses out
692	return(0)
693	endif
694
695	if(instrument == 7 && dist != 1531)
696	lens = 0
697	CheckBox checkLens,value=0
698	rw[28]=0 //flag for lenses out
699	return(0)
700	endif
701
702	// instrument specific wavelength requirements
703	if(instrument == 3 && !(lam == 8.4 \|\| lam == 17.2) )
704	lens = 0
705	CheckBox checkLens,value=0
706	rw[28]=0 //flag for lenses out
707	return(0)
708	endif
709
710	if(instrument == 7 && lam != 8.09 )
711	lens = 0
712	CheckBox checkLens,value=0
713	rw[28]=0 //flag for lenses out
714	return(0)
715	endif
716
717	endif
718
719	return(1) //return value not used
720	End
721
722	//simulation controls as a control bar that toggles on/off to the right
723	//
724	// depending on the state of the 2D flag, open the 1d or 2d control panel
725	Function SimCheckProc(ctrlName,checked) : CheckBoxControl
726	String ctrlName
727	Variable checked
728
729	if(checked)
730	NVAR do2D = root:Packages:NIST:SAS:gDoMonteCarlo
731
732	if(do2D)
733	DoWindow/F MC_SASCALC
734	if(V_flag==0)
735	Execute "MC_SASCALC()" //sets the variable
736	AutoPositionWindow/M=1/R=SASCALC MC_SASCALC
737	endif
738	else
739	//draw 1D panel
740	DoWindow/F Sim_1D_Panel
741	if(V_flag==0)
742	Execute "Sim_1D_Panel()"
743	AutoPositionWindow/M=1/R=SASCALC Sim_1D_Panel
744	endif
745	endif
746	NVAR doSim = root:Packages:NIST:SAS:doSimulation
747	doSim=1
748	else
749	//get rid of the controls (just try both)
750	DoWindow MC_SASCALC
751	if(V_flag!=0)
752	NVAR do2D = root:Packages:NIST:SAS:gDoMonteCarlo
753	do2D = 0
754	KillWindow MC_SASCALC
755	endif
756	DoWindow Sim_1D_Panel
757	if(V_flag!=0)
758	KillWindow Sim_1D_Panel
759	endif
760
761	NVAR doSim = root:Packages:NIST:SAS:doSimulation
762	doSim=0
763	endif
764
765	return(0)
766	End
767
768	// change the source aperture
769	//
770	Function SourceAperturePopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
771	String ctrlName
772	Variable popNum
773	String popStr
774
775	Variable a1 = sourceApertureDiam() //sets the new value in the wave
776
777	// if LensCheckProc is calling this, don't call it again.
778	// if this is called by a direct pop, then check
779	if(cmpstr(ctrlName,"") != 0)
780	LensCheckProc("",2) //make sure lenses are only selected for valid configurations
781	endif
782
783	ReCalculateInten(popnum) //skip the recalculation if I pass in a zero
784	End
785
786	// set sample aperture
787	//
788	Function SampleAperturePopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
789	String ctrlName
790	Variable popNum
791	String popStr
792
793	sampleApertureDiam()
794	ReCalculateInten(1)
795	End
796
797	// set sample to detector distance
798	//
799	Function SDDSetVarProc(ctrlName,varNum,varStr,varName) : SetVariableControl
800	String ctrlName
801	Variable varNum
802	String varStr
803	String varName
804
805	sampleToDetectorDist()
806	if(cmpstr(ctrlName,"") != 0)
807	LensCheckProc("",2) //make sure lenses are only selected for valid configurations
808	endif
809	ReCalculateInten(1)
810	End
811
812	// set offset
813	// -- also changes the beamcenter (x,y)
814	//
815	Function OffsetSetVarProc(ctrlName,varNum,varStr,varName) : SetVariableControl
816	String ctrlName
817	Variable varNum
818	String varStr
819	String varName
820
821	detectorOffset() //sets the offset in the wave and also the new (x,y) beamcenter
822
823	ReCalculateInten(1)
824	return(0)
825	End
826
827	// change the wavelength
828	Function LambdaSetVarProc(ctrlName,varNum,varStr,varName) : SetVariableControl
829	String ctrlName
830	Variable varNum
831	String varStr
832	String varName
833
834	WAVE rw=root:Packages:NIST:SAS:realsRead
835	Variable recalc=0
836
837	rw[26] = str2num(varStr)
838	if(cmpstr(ctrlName,"") != 0)
839	recalc=1
840	LensCheckProc("",2) //make sure lenses are only selected for valid configurations
841	endif
842	ReCalculateInten(recalc)
843	return(0)
844	End
845
846
847	// change the wavelength spread
848	Function DeltaLambdaPopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
849	String ctrlName
850	Variable popNum
851	String popStr
852
853	NVAR dl=root:Packages:NIST:SAS:gDeltaLambda
854	dl=str2num(popStr)
855	WAVE rw=root:Packages:NIST:SAS:realsRead
856	rw[27] = dl
857	ReCalculateInten(popnum) //skip the calculation if I pass in zero
858	return(0)
859	End
860
861
862	// calculate the intensity
863	// - either do MC or the straight calculation.
864	//
865	// *****currently the smeared calculation is turned off, by not prepending "fSmeared" to the FUNCREFs)
866	// --- the random deviate can't be calculated for the smeared model since I don't know the resolution
867	// function over an infinite q-range, just the detector. Maybe the interpolation is OK, but I
868	// don't really have a good way of testing this. Also, the resolution calculation explicitly multiplies
869	// by fShad, and this wrecks the random deviate calculation. The 2D looks great, but the probabilities
870	// are all wrong. fShad is only appropriate post-simulation.
871	//
872	//
873	// (--NO--) ALWAYS DOES THE RESOLUTION SMEARED CALCULATION
874	// even though the unsmeared model is plotted. this is more realistic to
875	// present to users that are planning base on what they will see in an experiment.
876	//
877	// some bits of the calculation are in a root:Simulation folder that are needed for the resolution smeared calculation
878	// all other bits (as possible) are in the SAS folder (a 2D work folder)
879	//
880	// passing in one does the calculation, "normal" or MC, depending on the global. Normal calculation is the default
881	// passing in zero from a control skips the calculation
882	//
883	Function ReCalculateInten(doIt)
884	Variable doIt
885
886	if(doIt==0)
887	return(0)
888	endif
889
890	// update the wave with the beamstop diameter here, since I don't know what
891	// combinations of parameters will change the BS - but anytime the curve is
892	// recalculated, or the text displayed, the right BS must be present
893	beamstopDiam()
894
895	// generate the resolution waves first, so they are present for a smearing calculation
896	// average the "fake" 2d data now to generate the smearing information
897	S_CircularAverageTo1D("SAS")
898	WAVE aveint=root:Packages:NIST:SAS:aveint
899	WAVE qval=root:Packages:NIST:SAS:qval
900	WAVE sigave=root:Packages:NIST:SAS:sigave
901	WAVE SigmaQ=root:Packages:NIST:SAS:sigmaQ
902	WAVE qbar=root:Packages:NIST:SAS:qbar
903	WAVE fSubS=root:Packages:NIST:SAS:fSubS
904
905	//generate a "fake" 1d data folder/set named "Simulation"
906	Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,"Simulation")
907
908	// do the simulation here, or not
909	Variable r1,xCtr,yCtr,sdd,pixSize,wavelength
910	String coefStr,abortStr,str
911
912	// now the cases are: simulation (0\|1), 1D (default) or 2D (hidden)
913	NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo // == 1 if 2D MonteCarlo set by hidden flag
914	NVAR doSimulation = root:Packages:NIST:SAS:doSimulation // == 1 if 1D simulated data, 0 if other from the checkbox
915	SVAR funcStr = root:Packages:NIST:SAS:gFuncStr //set by the popup
916
917	if(doSimulation == 1)
918	if(doMonteCarlo == 1)
919	//2D simulation (in MultiScatter_MonteCarlo_2D.ipf)
920
921	Simulate_2D_MC(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
922
923	//end 2D simulation
924	else
925	//1D simulation
926
927	if(exists(funcStr) != 0)
928	FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr) //a wrapper for the structure version
929	FUNCREF SANSModelAAO_MCproto funcUnsmeared=$(funcStr) //unsmeared
930	coefStr = MC_getFunctionCoef(funcStr)
931
932	if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
933	Abort "Function and coefficients do not match. You must plot the unsmeared function before simulation."
934	endif
935
936	Wave inten=$"root:Simulation:Simulation_i" // this will exist and send the smeared calculation to the corect DF
937
938	// the resolution-smeared intensity is calculated, including the incoherent background
939	func($coefStr,inten,qval)
940
941	NVAR imon = root:Packages:NIST:SAS:gImon
942	NVAR ctTime = root:Packages:NIST:SAS:gCntTime
943	NVAR thick = root:Packages:NIST:SAS:gThick
944	NVAR trans = root:Packages:NIST:SAS:gSamTrans
945	NVAR SimDetCts = root:Packages:NIST:SAS:g_1DTotCts //summed counts (simulated)
946	NVAR estDetCR = root:Packages:NIST:SAS:g_1DEstDetCR // estimated detector count rate
947	NVAR fracScat = root:Packages:NIST:SAS:g_1DFracScatt // fraction of beam captured on detector
948	NVAR estTrans = root:Packages:NIST:SAS:g_1DEstTrans // estimated transmission of sample
949	NVAR SimCountTime = root:Packages:NIST:SAS:gCntTime //counting time used for simulation
950
951	WAVE rw=root:Packages:NIST:SAS:realsRead
952	WAVE nCells=root:Packages:NIST:SAS:nCells
953
954	pixSize = rw[10]/10 // convert pix size in mm to cm
955	sdd = rw[18]*100 //convert header of [m] to [cm]
956	wavelength = rw[26] // in 1/A
957
958	imon = beamIntensity()
959
960	// calculate the scattering cross section simply to be able to estimate the transmission
961	Variable sig_sas=0
962
963	// remember that the random deviate is the coherent portion ONLY - the incoherent background is
964	// subtracted before the calculation.
965	CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
966
967	// if(sig_sas > 100)
968	// sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas
969	// endif
970	estTrans = exp(-1thicksig_sas) //thickness and sigma both in units of cm
971	Print "Sig_sas = ",sig_sas
972
973	Duplicate/O qval prob_i,countsInAnnulus
974
975	// not needed - nCells takes care of this when the error is correctly calculated
976	// Duplicate/O qval circle_fraction,rval,nCells_expected
977	// rval = sddtan(2asin(qval*wavelength/4/pi)) //radial distance in cm
978	// nCells_expected = 2pirval/pixSize //does this need to be an integer?
979	// circle_fraction = nCells / nCells_expected
980
981
982	// prob_i = transthicknCells(pixSize/sdd)^2inten //probability of a neutron in q-bin(i) that has nCells
983	prob_i = transthick(pixSize/sdd)^2*inten //probability of a neutron in q-bin(i)
984
985	Variable P_on = sum(prob_i,-inf,inf)
986	Print "P_on = ",P_on
987
988	// fracScat = P_on
989	fracScat = 1-estTrans
990
991	// aveint = (ImonctTime)prob_i / circle_fraction / nCells_expected
992	aveint = (ImonctTime)prob_i
993
994	countsInAnnulus = aveint*nCells
995	SimDetCts = sum(countsInAnnulus,-inf,inf)
996	estDetCR = SimDetCts/SimCountTime
997
998
999	NVAR doABS = root:Packages:NIST:SAS:g_1D_DoABS
1000	NVAR addNoise = root:Packages:NIST:SAS:g_1D_AddNoise
1001
1002	sigave = sqrt(aveint/nCells) // corrected based on John's memo, from 8/9/99
1003
1004	// add in random error in aveint based on the sigave
1005	if(addNoise)
1006	aveint += gnoise(sigave)
1007	endif
1008
1009	// convert to absolute scale
1010	if(doABS)
1011	Variable kappa = thick(pixSize/sdd)^2transiMonctTime
1012	aveint /= kappa
1013	sigave /= kappa
1014	endif
1015
1016	else
1017	//no function plotted, no simulation can be done
1018	DoAlert 0,"No function is selected or plotted, so no simulation is done. The default Debye function is used."
1019
1020	aveint = S_Debye(1000,100,0.0,qval)
1021	aveint *= fSubS // multiply either estimate by beamstop shadowing
1022	sigave = 0 //reset for model calculation
1023	endif
1024
1025	endif // end 1D simulation
1026	else
1027	//no simulation
1028
1029	aveint = S_Debye(1000,100,0.0,qval)
1030	aveint *= fSubS // multiply either estimate by beamstop shadowing
1031	sigave = 0 //reset for model calculation
1032
1033	//end no simulation
1034	endif
1035
1036
1037	//display the configuration text in a separate notebook
1038	DisplayConfigurationText()
1039
1040	return(0)
1041	End
1042
1043	//freezes the current configuration
1044	// -1- duplicates the trace on the graph
1045	// -2- copies the configuration text to a second notebook window for printing
1046	Function FreezeButtonProc(ctrlName) : ButtonControl
1047	String ctrlName
1048
1049	String str=""
1050	NVAR ct=root:Packages:NIST:SAS:gFreezeCount
1051
1052
1053	SetDataFolder root:Packages:NIST:SAS
1054
1055	Duplicate/O aveint,$("aveint_"+num2str(ct))
1056	Duplicate/O qval,$("qval_"+num2str(ct))
1057	Duplicate/O sigave,$("sigave_"+num2str(ct))
1058	Appendtograph $("aveint_"+num2str(ct)) vs $("qval_"+num2str(ct))
1059	ModifyGraph mode($("aveint_"+num2str(ct)))=3
1060	ModifyGraph marker($("aveint_"+num2str(ct)))=19
1061	ModifyGraph msize($("aveint_"+num2str(ct)))=2
1062	ErrorBars/T=0 $("aveint_"+num2str(ct)) Y,wave=($("sigave_"+num2str(ct)),$("sigave_"+num2str(ct)))
1063
1064	switch(mod(ct,10)) // 10 different colors - black is the unfrozen color
1065	case 0:
1066	ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,16385,16385)
1067	break
1068	case 1:
1069	ModifyGraph rgb($("aveint_"+num2str(ct)))=(2,39321,1)
1070	break
1071	case 2:
1072	ModifyGraph rgb($("aveint_"+num2str(ct)))=(0,0,65535)
1073	break
1074	case 3:
1075	ModifyGraph rgb($("aveint_"+num2str(ct)))=(39321,1,31457)
1076	break
1077	case 4:
1078	ModifyGraph rgb($("aveint_"+num2str(ct)))=(48059,48059,48059)
1079	break
1080	case 5:
1081	ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,32768,32768)
1082	break
1083	case 6:
1084	ModifyGraph rgb($("aveint_"+num2str(ct)))=(0,65535,0)
1085	break
1086	case 7:
1087	ModifyGraph rgb($("aveint_"+num2str(ct)))=(16385,65535,65535)
1088	break
1089	case 8:
1090	ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,32768,58981)
1091	break
1092	case 9:
1093	ModifyGraph rgb($("aveint_"+num2str(ct)))=(36873,14755,58982)
1094	break
1095	endswitch
1096
1097	NVAR doTraceOffset = root:Packages:NIST:SAS:gDoTraceOffset
1098	NVAR offset = root:Packages:NIST:SAS:gModelOffsetFactor
1099	if(doTraceOffset)
1100	offset = 2^ct
1101	//multiply by current offset (>=1)
1102	Wave inten = $("aveint_"+num2str(ct))
1103	inten *= offset
1104	// Print "new offset = ",offset
1105	endif
1106
1107	ct +=1
1108	SetDataFolder root:
1109
1110	// create or append the configuration to a new notebook
1111	if(WinType("Saved_Configurations")==0)
1112	NewNotebook/F=1/N=Saved_Configurations /W=(480,400,880,725)
1113	DoWindow/F SASCALC //return focus to SASCALC window
1114	endif
1115	//append the text
1116	sprintf str,"\rConfiguration #%d\r",ct-1
1117	Notebook Saved_Configurations showRuler=0,defaultTab=20,selection={endOfFile, endOfFile}
1118	Notebook Saved_Configurations font="Monaco",fSize=10,fstyle=1,text=str //bold
1119	Notebook Saved_Configurations font="Monaco",fSize=10,fstyle=0,text=SetConfigurationText()
1120	return(0)
1121	End
1122
1123	//clears the frozen traces on the graph, asks if you want to clear the saved text
1124	//
1125	Function S_ClearButtonProc(ctrlName) : ButtonControl
1126	String ctrlName
1127
1128	NVAR ct=root:Packages:NIST:SAS:gFreezeCount
1129	Variable ii
1130	Setdatafolder root:Packages:NIST:SAS
1131	for(ii=ct-1;ii>=1;ii-=1)
1132	//remove all traces, replace aveint
1133	// kill all waves _ct
1134	RemoveFromGraph $("aveint_"+num2str(ii))
1135	Killwaves/Z $("aveint_"+num2str(ii)),$("qval_"+num2str(ii))
1136	endfor
1137	ct=1
1138	setdatafolder root:
1139
1140	DoAlert 1,"Do you also want to clear the \"Saved Configurations\" text?"
1141	if(V_flag == 1) // yes
1142	DoWindow/K Saved_Configurations
1143	endif
1144
1145	//reset offset value
1146	NVAR offset = root:Packages:NIST:SAS:gModelOffsetFactor
1147	offset = 1
1148	ReCalculateInten(1)
1149	return(0)
1150	End
1151
1152
1153	///////////////////////////////////////////////////////////
1154	// 19MAR07 uses correction for beamstop diameter projection to get shadow factor correct
1155	//
1156	Function S_CircularAverageTo1D(type)
1157	String type
1158
1159	Variable isCircular = 1
1160
1161	//type is the data type to do the averaging on, and will be set as the current folder
1162	//get the current displayed data (so the correct folder is used)
1163	String destPath = "root:Packages:NIST:"+type
1164
1165	//
1166	Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist,dr,ddr
1167	Variable lambda,trans
1168	WAVE reals = $(destPath + ":RealsRead")
1169	// WAVE/T textread = $(destPath + ":TextRead")
1170	// String fileStr = textread[3]
1171
1172	// center of detector, for non-linear corrections
1173	Variable pixelsX=128,pixelsY=128
1174
1175	xcenter = pixelsX/2 + 0.5 // == 64.5 for 128x128 Ordela
1176	ycenter = pixelsY/2 + 0.5 // == 64.5 for 128x128 Ordela
1177
1178	// beam center, in pixels
1179	x0 = reals[16]
1180	y0 = reals[17]
1181	//detector calibration constants
1182	sx = reals[10] //mm/pixel (x)
1183	sx3 = reals[11] //nonlinear coeff
1184	sy = reals[13] //mm/pixel (y)
1185	sy3 = reals[14] //nonlinear coeff
1186
1187	dtsize = 10*reals[20] //det size in mm
1188	dtdist = 1000*reals[18] // det distance in mm
1189
1190	NVAR binWidth=root:Packages:NIST:SAS:gBinWidth
1191	// Variable binWidth = 1
1192
1193	dr = binWidth // *********annulus width set by user, default is one*********
1194	ddr = dr*sx //step size, in mm (this value should be passed to the resolution calculation, not dr 18NOV03)
1195
1196	Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,dxbm,dybm,ii
1197	Variable phi_rad,dphi_rad,phi_x,phi_y
1198	Variable forward,mirror
1199
1200	String side = "both"
1201	// side = StringByKey("SIDE",keyListStr,"=",";")
1202	// Print "side = ",side
1203
1204	/// data wave is data in the current folder which was set at the top of the function
1205	WAVE data=$(destPath + ":data")
1206
1207	// fake mask that uses all of the detector
1208	Make/D/O/N=(pixelsX,pixelsY) $(destPath + ":mask")
1209	Wave mask = $(destPath + ":mask")
1210	mask = 0
1211	//two pixels all around
1212	mask[0,1][] = 1
1213	mask[126,127][] = 1
1214	mask[][0,1] = 1
1215	mask[][126,127] = 1
1216	//
1217	//pixels within rcentr of beam center are broken into 9 parts (units of mm)
1218	rcentr = 100 //original
1219	// rcentr = 0
1220	// values for error if unable to estimate value
1221	//large_num = 1e10
1222	large_num = 1 //1e10 value (typically sig of last data point) plots poorly, arb set to 1
1223	small_num = 1e-10
1224
1225	// output wave are expected to exist (?) initialized to zero, what length?
1226	// 200 points on VAX --- use 300 here, or more if SAXS data is used with 1024x1024 detector (1000 pts seems good)
1227	Variable defWavePts=500
1228	Make/O/D/N=(defWavePts) $(destPath + ":qval"),$(destPath + ":aveint")
1229	Make/O/D/N=(defWavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave")
1230	Make/O/D/N=(defWavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar")
1231
1232	WAVE qval = $(destPath + ":qval")
1233	WAVE aveint = $(destPath + ":aveint")
1234	WAVE ncells = $(destPath + ":ncells")
1235	WAVE dsq = $(destPath + ":dsq")
1236	WAVE sigave = $(destPath + ":sigave")
1237	WAVE qbar = $(destPath + ":QBar")
1238	WAVE sigmaq = $(destPath + ":SigmaQ")
1239	WAVE fsubs = $(destPath + ":fSubS")
1240
1241	qval = 0
1242	aveint = 0
1243	ncells = 0
1244	dsq = 0
1245	sigave = 0
1246	qbar = 0
1247	sigmaq = 0
1248	fsubs = 0
1249
1250	dtdis2 = dtdist^2
1251	nq = 1
1252	xoffst=0
1253	//distance of beam center from detector center
1254	dxbm = S_FX(x0,sx3,xcenter,sx)
1255	dybm = S_FY(y0,sy3,ycenter,sy)
1256
1257	//BEGIN AVERAGE **********
1258	Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1
1259	Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p
1260
1261	// IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
1262	// loop index corresponds to FORTRAN (old code)
1263	// and the IGOR array indices must be adjusted (-1) to the correct address
1264	ii=1
1265	do
1266	xi = ii
1267	dxi = S_FX(xi,sx3,xcenter,sx)
1268	dx = dxi-dxbm //dx and dy are in mm
1269
1270	jj = 1
1271	do
1272	data_pixel = data[ii-1][jj-1] //assign to local variable
1273	yj = jj
1274	dyj = S_FY(yj,sy3,ycenter,sy)
1275	dy = dyj - dybm
1276	if(!(mask[ii-1][jj-1])) //masked pixels = 1, skip if masked (this way works...)
1277	dr2 = (dx^2 + dy^2)^(0.5) //distance from beam center NOTE dr2 used here - dr used above
1278	if(dr2>rcentr) //keep pixel whole
1279	nd = 1
1280	fd = 1
1281	else //break pixel into 9 equal parts
1282	nd = 3
1283	fd = 2
1284	endif
1285	nd2 = nd^2
1286	ll = 1 //"el-el" loop index
1287	do
1288	dxx = dx + (ll - fd)*sx/3
1289	kk = 1
1290	do
1291	dyy = dy + (kk - fd)*sy/3
1292	if(isCircular)
1293	//circular average, use all pixels
1294	//(increment)
1295	nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1296	else
1297	//a sector average - determine azimuthal angle
1298	dphi_p = S_dphi_pixel(dxx,dyy,phi_x,phi_y)
1299	if(dphi_p < dphi_rad)
1300	forward = 1 //within forward sector
1301	else
1302	forward = 0
1303	Endif
1304	if((Pi - dphi_p) < dphi_rad)
1305	mirror = 1 //within mirror sector
1306	else
1307	mirror = 0
1308	Endif
1309	//check if pixel lies within allowed sector(s)
1310	if(cmpstr(side,"both")==0) //both sectors
1311	if ( mirror \|\| forward)
1312	//increment
1313	nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1314	Endif
1315	else
1316	if(cmpstr(side,"right")==0) //forward sector only
1317	if(forward)
1318	//increment
1319	nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1320	Endif
1321	else //mirror sector only
1322	if(mirror)
1323	//increment
1324	nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1325	Endif
1326	Endif
1327	Endif //allowable sectors
1328	Endif //circular or sector check
1329	kk+=1
1330	while(kk<=nd)
1331	ll += 1
1332	while(ll<=nd)
1333	Endif //masked pixel check
1334	jj += 1
1335	while (jj<=pixelsY)
1336	ii += 1
1337	while(ii<=pixelsX) //end of the averaging
1338
1339	//compute q-values and errors
1340	Variable ntotal,rr,theta,avesq,aveisq,var
1341
1342	lambda = reals[26]
1343	ntotal = 0
1344	kk = 1
1345	do
1346	rr = (2kk-1)ddr/2
1347	theta = 0.5*atan(rr/dtdist)
1348	qval[kk-1] = (4Pi/lambda)sin(theta)
1349	if(ncells[kk-1] == 0)
1350	//no pixels in annuli, data unknown
1351	aveint[kk-1] = 0
1352	sigave[kk-1] = large_num
1353	else
1354	if(ncells[kk-1] <= 1)
1355	//need more than one pixel to determine error
1356	aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
1357	sigave[kk-1] = large_num
1358	else
1359	//assume that the intensity in each pixel in annuli is normally
1360	// distributed about mean...
1361	aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
1362	avesq = aveint[kk-1]^2
1363	aveisq = dsq[kk-1]/ncells[kk-1]
1364	var = aveisq-avesq
1365	if(var<=0)
1366	sigave[kk-1] = small_num
1367	else
1368	sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1))
1369	endif
1370	endif
1371	endif
1372	ntotal += ncells[kk-1]
1373	kk+=1
1374	while(kk<=nq)
1375
1376	//Print "NQ = ",nq
1377	// data waves were defined as 300 points (=defWavePts), but now have less than that (nq) points
1378	// use DeletePoints to remove junk from end of waves
1379	//WaveStats would be a more foolproof implementation, to get the # points in the wave
1380	Variable startElement,numElements
1381	startElement = nq
1382	numElements = defWavePts - startElement
1383	DeletePoints startElement,numElements, qval,aveint,ncells,dsq,sigave
1384
1385	//////////////end of VAX sector_ave()
1386
1387
1388
1389	// ***************************************************************
1390	//
1391	// Do the extra 3 columns of resolution calculations starting here.
1392	//
1393	// ***************************************************************
1394
1395	Variable L2 = reals[18]
1396	Variable BS = reals[21] //this the diameter is stored in mm
1397	// SRK - why was I using the projected diameter of the beam stop?? I added a step at the beginning of every recalculation
1398	// of the intensity to get the right beamstop diameter into RealsRead...
1399	// Variable BS = beamstopDiamProjection(1) * 10 //calculated projection in cm *10 = mm
1400	Variable S1 = reals[23]
1401	Variable S2 = reals[24]
1402	Variable L1 = reals[25]
1403	lambda = reals[26]
1404	Variable lambdaWidth = reals[27]
1405
1406	Variable DDet, apOff
1407	//typical value for NG3 and NG7 - distance between sample aperture and sample in (cm)
1408	apOff=5.0
1409	// hard wire value for Ordela detectors
1410	DDet = 0.5 // resolution in cm
1411	// String detStr=textRead[9]
1412	// DDet = DetectorPixelResolution(fileStr,detStr) //needs detector type and beamline
1413
1414	//Go from 0 to nq doing the calc for all three values at
1415	//every Q value
1416
1417	ii=0
1418	Variable ret1,ret2,ret3
1419	do
1420	S_getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,ret1,ret2,ret3)
1421	sigmaq[ii] = ret1 //res_wave[0]
1422	qbar[ii] = ret2 //res_wave[1]
1423	fsubs[ii] = ret3 //res_wave[2]
1424	ii+=1
1425	while(ii<nq)
1426	DeletePoints startElement,numElements, sigmaq, qbar, fsubs
1427
1428	fsubs += 1e-8 //keep the values from being too small
1429
1430	// End of resolution calculations
1431	// ***************************************************************
1432
1433	//get rid of the default mask, if one was created (it is in the current folder)
1434	//don't just kill "mask" since it might be pointing to the one in the MSK folder
1435	Killwaves/Z $(destPath+":mask")
1436
1437	//return to root folder (redundant)
1438	SetDataFolder root:
1439
1440	Return 0
1441	End
1442
1443	//returns nq, new number of q-values
1444	//arrays aveint,dsq,ncells are also changed by this function
1445	//
1446	Function S_IncrementPixel(dataPixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1447	Variable dataPixel,ddr,dxx,dyy
1448	Wave aveint,dsq,ncells
1449	Variable nq,nd2
1450
1451	Variable ir
1452
1453	ir = trunc(sqrt(dxxdxx+dyydyy)/ddr)+1
1454	if (ir>nq)
1455	nq = ir //resets maximum number of q-values
1456	endif
1457	aveint[ir-1] += dataPixel/nd2 //ir-1 must be used, since ir is physical
1458	dsq[ir-1] += dataPixel*dataPixel/nd2
1459	ncells[ir-1] += 1/nd2
1460
1461	Return nq
1462	End
1463
1464	//function determines azimuthal angle dphi that a vector connecting
1465	//center of detector to pixel makes with respect to vector
1466	//at chosen azimuthal angle phi -> [cos(phi),sin(phi)] = [phi_x,phi_y]
1467	//dphi is always positive, varying from 0 to Pi
1468	//
1469	Function S_dphi_pixel(dxx,dyy,phi_x,phi_y)
1470	Variable dxx,dyy,phi_x,phi_y
1471
1472	Variable val,rr,dot_prod
1473
1474	rr = sqrt(dxx^2 + dyy^2)
1475	dot_prod = (dxxphi_x + dyyphi_y)/rr
1476	//? correct for roundoff error? - is this necessary in IGOR, w/ double precision?
1477	if(dot_prod > 1)
1478	dot_prod =1
1479	Endif
1480	if(dot_prod < -1)
1481	dot_prod = -1
1482	Endif
1483
1484	val = acos(dot_prod)
1485
1486	return val
1487
1488	End
1489
1490	//calculates the x distance from the center of the detector, w/nonlinear corrections
1491	//
1492	Function S_FX(xx,sx3,xcenter,sx)
1493	Variable xx,sx3,xcenter,sx
1494
1495	Variable retval
1496
1497	retval = sx3tan((xx-xcenter)sx/sx3)
1498	Return retval
1499	End
1500
1501	//calculates the y distance from the center of the detector, w/nonlinear corrections
1502	//
1503	Function S_FY(yy,sy3,ycenter,sy)
1504	Variable yy,sy3,ycenter,sy
1505
1506	Variable retval
1507
1508	retval = sy3tan((yy-ycenter)sy/sy3)
1509	Return retval
1510	End
1511
1512	//**********************
1513	// Resolution calculation - used by the averaging routines
1514	// to calculate the resolution function at each q-value
1515	// - the return value is not used
1516	//
1517	// equivalent to John's routine on the VAX Q_SIGMA_AVE.FOR
1518	// Incorporates eqn. 3-15 from J. Appl. Cryst. (1995) v. 28 p105-114
1519	//
1520	Function/S S_getResolution(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,SigmaQ,QBar,fSubS)
1521	Variable inQ, lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r
1522	Variable &fSubS, &QBar, &SigmaQ //these are the output quantities at the input Q value
1523
1524	//lots of calculation variables
1525	Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g
1526	Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
1527
1528	//Constants
1529	//Variable del_r = .1
1530	Variable vz_1 = 3.956e5 //velocity [cm/s] of 1 A neutron
1531	Variable g = 981.0 //gravity acceleration [cm/s^2]
1532
1533	String results
1534	results ="Failure"
1535
1536	NVAR usingLenses = root:Packages:NIST:SAS:gUsingLenses
1537
1538	//rename for working variables, these must be gotten from global
1539	//variables
1540
1541	// Variable wLam, wLW, wL1, wL2, wS1, wS2
1542	// Variable wBS, wDDet, wApOff
1543	// wLam = lambda
1544	// wLW = lambdaWidth
1545	// wDDet = DDet
1546	// wApOff = apOff
1547	S1 = 0.50.1 //convert to radius and [cm]
1548	S2 = 0.50.1
1549
1550	L1 *= 100.0 // [cm]
1551	L1 -= apOff //correct the distance
1552
1553	L2 *= 100.0
1554	L2 += apOff
1555
1556	BS = 0.50.1 //convert to radius and [cm]
1557	del_r *= 0.1 //width of annulus, convert mm to [cm]
1558
1559	//Start resolution calculation
1560	a2 = S1L2/L1 + S2(L1+L2)/L1
1561	q_small = 2.0Pi(BS-a2)(1.0-lambdaWidth)/(lambdaL2)
1562	lp = 1.0/( 1.0/L1 + 1.0/L2)
1563
1564	v_lambda = lambdaWidth^2/6.0
1565
1566	if(usingLenses==1) //SRK 2007
1567	v_b = 0.25(S1L2/L1)^2 +0.25(2/3)(lambdaWidth/lambda)^2(S2L2/lp)^2 //correction to 2nd term
1568	else
1569	v_b = 0.25(S1L2/L1)^2 +0.25(S2L2/lp)^2 //original form
1570	endif
1571
1572	v_d = (DDet/2.3548)^2 + del_r^2/12.0
1573	vz = vz_1 / lambda
1574	yg = 0.5gL2*(L1+L2)/vz^2
1575	v_g = 2.0(2.0yg^2*v_lambda) //factor of 2 correction, B. Hammouda, 2007
1576
1577	r0 = L2tan(2.0asin(lambdainQ/(4.0Pi) ))
1578	delta = 0.5*(BS - r0)^2/v_d
1579
1580	if (r0 < BS)
1581	inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
1582	else
1583	inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
1584	endif
1585
1586	fSubS = 0.5(1.0+erf( (r0-BS)/sqrt(2.0v_d) ) )
1587	if (fSubS <= 0.0)
1588	fSubS = 1.e-10
1589	endif
1590	fr = 1.0 + sqrt(v_d)exp(-1.0delta) /(r0fSubSsqrt(2.0*Pi))
1591	fv = inc_gamma/(fSubSsqrt(Pi)) - r0^2(fr-1.0)^2/v_d
1592
1593	rmd = fr*r0
1594	v_r1 = v_b + fv*v_d +v_g
1595
1596	rm = rmd + 0.5*v_r1/rmd
1597	v_r = v_r1 - 0.5*(v_r1/rmd)^2
1598	if (v_r < 0.0)
1599	v_r = 0.0
1600	endif
1601	QBar = (4.0Pi/lambda)sin(0.5*atan(rm/L2))
1602	SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda)
1603
1604	results = "success"
1605	Return results
1606	End
1607
1608	Function S_Debye(scale,rg,bkg,x)
1609	Variable scale,rg,bkg
1610	Variable x
1611
1612	// variables are:
1613	//[0] scale factor
1614	//[1] radius of gyration [A]
1615	//[2] background [cm-1]
1616
1617	// calculates (scale*debye)+bkg
1618	Variable Pq,qr2
1619
1620	qr2=(x*rg)^2
1621	Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
1622
1623	//scale
1624	Pq *= scale
1625	// then add in the background
1626	return (Pq+bkg)
1627	End
1628
1629
1630	Function S_SphereForm(scale,radius,delrho,bkg,x)
1631	Variable scale,radius,delrho,bkg
1632	Variable x
1633
1634	// variables are:
1635	//[0] scale
1636	//[1] radius (A)
1637	//[2] delrho (A-2)
1638	//[3] background (cm-1)
1639
1640	// Variable scale,radius,delrho,bkg
1641	// scale = w[0]
1642	// radius = w[1]
1643	// delrho = w[2]
1644	// bkg = w[3]
1645
1646
1647	// calculates scale * f^2/Vol where f=Vol3delrho*((sin(qr)-qrcos(qr))/qr^3
1648	// and is rescaled to give [=] cm^-1
1649
1650	Variable bes,f,vol,f2
1651	//
1652	//handle q==0 separately
1653	If(x==0)
1654	f = 4/3piradius^3delrhodelrhoscale1e8 + bkg
1655	return(f)
1656	Endif
1657
1658	bes = 3(sin(xradius)-xradiuscos(x*radius))/x^3/radius^3
1659	vol = 4pi/3radius^3
1660	f = volbesdelrho // [=] A
1661	// normalize to single particle volume, convert to 1/cm
1662	f2 = f * f / vol * 1.0e8 // [=] 1/cm
1663
1664	return (scale*f2+bkg) // Scale, then add in the background
1665
1666	End
1667
1668	Function/S SetConfigurationText()
1669
1670	String str="",temp
1671
1672	SetDataFolder root:Packages:NIST:SAS
1673
1674	NVAR numberOfGuides=gNg
1675	NVAR gTable=gTable //2=chamber, 1=table
1676	NVAR wavelength=gLambda
1677	NVAR lambdaWidth=gDeltaLambda
1678	NVAR instrument = instrument
1679	NVAR L2diff = L2diff
1680	NVAR lens = root:Packages:NIST:SAS:gUsingLenses
1681	SVAR aStr = root:myGlobals:gAngstStr
1682
1683	sprintf temp,"Source Aperture Diameter =\t\t%6.2f cm\r",sourceApertureDiam()
1684	str += temp
1685	sprintf temp,"Source to Sample =\t\t\t\t%6.0f cm\r",sourceToSampleDist()
1686	str += temp
1687	sprintf temp,"Sample Aperture to Detector =\t%6.0f cm\r",sampleToDetectorDist() + L2diff
1688	str += temp
1689	sprintf temp,"Beam diameter =\t\t\t\t\t%6.2f cm\r",beamDiameter("maximum")
1690	str += temp
1691	sprintf temp,"Beamstop diameter =\t\t\t\t%6.2f inches\r",beamstopDiam()/2.54
1692	str += temp
1693	sprintf temp,"Minimum Q-value =\t\t\t\t%6.4f 1/%s (sigQ/Q = %3.1f %%)\r",qMin(),aStr,deltaQ(qMin())
1694	str += temp
1695	sprintf temp,"Maximum Horizontal Q-value =\t%6.4f 1/%s\r",qMaxHoriz(),aStr
1696	str += temp
1697	sprintf temp,"Maximum Vertical Q-value =\t\t%6.4f 1/%s\r",qMaxVert(),aStr
1698	str += temp
1699	sprintf temp,"Maximum Q-value =\t\t\t\t%6.4f 1/%s (sigQ/Q = %3.1f %%)\r",qMaxCorner(),aStr,deltaQ(qMaxCorner())
1700	str += temp
1701	sprintf temp,"Beam Intensity =\t\t\t\t%.0f counts/s\r",beamIntensity()
1702	str += temp
1703	sprintf temp,"Figure of Merit =\t\t\t\t%3.3g %s^2/s\r",figureOfMerit(),aStr
1704	str += temp
1705	sprintf temp,"Attenuator transmission =\t\t%3.3g = Atten # %d\r",attenuatorTransmission(),attenuatorNumber()
1706	str += temp
1707	//
1708	// // add text of the user-edited values
1709	// //
1710	sprintf temp,"*************** NG %d ***************\r",instrument
1711	str += temp
1712	sprintf temp,"Sample Aperture Diameter =\t\t\t\t%.2f cm\r",sampleApertureDiam()
1713	str += temp
1714	sprintf temp,"Number of Guides =\t\t\t\t\t\t%d \r", numberOfGuides
1715	str += temp
1716	sprintf temp,"Sample Chamber to Detector =\t\t\t%.1f cm\r", chamberToDetectorDist()
1717	str += temp
1718	if(gTable==1)
1719	sprintf temp,"Sample Position is \t\t\t\t\t\tHuber\r"
1720	else
1721	sprintf temp,"Sample Position is \t\t\t\t\t\tChamber\r"
1722	endif
1723	str += temp
1724	sprintf temp,"Detector Offset =\t\t\t\t\t\t%.1f cm\r", detectorOffset()
1725	str += temp
1726	sprintf temp,"Neutron Wavelength =\t\t\t\t\t%.2f %s\r", wavelength,aStr
1727	str += temp
1728	sprintf temp,"Wavelength Spread, FWHM =\t\t\t\t%.3f\r", lambdaWidth
1729	str += temp
1730	sprintf temp,"Sample Aperture to Sample Position =\t%.2f cm\r", L2Diff
1731	str += temp
1732	if(lens==1)
1733	sprintf temp,"Lenses are IN\r"
1734	else
1735	sprintf temp,"Lenses are OUT\r"
1736	endif
1737	str += temp
1738
1739	setDataFolder root:
1740	return str
1741	End
1742
1743	Function DisplayConfigurationText()
1744
1745	if(WinType("Trial_Configuration")==0)
1746	NewNotebook/F=0/K=1/N=Trial_Configuration /W=(480,44,880,369)
1747	endif
1748	//replace the text
1749	Notebook Trial_Configuration selection={startOfFile, endOfFile}
1750	Notebook Trial_Configuration font="Monaco",fSize=10,text=SetConfigurationText()
1751	return(0)
1752	end
1753
1754	//parses the control for A1 diam
1755	// updates the wave
1756	Function sourceApertureDiam()
1757	ControlInfo/W=SASCALC popup0
1758	Variable diam
1759	sscanf S_Value, "%g cm", diam
1760	WAVE rw=root:Packages:NIST:SAS:realsRead
1761	rw[23] = diam*10 //sample aperture diameter in mm
1762	return(diam)
1763	end
1764
1765	// change the sample aperture to a non-standard value
1766	Function SampleApOtherSetVarProc(ctrlName,varNum,varStr,varName) : SetVariableControl
1767	String ctrlName
1768	Variable varNum
1769	String varStr
1770	String varName
1771
1772	WAVE rw=root:Packages:NIST:SAS:realsRead
1773	rw[24] = varNum //sample aperture diameter in mm
1774	ReCalculateInten(1)
1775	return(0)
1776	End
1777
1778	//parses the control for A2 diam
1779	// updates the wave and global
1780	// returns a2 in cm
1781	Function sampleApertureDiam()
1782	ControlInfo/W=SASCALC popup0_1
1783
1784	// Print "In sampleApertureDiam()"
1785	//set the global
1786	NVAR a2=root:Packages:NIST:SAS:gSamAp
1787
1788	if(cmpstr(S_Value,"other") == 0) // "other" selected
1789	//enable the setvar, diameter in mm!
1790	SetVariable setvar0_3 disable=0
1791	// read its value (a global)
1792	NVAR a2other = root:Packages:NIST:SAS:gSamApOther
1793	a2=a2other/10 //a2 in cm
1794	else
1795	SetVariable setvar0_3 disable=1
1796	//1st item is 1/16", popup steps by 1/16"
1797	a2 = 2.54/16.0 * (V_Value) //convert to cm
1798	endif
1799	WAVE rw=root:Packages:NIST:SAS:realsRead
1800	rw[24] = a2*10 //sample aperture diameter in mm
1801
1802	return(a2)
1803	end
1804
1805	//1=huber 2=chamber
1806	Function tableposition()
1807	ControlInfo/W=SASCALC checkHuber
1808	if(V_Value)
1809	return(1) //huber
1810	else
1811	return(2) //chamber
1812	endif
1813	End
1814
1815	//compute SSD and update both the global and the wave
1816	//
1817	Function sourceToSampleDist()
1818
1819	NVAR NG=root:Packages:NIST:SAS:gNg
1820	NVAR S12 = root:Packages:NIST:SAS:S12
1821	NVAR L2Diff = root:Packages:NIST:SAS:L2Diff
1822	NVAR SSD = root:Packages:NIST:SAS:gSSD
1823
1824	SSD = 1632 - 155NG - s12(2-tableposition()) - L2Diff
1825
1826	WAVE rw=root:Packages:NIST:SAS:realsRead
1827	rw[25] = SSD/100 // in meters
1828	return(SSD)
1829	End
1830
1831
1832	// not part of SASCALC, but can be used to convert the SSD to number of guides
1833	//
1834	// SSD in meters
1835	Function numGuides(SSD)
1836	variable SSD
1837
1838	Variable Ng
1839	Ng = SSD*100 + 5 - 1632
1840	Ng /= -155
1841
1842	Ng = round(Ng)
1843	return(Ng)
1844	End
1845
1846
1847	//returns the offset value
1848	// slider and setVar are linked to the same global
1849	// updates the wave and changes the beamcenter (x,y) in the wave
1850	Function detectorOffset()
1851
1852	WAVE rw=root:Packages:NIST:SAS:RealsRead
1853	NVAR val = root:Packages:NIST:SAS:gOffset
1854	rw[19] = val // already in cm
1855	//move the beamcenter, make it an integer value for the MC simulation
1856	rw[16] = 64 + round(2*rw[19]) //approximate beam X is 64 w/no offset, 114 w/25 cm offset
1857	rw[17] = 64 //typical value
1858
1859	return(val)
1860	end
1861
1862	//returns the detector distance (slider and setVar are linked by the global)
1863	//
1864	Function chamberToDetectorDist()
1865
1866	NVAR val = root:Packages:NIST:SAS:gDetDist
1867	return(val)
1868	End
1869
1870	//sets the SDD (slider and setVar are linked by the global and is the detector position
1871	// relative to the chamber)
1872	// updates the wave
1873	Function sampleToDetectorDist()
1874
1875	NVAR detDist = root:Packages:NIST:SAS:gDetDist
1876	NVAR S12 = root:Packages:NIST:SAS:S12
1877	WAVE rw=root:Packages:NIST:SAS:RealsRead
1878	Variable SDD
1879
1880	SDD = detDist + s12*(2-tableposition())
1881	rw[18] = SDD/100 // convert to meters for header
1882	return(SDD)
1883	End
1884
1885	//direction = one of "vertical;horizontal;maximum;"
1886	// all of this is bypassed if the lenses are in
1887	//
1888	Function beamDiameter(direction)
1889	String direction
1890
1891	NVAR lens = root:Packages:NIST:SAS:gUsingLenses
1892	if(lens)
1893	return sourceApertureDiam()
1894	endif
1895
1896	Variable l1 = sourceToSampleDist()
1897	Variable l2 //= sampleAperToDetDist()
1898	Variable d1,d2,bh,bv,bm,umbra,a1,a2
1899
1900	NVAR L2diff = root:Packages:NIST:SAS:L2diff
1901	NVAR lambda = root:Packages:NIST:SAS:gLambda
1902	NVAR lambda_width = root:Packages:NIST:SAS:gDeltaLambda
1903	NVAR bs_factor = root:Packages:NIST:SAS:bs_factor
1904
1905	l2 = sampleToDetectorDist() + L2diff
1906	a1 = sourceApertureDiam()
1907	a2 = sampleApertureDiam()
1908
1909	d1 = a1*l2/l1
1910	d2 = a2*(l1+l2)/l1
1911	bh = d1+d2 //beam size in horizontal direction
1912	umbra = abs(d1-d2)
1913	//vertical spreading due to gravity
1914	bv = bh + 1.25e-8(l1+l2)l2lambdalambda*lambda_width
1915	bm = (bs_factorbh > bv) ? bs_factorbh : bv //use the larger of horiz*safety or vertical
1916
1917	strswitch(direction) // string switch
1918	case "vertical": // execute if case matches expression
1919	return(bv)
1920	break // exit from switch
1921	case "horizontal": // execute if case matches expression
1922	return(bh)
1923	break
1924	case "maximum": // execute if case matches expression
1925	return(bm)
1926	break
1927	default: // optional default expression executed
1928	return(bm) // when no case matches
1929	endswitch
1930	End
1931
1932	//on NG3 and NG7, allowable sizes are 1,2,3,4" diameter
1933	//will return values larger than 4.0*2.54 if a larger beam is needed
1934	//
1935	// - in an approximate way, account for lenses
1936	Function beamstopDiam()
1937
1938	NVAR yesLens = root:Packages:NIST:SAS:gUsingLenses
1939	Variable bm=0
1940	Variable bs=0.0
1941
1942	if(yesLens)
1943	//bm = sourceApertureDiam() //ideal result, not needed
1944	bs = 1 //force the diameter to 1"
1945	else
1946	bm = beamDiameter("maximum")
1947	do
1948	bs += 1
1949	while ( (bs*2.54 < bm) \|\| (bs > 30.0)) //30 = ridiculous limit to avoid inf loop
1950	endif
1951
1952	//update the wave
1953	WAVE rw=root:Packages:NIST:SAS:realsRead
1954	rw[21] = bs*25.4 //store the BS diameter in mm
1955
1956	return (bs*2.54) //return diameter in cm, not inches for txt
1957	End
1958
1959	//returns the projected diameter of the beamstop at the anode plane.
1960	// most noticeable at short SDD
1961	//if flag == 0 use conservative estimate = largest diameter (for SASCALC, default)
1962	//if flag != 0 use point aperture = average diameter (for resolution calculation)
1963	Function beamstopDiamProjection(flag)
1964	Variable flag
1965
1966	NVAR L2diff = root:Packages:NIST:SAS:L2diff
1967	Variable a2 = sampleApertureDiam()
1968	Variable bs = beamstopDiam()
1969	Variable l2, LB, BS_P
1970
1971	l2 = sampleToDetectorDist() + L2diff
1972	LB = 20.1 + 1.61*BS //distance in cm from beamstop to anode plane (empirical)
1973	if(flag==0)
1974	BS_P = bs + (bs+a2)*lb/(l2-lb) //diameter of shadow from parallax
1975	else
1976	BS_P = bs + bs*lb/(l2-lb) //diameter of shadow, point A2
1977	endif
1978	return (bs_p) //return projected diameter in cm
1979	End
1980
1981	// 19MAR07 - using correction from John for an estimate of the shadow of the beamstop
1982	// at the detection plane. This is a noticeable effect at short SDD, where the projected
1983	// diameter of the beamstop is much larger than the physical diameter.
1984	Function qMin()
1985
1986	Variable l2s = sampleToDetectorDist() //distance from sample to detector in cm
1987	// Variable bs = beamstopDiam() //beamstop diameter in cm
1988	Variable bs_p = beamstopDiamProjection(0) //projected beamstop diameter in cm
1989	NVAR lambda = root:Packages:NIST:SAS:gLambda
1990	NVAR d_det = root:Packages:NIST:SAS:d_det //cm
1991	NVAR a_pixel = root:Packages:NIST:SAS:a_pixel //cm
1992
1993	return( (pi/lambda)*(bs_p + d_det + a_pixel)/l2s ) //use bs_p rather than bs
1994	// return( (pi/lambda)*(bs + d_det + a_pixel)/l2s ) //use bs (incorrect)
1995	End
1996
1997	Function qMaxVert()
1998
1999	Variable theta
2000	Variable l2s = sampleToDetectorDist() //distance from sample to detector
2001	NVAR lambda = root:Packages:NIST:SAS:gLambda
2002	NVAR det_width = root:Packages:NIST:SAS:det_width
2003
2004	theta = atan( (det_width/2.0)/l2s )
2005
2006	return ( 4.0pi/lambda sin(theta/2.0) )
2007	end
2008
2009	Function qMaxCorner()
2010
2011	Variable l2s = sampleToDetectorDist() //distance from sample to detector
2012	Variable radial
2013	NVAR lambda = root:Packages:NIST:SAS:gLambda
2014	NVAR det_off = root:Packages:NIST:SAS:gOffset
2015	NVAR det_width = root:Packages:NIST:SAS:det_width
2016
2017	radial=sqrt( (0.5det_width)(0.5det_width)+(0.5det_width+det_off)(0.5det_width+det_off) )
2018
2019	return ( 4pi/lambdasin(0.5*atan(radial/l2s)) )
2020	End
2021
2022	Function qMaxHoriz()
2023
2024	Variable theta
2025	Variable l2s = sampleToDetectorDist() //distance from sample to detector
2026	NVAR lambda = root:Packages:NIST:SAS:gLambda
2027	NVAR det_off = root:Packages:NIST:SAS:gOffset
2028	NVAR det_width = root:Packages:NIST:SAS:det_width
2029
2030	theta = atan( ((det_width/2.0)+det_off)/l2s ) //from the instance variables
2031
2032	return ( 4.0pi/lambda sin(theta/2.0) )
2033	End
2034
2035	// calculate sigma for the resolution function at either limit of q-range
2036	Function deltaQ(atQ)
2037	Variable atQ
2038
2039	Variable k02,lp,l1,l2,sig_02,sigQ2,a1,a2
2040	NVAR l2Diff = root:Packages:NIST:SAS:L2diff
2041	NVAR lambda = root:Packages:NIST:SAS:gLambda
2042	NVAR lambda_width = root:Packages:NIST:SAS:gDeltaLambda
2043	NVAR d_det = root:Packages:NIST:SAS:d_det
2044	NVAR del_r = root:Packages:NIST:SAS:del_r
2045
2046
2047	l1 = sourceToSampleDist()
2048	l2 = sampleToDetectorDist() + L2diff
2049	a1 = sourceApertureDiam()
2050	a2 = sampleApertureDiam()
2051
2052	k02 = (6.2832/lambda)*(6.2832/lambda)
2053	lp = 1/(1/l1 + 1/l2)
2054
2055	sig_02 = (0.25a1/l1)(0.25*a1/l1)
2056	sig_02 += (0.25a2/lp)(0.25*a2/lp)
2057	sig_02 += (d_det/(2.355l2))(d_det/(2.355*l2))
2058	sig_02 += (del_r/l2)*(del_r/l2)/12
2059	sig_02 *= k02
2060
2061	sigQ2 = sig_02 + (atQlambda_width)(atQ*lambda_width)/6
2062
2063	return(100*sqrt(sigQ2)/atQ)
2064	End
2065
2066
2067	Function beamIntensity()
2068
2069	Variable alpha,f,t,t4,t5,t6,as,solid_angle,l1,d2_phi
2070	Variable a1,a2,retVal
2071	SetDataFolder root:Packages:NIST:SAS
2072	NVAR l_gap=l_gap,guide_width =guide_width,ng = gNg
2073	NVAR lambda_t=lambda_t,b=b,c=c
2074	NVAR lambda=gLambda,t1=t1,t2=t2,t3=t3,phi_0=phi_0
2075	NVAR lambda_width=gDeltaLambda
2076
2077	l1 = sourceToSampleDist()
2078	a1 = sourceApertureDiam()
2079	a2 = sampleApertureDiam()
2080
2081
2082	alpha = (a1+a2)/(2*l1) //angular divergence of beam
2083	f = l_gapalpha/(2guide_width)
2084	t4 = (1-f)*(1-f)
2085	t5 = exp(ng*ln(0.96)) // trans losses of guides in pre-sample flight
2086	t6 = 1 - lambda(b-(ng/8)(b-c)) //experimental correction factor
2087	t = t1t2t3t4t5*t6
2088
2089	as = pi/4a2a2 //area of sample in the beam
2090	d2_phi = phi_0/(2*pi)
2091	d2_phi = exp(4ln(lambda_t/lambda))
2092	d2_phi = exp(-1(lambda_t*lambda_t/lambda/lambda))
2093
2094	solid_angle = pi/4* (a1/l1)*(a1/l1)
2095
2096	retVal = as * d2_phi * lambda_width * solid_angle * t
2097	SetDataFolder root:
2098	return (retVal)
2099	end
2100
2101	Function figureOfMerit()
2102
2103	Variable bi = beamIntensity()
2104	NVAR lambda = root:Packages:NIST:SAS:gLambda
2105
2106	return (lambdalambdabi)
2107	End
2108
2109	//estimate the number of pixels in the beam, and enforce the maximum countrate per pixel (idmax)
2110	Function attenuatorTransmission()
2111
2112	Variable num_pixels,i_pix //i_pix = id in John's notation
2113	Variable bDiam = beamDiameter("horizontal") //!! note that prev calculations used bh (horizontal)
2114	Variable atten,a2
2115	SetDataFolder root:Packages:NIST:SAS
2116	NVAR a_pixel=a_pixel,idmax=idmax
2117
2118	a2 = sampleApertureDiam()
2119
2120	num_pixels = pi/4(0.5(a2+bDiam))(0.5(a2+bDiam))/a_pixel/a_pixel
2121	i_pix = ( beamIntensity() )/num_pixels
2122
2123	atten = (i_pix < idmax) ? 1.0 : idmax/i_pix
2124	SetDataFolder root:
2125	return(atten)
2126	End
2127
2128	Function attenuatorNumber()
2129
2130	Variable atten = attenuatorTransmission()
2131	Variable af,nf,numAtten
2132	SetDataFolder root:Packages:NIST:SAS
2133	NVAR lambda=gLambda
2134
2135	af = 0.498 + 0.0792lambda - 1.66e-3lambda*lambda
2136	nf = -ln(atten)/af //floating point
2137
2138	numAtten = trunc(nf) + 1 //in c, (int)nf
2139	//correct for larger step thickness at n > 6
2140	if(numAtten > 6)
2141	numAtten = 7 + trunc( (numAtten-6)/2 ) //in c, numAtten = 7 + (int)( (numAtten-6)/2 )
2142	endif
2143
2144	SetDatafolder root:
2145	return (numAtten)
2146	End

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

Download in other formats:

Original Format