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

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

Procedures for mostly my use that enable me to (crudely) script the MC simulations so that they can be done in a big batch overnight.

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