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

Last change on this file since 597 was 597, checked in by srkline, 13 years ago

Updated SASCALC with newly measured flux data for both NG3 and NG7.

-At 6A on NG7, now SASCALC is within +/- 5% of the measured value.
-At 6A on NG3, now SASCALC is +15% compared to the measured value.

The average detector efficiency is now included in the simulation (e=0.75) since the cts/s from SASCALC is flux on sample, not as would be measured at the detector.

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