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

Last change on this file since 546 was 546, checked in by srkline, 14 years ago

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