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

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

added a more realistic shadow correction for resolution

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