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

Last change on this file since 842 was 842, checked in by srkline, 10 years ago

Additions to SASCALC for the new sample apertures in boxes 6 and 7 @ NG3. Now 6 guides is an allowed configuration. A polarizer has been added to NG7 in box 2, but no changes have been made to SASCALC, as no guide configurations are excluded.

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