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

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

A variety of fixes to make some procedures compatible with the new Packages:NIST subfolder structure

Also, some additional changes to include the proper procedures when loading overall packages

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