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

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

Added a "signature" to the simulated 2D data set, setting a few pixels in the lower-left corner to known values. These points are behind the default mask.

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