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

Last change on this file since 779 was 764, checked in by srkline, 12 years ago
  • fixed a few bugs with GenCurveFit?, and how the reports are generated


  • the DEFAULT.MASK file is automatically loaded when the path is set, if it can be found. this only happens from the main (yellow) panel. Nothing happens if that exact file is not found.


  • a "Sector_PlusMinus" averaging option is added. This defines the LEFT sector as being "negative" q-values. Everything else behaves as a normal sector average. This is from Lionel, a very old ticket #31


  • if sectors or annular regions are drawn on RAW data files, the drawn lines are re-drawn correctly as the data is scrolled using the < and > buttons.


  • a super secret option for a "histogram pair" has been added. May be useful for alignment, may ditch if Jeff and Cedric don't like it. To do this, put cursor A on the 2D image at the center of where you want the vertical and horizontal swath to be. +-5 pixels is hard-wired. draw any marquee(size, location is ignored) and select SANS Histogram, and you get the pair. If cursor A is not on the graph, you get the normal histogram as defined by the marquee.


  • arrow buttons on RAW 2D data display now search +- 10 data files for "next", in case there are missing file numbers.


  • Incorporated Lionel's changes to ILL_* files for his generation of a "DIV" file for D22


  • Added the offset traces checkbox back to the SASCALC panel. previously it was hidden on the simulation panels.


  • loosened the tolerance for SDD matching onn the MRED panel to 0.1 m



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