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

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

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