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

Last change on this file since 472 was 472, checked in by srkline, 14 years ago

A large number of changes and fixes:

--168/169/170: panels and windows are now at least on-screen for all packages. Tested
with 1024x768 resolution.
-- closed ticket 176 which was a question about resampling data to generate error estimates
on fitted parameters. Useful for reflectometry, not needed for SANS.
--157: bug of low Q power law extrapolation in Invariant fixed by avoiding q==0
--178/180: Tr/Tw? notification in USANS. log/lin checkbox for display.
--167: saveData checkbox for USANS not behaving well. turns off/on better now.
--197: changed all (?) 1D writing routines to enforce 26 characters as the maximum length
to make sure that file loading will never exceed 31 characters

-- lots of changes to MonteCarlo? and SASCALC

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