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

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

modified the UCALC panel to have the control bar on the left, so that it wll better fit on-screen.

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