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

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

Added a calculation of the total fraction of multiple coherent scattering. This value is now displayed on the 1D sim panel and on the UCALC panel. If the total fraction of multiple coherent scattering (relative to the fraction coherently scattered) is greater than 0.10, then the field turns red as a warning.

Changed the numbering and suffix scheme for saving simulated VAX raw data files.

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