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

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

Added a new template "TotalTemplate?.pxt" that includes the package loader as part of the macros menu. from there, everything can be loaded (but nothing unloaded).

Many changes to SASCALC and MultiScatter? to make them correct, and play nice together. An XOP version is functional, but has not yet been added, and won't be until the MC simulation is checked for correctness. Without incoherent scattering, it looks good.

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