source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/USANS/U_CALC.ipf @ 550

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

sample aperture popup now properly updates the displayed empty cell data

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