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

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

Made preferences a common panel (moved to PlotUtilsMacro?.ipf and globals to root:Packages:NIST:) and added menu items for all packages. Many files had to be modified so that the preferences could be properly accessed

File Open dialog now is set to "All files" so that XML can be selected. I think that all open access that doesn't already have the full path go through this common function.

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