Context Navigation

SASCALC.ipf @ 598

Last change on this file since 598 was 598, checked in by srkline, 13 years ago
Procedures for mostly my use that enable me to (crudely) script the MC simulations so that they can be done in a big batch overnight.
File size: 62.5 KB

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

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

Download in other formats:

Original Format