source: sans/SASCalc/trunk/SASCALC.ipf @ 233

Last change on this file since 233 was 233, checked in by srkline, 15 years ago

removed thumb color from the sliders to get the native GUI appearance

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