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

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

Added a procedure file to allow "scripting" of the MonteCarlo? simulation: MC_SimulationScripting.ipf. See the beginning of the file for instructions, and then the examples for numbered instructions of what needs to be edited to write your own script. Not casual-user-friendly, but makes it easire for more advanced users to simulate an entire experiment.

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