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

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

Change (1):
In preparation for release, updated pragma IgorVersion?=6.1 in all procedures

Change (2):
As a side benefit of requiring 6.1, we can use the MultiThread? keyword to thread any model function we like. The speed benefit is only noticeable on functions that require at least one integration and at least 100 points (resolution smearing is NOT threaded, too many threadSafe issues, too little benefit). I have chosen to use the MultiThread? only on the XOP assignment. In the Igor code there are too many functions that are not explicitly declared threadsafe, making for a mess.

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