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

Last change on this file since 752 was 752, checked in by srkline, 12 years ago

explicitly named the window for the control actions (win=SASCALC) to prevent the controls from being randomly drawn on other panels. This seems to happen frequently to me with the magic mouse. Inadvertent scrolls move controls on bacground windows, drawing them on the top window.

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