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

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

A bunch of random changes:
-SASCALC - added the 2.5 cm aperture NG3/7g/polarizer back in, but made the 5 cm the default
-Some fiddling with 2D functionality to enable 2D resolution smearing. Still a work in progress
-Added two new model functions, and added them to the list.
-Changes to the wrapper for the new release to generate kw=val strings as needed for each function
and its coefficients and parameters. This behavior has been changed in the new release, so this
fix should keep old analysis experiments compatible with the new version.

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