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

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

more changes to UCALC, display vs. countrate and add a display of experimental empty measurements along with the simulation. empty beam and empty call measurements are so similar - they have not been differentiated, especially since they are only for comparison.

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