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

Last change on this file since 139 was 139, checked in by srkline, 16 years ago

Added the option for an arbitrary diameter of sample aperture

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