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

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

More tweaking of the MonteCarlo?, trying in vain to thread the XOP. As is, works fine with the XOP and single threading.

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