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

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

Fixed logic error in writing the XY box count values to the empty beam header that could result in INF transmission.

Updated critera for "RAW" data files to include "SIM" as well.

Added NCNR_Utils to the includes list for Analysis.

File size: 65.2 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                                // this is where the number of cells comes in - the calculation of the error bars
1003                                // sigma[i] = SUM(sigma[ij]^2) / nCells^2
1004                                // and since in the simulation, SUM(sigma[ij]^2) = nCells*sigma[ij]^2 = nCells*Inten
1005                                // then...
1006                                sigave = sqrt(aveint/nCells)            // corrected based on John's memo, from 8/9/99
1007                               
1008                                // add in random error in aveint based on the sigave
1009                                if(addNoise)
1010                                        aveint += gnoise(sigave)
1011                                endif
1012
1013                                // convert to absolute scale
1014                                if(doABS)
1015                                        Variable kappa = thick*(pixSize/sdd)^2*trans*iMon*ctTime
1016                                        aveint /= kappa
1017                                        sigave /= kappa
1018                                endif
1019                               
1020                        else
1021                                //no function plotted, no simulation can be done
1022                                DoAlert 0,"No function is selected or plotted, so no simulation is done. The default Debye function is used."
1023               
1024                                aveint = S_Debye(1000,100,0.0,qval)
1025                                aveint *= fSubS         // multiply either estimate by beamstop shadowing
1026                                sigave = 0              //reset for model calculation
1027                        endif
1028                       
1029                endif // end 1D simulation
1030        else
1031                //no simulation
1032               
1033                aveint = S_Debye(1000,100,0.0,qval)
1034                aveint *= fSubS         // multiply either estimate by beamstop shadowing
1035                sigave = 0              //reset for model calculation
1036               
1037                //end no simulation     
1038        endif
1039       
1040
1041        //display the configuration text in a separate notebook
1042        DisplayConfigurationText()
1043       
1044        return(0)
1045End
1046
1047//freezes the current configuration
1048// -1- duplicates the trace on the graph
1049// -2- copies the configuration text to a second notebook window for printing
1050Function FreezeButtonProc(ctrlName) : ButtonControl
1051        String ctrlName
1052
1053        String str=""
1054        NVAR ct=root:Packages:NIST:SAS:gFreezeCount
1055       
1056       
1057        SetDataFolder root:Packages:NIST:SAS
1058       
1059        Duplicate/O aveint,$("aveint_"+num2str(ct))
1060        Duplicate/O qval,$("qval_"+num2str(ct))
1061        Duplicate/O sigave,$("sigave_"+num2str(ct))
1062        Appendtograph $("aveint_"+num2str(ct)) vs $("qval_"+num2str(ct))
1063        ModifyGraph mode($("aveint_"+num2str(ct)))=3
1064        ModifyGraph marker($("aveint_"+num2str(ct)))=19
1065        ModifyGraph msize($("aveint_"+num2str(ct)))=2
1066        ErrorBars/T=0 $("aveint_"+num2str(ct)) Y,wave=($("sigave_"+num2str(ct)),$("sigave_"+num2str(ct)))
1067       
1068        switch(mod(ct,10))      // 10 different colors - black is the unfrozen color
1069                case 0:
1070                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,16385,16385)
1071                        break
1072                case 1:
1073                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(2,39321,1)
1074                        break
1075                case 2:
1076                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(0,0,65535)
1077                        break
1078                case 3:
1079                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(39321,1,31457)
1080                        break
1081                case 4:
1082                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(48059,48059,48059)
1083                        break
1084                case 5:
1085                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,32768,32768)
1086                        break
1087                case 6:
1088                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(0,65535,0)
1089                        break
1090                case 7:
1091                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(16385,65535,65535)
1092                        break
1093                case 8:
1094                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,32768,58981)
1095                        break
1096                case 9:
1097                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(36873,14755,58982)
1098                        break
1099        endswitch
1100       
1101        NVAR doTraceOffset = root:Packages:NIST:SAS:gDoTraceOffset
1102        NVAR offset = root:Packages:NIST:SAS:gModelOffsetFactor
1103        if(doTraceOffset)
1104                offset = 2^ct
1105                //multiply by current offset (>=1)
1106                Wave inten = $("aveint_"+num2str(ct))
1107                inten *= offset
1108                //      Print "new offset = ",offset
1109        endif
1110       
1111        ct +=1
1112        SetDataFolder root:
1113       
1114        // create or append the configuration to a new notebook
1115        if(WinType("Saved_Configurations")==0)
1116                NewNotebook/F=1/N=Saved_Configurations /W=(480,400,880,725)
1117                DoWindow/F SASCALC              //return focus to SASCALC window
1118        endif
1119        //append the text
1120        sprintf str,"\rConfiguration #%d\r",ct-1
1121        Notebook Saved_Configurations showRuler=0,defaultTab=20,selection={endOfFile, endOfFile}
1122        Notebook Saved_Configurations font="Monaco",fSize=10,fstyle=1,text=str          //bold
1123        Notebook Saved_Configurations font="Monaco",fSize=10,fstyle=0,text=SetConfigurationText()
1124        return(0)
1125End
1126
1127//clears the frozen traces on the graph, asks if you want to clear the saved text
1128//
1129Function S_ClearButtonProc(ctrlName) : ButtonControl
1130        String ctrlName
1131
1132        NVAR ct=root:Packages:NIST:SAS:gFreezeCount
1133        Variable ii
1134        Setdatafolder root:Packages:NIST:SAS
1135        for(ii=ct-1;ii>=1;ii-=1)
1136        //remove all traces, replace aveint
1137        // kill all waves _ct
1138                RemoveFromGraph $("aveint_"+num2str(ii))
1139                Killwaves/Z $("aveint_"+num2str(ii)),$("qval_"+num2str(ii))
1140        endfor
1141        ct=1
1142        setdatafolder root:
1143       
1144        DoAlert 1,"Do you also want to clear the \"Saved Configurations\" text?"
1145        if(V_flag == 1)         // yes
1146                DoWindow/K Saved_Configurations
1147        endif
1148       
1149        //reset offset value
1150        NVAR offset = root:Packages:NIST:SAS:gModelOffsetFactor
1151        offset = 1
1152        ReCalculateInten(1)
1153        return(0)
1154End
1155
1156
1157///////////////////////////////////////////////////////////
1158// 19MAR07 uses correction for beamstop diameter projection to get shadow factor correct
1159//
1160Function S_CircularAverageTo1D(type)
1161        String type
1162       
1163        Variable isCircular = 1
1164       
1165        //type is the data type to do the averaging on, and will be set as the current folder
1166        //get the current displayed data (so the correct folder is used)
1167        String destPath = "root:Packages:NIST:"+type
1168       
1169        //
1170        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist,dr,ddr
1171        Variable lambda,trans
1172        WAVE reals = $(destPath + ":RealsRead")
1173//      WAVE/T textread = $(destPath + ":TextRead")
1174//      String fileStr = textread[3]
1175       
1176        // center of detector, for non-linear corrections
1177        Variable pixelsX=128,pixelsY=128
1178       
1179        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
1180        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
1181       
1182        // beam center, in pixels
1183        x0 = reals[16]
1184        y0 = reals[17]
1185        //detector calibration constants
1186        sx = reals[10]          //mm/pixel (x)
1187        sx3 = reals[11]         //nonlinear coeff
1188        sy = reals[13]          //mm/pixel (y)
1189        sy3 = reals[14]         //nonlinear coeff
1190       
1191        dtsize = 10*reals[20]           //det size in mm
1192        dtdist = 1000*reals[18]         // det distance in mm
1193       
1194        NVAR binWidth=root:Packages:NIST:SAS:gBinWidth
1195//      Variable binWidth = 1
1196       
1197        dr = binWidth           // ***********annulus width set by user, default is one***********
1198        ddr = dr*sx             //step size, in mm (this value should be passed to the resolution calculation, not dr 18NOV03)
1199               
1200        Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,dxbm,dybm,ii
1201        Variable phi_rad,dphi_rad,phi_x,phi_y
1202        Variable forward,mirror
1203       
1204        String side = "both"
1205//      side = StringByKey("SIDE",keyListStr,"=",";")
1206//      Print "side = ",side
1207       
1208        /// data wave is data in the current folder which was set at the top of the function
1209        WAVE data=$(destPath + ":data")
1210
1211// fake mask that uses all of the detector
1212        Make/D/O/N=(pixelsX,pixelsY) $(destPath + ":mask")
1213        Wave mask = $(destPath + ":mask")
1214        mask = 0
1215        //two pixels all around
1216        mask[0,1][] = 1
1217        mask[126,127][] = 1
1218        mask[][0,1] = 1
1219        mask[][126,127] = 1
1220        //
1221        //pixels within rcentr of beam center are broken into 9 parts (units of mm)
1222        rcentr = 100            //original
1223//      rcentr = 0
1224        // values for error if unable to estimate value
1225        //large_num = 1e10
1226        large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1
1227        small_num = 1e-10
1228       
1229        // output wave are expected to exist (?) initialized to zero, what length?
1230        // 200 points on VAX --- use 300 here, or more if SAXS data is used with 1024x1024 detector (1000 pts seems good)
1231        Variable defWavePts=500
1232        Make/O/D/N=(defWavePts) $(destPath + ":qval"),$(destPath + ":aveint")
1233        Make/O/D/N=(defWavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave")
1234        Make/O/D/N=(defWavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar")
1235
1236        WAVE qval = $(destPath + ":qval")
1237        WAVE aveint = $(destPath + ":aveint")
1238        WAVE ncells = $(destPath + ":ncells")
1239        WAVE dsq = $(destPath + ":dsq")
1240        WAVE sigave = $(destPath + ":sigave")
1241        WAVE qbar = $(destPath + ":QBar")
1242        WAVE sigmaq = $(destPath + ":SigmaQ")
1243        WAVE fsubs = $(destPath + ":fSubS")
1244       
1245        qval = 0
1246        aveint = 0
1247        ncells = 0
1248        dsq = 0
1249        sigave = 0
1250        qbar = 0
1251        sigmaq = 0
1252        fsubs = 0
1253
1254        dtdis2 = dtdist^2
1255        nq = 1
1256        xoffst=0
1257        //distance of beam center from detector center
1258        dxbm = S_FX(x0,sx3,xcenter,sx)
1259        dybm = S_FY(y0,sy3,ycenter,sy)
1260               
1261        //BEGIN AVERAGE **********
1262        Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1
1263        Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p
1264       
1265        // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
1266        // loop index corresponds to FORTRAN (old code)
1267        // and the IGOR array indices must be adjusted (-1) to the correct address
1268        ii=1
1269        do
1270                xi = ii
1271                dxi = S_FX(xi,sx3,xcenter,sx)
1272                dx = dxi-dxbm           //dx and dy are in mm
1273               
1274                jj = 1
1275                do
1276                        data_pixel = data[ii-1][jj-1]           //assign to local variable
1277                        yj = jj
1278                        dyj = S_FY(yj,sy3,ycenter,sy)
1279                        dy = dyj - dybm
1280                        if(!(mask[ii-1][jj-1]))                 //masked pixels = 1, skip if masked (this way works...)
1281                                dr2 = (dx^2 + dy^2)^(0.5)               //distance from beam center NOTE dr2 used here - dr used above
1282                                if(dr2>rcentr)          //keep pixel whole
1283                                        nd = 1
1284                                        fd = 1
1285                                else                            //break pixel into 9 equal parts
1286                                        nd = 3
1287                                        fd = 2
1288                                endif
1289                                nd2 = nd^2
1290                                ll = 1          //"el-el" loop index
1291                                do
1292                                        dxx = dx + (ll - fd)*sx/3
1293                                        kk = 1
1294                                        do
1295                                                dyy = dy + (kk - fd)*sy/3
1296                                                if(isCircular)
1297                                                        //circular average, use all pixels
1298                                                        //(increment)
1299                                                        nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1300                                                else
1301                                                        //a sector average - determine azimuthal angle
1302                                                        dphi_p = S_dphi_pixel(dxx,dyy,phi_x,phi_y)
1303                                                        if(dphi_p < dphi_rad)
1304                                                                forward = 1                     //within forward sector
1305                                                        else
1306                                                                forward = 0
1307                                                        Endif
1308                                                        if((Pi - dphi_p) < dphi_rad)
1309                                                                mirror = 1              //within mirror sector
1310                                                        else
1311                                                                mirror = 0
1312                                                        Endif
1313                                                        //check if pixel lies within allowed sector(s)
1314                                                        if(cmpstr(side,"both")==0)              //both sectors
1315                                                                if ( mirror || forward)
1316                                                                        //increment
1317                                                                        nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1318                                                                Endif
1319                                                        else
1320                                                                if(cmpstr(side,"right")==0)             //forward sector only
1321                                                                        if(forward)
1322                                                                                //increment
1323                                                                                nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1324                                                                        Endif
1325                                                                else                    //mirror sector only
1326                                                                        if(mirror)
1327                                                                                //increment
1328                                                                                nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1329                                                                        Endif
1330                                                                Endif
1331                                                        Endif           //allowable sectors
1332                                                Endif   //circular or sector check
1333                                                kk+=1
1334                                        while(kk<=nd)
1335                                        ll += 1
1336                                while(ll<=nd)
1337                        Endif           //masked pixel check
1338                        jj += 1
1339                while (jj<=pixelsY)
1340                ii += 1
1341        while(ii<=pixelsX)              //end of the averaging
1342               
1343        //compute q-values and errors
1344        Variable ntotal,rr,theta,avesq,aveisq,var
1345       
1346        lambda = reals[26]
1347        ntotal = 0
1348        kk = 1
1349        do
1350                rr = (2*kk-1)*ddr/2
1351                theta = 0.5*atan(rr/dtdist)
1352                qval[kk-1] = (4*Pi/lambda)*sin(theta)
1353                if(ncells[kk-1] == 0)
1354                        //no pixels in annuli, data unknown
1355                        aveint[kk-1] = 0
1356                        sigave[kk-1] = large_num
1357                else
1358                        if(ncells[kk-1] <= 1)
1359                                //need more than one pixel to determine error
1360                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
1361                                sigave[kk-1] = large_num
1362                        else
1363                                //assume that the intensity in each pixel in annuli is normally
1364                                // distributed about mean...
1365                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
1366                                avesq = aveint[kk-1]^2
1367                                aveisq = dsq[kk-1]/ncells[kk-1]
1368                                var = aveisq-avesq
1369                                if(var<=0)
1370                                        sigave[kk-1] = small_num
1371                                else
1372                                        sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1))
1373                                endif
1374                        endif
1375                endif
1376                ntotal += ncells[kk-1]
1377                kk+=1
1378        while(kk<=nq)
1379       
1380        //Print "NQ = ",nq
1381        // data waves were defined as 300 points (=defWavePts), but now have less than that (nq) points
1382        // use DeletePoints to remove junk from end of waves
1383        //WaveStats would be a more foolproof implementation, to get the # points in the wave
1384        Variable startElement,numElements
1385        startElement = nq
1386        numElements = defWavePts - startElement
1387        DeletePoints startElement,numElements, qval,aveint,ncells,dsq,sigave
1388       
1389        //////////////end of VAX sector_ave()
1390               
1391
1392
1393// ***************************************************************
1394//
1395// Do the extra 3 columns of resolution calculations starting here.
1396//
1397// ***************************************************************
1398
1399        Variable L2 = reals[18]
1400        Variable BS = reals[21]         //this the diameter is stored in mm
1401//  SRK - why was I using the projected diameter of the beam stop?? I added a step at the beginning of every recalculation
1402// of the intensity to get the right beamstop diameter into RealsRead...
1403//      Variable BS = beamstopDiamProjection(1) * 10            //calculated projection in cm *10 = mm
1404        Variable S1 = reals[23]
1405        Variable S2 = reals[24]
1406        Variable L1 = reals[25]
1407        lambda = reals[26]
1408        Variable lambdaWidth = reals[27]
1409
1410        Variable DDet, apOff
1411        //typical value for NG3 and NG7 - distance between sample aperture and sample in (cm)
1412        apOff=5.0
1413        // hard wire value for Ordela detectors
1414        DDet = 0.5              // resolution in cm
1415        //      String detStr=textRead[9]
1416        //      DDet = DetectorPixelResolution(fileStr,detStr)          //needs detector type and beamline
1417
1418        //Go from 0 to nq doing the calc for all three values at
1419        //every Q value
1420
1421        ii=0
1422        Variable ret1,ret2,ret3
1423        do
1424                S_getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,ret1,ret2,ret3)
1425                sigmaq[ii] = ret1       //res_wave[0]
1426                qbar[ii] = ret2         //res_wave[1]
1427                fsubs[ii] = ret3                //res_wave[2]
1428                ii+=1
1429        while(ii<nq)
1430        DeletePoints startElement,numElements, sigmaq, qbar, fsubs
1431
1432        fsubs += 1e-8           //keep the values from being too small
1433
1434// End of resolution calculations
1435// ***************************************************************
1436       
1437        //get rid of the default mask, if one was created (it is in the current folder)
1438        //don't just kill "mask" since it might be pointing to the one in the MSK folder
1439        Killwaves/Z $(destPath+":mask")
1440       
1441        //return to root folder (redundant)
1442        SetDataFolder root:
1443       
1444        Return 0
1445End
1446
1447//returns nq, new number of q-values
1448//arrays aveint,dsq,ncells are also changed by this function
1449//
1450Function S_IncrementPixel(dataPixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1451        Variable dataPixel,ddr,dxx,dyy
1452        Wave aveint,dsq,ncells
1453        Variable nq,nd2
1454       
1455        Variable ir
1456       
1457        ir = trunc(sqrt(dxx*dxx+dyy*dyy)/ddr)+1
1458        if (ir>nq)
1459                nq = ir         //resets maximum number of q-values
1460        endif
1461        aveint[ir-1] += dataPixel/nd2           //ir-1 must be used, since ir is physical
1462        dsq[ir-1] += dataPixel*dataPixel/nd2
1463        ncells[ir-1] += 1/nd2
1464       
1465        Return nq
1466End
1467
1468//function determines azimuthal angle dphi that a vector connecting
1469//center of detector to pixel makes with respect to vector
1470//at chosen azimuthal angle phi -> [cos(phi),sin(phi)] = [phi_x,phi_y]
1471//dphi is always positive, varying from 0 to Pi
1472//
1473Function S_dphi_pixel(dxx,dyy,phi_x,phi_y)
1474        Variable dxx,dyy,phi_x,phi_y
1475       
1476        Variable val,rr,dot_prod
1477       
1478        rr = sqrt(dxx^2 + dyy^2)
1479        dot_prod = (dxx*phi_x + dyy*phi_y)/rr
1480        //? correct for roundoff error? - is this necessary in IGOR, w/ double precision?
1481        if(dot_prod > 1)
1482                dot_prod =1
1483        Endif
1484        if(dot_prod < -1)
1485                dot_prod = -1
1486        Endif
1487       
1488        val = acos(dot_prod)
1489       
1490        return val
1491
1492End
1493
1494//calculates the x distance from the center of the detector, w/nonlinear corrections
1495//
1496Function S_FX(xx,sx3,xcenter,sx)               
1497        Variable xx,sx3,xcenter,sx
1498       
1499        Variable retval
1500       
1501        retval = sx3*tan((xx-xcenter)*sx/sx3)
1502        Return retval
1503End
1504
1505//calculates the y distance from the center of the detector, w/nonlinear corrections
1506//
1507Function S_FY(yy,sy3,ycenter,sy)               
1508        Variable yy,sy3,ycenter,sy
1509       
1510        Variable retval
1511       
1512        retval = sy3*tan((yy-ycenter)*sy/sy3)
1513        Return retval
1514End
1515
1516//**********************
1517// Resolution calculation - used by the averaging routines
1518// to calculate the resolution function at each q-value
1519// - the return value is not used
1520//
1521// equivalent to John's routine on the VAX Q_SIGMA_AVE.FOR
1522// Incorporates eqn. 3-15 from J. Appl. Cryst. (1995) v. 28 p105-114
1523//
1524Function/S S_getResolution(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,SigmaQ,QBar,fSubS)
1525        Variable inQ, lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r
1526        Variable &fSubS, &QBar, &SigmaQ         //these are the output quantities at the input Q value
1527       
1528        //lots of calculation variables
1529        Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g
1530        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
1531
1532        //Constants
1533        //Variable del_r = .1
1534        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
1535        Variable g = 981.0                              //gravity acceleration [cm/s^2]
1536
1537        String results
1538        results ="Failure"
1539       
1540        NVAR usingLenses = root:Packages:NIST:SAS:gUsingLenses
1541
1542        //rename for working variables,  these must be gotten from global
1543        //variables
1544
1545//      Variable wLam, wLW, wL1, wL2, wS1, wS2
1546//      Variable wBS, wDDet, wApOff
1547//      wLam = lambda
1548//      wLW = lambdaWidth
1549//      wDDet = DDet
1550//      wApOff = apOff
1551        S1 *= 0.5*0.1                   //convert to radius and [cm]
1552        S2 *= 0.5*0.1
1553
1554        L1 *= 100.0                     // [cm]
1555        L1 -= apOff                             //correct the distance
1556
1557        L2 *= 100.0
1558        L2 += apOff
1559
1560        BS *= 0.5*0.1                   //convert to radius and [cm]
1561        del_r *= 0.1                            //width of annulus, convert mm to [cm]
1562       
1563        //Start resolution calculation
1564        a2 = S1*L2/L1 + S2*(L1+L2)/L1
1565        q_small = 2.0*Pi*(BS-a2)*(1.0-lambdaWidth)/(lambda*L2)
1566        lp = 1.0/( 1.0/L1 + 1.0/L2)
1567
1568        v_lambda = lambdaWidth^2/6.0
1569       
1570        if(usingLenses==1)                      //SRK 2007
1571                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term
1572        else
1573                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form
1574        endif
1575       
1576        v_d = (DDet/2.3548)^2 + del_r^2/12.0
1577        vz = vz_1 / lambda
1578        yg = 0.5*g*L2*(L1+L2)/vz^2
1579        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007
1580
1581        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
1582        delta = 0.5*(BS - r0)^2/v_d
1583
1584        if (r0 < BS)
1585                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
1586        else
1587                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
1588        endif
1589
1590        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
1591        if (fSubS <= 0.0)
1592                fSubS = 1.e-10
1593        endif
1594        fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
1595        fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
1596
1597        rmd = fr*r0
1598        v_r1 = v_b + fv*v_d +v_g
1599
1600        rm = rmd + 0.5*v_r1/rmd
1601        v_r = v_r1 - 0.5*(v_r1/rmd)^2
1602        if (v_r < 0.0)
1603                v_r = 0.0
1604        endif
1605        QBar = (4.0*Pi/lambda)*sin(0.5*atan(rm/L2))
1606        SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda)
1607       
1608        results = "success"
1609        Return results
1610End
1611
1612Function S_Debye(scale,rg,bkg,x)
1613        Variable scale,rg,bkg
1614        Variable x
1615       
1616        // variables are:
1617        //[0] scale factor
1618        //[1] radius of gyration [A]
1619        //[2] background        [cm-1]
1620       
1621        // calculates (scale*debye)+bkg
1622        Variable Pq,qr2
1623       
1624        qr2=(x*rg)^2
1625        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
1626       
1627        //scale
1628        Pq *= scale
1629        // then add in the background
1630        return (Pq+bkg)
1631End
1632
1633
1634Function S_SphereForm(scale,radius,delrho,bkg,x)                               
1635        Variable scale,radius,delrho,bkg
1636        Variable x
1637       
1638        // variables are:                                                       
1639        //[0] scale
1640        //[1] radius (A)
1641        //[2] delrho (A-2)
1642        //[3] background (cm-1)
1643       
1644//      Variable scale,radius,delrho,bkg                               
1645//      scale = w[0]
1646//      radius = w[1]
1647//      delrho = w[2]
1648//      bkg = w[3]
1649       
1650       
1651        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
1652        // and is rescaled to give [=] cm^-1
1653       
1654        Variable bes,f,vol,f2
1655        //
1656        //handle q==0 separately
1657        If(x==0)
1658                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
1659                return(f)
1660        Endif
1661       
1662        bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
1663        vol = 4*pi/3*radius^3
1664        f = vol*bes*delrho              // [=] A
1665        // normalize to single particle volume, convert to 1/cm
1666        f2 = f * f / vol * 1.0e8                // [=] 1/cm
1667       
1668        return (scale*f2+bkg)   // Scale, then add in the background
1669       
1670End
1671
1672Function/S SetConfigurationText()
1673
1674        String str="",temp
1675
1676        SetDataFolder root:Packages:NIST:SAS
1677       
1678        NVAR numberOfGuides=gNg
1679        NVAR gTable=gTable              //2=chamber, 1=table
1680        NVAR wavelength=gLambda
1681        NVAR lambdaWidth=gDeltaLambda
1682        NVAR instrument = instrument
1683        NVAR L2diff = L2diff
1684        NVAR lens = root:Packages:NIST:SAS:gUsingLenses
1685        SVAR aStr = root:myGlobals:gAngstStr
1686       
1687        sprintf temp,"Source Aperture Diameter =\t\t%6.2f cm\r",sourceApertureDiam()
1688        str += temp
1689        sprintf temp,"Source to Sample =\t\t\t\t%6.0f cm\r",sourceToSampleDist()
1690        str += temp
1691        sprintf temp,"Sample Aperture to Detector =\t%6.0f cm\r",sampleToDetectorDist() + L2diff
1692        str += temp
1693        sprintf temp,"Beam diameter =\t\t\t\t\t%6.2f cm\r",beamDiameter("maximum")
1694        str += temp
1695        sprintf temp,"Beamstop diameter =\t\t\t\t%6.2f inches\r",beamstopDiam()/2.54
1696        str += temp
1697        sprintf temp,"Minimum Q-value =\t\t\t\t%6.4f 1/%s (sigQ/Q = %3.1f %%)\r",qMin(),aStr,deltaQ(qMin())
1698        str += temp
1699        sprintf temp,"Maximum Horizontal Q-value =\t%6.4f 1/%s\r",qMaxHoriz(),aStr
1700        str += temp
1701        sprintf temp,"Maximum Vertical Q-value =\t\t%6.4f 1/%s\r",qMaxVert(),aStr
1702        str += temp
1703        sprintf temp,"Maximum Q-value =\t\t\t\t%6.4f 1/%s (sigQ/Q = %3.1f %%)\r",qMaxCorner(),aStr,deltaQ(qMaxCorner())
1704        str += temp
1705        sprintf temp,"Beam Intensity =\t\t\t\t%.0f counts/s\r",beamIntensity()
1706        str += temp
1707        sprintf temp,"Figure of Merit =\t\t\t\t%3.3g %s^2/s\r",figureOfMerit(),aStr
1708        str += temp
1709        sprintf temp,"Attenuator transmission =\t\t%3.3g = Atten # %d\r",attenuatorTransmission(),attenuatorNumber()
1710        str += temp
1711//     
1712//      // add text of the user-edited values
1713//      //
1714        sprintf temp,"***************** NG %d *****************\r",instrument
1715        str += temp
1716        sprintf temp,"Sample Aperture Diameter =\t\t\t\t%.2f cm\r",sampleApertureDiam()
1717        str += temp
1718        sprintf temp,"Number of Guides =\t\t\t\t\t\t%d \r", numberOfGuides
1719        str += temp
1720        sprintf temp,"Sample Chamber to Detector =\t\t\t%.1f cm\r", chamberToDetectorDist()
1721        str += temp
1722        if(gTable==1)
1723                sprintf temp,"Sample Position is \t\t\t\t\t\tHuber\r"
1724        else
1725                sprintf temp,"Sample Position is \t\t\t\t\t\tChamber\r"
1726        endif
1727        str += temp
1728        sprintf temp,"Detector Offset =\t\t\t\t\t\t%.1f cm\r", detectorOffset()
1729        str += temp
1730        sprintf temp,"Neutron Wavelength =\t\t\t\t\t%.2f %s\r", wavelength,aStr
1731        str += temp
1732        sprintf temp,"Wavelength Spread, FWHM =\t\t\t\t%.3f\r", lambdaWidth
1733        str += temp
1734        sprintf temp,"Sample Aperture to Sample Position =\t%.2f cm\r", L2Diff
1735        str += temp
1736        if(lens==1)
1737                sprintf temp,"Lenses are IN\r"
1738        else
1739                sprintf temp,"Lenses are OUT\r"
1740        endif
1741        str += temp
1742       
1743   setDataFolder root:
1744   return str                   
1745End
1746
1747Function DisplayConfigurationText()
1748
1749        if(WinType("Trial_Configuration")==0)
1750                NewNotebook/F=0/K=1/N=Trial_Configuration /W=(480,44,880,369)
1751        endif
1752        //replace the text
1753        Notebook Trial_Configuration selection={startOfFile, endOfFile}
1754        Notebook Trial_Configuration font="Monaco",fSize=10,text=SetConfigurationText()
1755        return(0)
1756end
1757
1758//parses the control for A1 diam
1759// updates the wave
1760Function sourceApertureDiam()
1761        ControlInfo/W=SASCALC popup0
1762        Variable diam
1763        sscanf S_Value, "%g cm", diam
1764        WAVE rw=root:Packages:NIST:SAS:realsRead
1765        rw[23] = diam*10                        //sample aperture diameter in mm
1766        return(diam)
1767end
1768
1769// change the sample aperture to a non-standard value
1770Function SampleApOtherSetVarProc(ctrlName,varNum,varStr,varName) : SetVariableControl
1771        String ctrlName
1772        Variable varNum
1773        String varStr
1774        String varName
1775               
1776        WAVE rw=root:Packages:NIST:SAS:realsRead
1777        rw[24] = varNum                 //sample aperture diameter in mm
1778        ReCalculateInten(1)
1779        return(0)
1780End
1781
1782//parses the control for A2 diam
1783// updates the wave and global
1784// returns a2 in cm
1785Function sampleApertureDiam()
1786        ControlInfo/W=SASCALC popup0_1
1787       
1788//      Print "In sampleApertureDiam()"
1789        //set the global
1790        NVAR a2=root:Packages:NIST:SAS:gSamAp
1791       
1792        if(cmpstr(S_Value,"other") == 0)                // "other" selected
1793                //enable the setvar, diameter in mm!
1794                SetVariable setvar0_3 disable=0
1795                // read its value (a global)
1796                NVAR a2other = root:Packages:NIST:SAS:gSamApOther
1797                a2=a2other/10                           //a2 in cm
1798        else
1799                SetVariable setvar0_3 disable=1
1800                //1st item is 1/16", popup steps by 1/16"
1801                a2 = 2.54/16.0 * (V_Value)                      //convert to cm         
1802        endif
1803        WAVE rw=root:Packages:NIST:SAS:realsRead
1804        rw[24] = a2*10                  //sample aperture diameter in mm
1805       
1806        return(a2)
1807end
1808
1809//1=huber 2=chamber
1810Function tableposition()
1811        ControlInfo/W=SASCALC checkHuber
1812        if(V_Value)
1813                return(1)               //huber
1814        else
1815                return(2)               //chamber
1816        endif
1817End
1818
1819//compute SSD and update both the global and the wave
1820//
1821Function sourceToSampleDist()
1822
1823        NVAR NG=root:Packages:NIST:SAS:gNg
1824        NVAR S12 = root:Packages:NIST:SAS:S12
1825        NVAR L2Diff = root:Packages:NIST:SAS:L2Diff
1826        NVAR SSD = root:Packages:NIST:SAS:gSSD
1827       
1828        SSD = 1632 - 155*NG - s12*(2-tableposition()) - L2Diff
1829       
1830        WAVE rw=root:Packages:NIST:SAS:realsRead
1831        rw[25] = SSD/100                // in meters
1832        return(SSD)
1833End
1834
1835
1836// not part of SASCALC, but can be used to convert the SSD to number of guides
1837//
1838// SSD in meters
1839Function numGuides(SSD)
1840        variable SSD
1841       
1842        Variable Ng
1843        Ng = SSD*100 + 5 - 1632
1844        Ng /= -155
1845       
1846        Ng = round(Ng)
1847        return(Ng)
1848End
1849
1850
1851//returns the offset value
1852// slider and setVar are linked to the same global
1853// updates the wave and changes the beamcenter (x,y) in the wave
1854Function detectorOffset()
1855       
1856        WAVE rw=root:Packages:NIST:SAS:RealsRead
1857        NVAR val = root:Packages:NIST:SAS:gOffset
1858        rw[19] = val            // already in cm
1859        //move the beamcenter, make it an integer value for the MC simulation
1860        rw[16] = 64 + round(2*rw[19])           //approximate beam X is 64 w/no offset, 114 w/25 cm offset
1861        rw[17] = 64             //typical value
1862       
1863        return(val)
1864end
1865
1866//returns the detector distance (slider and setVar are linked by the global)
1867//
1868Function chamberToDetectorDist()
1869
1870        NVAR val = root:Packages:NIST:SAS:gDetDist
1871        return(val)
1872End
1873
1874//sets the SDD (slider and setVar are linked by the global and is the detector position
1875//  relative to the chamber)
1876// updates the wave
1877Function sampleToDetectorDist()
1878
1879        NVAR detDist = root:Packages:NIST:SAS:gDetDist
1880        NVAR S12 = root:Packages:NIST:SAS:S12
1881        WAVE rw=root:Packages:NIST:SAS:RealsRead       
1882        Variable SDD   
1883       
1884        SDD = detDist + s12*(2-tableposition())
1885        rw[18] = SDD/100                // convert to meters for header
1886        return(SDD)
1887End
1888
1889//direction = one of "vertical;horizontal;maximum;"
1890// all of this is bypassed if the lenses are in
1891//
1892Function beamDiameter(direction)
1893        String direction
1894
1895        NVAR lens = root:Packages:NIST:SAS:gUsingLenses
1896        if(lens)
1897                return sourceApertureDiam()
1898        endif
1899       
1900    Variable l1 = sourceToSampleDist()
1901    Variable l2 //= sampleAperToDetDist()
1902    Variable d1,d2,bh,bv,bm,umbra,a1,a2
1903   
1904    NVAR L2diff = root:Packages:NIST:SAS:L2diff
1905    NVAR lambda = root:Packages:NIST:SAS:gLambda
1906    NVAR lambda_width = root:Packages:NIST:SAS:gDeltaLambda
1907    NVAR bs_factor = root:Packages:NIST:SAS:bs_factor
1908   
1909    l2 = sampleToDetectorDist() + L2diff
1910    a1 = sourceApertureDiam()
1911    a2 = sampleApertureDiam()
1912   
1913    d1 = a1*l2/l1
1914    d2 = a2*(l1+l2)/l1
1915    bh = d1+d2          //beam size in horizontal direction
1916    umbra = abs(d1-d2)
1917    //vertical spreading due to gravity
1918    bv = bh + 1.25e-8*(l1+l2)*l2*lambda*lambda*lambda_width
1919    bm = (bs_factor*bh > bv) ? bs_factor*bh : bv //use the larger of horiz*safety or vertical
1920   
1921    strswitch(direction)        // string switch
1922        case "vertical":                // execute if case matches expression
1923                return(bv)
1924                break                                           // exit from switch
1925        case "horizontal":              // execute if case matches expression
1926                return(bh)
1927                break
1928        case "maximum":         // execute if case matches expression
1929                return(bm)
1930                break
1931        default:                                                        // optional default expression executed
1932                return(bm)                                              // when no case matches
1933    endswitch
1934End
1935
1936//on NG3 and NG7, allowable sizes are 1,2,3,4" diameter
1937//will return values larger than 4.0*2.54 if a larger beam is needed
1938//
1939// - in an approximate way, account for lenses
1940Function beamstopDiam()
1941
1942        NVAR yesLens = root:Packages:NIST:SAS:gUsingLenses
1943        Variable bm=0
1944        Variable bs=0.0
1945   
1946        if(yesLens)
1947                //bm = sourceApertureDiam()             //ideal result, not needed
1948                bs = 1                                                          //force the diameter to 1"
1949        else
1950                bm = beamDiameter("maximum")
1951                do
1952                bs += 1
1953           while ( (bs*2.54 < bm) || (bs > 30.0))                       //30 = ridiculous limit to avoid inf loop
1954        endif
1955
1956        //update the wave
1957        WAVE rw=root:Packages:NIST:SAS:realsRead
1958        rw[21] = bs*25.4                //store the BS diameter in mm
1959       
1960    return (bs*2.54)            //return diameter in cm, not inches for txt
1961End
1962
1963//returns the projected diameter of the beamstop at the anode plane.
1964// most noticeable at short SDD
1965//if flag == 0 use conservative estimate = largest diameter (for SASCALC, default)
1966//if flag != 0 use point aperture = average diameter (for resolution calculation)
1967Function beamstopDiamProjection(flag)
1968        Variable flag
1969       
1970        NVAR L2diff = root:Packages:NIST:SAS:L2diff
1971        Variable a2 = sampleApertureDiam()
1972        Variable bs = beamstopDiam()
1973        Variable l2, LB, BS_P
1974   
1975        l2 = sampleToDetectorDist() + L2diff
1976        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical)
1977        if(flag==0)
1978                BS_P = bs + (bs+a2)*lb/(l2-lb)          //diameter of shadow from parallax
1979        else
1980                BS_P = bs + bs*lb/(l2-lb)               //diameter of shadow, point A2
1981        endif
1982        return (bs_p)           //return projected diameter in cm
1983End
1984
1985// 19MAR07 - using correction from John for an estimate of the shadow of the beamstop
1986// at the detection plane. This is a noticeable effect at short SDD, where the projected
1987// diameter of the beamstop is much larger than the physical diameter.
1988Function qMin()
1989
1990    Variable l2s = sampleToDetectorDist()       //distance from sample to detector in cm
1991//    Variable bs = beamstopDiam()              //beamstop diameter in cm
1992    Variable bs_p = beamstopDiamProjection(0)           //projected beamstop diameter in cm
1993    NVAR lambda = root:Packages:NIST:SAS:gLambda
1994    NVAR d_det = root:Packages:NIST:SAS:d_det           //cm
1995    NVAR a_pixel = root:Packages:NIST:SAS:a_pixel       //cm
1996
1997    return( (pi/lambda)*(bs_p + d_det + a_pixel)/l2s )          //use bs_p rather than bs
1998   // return( (pi/lambda)*(bs + d_det + a_pixel)/l2s )          //use bs (incorrect)
1999End
2000
2001Function qMaxVert()
2002
2003    Variable theta
2004    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
2005    NVAR lambda = root:Packages:NIST:SAS:gLambda
2006        NVAR det_width = root:Packages:NIST:SAS:det_width
2007       
2008    theta = atan( (det_width/2.0)/l2s )
2009   
2010    return ( 4.0*pi/lambda * sin(theta/2.0) )
2011end
2012
2013Function qMaxCorner()
2014
2015    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
2016    Variable radial
2017    NVAR lambda = root:Packages:NIST:SAS:gLambda
2018        NVAR det_off = root:Packages:NIST:SAS:gOffset
2019        NVAR det_width = root:Packages:NIST:SAS:det_width
2020
2021    radial=sqrt( (0.5*det_width)*(0.5*det_width)+(0.5*det_width+det_off)*(0.5*det_width+det_off) )
2022   
2023    return ( 4*pi/lambda*sin(0.5*atan(radial/l2s)) )
2024End
2025
2026Function qMaxHoriz()
2027
2028    Variable theta
2029    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
2030    NVAR lambda = root:Packages:NIST:SAS:gLambda
2031        NVAR det_off = root:Packages:NIST:SAS:gOffset
2032        NVAR det_width = root:Packages:NIST:SAS:det_width
2033
2034    theta = atan( ((det_width/2.0)+det_off)/l2s )       //from the instance variables
2035   
2036    return ( 4.0*pi/lambda * sin(theta/2.0) )
2037End
2038
2039// calculate sigma for the resolution function at either limit of q-range
2040Function deltaQ(atQ)
2041        Variable atQ
2042       
2043    Variable k02,lp,l1,l2,sig_02,sigQ2,a1,a2
2044    NVAR l2Diff = root:Packages:NIST:SAS:L2diff
2045        NVAR lambda = root:Packages:NIST:SAS:gLambda
2046        NVAR lambda_width = root:Packages:NIST:SAS:gDeltaLambda
2047        NVAR d_det = root:Packages:NIST:SAS:d_det
2048        NVAR del_r = root:Packages:NIST:SAS:del_r
2049       
2050       
2051    l1 = sourceToSampleDist()
2052    l2 = sampleToDetectorDist() + L2diff
2053    a1 = sourceApertureDiam()
2054    a2 = sampleApertureDiam()
2055   
2056    k02 = (6.2832/lambda)*(6.2832/lambda)
2057    lp = 1/(1/l1 + 1/l2)
2058   
2059    sig_02 = (0.25*a1/l1)*(0.25*a1/l1)
2060    sig_02 += (0.25*a2/lp)*(0.25*a2/lp)
2061    sig_02 += (d_det/(2.355*l2))*(d_det/(2.355*l2))
2062    sig_02 += (del_r/l2)*(del_r/l2)/12
2063    sig_02 *= k02
2064   
2065    sigQ2 = sig_02 + (atQ*lambda_width)*(atQ*lambda_width)/6
2066
2067    return(100*sqrt(sigQ2)/atQ)
2068End
2069
2070
2071Function beamIntensity()
2072
2073    Variable alpha,f,t,t4,t5,t6,as,solid_angle,l1,d2_phi
2074    Variable a1,a2,retVal
2075    SetDataFolder root:Packages:NIST:SAS
2076    NVAR l_gap=l_gap,guide_width =guide_width,ng = gNg
2077    NVAR lambda_t=lambda_t,b=b,c=c
2078    NVAR lambda=gLambda,t1=t1,t2=t2,t3=t3,phi_0=phi_0
2079    NVAR lambda_width=gDeltaLambda
2080   
2081    l1 = sourceToSampleDist()
2082    a1 = sourceApertureDiam()
2083    a2 = sampleApertureDiam()
2084   
2085   
2086    alpha = (a1+a2)/(2*l1)      //angular divergence of beam
2087    f = l_gap*alpha/(2*guide_width)
2088    t4 = (1-f)*(1-f)
2089    t5 = exp(ng*ln(0.96))       // trans losses of guides in pre-sample flight
2090    t6 = 1 - lambda*(b-(ng/8)*(b-c))            //experimental correction factor
2091    t = t1*t2*t3*t4*t5*t6
2092   
2093    as = pi/4*a2*a2             //area of sample in the beam
2094    d2_phi = phi_0/(2*pi)
2095    d2_phi *= exp(4*ln(lambda_t/lambda))
2096    d2_phi *= exp(-1*(lambda_t*lambda_t/lambda/lambda))
2097
2098    solid_angle = pi/4* (a1/l1)*(a1/l1)
2099
2100    retVal = as * d2_phi * lambda_width * solid_angle * t
2101     SetDataFolder root:
2102    return (retVal)
2103end
2104
2105Function figureOfMerit()
2106
2107        Variable bi = beamIntensity()
2108        NVAR lambda = root:Packages:NIST:SAS:gLambda
2109       
2110    return (lambda*lambda*bi)
2111End
2112
2113//estimate the number of pixels in the beam, and enforce the maximum countrate per pixel (idmax)
2114Function attenuatorTransmission()
2115
2116    Variable num_pixels,i_pix           //i_pix = id in John's notation
2117    Variable bDiam = beamDiameter("horizontal") //!! note that prev calculations used bh (horizontal)
2118    Variable atten,a2
2119    SetDataFolder root:Packages:NIST:SAS
2120    NVAR a_pixel=a_pixel,idmax=idmax
2121   
2122    a2 = sampleApertureDiam()
2123   
2124    num_pixels = pi/4*(0.5*(a2+bDiam))*(0.5*(a2+bDiam))/a_pixel/a_pixel
2125    i_pix = ( beamIntensity() )/num_pixels
2126   
2127    atten = (i_pix < idmax) ? 1.0 : idmax/i_pix
2128    SetDataFolder root:
2129    return(atten)
2130End
2131
2132Function attenuatorNumber()
2133
2134    Variable atten = attenuatorTransmission()
2135    Variable af,nf,numAtten
2136    SetDataFolder root:Packages:NIST:SAS
2137    NVAR lambda=gLambda
2138   
2139    af = 0.498 + 0.0792*lambda - 1.66e-3*lambda*lambda
2140    nf = -ln(atten)/af          //floating point
2141   
2142    numAtten = trunc(nf) + 1                    //in c, (int)nf
2143    //correct for larger step thickness at n > 6
2144    if(numAtten > 6)
2145        numAtten = 7 + trunc( (numAtten-6)/2 )          //in c, numAtten = 7 + (int)( (numAtten-6)/2 )
2146    endif
2147   
2148    SetDatafolder root:
2149    return (numAtten)
2150End
Note: See TracBrowser for help on using the repository browser.