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

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

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