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

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

Added UCALC window to list of windows is package loader

Box Sum (from the marquee) asks for normalized (SAM) or RAW data

SASCALC uses real beamstop size rather than projected sizewhen doing simulation

Fake1DDataFolder() in MultiScatter_MonteCarlo ipf bug fixed where a bogus beamstop diameter (not 1"-4") would cause oddities in the resolution wave, and an (obviously) incorrect smearing of the simulation.

Fixed range 7 defaults in UCALC panel.

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