source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/SASCALC.ipf @ 868

Last change on this file since 868 was 868, checked in by srkline, 10 years ago

Modified SASCALC to be *almost* rady for the 10m SANS instrument. the radio button is currently hidden, but can be easily turned back on. Most all of the functionality works, but there needs to be a lot of instrument-specific dimensions entered, and instrument-specific functionality. But it should be a relatively quick fix to enter in the necessary numbers.

Made a few changes to the NCNR_Utils to prepare for the 10m SANS. Here, mostly initialization constants, and the attenuator table. Everything else that the reduction needs is required to be in the raw data header.

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