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

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

changes to how a fake VAX raw data file is written out. now call accepts optional parameters so that the index and label can be supplied, and the MC simulation can be very crudely scripted for unattended operation.

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