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

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

Lots of little changes:

changed debye function to a limiting value at low qRg to avoid numerical errors. There is no XOP, so no changes needed there.

Added MonteCarlo? simulation routine (1st draft) to the SANS Reduction, including changes to SASCALC to implement. Still very alpha with lots of bits to add. See the top of the monte carlo file for the list.

added MonteCarlo? file to SANS includes
added SAS folder to 2D display list to view MonteCarlo? results
removed more "non-functions" from the function popup in analysis. these were a result of motofit/genFit, and simultaneous loading of SANS/USANS reduction. More will follow.

File size: 53.2 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
650       
651//      Variable imon,thick,r2,sig_incoh
652//      String funcStr
653//      imon = 10000
654//      thick = 0.1
655//      sig_incoh = 0.1
656//      funcStr = "SphereForm"
657//      r2 = 2.54/2                     //typical 1" diameter sample, convert to radius in cm
658
659        NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo                // == 1 if MC, 0 if other
660        SVAR funcStr = root:Packages:NIST:SAS:gFuncStr
661
662        if(doMonteCarlo == 1)
663                WAVE rw=root:Packages:NIST:SAS:realsRead
664               
665                NVAR imon = root:Packages:NIST:SAS:gImon
666                NVAR thick = root:Packages:NIST:SAS:gThick
667                NVAR sig_incoh = root:Packages:NIST:SAS:gSig_incoh
668                NVAR r2 = root:Packages:NIST:SAS:gR2
669       
670                r1 = rw[24]/2/10                // sample diameter convert diam in [mm] to radius in cm
671                xCtr = rw[16]
672                yCtr = rw[17]
673                sdd = rw[18]*100                //conver header of [m] to [cm]
674                pixSize = rw[10]/10             // convert pix size in mm to cm
675                wavelength = rw[26]
676                coefStr = MC_getFunctionCoef(funcStr)
677               
678                if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
679                        Abort "The coefficients and function type do not match. Please correct the selections in the popup menus."
680                endif
681               
682                FUNCREF SANSModelAAO_MCproto func=$funcStr
683                WAVE results = root:Packages:NIST:SAS:results
684                results = 0
685               
686                Monte_SANS(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,func,$coefStr,results)
687               
688                // convert to absolute scale
689                Variable kappa,beaminten = beamIntensity()
690                // results[6] is the fraction transmitted
691//              kappa = beamInten*pi*r1*r1*thick*(pixSize/sdd)^2*results[6]*(iMon/beaminten)
692                kappa = thick*(pixSize/sdd)^2*results[6]*iMon *2        //why the factor of 2?
693               
694                WAVE linear_data = root:Packages:NIST:SAS:linear_data
695                WAVE data = root:Packages:NIST:SAS:data
696                linear_data = linear_data / kappa
697                data = linear_data
698       
699        endif
700       
701        // update the wave with the beamstop diameter here, since I don't know what
702        // combinations of parameters will change the BS - but anytime the curve is
703        // recalculated, or the text displayed, the right BS must be present
704        beamstopDiam()
705       
706        S_CircularAverageTo1D("SAS")
707        WAVE aveint=root:Packages:NIST:SAS:aveint
708        WAVE qval=root:Packages:NIST:SAS:qval
709        WAVE fSubS=root:Packages:NIST:SAS:fSubS
710       
711        //multiply by sphere FF (S_SphereForm(scale,radius,delrho,bkg,x))
712        // or Debye Function S_Debye(scale,rg,bkg,x)
713       
714        //aveint = S_SphereForm(1,80,1e-6,0,qval)
715        if(doMonteCarlo != 1)
716                if(exists(funcStr) != 0)
717                        FUNCREF SANSModelAAO_MCproto func=$funcStr
718                        coefStr = MC_getFunctionCoef(funcStr)
719                       
720                        if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
721                                Abort "The coefficients and function type do not match. Please correct the selections in the popup menus."
722                        endif
723                        func($coefStr,aveint,qval)
724                else
725                        aveint = S_Debye(1000,100,0.0,qval)
726                endif
727        endif
728       
729        // multiply either estimate by beamstop shadowing
730        aveint *= fSubS
731       
732        //display the configuration text in a separate notebook
733        DisplayConfigurationText()
734       
735        return(0)
736End
737
738//freezes the current configuration
739// -1- duplicates the trace on the graph
740// -2- copies the configuration text to a second notebook window for printing
741Function FreezeButtonProc(ctrlName) : ButtonControl
742        String ctrlName
743
744        String str=""
745        NVAR ct=root:Packages:NIST:SAS:gFreezeCount
746       
747       
748        SetDataFolder root:Packages:NIST:SAS
749       
750        Duplicate/O aveint,$("aveint_"+num2str(ct))
751        Duplicate/O qval,$("qval_"+num2str(ct))
752        Appendtograph $("aveint_"+num2str(ct)) vs $("qval_"+num2str(ct))
753        ModifyGraph mode=3
754        ModifyGraph marker=19
755        switch(mod(ct,10))      // 10 different colors - black is the unfrozen color
756                case 0:
757                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,16385,16385)
758                        break
759                case 1:
760                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(2,39321,1)
761                        break
762                case 2:
763                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(0,0,65535)
764                        break
765                case 3:
766                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(39321,1,31457)
767                        break
768                case 4:
769                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(48059,48059,48059)
770                        break
771                case 5:
772                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,32768,32768)
773                        break
774                case 6:
775                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(0,65535,0)
776                        break
777                case 7:
778                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(16385,65535,65535)
779                        break
780                case 8:
781                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,32768,58981)
782                        break
783                case 9:
784                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(36873,14755,58982)
785                        break
786        endswitch
787       
788        NVAR offset = root:Packages:NIST:SAS:gModelOffsetFactor
789        offset = 2^ct
790        //multiply by current offset (>=1)
791        Wave inten = $("aveint_"+num2str(ct))
792        inten *= offset
793//      Print "new offset = ",offset
794       
795        ct +=1
796        SetDataFolder root:
797       
798        // create or append the configuration to a new notebook
799        if(WinType("Saved_Configurations")==0)
800                NewNotebook/F=1/N=Saved_Configurations /W=(480,400,880,725)
801                DoWindow/F SASCALC              //return focus to SASCALC window
802        endif
803        //append the text
804        sprintf str,"\rConfiguration #%d\r",ct-1
805        Notebook Saved_Configurations showRuler=0,defaultTab=20,selection={endOfFile, endOfFile}
806        Notebook Saved_Configurations font="Monaco",fSize=10,fstyle=1,text=str          //bold
807        Notebook Saved_Configurations font="Monaco",fSize=10,fstyle=0,text=SetConfigurationText()
808        return(0)
809End
810
811//clears the frozen traces on the graph, asks if you want to clear the saved text
812//
813Function S_ClearButtonProc(ctrlName) : ButtonControl
814        String ctrlName
815
816        NVAR ct=root:Packages:NIST:SAS:gFreezeCount
817        Variable ii
818        Setdatafolder root:Packages:NIST:SAS
819        for(ii=ct-1;ii>=1;ii-=1)
820        //remove all traces, replace aveint
821        // kill all waves _ct
822                RemoveFromGraph $("aveint_"+num2str(ii))
823                Killwaves/Z $("aveint_"+num2str(ii)),$("qval_"+num2str(ii))
824        endfor
825        ct=1
826        setdatafolder root:
827       
828        DoAlert 1,"Do you also want to clear the \"Saved Configurations\" text?"
829        if(V_flag == 1)         // yes
830                DoWindow/K Saved_Configurations
831        endif
832       
833        //reset offset value
834        NVAR offset = root:Packages:NIST:SAS:gModelOffsetFactor
835        offset = 1
836        ReCalculateInten(1)
837        return(0)
838End
839
840
841///////////////////////////////////////////////////////////
842// 19MAR07 uses correction for beamstop diameter projection to get shadow factor correct
843//
844Function S_CircularAverageTo1D(type)
845        String type
846       
847        Variable isCircular = 1
848       
849        //type is the data type to do the averaging on, and will be set as the current folder
850        //get the current displayed data (so the correct folder is used)
851        String destPath = "root:Packages:NIST:"+type
852       
853        //
854        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist,dr,ddr
855        Variable lambda,trans
856        WAVE reals = $(destPath + ":RealsRead")
857//      WAVE/T textread = $(destPath + ":TextRead")
858//      String fileStr = textread[3]
859       
860        // center of detector, for non-linear corrections
861        Variable pixelsX=128,pixelsY=128
862       
863        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
864        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
865       
866        // beam center, in pixels
867        x0 = reals[16]
868        y0 = reals[17]
869        //detector calibration constants
870        sx = reals[10]          //mm/pixel (x)
871        sx3 = reals[11]         //nonlinear coeff
872        sy = reals[13]          //mm/pixel (y)
873        sy3 = reals[14]         //nonlinear coeff
874       
875        dtsize = 10*reals[20]           //det size in mm
876        dtdist = 1000*reals[18]         // det distance in mm
877       
878        NVAR binWidth=root:Packages:NIST:SAS:gBinWidth
879//      Variable binWidth = 1
880       
881        dr = binWidth           // ***********annulus width set by user, default is one***********
882        ddr = dr*sx             //step size, in mm (this value should be passed to the resolution calculation, not dr 18NOV03)
883               
884        Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,dxbm,dybm,ii
885        Variable phi_rad,dphi_rad,phi_x,phi_y
886        Variable forward,mirror
887       
888        String side = "both"
889//      side = StringByKey("SIDE",keyListStr,"=",";")
890//      Print "side = ",side
891       
892        /// data wave is data in the current folder which was set at the top of the function
893        WAVE data=$(destPath + ":data")
894
895// fake mask that uses all of the detector
896        Make/O/N=(pixelsX,pixelsY) $(destPath + ":mask")
897        Wave mask = $(destPath + ":mask")
898        mask = 0
899        //two pixels all around
900        mask[0,1][] = 1
901        mask[126,127][] = 1
902        mask[][0,1] = 1
903        mask[][126,127] = 1
904        //
905        //pixels within rcentr of beam center are broken into 9 parts (units of mm)
906        rcentr = 100            //original
907//      rcentr = 0
908        // values for error if unable to estimate value
909        //large_num = 1e10
910        large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1
911        small_num = 1e-10
912       
913        // output wave are expected to exist (?) initialized to zero, what length?
914        // 200 points on VAX --- use 300 here, or more if SAXS data is used with 1024x1024 detector (1000 pts seems good)
915        Variable defWavePts=500
916        Make/O/N=(defWavePts) $(destPath + ":qval"),$(destPath + ":aveint")
917        Make/O/N=(defWavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave")
918        Make/O/N=(defWavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar")
919
920        WAVE qval = $(destPath + ":qval")
921        WAVE aveint = $(destPath + ":aveint")
922        WAVE ncells = $(destPath + ":ncells")
923        WAVE dsq = $(destPath + ":dsq")
924        WAVE sigave = $(destPath + ":sigave")
925        WAVE qbar = $(destPath + ":QBar")
926        WAVE sigmaq = $(destPath + ":SigmaQ")
927        WAVE fsubs = $(destPath + ":fSubS")
928       
929        qval = 0
930        aveint = 0
931        ncells = 0
932        dsq = 0
933        sigave = 0
934        qbar = 0
935        sigmaq = 0
936        fsubs = 0
937
938        dtdis2 = dtdist^2
939        nq = 1
940        xoffst=0
941        //distance of beam center from detector center
942        dxbm = S_FX(x0,sx3,xcenter,sx)
943        dybm = S_FY(y0,sy3,ycenter,sy)
944               
945        //BEGIN AVERAGE **********
946        Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1
947        Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p
948       
949        // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
950        // loop index corresponds to FORTRAN (old code)
951        // and the IGOR array indices must be adjusted (-1) to the correct address
952        ii=1
953        do
954                xi = ii
955                dxi = S_FX(xi,sx3,xcenter,sx)
956                dx = dxi-dxbm           //dx and dy are in mm
957               
958                jj = 1
959                do
960                        data_pixel = data[ii-1][jj-1]           //assign to local variable
961                        yj = jj
962                        dyj = S_FY(yj,sy3,ycenter,sy)
963                        dy = dyj - dybm
964                        if(!(mask[ii-1][jj-1]))                 //masked pixels = 1, skip if masked (this way works...)
965                                dr2 = (dx^2 + dy^2)^(0.5)               //distance from beam center NOTE dr2 used here - dr used above
966                                if(dr2>rcentr)          //keep pixel whole
967                                        nd = 1
968                                        fd = 1
969                                else                            //break pixel into 9 equal parts
970                                        nd = 3
971                                        fd = 2
972                                endif
973                                nd2 = nd^2
974                                ll = 1          //"el-el" loop index
975                                do
976                                        dxx = dx + (ll - fd)*sx/3
977                                        kk = 1
978                                        do
979                                                dyy = dy + (kk - fd)*sy/3
980                                                if(isCircular)
981                                                        //circular average, use all pixels
982                                                        //(increment)
983                                                        nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
984                                                else
985                                                        //a sector average - determine azimuthal angle
986                                                        dphi_p = S_dphi_pixel(dxx,dyy,phi_x,phi_y)
987                                                        if(dphi_p < dphi_rad)
988                                                                forward = 1                     //within forward sector
989                                                        else
990                                                                forward = 0
991                                                        Endif
992                                                        if((Pi - dphi_p) < dphi_rad)
993                                                                mirror = 1              //within mirror sector
994                                                        else
995                                                                mirror = 0
996                                                        Endif
997                                                        //check if pixel lies within allowed sector(s)
998                                                        if(cmpstr(side,"both")==0)              //both sectors
999                                                                if ( mirror || forward)
1000                                                                        //increment
1001                                                                        nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1002                                                                Endif
1003                                                        else
1004                                                                if(cmpstr(side,"right")==0)             //forward sector only
1005                                                                        if(forward)
1006                                                                                //increment
1007                                                                                nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1008                                                                        Endif
1009                                                                else                    //mirror sector only
1010                                                                        if(mirror)
1011                                                                                //increment
1012                                                                                nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1013                                                                        Endif
1014                                                                Endif
1015                                                        Endif           //allowable sectors
1016                                                Endif   //circular or sector check
1017                                                kk+=1
1018                                        while(kk<=nd)
1019                                        ll += 1
1020                                while(ll<=nd)
1021                        Endif           //masked pixel check
1022                        jj += 1
1023                while (jj<=pixelsY)
1024                ii += 1
1025        while(ii<=pixelsX)              //end of the averaging
1026               
1027        //compute q-values and errors
1028        Variable ntotal,rr,theta,avesq,aveisq,var
1029       
1030        lambda = reals[26]
1031        ntotal = 0
1032        kk = 1
1033        do
1034                rr = (2*kk-1)*ddr/2
1035                theta = 0.5*atan(rr/dtdist)
1036                qval[kk-1] = (4*Pi/lambda)*sin(theta)
1037                if(ncells[kk-1] == 0)
1038                        //no pixels in annuli, data unknown
1039                        aveint[kk-1] = 0
1040                        sigave[kk-1] = large_num
1041                else
1042                        if(ncells[kk-1] <= 1)
1043                                //need more than one pixel to determine error
1044                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
1045                                sigave[kk-1] = large_num
1046                        else
1047                                //assume that the intensity in each pixel in annuli is normally
1048                                // distributed about mean...
1049                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
1050                                avesq = aveint[kk-1]^2
1051                                aveisq = dsq[kk-1]/ncells[kk-1]
1052                                var = aveisq-avesq
1053                                if(var<=0)
1054                                        sigave[kk-1] = small_num
1055                                else
1056                                        sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1))
1057                                endif
1058                        endif
1059                endif
1060                ntotal += ncells[kk-1]
1061                kk+=1
1062        while(kk<=nq)
1063       
1064        //Print "NQ = ",nq
1065        // data waves were defined as 300 points (=defWavePts), but now have less than that (nq) points
1066        // use DeletePoints to remove junk from end of waves
1067        //WaveStats would be a more foolproof implementation, to get the # points in the wave
1068        Variable startElement,numElements
1069        startElement = nq
1070        numElements = defWavePts - startElement
1071        DeletePoints startElement,numElements, qval,aveint,ncells,dsq,sigave
1072       
1073        //////////////end of VAX sector_ave()
1074               
1075
1076
1077// ***************************************************************
1078//
1079// Do the extra 3 columns of resolution calculations starting here.
1080//
1081// ***************************************************************
1082
1083        Variable L2 = reals[18]
1084//      Variable BS = reals[21]         //this the diameter is stored in mm
1085        Variable BS = beamstopDiamProjection(1) * 10            //calculated projection in cm *10 = mm
1086        Variable S1 = reals[23]
1087        Variable S2 = reals[24]
1088        Variable L1 = reals[25]
1089        lambda = reals[26]
1090        Variable lambdaWidth = reals[27]
1091
1092        Variable DDet, apOff
1093        //typical value for NG3 and NG7 - distance between sample aperture and sample in (cm)
1094        apOff=5.0
1095        // hard wire value for Ordela detectors
1096        DDet = 0.5              // resolution in cm
1097        //      String detStr=textRead[9]
1098        //      DDet = DetectorPixelResolution(fileStr,detStr)          //needs detector type and beamline
1099
1100        //Go from 0 to nq doing the calc for all three values at
1101        //every Q value
1102
1103        ii=0
1104        Variable ret1,ret2,ret3
1105        do
1106                S_getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,ret1,ret2,ret3)
1107                sigmaq[ii] = ret1       //res_wave[0]
1108                qbar[ii] = ret2         //res_wave[1]
1109                fsubs[ii] = ret3                //res_wave[2]
1110                ii+=1
1111        while(ii<nq)
1112        DeletePoints startElement,numElements, sigmaq, qbar, fsubs
1113
1114        fsubs += 1e-8           //keep the values from being too small
1115
1116// End of resolution calculations
1117// ***************************************************************
1118       
1119        //get rid of the default mask, if one was created (it is in the current folder)
1120        //don't just kill "mask" since it might be pointing to the one in the MSK folder
1121        Killwaves/Z $(destPath+":mask")
1122       
1123        //return to root folder (redundant)
1124        SetDataFolder root:
1125       
1126        Return 0
1127End
1128
1129//returns nq, new number of q-values
1130//arrays aveint,dsq,ncells are also changed by this function
1131//
1132Function S_IncrementPixel(dataPixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1133        Variable dataPixel,ddr,dxx,dyy
1134        Wave aveint,dsq,ncells
1135        Variable nq,nd2
1136       
1137        Variable ir
1138       
1139        ir = trunc(sqrt(dxx*dxx+dyy*dyy)/ddr)+1
1140        if (ir>nq)
1141                nq = ir         //resets maximum number of q-values
1142        endif
1143        aveint[ir-1] += dataPixel/nd2           //ir-1 must be used, since ir is physical
1144        dsq[ir-1] += dataPixel*dataPixel/nd2
1145        ncells[ir-1] += 1/nd2
1146       
1147        Return nq
1148End
1149
1150//function determines azimuthal angle dphi that a vector connecting
1151//center of detector to pixel makes with respect to vector
1152//at chosen azimuthal angle phi -> [cos(phi),sin(phi)] = [phi_x,phi_y]
1153//dphi is always positive, varying from 0 to Pi
1154//
1155Function S_dphi_pixel(dxx,dyy,phi_x,phi_y)
1156        Variable dxx,dyy,phi_x,phi_y
1157       
1158        Variable val,rr,dot_prod
1159       
1160        rr = sqrt(dxx^2 + dyy^2)
1161        dot_prod = (dxx*phi_x + dyy*phi_y)/rr
1162        //? correct for roundoff error? - is this necessary in IGOR, w/ double precision?
1163        if(dot_prod > 1)
1164                dot_prod =1
1165        Endif
1166        if(dot_prod < -1)
1167                dot_prod = -1
1168        Endif
1169       
1170        val = acos(dot_prod)
1171       
1172        return val
1173
1174End
1175
1176//calculates the x distance from the center of the detector, w/nonlinear corrections
1177//
1178Function S_FX(xx,sx3,xcenter,sx)               
1179        Variable xx,sx3,xcenter,sx
1180       
1181        Variable retval
1182       
1183        retval = sx3*tan((xx-xcenter)*sx/sx3)
1184        Return retval
1185End
1186
1187//calculates the y distance from the center of the detector, w/nonlinear corrections
1188//
1189Function S_FY(yy,sy3,ycenter,sy)               
1190        Variable yy,sy3,ycenter,sy
1191       
1192        Variable retval
1193       
1194        retval = sy3*tan((yy-ycenter)*sy/sy3)
1195        Return retval
1196End
1197
1198//**********************
1199// Resolution calculation - used by the averaging routines
1200// to calculate the resolution function at each q-value
1201// - the return value is not used
1202//
1203// equivalent to John's routine on the VAX Q_SIGMA_AVE.FOR
1204// Incorporates eqn. 3-15 from J. Appl. Cryst. (1995) v. 28 p105-114
1205//
1206Function/S S_getResolution(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,SigmaQ,QBar,fSubS)
1207        Variable inQ, lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r
1208        Variable &fSubS, &QBar, &SigmaQ         //these are the output quantities at the input Q value
1209       
1210        //lots of calculation variables
1211        Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g
1212        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
1213
1214        //Constants
1215        //Variable del_r = .1
1216        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
1217        Variable g = 981.0                              //gravity acceleration [cm/s^2]
1218
1219        String results
1220        results ="Failure"
1221       
1222        NVAR usingLenses = root:Packages:NIST:SAS:gUsingLenses
1223
1224        //rename for working variables,  these must be gotten from global
1225        //variables
1226
1227//      Variable wLam, wLW, wL1, wL2, wS1, wS2
1228//      Variable wBS, wDDet, wApOff
1229//      wLam = lambda
1230//      wLW = lambdaWidth
1231//      wDDet = DDet
1232//      wApOff = apOff
1233        S1 *= 0.5*0.1                   //convert to radius and [cm]
1234        S2 *= 0.5*0.1
1235
1236        L1 *= 100.0                     // [cm]
1237        L1 -= apOff                             //correct the distance
1238
1239        L2 *= 100.0
1240        L2 += apOff
1241
1242        BS *= 0.5*0.1                   //convert to radius and [cm]
1243        del_r *= 0.1                            //width of annulus, convert mm to [cm]
1244       
1245        //Start resolution calculation
1246        a2 = S1*L2/L1 + S2*(L1+L2)/L1
1247        q_small = 2.0*Pi*(BS-a2)*(1.0-lambdaWidth)/(lambda*L2)
1248        lp = 1.0/( 1.0/L1 + 1.0/L2)
1249
1250        v_lambda = lambdaWidth^2/6.0
1251       
1252        if(usingLenses==1)                      //SRK 2007
1253                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term
1254        else
1255                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form
1256        endif
1257       
1258        v_d = (DDet/2.3548)^2 + del_r^2/12.0
1259        vz = vz_1 / lambda
1260        yg = 0.5*g*L2*(L1+L2)/vz^2
1261        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007
1262
1263        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
1264        delta = 0.5*(BS - r0)^2/v_d
1265
1266        if (r0 < BS)
1267                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
1268        else
1269                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
1270        endif
1271
1272        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
1273        if (fSubS <= 0.0)
1274                fSubS = 1.e-10
1275        endif
1276        fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
1277        fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
1278
1279        rmd = fr*r0
1280        v_r1 = v_b + fv*v_d +v_g
1281
1282        rm = rmd + 0.5*v_r1/rmd
1283        v_r = v_r1 - 0.5*(v_r1/rmd)^2
1284        if (v_r < 0.0)
1285                v_r = 0.0
1286        endif
1287        QBar = (4.0*Pi/lambda)*sin(0.5*atan(rm/L2))
1288        SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda)
1289
1290        results = "success"
1291        Return results
1292End
1293
1294Function S_Debye(scale,rg,bkg,x)
1295        Variable scale,rg,bkg
1296        Variable x
1297       
1298        // variables are:
1299        //[0] scale factor
1300        //[1] radius of gyration [A]
1301        //[2] background        [cm-1]
1302       
1303        // calculates (scale*debye)+bkg
1304        Variable Pq,qr2
1305       
1306        qr2=(x*rg)^2
1307        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
1308       
1309        //scale
1310        Pq *= scale
1311        // then add in the background
1312        return (Pq+bkg)
1313End
1314
1315
1316Function S_SphereForm(scale,radius,delrho,bkg,x)                               
1317        Variable scale,radius,delrho,bkg
1318        Variable x
1319       
1320        // variables are:                                                       
1321        //[0] scale
1322        //[1] radius (A)
1323        //[2] delrho (A-2)
1324        //[3] background (cm-1)
1325       
1326//      Variable scale,radius,delrho,bkg                               
1327//      scale = w[0]
1328//      radius = w[1]
1329//      delrho = w[2]
1330//      bkg = w[3]
1331       
1332       
1333        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
1334        // and is rescaled to give [=] cm^-1
1335       
1336        Variable bes,f,vol,f2
1337        //
1338        //handle q==0 separately
1339        If(x==0)
1340                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
1341                return(f)
1342        Endif
1343       
1344        bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
1345        vol = 4*pi/3*radius^3
1346        f = vol*bes*delrho              // [=] A
1347        // normalize to single particle volume, convert to 1/cm
1348        f2 = f * f / vol * 1.0e8                // [=] 1/cm
1349       
1350        return (scale*f2+bkg)   // Scale, then add in the background
1351       
1352End
1353
1354Function/S SetConfigurationText()
1355
1356        String str="",temp
1357
1358        SetDataFolder root:Packages:NIST:SAS
1359       
1360        NVAR numberOfGuides=gNg
1361        NVAR gTable=gTable              //2=chamber, 1=table
1362        NVAR wavelength=gLambda
1363        NVAR lambdaWidth=gDeltaLambda
1364        NVAR instrument = instrument
1365        NVAR L2diff = L2diff
1366        NVAR lens = root:Packages:NIST:SAS:gUsingLenses
1367        SVAR aStr = root:myGlobals:gAngstStr
1368       
1369        sprintf temp,"Source Aperture Diameter =\t\t%6.2f cm\r",sourceApertureDiam()
1370        str += temp
1371        sprintf temp,"Source to Sample =\t\t\t\t%6.0f cm\r",sourceToSampleDist()
1372        str += temp
1373        sprintf temp,"Sample Aperture to Detector =\t%6.0f cm\r",sampleToDetectorDist() + L2diff
1374        str += temp
1375        sprintf temp,"Beam diameter =\t\t\t\t\t%6.2f cm\r",beamDiameter("maximum")
1376        str += temp
1377        sprintf temp,"Beamstop diameter =\t\t\t\t%6.2f inches\r",beamstopDiam()/2.54
1378        str += temp
1379        sprintf temp,"Minimum Q-value =\t\t\t\t%6.4f 1/%s (sigQ/Q = %3.1f %%)\r",qMin(),aStr,deltaQ(qMin())
1380        str += temp
1381        sprintf temp,"Maximum Horizontal Q-value =\t%6.4f 1/%s\r",qMaxHoriz(),aStr
1382        str += temp
1383        sprintf temp,"Maximum Vertical Q-value =\t\t%6.4f 1/%s\r",qMaxVert(),aStr
1384        str += temp
1385        sprintf temp,"Maximum Q-value =\t\t\t\t%6.4f 1/%s (sigQ/Q = %3.1f %%)\r",qMaxCorner(),aStr,deltaQ(qMaxCorner())
1386        str += temp
1387        sprintf temp,"Beam Intensity =\t\t\t\t%.0f counts/s\r",beamIntensity()
1388        str += temp
1389        sprintf temp,"Figure of Merit =\t\t\t\t%3.3g %s^2/s\r",figureOfMerit(),aStr
1390        str += temp
1391        sprintf temp,"Attenuator transmission =\t\t%3.3g = Atten # %d\r",attenuatorTransmission(),attenuatorNumber()
1392        str += temp
1393//     
1394//      // add text of the user-edited values
1395//      //
1396        sprintf temp,"***************** NG %d *****************\r",instrument
1397        str += temp
1398        sprintf temp,"Sample Aperture Diameter =\t\t\t\t%.2f cm\r",sampleApertureDiam()
1399        str += temp
1400        sprintf temp,"Number of Guides =\t\t\t\t\t\t%d \r", numberOfGuides
1401        str += temp
1402        sprintf temp,"Sample Chamber to Detector =\t\t\t%.1f cm\r", chamberToDetectorDist()
1403        str += temp
1404        if(gTable==1)
1405                sprintf temp,"Sample Position is \t\t\t\t\t\tHuber\r"
1406        else
1407                sprintf temp,"Sample Position is \t\t\t\t\t\tChamber\r"
1408        endif
1409        str += temp
1410        sprintf temp,"Detector Offset =\t\t\t\t\t\t%.1f cm\r", detectorOffset()
1411        str += temp
1412        sprintf temp,"Neutron Wavelength =\t\t\t\t\t%.2f %s\r", wavelength,aStr
1413        str += temp
1414        sprintf temp,"Wavelength Spread, FWHM =\t\t\t\t%.3f\r", lambdaWidth
1415        str += temp
1416        sprintf temp,"Sample Aperture to Sample Position =\t%.2f cm\r", L2Diff
1417        str += temp
1418        if(lens==1)
1419                sprintf temp,"Lenses are IN\r"
1420        else
1421                sprintf temp,"Lenses are OUT\r"
1422        endif
1423        str += temp
1424       
1425   setDataFolder root:
1426   return str                   
1427End
1428
1429Function DisplayConfigurationText()
1430
1431        if(WinType("Trial_Configuration")==0)
1432                NewNotebook/F=0/K=1/N=Trial_Configuration /W=(480,44,880,369)
1433        endif
1434        //replace the text
1435        Notebook Trial_Configuration selection={startOfFile, endOfFile}
1436        Notebook Trial_Configuration font="Monaco",fSize=10,text=SetConfigurationText()
1437        return(0)
1438end
1439
1440//parses the control for A1 diam
1441// updates the wave
1442Function sourceApertureDiam()
1443        ControlInfo/W=SASCALC popup0
1444        Variable diam
1445        sscanf S_Value, "%g cm", diam
1446        WAVE rw=root:Packages:NIST:SAS:realsRead
1447        rw[23] = diam*10                        //sample aperture diameter in mm
1448        return(diam)
1449end
1450
1451// change the sample aperture to a non-standard value
1452Function SampleApOtherSetVarProc(ctrlName,varNum,varStr,varName) : SetVariableControl
1453        String ctrlName
1454        Variable varNum
1455        String varStr
1456        String varName
1457               
1458        WAVE rw=root:Packages:NIST:SAS:realsRead
1459        rw[24] = varNum                 //sample aperture diameter in mm
1460        ReCalculateInten(1)
1461        return(0)
1462End
1463
1464//parses the control for A2 diam
1465// updates the wave and global
1466// returns a2 in cm
1467Function sampleApertureDiam()
1468        ControlInfo/W=SASCALC popup0_1
1469       
1470//      Print "In sampleApertureDiam()"
1471        //set the global
1472        NVAR a2=root:Packages:NIST:SAS:gSamAp
1473       
1474        if(cmpstr(S_Value,"other") == 0)                // "other" selected
1475                //enable the setvar, diameter in mm!
1476                SetVariable setvar0_3 disable=0
1477                // read its value (a global)
1478                NVAR a2other = root:Packages:NIST:SAS:gSamApOther
1479                a2=a2other/10                           //a2 in cm
1480        else
1481                SetVariable setvar0_3 disable=1
1482                //1st item is 1/16", popup steps by 1/16"
1483                a2 = 2.54/16.0 * (V_Value)                      //convert to cm         
1484        endif
1485        WAVE rw=root:Packages:NIST:SAS:realsRead
1486        rw[24] = a2*10                  //sample aperture diameter in mm
1487       
1488        return(a2)
1489end
1490
1491//1=huber 2=chamber
1492Function tableposition()
1493        ControlInfo/W=SASCALC checkHuber
1494        if(V_Value)
1495                return(1)               //huber
1496        else
1497                return(2)               //chamber
1498        endif
1499End
1500
1501//compute SSD and update both the global and the wave
1502//
1503Function sourceToSampleDist()
1504
1505        NVAR NG=root:Packages:NIST:SAS:gNg
1506        NVAR S12 = root:Packages:NIST:SAS:S12
1507        NVAR L2Diff = root:Packages:NIST:SAS:L2Diff
1508        NVAR SSD = root:Packages:NIST:SAS:gSSD
1509       
1510        SSD = 1632 - 155*NG - s12*(2-tableposition()) - L2Diff
1511       
1512        WAVE rw=root:Packages:NIST:SAS:realsRead
1513        rw[25] = SSD/100                // in meters
1514        return(SSD)
1515End
1516
1517//returns the offset value
1518// slider and setVar are linked to the same global
1519// updates the wave and changes the beamcenter (x,y) in the wave
1520Function detectorOffset()
1521       
1522        WAVE rw=root:Packages:NIST:SAS:RealsRead
1523        NVAR val = root:Packages:NIST:SAS:gOffset
1524        rw[19] = val            // already in cm
1525        //move the beamcenter
1526        rw[16] = 64 + 2*rw[19]          //approximate beam X is 64 w/no offset, 114 w/25 cm offset
1527        rw[17] = 64             //typical value
1528       
1529        return(val)
1530end
1531
1532//returns the detector distance (slider and setVar are linked by the global)
1533//
1534Function chamberToDetectorDist()
1535
1536        NVAR val = root:Packages:NIST:SAS:gDetDist
1537        return(val)
1538End
1539
1540//sets the SDD (slider and setVar are linked by the global and is the detector position
1541//  relative to the chamber)
1542// updates the wave
1543Function sampleToDetectorDist()
1544
1545        NVAR detDist = root:Packages:NIST:SAS:gDetDist
1546        NVAR S12 = root:Packages:NIST:SAS:S12
1547        WAVE rw=root:Packages:NIST:SAS:RealsRead       
1548        Variable SDD   
1549       
1550        SDD = detDist + s12*(2-tableposition())
1551        rw[18] = SDD/100                // convert to meters for header
1552        return(SDD)
1553End
1554
1555//direction = one of "vertical;horizontal;maximum;"
1556// all of this is bypassed if the lenses are in
1557//
1558Function beamDiameter(direction)
1559        String direction
1560
1561        NVAR lens = root:Packages:NIST:SAS:gUsingLenses
1562        if(lens)
1563                return sourceApertureDiam()
1564        endif
1565       
1566    Variable l1 = sourceToSampleDist()
1567    Variable l2 //= sampleAperToDetDist()
1568    Variable d1,d2,bh,bv,bm,umbra,a1,a2
1569   
1570    NVAR L2diff = root:Packages:NIST:SAS:L2diff
1571    NVAR lambda = root:Packages:NIST:SAS:gLambda
1572    NVAR lambda_width = root:Packages:NIST:SAS:gDeltaLambda
1573    NVAR bs_factor = root:Packages:NIST:SAS:bs_factor
1574   
1575    l2 = sampleToDetectorDist() + L2diff
1576    a1 = sourceApertureDiam()
1577    a2 = sampleApertureDiam()
1578   
1579    d1 = a1*l2/l1
1580    d2 = a2*(l1+l2)/l1
1581    bh = d1+d2          //beam size in horizontal direction
1582    umbra = abs(d1-d2)
1583    //vertical spreading due to gravity
1584    bv = bh + 1.25e-8*(l1+l2)*l2*lambda*lambda*lambda_width
1585    bm = (bs_factor*bh > bv) ? bs_factor*bh : bv //use the larger of horiz*safety or vertical
1586   
1587    strswitch(direction)        // string switch
1588        case "vertical":                // execute if case matches expression
1589                return(bv)
1590                break                                           // exit from switch
1591        case "horizontal":              // execute if case matches expression
1592                return(bh)
1593                break
1594        case "maximum":         // execute if case matches expression
1595                return(bm)
1596                break
1597        default:                                                        // optional default expression executed
1598                return(bm)                                              // when no case matches
1599    endswitch
1600End
1601
1602//on NG3 and NG7, allowable sizes are 1,2,3,4" diameter
1603//will return values larger than 4.0*2.54 if a larger beam is needed
1604//
1605// - in an approximate way, account for lenses
1606Function beamstopDiam()
1607
1608        NVAR yesLens = root:Packages:NIST:SAS:gUsingLenses
1609        Variable bm=0
1610        Variable bs=0.0
1611   
1612        if(yesLens)
1613                //bm = sourceApertureDiam()             //ideal result, not needed
1614                bs = 1                                                          //force the diameter to 1"
1615        else
1616                bm = beamDiameter("maximum")
1617                do
1618                bs += 1
1619           while ( (bs*2.54 < bm) || (bs > 30.0))                       //30 = ridiculous limit to avoid inf loop
1620        endif
1621
1622        //update the wave
1623        WAVE rw=root:Packages:NIST:SAS:realsRead
1624        rw[21] = bs*25.4                //store the BS diameter in mm
1625       
1626    return (bs*2.54)            //return diameter in cm, not inches for txt
1627End
1628
1629//returns the projected diameter of the beamstop at the anode plane.
1630// most noticeable at short SDD
1631//if flag == 0 use conservative estimate = largest diameter (for SASCALC, default)
1632//if flag != 0 use point aperture = average diameter (for resolution calculation)
1633Function beamstopDiamProjection(flag)
1634        Variable flag
1635       
1636        NVAR L2diff = root:Packages:NIST:SAS:L2diff
1637        Variable a2 = sampleApertureDiam()
1638        Variable bs = beamstopDiam()
1639        Variable l2, LB, BS_P
1640   
1641        l2 = sampleToDetectorDist() + L2diff
1642        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical)
1643        if(flag==0)
1644                BS_P = bs + (bs+a2)*lb/(l2-lb)          //diameter of shadow from parallax
1645        else
1646                BS_P = bs + bs*lb/(l2-lb)               //diameter of shadow, point A2
1647        endif
1648        return (bs_p)           //return projected diameter in cm
1649End
1650
1651// 19MAR07 - using correction from John for an estimate of the shadow of the beamstop
1652// at the detection plane. This is a noticeable effect at short SDD, where the projected
1653// diameter of the beamstop is much larger than the physical diameter.
1654Function qMin()
1655
1656    Variable l2s = sampleToDetectorDist()       //distance from sample to detector in cm
1657//    Variable bs = beamstopDiam()              //beamstop diameter in cm
1658    Variable bs_p = beamstopDiamProjection(0)           //projected beamstop diameter in cm
1659    NVAR lambda = root:Packages:NIST:SAS:gLambda
1660    NVAR d_det = root:Packages:NIST:SAS:d_det           //cm
1661    NVAR a_pixel = root:Packages:NIST:SAS:a_pixel       //cm
1662
1663    return( (pi/lambda)*(bs_p + d_det + a_pixel)/l2s )          //use bs_p rather than bs
1664   // return( (pi/lambda)*(bs + d_det + a_pixel)/l2s )          //use bs (incorrect)
1665End
1666
1667Function qMaxVert()
1668
1669    Variable theta
1670    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
1671    NVAR lambda = root:Packages:NIST:SAS:gLambda
1672        NVAR det_width = root:Packages:NIST:SAS:det_width
1673       
1674    theta = atan( (det_width/2.0)/l2s )
1675   
1676    return ( 4.0*pi/lambda * sin(theta/2.0) )
1677end
1678
1679Function qMaxCorner()
1680
1681    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
1682    Variable radial
1683    NVAR lambda = root:Packages:NIST:SAS:gLambda
1684        NVAR det_off = root:Packages:NIST:SAS:gOffset
1685        NVAR det_width = root:Packages:NIST:SAS:det_width
1686
1687    radial=sqrt( (0.5*det_width)*(0.5*det_width)+(0.5*det_width+det_off)*(0.5*det_width+det_off) )
1688   
1689    return ( 4*pi/lambda*sin(0.5*atan(radial/l2s)) )
1690End
1691
1692Function qMaxHoriz()
1693
1694    Variable theta
1695    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
1696    NVAR lambda = root:Packages:NIST:SAS:gLambda
1697        NVAR det_off = root:Packages:NIST:SAS:gOffset
1698        NVAR det_width = root:Packages:NIST:SAS:det_width
1699
1700    theta = atan( ((det_width/2.0)+det_off)/l2s )       //from the instance variables
1701   
1702    return ( 4.0*pi/lambda * sin(theta/2.0) )
1703End
1704
1705// calculate sigma for the resolution function at either limit of q-range
1706Function deltaQ(atQ)
1707        Variable atQ
1708       
1709    Variable k02,lp,l1,l2,sig_02,sigQ2,a1,a2
1710    NVAR l2Diff = root:Packages:NIST:SAS:L2diff
1711        NVAR lambda = root:Packages:NIST:SAS:gLambda
1712        NVAR lambda_width = root:Packages:NIST:SAS:gDeltaLambda
1713        NVAR d_det = root:Packages:NIST:SAS:d_det
1714        NVAR del_r = root:Packages:NIST:SAS:del_r
1715       
1716       
1717    l1 = sourceToSampleDist()
1718    l2 = sampleToDetectorDist() + L2diff
1719    a1 = sourceApertureDiam()
1720    a2 = sampleApertureDiam()
1721   
1722    k02 = (6.2832/lambda)*(6.2832/lambda)
1723    lp = 1/(1/l1 + 1/l2)
1724   
1725    sig_02 = (0.25*a1/l1)*(0.25*a1/l1)
1726    sig_02 += (0.25*a2/lp)*(0.25*a2/lp)
1727    sig_02 += (d_det/(2.355*l2))*(d_det/(2.355*l2))
1728    sig_02 += (del_r/l2)*(del_r/l2)/12
1729    sig_02 *= k02
1730   
1731    sigQ2 = sig_02 + (atQ*lambda_width)*(atQ*lambda_width)/6
1732
1733    return(100*sqrt(sigQ2)/atQ)
1734End
1735
1736
1737Function beamIntensity()
1738
1739    Variable alpha,f,t,t4,t5,t6,as,solid_angle,l1,d2_phi
1740    Variable a1,a2,retVal
1741    SetDataFolder root:Packages:NIST:SAS
1742    NVAR l_gap=l_gap,guide_width =guide_width,ng = gNg
1743    NVAR lambda_t=lambda_t,b=b,c=c
1744    NVAR lambda=gLambda,t1=t1,t2=t2,t3=t3,phi_0=phi_0
1745    NVAR lambda_width=gDeltaLambda
1746   
1747    l1 = sourceToSampleDist()
1748    a1 = sourceApertureDiam()
1749    a2 = sampleApertureDiam()
1750   
1751   
1752    alpha = (a1+a2)/(2*l1)      //angular divergence of beam
1753    f = l_gap*alpha/(2*guide_width)
1754    t4 = (1-f)*(1-f)
1755    t5 = exp(ng*ln(0.96))       // trans losses of guides in pre-sample flight
1756    t6 = 1 - lambda*(b-(ng/8)*(b-c))            //experimental correction factor
1757    t = t1*t2*t3*t4*t5*t6
1758   
1759    as = pi/4*a2*a2             //area of sample in the beam
1760    d2_phi = phi_0/(2*pi)
1761    d2_phi *= exp(4*ln(lambda_t/lambda))
1762    d2_phi *= exp(-1*(lambda_t*lambda_t/lambda/lambda))
1763
1764    solid_angle = pi/4* (a1/l1)*(a1/l1)
1765
1766    retVal = as * d2_phi * lambda_width * solid_angle * t
1767     SetDataFolder root:
1768    return (retVal)
1769end
1770
1771Function figureOfMerit()
1772
1773        Variable bi = beamIntensity()
1774        NVAR lambda = root:Packages:NIST:SAS:gLambda
1775       
1776    return (lambda*lambda*bi)
1777End
1778
1779//estimate the number of pixels in the beam, and enforce the maximum countrate per pixel (idmax)
1780Function attenuatorTransmission()
1781
1782    Variable num_pixels,i_pix           //i_pix = id in John's notation
1783    Variable bDiam = beamDiameter("horizontal") //!! note that prev calculations used bh (horizontal)
1784    Variable atten,a2
1785    SetDataFolder root:Packages:NIST:SAS
1786    NVAR a_pixel=a_pixel,idmax=idmax
1787   
1788    a2 = sampleApertureDiam()
1789   
1790    num_pixels = pi/4*(0.5*(a2+bDiam))*(0.5*(a2+bDiam))/a_pixel/a_pixel
1791    i_pix = ( beamIntensity() )/num_pixels
1792   
1793    atten = (i_pix < idmax) ? 1.0 : idmax/i_pix
1794    SetDataFolder root:
1795    return(atten)
1796End
1797
1798Function attenuatorNumber()
1799
1800    Variable atten = attenuatorTransmission()
1801    Variable af,nf,numAtten
1802    SetDataFolder root:Packages:NIST:SAS
1803    NVAR lambda=gLambda
1804   
1805    af = 0.498 + 0.0792*lambda - 1.66e-3*lambda*lambda
1806    nf = -ln(atten)/af          //floating point
1807   
1808    numAtten = trunc(nf) + 1                    //in c, (int)nf
1809    //correct for larger step thickness at n > 6
1810    if(numAtten > 6)
1811        numAtten = 7 + trunc( (numAtten-6)/2 )          //in c, numAtten = 7 + (int)( (numAtten-6)/2 )
1812    endif
1813   
1814    SetDatafolder root:
1815    return (numAtten)
1816End
Note: See TracBrowser for help on using the repository browser.