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

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

Put I(0) errors in most of the places that I could make sense to put them. #215

Also removed FIT_Ops entirely by copying FITRPA functions to LinearizedFIts, and simply pointing FIT_Ops to include LinearizedFits? (this was done for backwards compatibility). part of ticket #275

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