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

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

More changes to the event mode processing, updating the panel and adding some methods for removing misencoded time point

  1. Ready for some more testing to figure out what is good and what is bad.

Added some more values in to the SASCALC for the 10m instrument (still hidden)

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