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

Last change on this file since 52 was 52, checked in by srkline, 16 years ago
  • changed pxt to include the procedure file
  • added freezing of configuration text for printing (RT#713)
  • colorized frozen traces
File size: 43.2 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//
11// calculate what q-values you get based on the instruments settings
12// that you would typically input to SASCALC
13// calculation is for (80A radius spheres) * (beam stop shadow factor)
14// or a Debye function (Rg=50A) * (beam stop shadow factor)
15// - NOT true intensity, not counts, just a display
16//
17// To Do:
18// - add in instrument conditions for lens/(lens+prism) configurations
19// - proper resolution calculation for lens/prism
20// - don't use the corners... misleading high Q limit (just max horiz?) - make a default mask
21
22//
23// Optional:
24// - freeze configurations with a user defined tag
25// - different model functions (+ change simple parameters)
26// - resolution smeared models
27// - "simulation" of data and error bars given a model and a total number of detector counts
28// - sequence of colors for frozen plots
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=0
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        endif
651        //append the text
652        sprintf str,"\rConfiguration #%d\r",ct-1
653        Notebook Saved_Configurations showRuler=0,defaultTab=20,selection={endOfFile, endOfFile}
654        Notebook Saved_Configurations font="Monaco",fSize=10,fstyle=1,text=str          //bold
655        Notebook Saved_Configurations font="Monaco",fSize=10,fstyle=0,text=SetConfigurationText()
656        return(0)
657End
658
659//clears the frozen traces on the graph, asks if you want to clear the saved text
660//
661Function ClearButtonProc(ctrlName) : ButtonControl
662        String ctrlName
663
664        NVAR ct=root:SAS:gFreezeCount
665        Variable ii
666        Setdatafolder root:SAS
667        for(ii=ct-1;ii>=0;ii-=1)
668        //remove all traces, replace aveint
669        // kill all waves _ct
670                RemoveFromGraph $("aveint_"+num2str(ii))
671                Killwaves/Z $("aveint_"+num2str(ii)),$("qval_"+num2str(ii))
672        endfor
673        ct=0
674        setdatafolder root:
675       
676        DoAlert 1,"Do you also want to clear the \"Saved Configurations\" text?"
677        if(V_flag == 1)         // yes
678                DoWindow/K Saved_Configurations
679        endif
680       
681        return(0)
682End
683
684
685///////////////////////////////////////////////////////////
686Function S_CircularAverageTo1D(type)
687        String type
688       
689        Variable isCircular = 1
690       
691        //type is the data type to do the averaging on, and will be set as the current folder
692        //get the current displayed data (so the correct folder is used)
693        String destPath = "root:"+type
694       
695        //
696        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist,dr,ddr
697        Variable lambda,trans
698        WAVE reals = $(destPath + ":RealsRead")
699//      WAVE/T textread = $(destPath + ":TextRead")
700//      String fileStr = textread[3]
701       
702        // center of detector, for non-linear corrections
703//      NVAR pixelsX = root:myGlobals:gNPixelsX
704//      NVAR pixelsY = root:myGlobals:gNPixelsY
705        Variable pixelsX=128,pixelsY=128
706       
707        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
708        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
709       
710        // beam center, in pixels
711        x0 = reals[16]
712        y0 = reals[17]
713        //detector calibration constants
714        sx = reals[10]          //mm/pixel (x)
715        sx3 = reals[11]         //nonlinear coeff
716        sy = reals[13]          //mm/pixel (y)
717        sy3 = reals[14]         //nonlinear coeff
718       
719        dtsize = 10*reals[20]           //det size in mm
720        dtdist = 1000*reals[18]         // det distance in mm
721       
722        NVAR binWidth=root:myGlobals:gBinWidth
723//      Variable binWidth = 1
724       
725        dr = binWidth           // ***********annulus width set by user, default is one***********
726        ddr = dr*sx             //step size, in mm (this value should be passed to the resolution calculation, not dr 18NOV03)
727               
728        Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,dxbm,dybm,ii
729        Variable phi_rad,dphi_rad,phi_x,phi_y
730        Variable forward,mirror
731       
732        String side = "both"
733//      side = StringByKey("SIDE",keyListStr,"=",";")
734//      Print "side = ",side
735       
736        /// data wave is data in the current folder which was set at the top of the function
737        WAVE data=$(destPath + ":data")
738
739// fake mask that uses all of the detector
740        Make/O/N=(pixelsX,pixelsY) $(destPath + ":mask")
741        Wave mask = $(destPath + ":mask")
742        mask = 0
743        //two pixels all around
744        mask[0,1][] = 1
745        mask[126,127][] = 1
746        mask[][0,1] = 1
747        mask[][126,127] = 1
748        //
749        //pixels within rcentr of beam center are broken into 9 parts (units of mm)
750        rcentr = 100            //original
751//      rcentr = 0
752        // values for error if unable to estimate value
753        //large_num = 1e10
754        large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1
755        small_num = 1e-10
756       
757        // output wave are expected to exist (?) initialized to zero, what length?
758        // 200 points on VAX --- use 300 here, or more if SAXS data is used with 1024x1024 detector (1000 pts seems good)
759        Variable defWavePts=500
760        Make/O/N=(defWavePts) $(destPath + ":qval"),$(destPath + ":aveint")
761        Make/O/N=(defWavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave")
762        Make/O/N=(defWavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar")
763
764        WAVE qval = $(destPath + ":qval")
765        WAVE aveint = $(destPath + ":aveint")
766        WAVE ncells = $(destPath + ":ncells")
767        WAVE dsq = $(destPath + ":dsq")
768        WAVE sigave = $(destPath + ":sigave")
769        WAVE qbar = $(destPath + ":QBar")
770        WAVE sigmaq = $(destPath + ":SigmaQ")
771        WAVE fsubs = $(destPath + ":fSubS")
772       
773        qval = 0
774        aveint = 0
775        ncells = 0
776        dsq = 0
777        sigave = 0
778        qbar = 0
779        sigmaq = 0
780        fsubs = 0
781
782        dtdis2 = dtdist^2
783        nq = 1
784        xoffst=0
785        //distance of beam center from detector center
786        dxbm = S_FX(x0,sx3,xcenter,sx)
787        dybm = S_FY(y0,sy3,ycenter,sy)
788               
789        //BEGIN AVERAGE **********
790        Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1
791        Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p
792       
793        // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
794        // loop index corresponds to FORTRAN (old code)
795        // and the IGOR array indices must be adjusted (-1) to the correct address
796        ii=1
797        do
798                xi = ii
799                dxi = S_FX(xi,sx3,xcenter,sx)
800                dx = dxi-dxbm           //dx and dy are in mm
801               
802                jj = 1
803                do
804                        data_pixel = data[ii-1][jj-1]           //assign to local variable
805                        yj = jj
806                        dyj = S_FY(yj,sy3,ycenter,sy)
807                        dy = dyj - dybm
808                        if(!(mask[ii-1][jj-1]))                 //masked pixels = 1, skip if masked (this way works...)
809                                dr2 = (dx^2 + dy^2)^(0.5)               //distance from beam center NOTE dr2 used here - dr used above
810                                if(dr2>rcentr)          //keep pixel whole
811                                        nd = 1
812                                        fd = 1
813                                else                            //break pixel into 9 equal parts
814                                        nd = 3
815                                        fd = 2
816                                endif
817                                nd2 = nd^2
818                                ll = 1          //"el-el" loop index
819                                do
820                                        dxx = dx + (ll - fd)*sx/3
821                                        kk = 1
822                                        do
823                                                dyy = dy + (kk - fd)*sy/3
824                                                if(isCircular)
825                                                        //circular average, use all pixels
826                                                        //(increment)
827                                                        nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
828                                                else
829                                                        //a sector average - determine azimuthal angle
830                                                        dphi_p = S_dphi_pixel(dxx,dyy,phi_x,phi_y)
831                                                        if(dphi_p < dphi_rad)
832                                                                forward = 1                     //within forward sector
833                                                        else
834                                                                forward = 0
835                                                        Endif
836                                                        if((Pi - dphi_p) < dphi_rad)
837                                                                mirror = 1              //within mirror sector
838                                                        else
839                                                                mirror = 0
840                                                        Endif
841                                                        //check if pixel lies within allowed sector(s)
842                                                        if(cmpstr(side,"both")==0)              //both sectors
843                                                                if ( mirror || forward)
844                                                                        //increment
845                                                                        nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
846                                                                Endif
847                                                        else
848                                                                if(cmpstr(side,"right")==0)             //forward sector only
849                                                                        if(forward)
850                                                                                //increment
851                                                                                nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
852                                                                        Endif
853                                                                else                    //mirror sector only
854                                                                        if(mirror)
855                                                                                //increment
856                                                                                nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
857                                                                        Endif
858                                                                Endif
859                                                        Endif           //allowable sectors
860                                                Endif   //circular or sector check
861                                                kk+=1
862                                        while(kk<=nd)
863                                        ll += 1
864                                while(ll<=nd)
865                        Endif           //masked pixel check
866                        jj += 1
867                while (jj<=pixelsY)
868                ii += 1
869        while(ii<=pixelsX)              //end of the averaging
870               
871        //compute q-values and errors
872        Variable ntotal,rr,theta,avesq,aveisq,var
873       
874        lambda = reals[26]
875        ntotal = 0
876        kk = 1
877        do
878                rr = (2*kk-1)*ddr/2
879                theta = 0.5*atan(rr/dtdist)
880                qval[kk-1] = (4*Pi/lambda)*sin(theta)
881                if(ncells[kk-1] == 0)
882                        //no pixels in annuli, data unknown
883                        aveint[kk-1] = 0
884                        sigave[kk-1] = large_num
885                else
886                        if(ncells[kk-1] <= 1)
887                                //need more than one pixel to determine error
888                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
889                                sigave[kk-1] = large_num
890                        else
891                                //assume that the intensity in each pixel in annuli is normally
892                                // distributed about mean...
893                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
894                                avesq = aveint[kk-1]^2
895                                aveisq = dsq[kk-1]/ncells[kk-1]
896                                var = aveisq-avesq
897                                if(var<=0)
898                                        sigave[kk-1] = small_num
899                                else
900                                        sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1))
901                                endif
902                        endif
903                endif
904                ntotal += ncells[kk-1]
905                kk+=1
906        while(kk<=nq)
907       
908        //Print "NQ = ",nq
909        // data waves were defined as 300 points (=defWavePts), but now have less than that (nq) points
910        // use DeletePoints to remove junk from end of waves
911        //WaveStats would be a more foolproof implementation, to get the # points in the wave
912        Variable startElement,numElements
913        startElement = nq
914        numElements = defWavePts - startElement
915        DeletePoints startElement,numElements, qval,aveint,ncells,dsq,sigave
916       
917        //////////////end of VAX sector_ave()
918               
919
920
921// ***************************************************************
922//
923// Do the extra 3 columns of resolution calculations starting here.
924//
925// ***************************************************************
926
927        Variable L2 = reals[18]
928        Variable BS = reals[21]
929        Variable S1 = reals[23]
930        Variable S2 = reals[24]
931        Variable L1 = reals[25]
932        lambda = reals[26]
933        Variable lambdaWidth = reals[27]
934
935        Variable DDet, apOff=0.0
936        // hard wire value for Ordela detectors
937        DDet = 0.5              // resoltion in cm
938        //      String detStr=textRead[9]
939        //      DDet = DetectorPixelResolution(fileStr,detStr)          //needs detector type and beamline
940
941        //Go from 0 to nq doing the calc for all three values at
942        //every Q value
943
944        ii=0
945//      String res_string="root:myGlobals:Res_vals"
946//      Make/O/D/N=3 $res_string
947//      Wave res_wave=$res_string
948        Variable ret1,ret2,ret3
949        do
950                S_getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,ret1,ret2,ret3)
951                sigmaq[ii] = ret1       //res_wave[0]
952                qbar[ii] = ret2         //res_wave[1]
953                fsubs[ii] = ret3                //res_wave[2]
954                ii+=1
955        while(ii<nq)
956        DeletePoints startElement,numElements, sigmaq, qbar, fsubs
957
958        fsubs += 1e-8           //keep the values from being too small
959
960// End of resolution calculations
961// ***************************************************************
962       
963        //get rid of the default mask, if one was created (it is in the current folder)
964        //don't just kill "mask" since it might be pointing to the one in the MSK folder
965        Killwaves/Z $(destPath+":mask")
966       
967        //return to root folder (redundant)
968        SetDataFolder root:
969       
970        Return 0
971End
972
973//returns nq, new number of q-values
974//arrays aveint,dsq,ncells are also changed by this function
975//
976Function S_IncrementPixel(dataPixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
977        Variable dataPixel,ddr,dxx,dyy
978        Wave aveint,dsq,ncells
979        Variable nq,nd2
980       
981        Variable ir
982       
983        ir = trunc(sqrt(dxx*dxx+dyy*dyy)/ddr)+1
984        if (ir>nq)
985                nq = ir         //resets maximum number of q-values
986        endif
987        aveint[ir-1] += dataPixel/nd2           //ir-1 must be used, since ir is physical
988        dsq[ir-1] += dataPixel*dataPixel/nd2
989        ncells[ir-1] += 1/nd2
990       
991        Return nq
992End
993
994//function determines azimuthal angle dphi that a vector connecting
995//center of detector to pixel makes with respect to vector
996//at chosen azimuthal angle phi -> [cos(phi),sin(phi)] = [phi_x,phi_y]
997//dphi is always positive, varying from 0 to Pi
998//
999Function S_dphi_pixel(dxx,dyy,phi_x,phi_y)
1000        Variable dxx,dyy,phi_x,phi_y
1001       
1002        Variable val,rr,dot_prod
1003       
1004        rr = sqrt(dxx^2 + dyy^2)
1005        dot_prod = (dxx*phi_x + dyy*phi_y)/rr
1006        //? correct for roundoff error? - is this necessary in IGOR, w/ double precision?
1007        if(dot_prod > 1)
1008                dot_prod =1
1009        Endif
1010        if(dot_prod < -1)
1011                dot_prod = -1
1012        Endif
1013       
1014        val = acos(dot_prod)
1015       
1016        return val
1017
1018End
1019
1020//calculates the x distance from the center of the detector, w/nonlinear corrections
1021//
1022Function S_FX(xx,sx3,xcenter,sx)               
1023        Variable xx,sx3,xcenter,sx
1024       
1025        Variable retval
1026       
1027        retval = sx3*tan((xx-xcenter)*sx/sx3)
1028        Return retval
1029End
1030
1031//calculates the y distance from the center of the detector, w/nonlinear corrections
1032//
1033Function S_FY(yy,sy3,ycenter,sy)               
1034        Variable yy,sy3,ycenter,sy
1035       
1036        Variable retval
1037       
1038        retval = sy3*tan((yy-ycenter)*sy/sy3)
1039        Return retval
1040End
1041
1042//**********************
1043// Resolution calculation - used by the averaging routines
1044// to calculate the resolution function at each q-value
1045// - the return value is not used
1046//
1047// equivalent to John's routine on the VAX Q_SIGMA_AVE.FOR
1048// Incorporates eqn. 3-15 from J. Appl. Cryst. (1995) v. 28 p105-114
1049//
1050Function/S S_getResolution(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,SigmaQ,QBar,fSubS)
1051        Variable inQ, lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r
1052        Variable &fSubS, &QBar, &SigmaQ         //these are the output quantities at the input Q value
1053       
1054        //lots of calculation variables
1055        Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g
1056        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
1057
1058        //Constants
1059        //Variable del_r = .1
1060        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
1061        Variable g = 981.0                              //gravity acceleration [cm/s^2]
1062
1063        String results
1064        results ="Failure"
1065
1066        //rename for working variables,  these must be gotten from global
1067        //variables
1068
1069//      Variable wLam, wLW, wL1, wL2, wS1, wS2
1070//      Variable wBS, wDDet, wApOff
1071//      wLam = lambda
1072//      wLW = lambdaWidth
1073//      wDDet = DDet
1074//      wApOff = apOff
1075        S1 *= 0.5*0.1                   //convert to radius and [cm]
1076        S2 *= 0.5*0.1
1077
1078        L1 *= 100.0                     // [cm]
1079        L1 -= apOff                             //correct the distance
1080
1081        L2 *= 100.0
1082        L2 += apOff
1083
1084        BS *= 0.5*0.1                   //convert to radius and [cm]
1085        del_r *= 0.1                            //width of annulus, convert mm to [cm]
1086       
1087        //Start resolution calculation
1088        a2 = S1*L2/L1 + S2*(L1+L2)/L1
1089        q_small = 2.0*Pi*(BS-a2)*(1.0-lambdaWidth)/(lambda*L2)
1090        lp = 1.0/( 1.0/L1 + 1.0/L2)
1091
1092        v_lambda = lambdaWidth^2/6.0
1093        v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2
1094        v_d = (DDet/2.3548)^2 + del_r^2/12.0
1095        vz = vz_1 / lambda
1096        yg = 0.5*g*L2*(L1+L2)/vz^2
1097        v_g = 2.0*yg^2*v_lambda
1098
1099        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
1100        delta = 0.5*(BS - r0)^2/v_d
1101
1102        if (r0 < BS)
1103                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
1104        else
1105                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
1106        endif
1107
1108        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
1109        if (fSubS <= 0.0)
1110                fSubS = 1.e-10
1111        endif
1112        fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
1113        fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
1114
1115        rmd = fr*r0
1116        v_r1 = v_b + fv*v_d +v_g
1117
1118        rm = rmd + 0.5*v_r1/rmd
1119        v_r = v_r1 - 0.5*(v_r1/rmd)^2
1120        if (v_r < 0.0)
1121                v_r = 0.0
1122        endif
1123        QBar = (4.0*Pi/lambda)*sin(0.5*atan(rm/L2))
1124        SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda)
1125
1126        results = "success"
1127        Return results
1128End
1129
1130Function S_Debye(scale,rg,bkg,x)
1131        Variable scale,rg,bkg
1132        Variable x
1133       
1134        // variables are:
1135        //[0] scale factor
1136        //[1] radius of gyration []
1137        //[2] background        [cm-1]
1138       
1139        // calculates (scale*debye)+bkg
1140        Variable Pq,qr2
1141       
1142        qr2=(x*rg)^2
1143        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
1144       
1145        //scale
1146        Pq *= scale
1147        // then add in the background
1148        return (Pq+bkg)
1149End
1150
1151
1152Function S_SphereForm(scale,radius,delrho,bkg,x)                               
1153        Variable scale,radius,delrho,bkg
1154        Variable x
1155       
1156        // variables are:                                                       
1157        //[0] scale
1158        //[1] radius ()
1159        //[2] delrho (-2)
1160        //[3] background (cm-1)
1161       
1162//      Variable scale,radius,delrho,bkg                               
1163//      scale = w[0]
1164//      radius = w[1]
1165//      delrho = w[2]
1166//      bkg = w[3]
1167       
1168       
1169        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
1170        // and is rescaled to give [=] cm^-1
1171       
1172        Variable bes,f,vol,f2
1173        //
1174        //handle q==0 separately
1175        If(x==0)
1176                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
1177                return(f)
1178        Endif
1179       
1180        bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
1181        vol = 4*pi/3*radius^3
1182        f = vol*bes*delrho              // [=] 
1183        // normalize to single particle volume, convert to 1/cm
1184        f2 = f * f / vol * 1.0e8                // [=] 1/cm
1185       
1186        return (scale*f2+bkg)   // Scale, then add in the background
1187       
1188End
1189
1190Function/S SetConfigurationText()
1191
1192        String str="",temp
1193
1194        SetDataFolder root:SAS
1195       
1196        NVAR numberOfGuides=gNg
1197        NVAR gTable=gTable              //2=chamber, 1=table
1198        NVAR wavelength=gLambda
1199        NVAR lambdaWidth=gDeltaLambda
1200        NVAR instrument = instrument
1201        NVAR L2diff = L2diff
1202       
1203       
1204        sprintf temp,"Source Aperture Diameter =\t\t%6.2f cm\r",sourceApertureDiam()
1205        str += temp
1206        sprintf temp,"Source to Sample =\t\t\t\t%6.0f cm\r",sourceToSampleDist()
1207        str += temp
1208        sprintf temp,"Sample Aperture to Detector =\t%6.0f cm\r",sampleToDetectorDist() + L2diff
1209        str += temp
1210        sprintf temp,"Beam diameter =\t\t\t\t\t%6.2f cm\r",beamDiameter("maximum")
1211        str += temp
1212        sprintf temp,"Beamstop diameter =\t\t\t\t%6.2f inches\r",beamstopDiam()/2.54
1213        str += temp
1214        sprintf temp,"Minimum Q-value =\t\t\t\t%6.4f 1/ (sigQ/Q = %3.1f %%)\r",qMin(),deltaQ(qMin())
1215        str += temp
1216        sprintf temp,"Maximum Horizontal Q-value =\t%6.4f 1/\r",qMaxHoriz()
1217        str += temp
1218        sprintf temp,"Maximum Vertical Q-value =\t\t%6.4f 1/\r",qMaxVert()
1219        str += temp
1220        sprintf temp,"Maximum Q-value =\t\t\t\t%6.4f 1/ (sigQ/Q = %3.1f %%)\r",qMaxCorner(),deltaQ(qMaxCorner())
1221        str += temp
1222        sprintf temp,"Beam Intensity =\t\t\t\t%.0f counts/s\r",beamIntensity()
1223        str += temp
1224        sprintf temp,"Figure of Merit =\t\t\t\t%3.3g ^2/s\r",figureOfMerit()
1225        str += temp
1226        sprintf temp,"Attenuator transmission =\t\t%3.3g = Atten # %d\r",attenuatorTransmission(),attenuatorNumber()
1227        str += temp
1228//     
1229//      // add text of the user-edited values
1230//      //
1231        sprintf temp,"***************** NG %d *****************\r",instrument
1232        str += temp
1233        sprintf temp,"Sample Aperture Diameter =\t\t\t\t%.2f cm\r",sampleApertureDiam()
1234        str += temp
1235        sprintf temp,"Number of Guides =\t\t\t\t\t\t%d \r", numberOfGuides
1236        str += temp
1237        sprintf temp,"Sample Chamber to Detector =\t\t\t%.1f cm\r", chamberToDetectorDist()
1238        str += temp
1239        if(gTable==1)
1240                sprintf temp,"Sample Position is \t\t\t\t\t\tHuber\r"
1241        else
1242                sprintf temp,"Sample Position is \t\t\t\t\t\tChamber\r"
1243        endif
1244        str += temp
1245        sprintf temp,"Detector Offset =\t\t\t\t\t\t%.1f cm\r", detectorOffset()
1246        str += temp
1247        sprintf temp,"Neutron Wavelength =\t\t\t\t\t%.2f \r", wavelength
1248        str += temp
1249        sprintf temp,"Wavelength Spread, FWHM =\t\t\t\t%.3f\r", lambdaWidth
1250        str += temp
1251        sprintf temp,"Sample Aperture to Sample Position =\t%.2f cm\r", L2Diff
1252        str += temp
1253       
1254   setDataFolder root:
1255   return str                   
1256End
1257
1258Function DisplayConfigurationText()
1259
1260        if(WinType("Trial_Configuration")==0)
1261                NewNotebook/F=0/K=1/N=Trial_Configuration /W=(480,44,880,369)
1262        endif
1263        //replace the text
1264        Notebook Trial_Configuration selection={startOfFile, endOfFile}
1265        Notebook Trial_Configuration font="Monaco",fSize=10,text=SetConfigurationText()
1266        return(0)
1267end
1268
1269//parses the control for A1 diam
1270// updates the wave
1271Function sourceApertureDiam()
1272        ControlInfo popup0
1273        Variable diam
1274        sscanf S_Value, "%g cm", diam
1275        WAVE rw=root:SAS:realsRead
1276        rw[23] = diam*10                        //sample aperture diameter in mm
1277        return(diam)
1278end
1279
1280//parses the control for A2 diam
1281// updates the wave and global
1282// returns a2 in cm
1283Function sampleApertureDiam()
1284        ControlInfo popup0_1
1285       
1286        //set the global
1287        NVAR a2=root:SAS:gSamAp
1288        //1st item is 1/16", popup steps by 1/16"
1289        a2 = 2.54/16.0 * (V_Value)
1290       
1291        WAVE rw=root:SAS:realsRead
1292        rw[24] = a2*10                  //sample aperture diameter in mm
1293       
1294        return(a2)
1295end
1296
1297//1=huber 2=chamber
1298Function tableposition()
1299        ControlInfo checkHuber
1300        if(V_Value)
1301                return(1)               //huber
1302        else
1303                return(2)               //chamber
1304        endif
1305End
1306
1307//compute SSD and update both the global and the wave
1308//
1309Function sourceToSampleDist()
1310
1311        NVAR NG=root:SAS:gNg
1312        NVAR S12 = root:SAS:S12
1313        NVAR L2Diff = root:SAS:L2Diff
1314        NVAR SSD = root:SAS:gSSD
1315       
1316        SSD = 1632 - 155*NG - s12*(2-tableposition()) - L2Diff
1317       
1318        WAVE rw=root:SAS:realsRead
1319        rw[25] = SSD/100                // in meters
1320        return(SSD)
1321End
1322
1323//returns the offset value
1324// slider and setVar are linked to the same global
1325// updates the wave and changes the beamcenter (x,y) in the wave
1326Function detectorOffset()
1327       
1328        WAVE rw=root:SAS:RealsRead
1329        NVAR val = root:SAS:gOffset
1330        rw[19] = val            // already in cm
1331        //move the beamcenter
1332        rw[16] = 64 + 2*rw[19]          //approximate beam X is 64 w/no offset, 114 w/25 cm offset
1333        rw[17] = 64             //typical value
1334       
1335        return(val)
1336end
1337
1338//returns the detector distance (slider and setVar are linked by the global)
1339//
1340Function chamberToDetectorDist()
1341
1342        NVAR val = root:SAS:gDetDist
1343        return(val)
1344End
1345
1346//sets the SDD (slider and setVar are linked by the global and is the detector position
1347//  relative to the chamber)
1348// updates the wave
1349Function sampleToDetectorDist()
1350
1351        NVAR detDist = root:SAS:gDetDist
1352        NVAR S12 = root:SAS:S12
1353        WAVE rw=root:SAS:RealsRead     
1354        Variable SDD   
1355       
1356        SDD = detDist + s12*(2-tableposition())
1357        rw[18] = SDD/100                // convert to meters
1358        return(SDD)
1359End
1360
1361//direction = one of "vertical;horizontal;maximum;"
1362Function beamDiameter(direction)
1363        String direction
1364
1365    Variable l1 = sourceToSampleDist()
1366    Variable l2 //= sampleAperToDetDist()
1367    Variable d1,d2,bh,bv,bm,umbra,a1,a2
1368   
1369    NVAR L2diff = root:SAS:L2diff
1370    NVAR lambda = root:SAS:gLambda
1371    NVAR lambda_width = root:SAS:gDeltaLambda
1372    NVAR bs_factor = root:SAS:bs_factor
1373   
1374    l2 = sampleToDetectorDist() + L2diff
1375    a1 = sourceApertureDiam()
1376    a2 = sampleApertureDiam()
1377   
1378    d1 = a1*l2/l1
1379    d2 = a2*(l1+l2)/l1
1380    bh = d1+d2          //beam size in horizontal direction
1381    umbra = abs(d1-d2)
1382    //vertical spreading due to gravity
1383    bv = bh + 1.25e-8*(l1+l2)*l2*lambda*lambda*lambda_width
1384    bm = (bs_factor*bh > bv) ? bs_factor*bh : bv //use the larger of horiz*safety or vertical
1385   
1386    strswitch(direction)        // string switch
1387        case "vertical":                // execute if case matches expression
1388                return(bv)
1389                break                                           // exit from switch
1390        case "horizontal":              // execute if case matches expression
1391                return(bh)
1392                break
1393        case "maximum":         // execute if case matches expression
1394                return(bm)
1395                break
1396        default:                                                        // optional default expression executed
1397                return(bm)                                              // when no case matches
1398    endswitch
1399End
1400
1401//on NG3 and NG7, allowable sizes are 1,2,3,4" diameter
1402//will return values larger than 4.0*2.54 if a larger beam is needed
1403Function beamstopDiam()
1404
1405    Variable bm = beamDiameter("maximum")
1406    Variable bs=0.0
1407   
1408    do
1409        bs += 1
1410    while ( (bs*2.54 < bm) || (bs > 30.0))                      //30 = ridiculous limit to avoid inf loop
1411 
1412        //update the wave
1413        WAVE rw=root:SAS:realsRead
1414        rw[21] = bs*25.4                //store the BS diameter in mm
1415       
1416    return (bs*2.54)            //return diameter in cm, not inches for txt
1417End
1418
1419
1420Function qMin()
1421
1422    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
1423    Variable bs = beamstopDiam()                //beamstop diameter
1424    NVAR lambda = root:SAS:gLambda
1425    NVAR d_det = root:SAS:d_det
1426    NVAR a_pixel = root:SAS:a_pixel
1427   
1428    return( (pi/lambda)*(bs + d_det + a_pixel)/l2s )
1429End
1430
1431Function qMaxVert()
1432
1433    Variable theta
1434    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
1435    NVAR lambda = root:SAS:gLambda
1436        NVAR det_width = root:SAS:det_width
1437       
1438    theta = atan( (det_width/2.0)/l2s )
1439   
1440    return ( 4.0*pi/lambda * sin(theta/2.0) )
1441end
1442
1443Function qMaxCorner()
1444
1445    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
1446    Variable radial
1447    NVAR lambda = root:SAS:gLambda
1448        NVAR det_off = root:SAS:gOffset
1449        NVAR det_width = root:SAS:det_width
1450
1451    radial=sqrt( (0.5*det_width)*(0.5*det_width)+(0.5*det_width+det_off)*(0.5*det_width+det_off) )
1452   
1453    return ( 4*pi/lambda*sin(0.5*atan(radial/l2s)) )
1454End
1455
1456Function qMaxHoriz()
1457
1458    Variable theta
1459    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
1460    NVAR lambda = root:SAS:gLambda
1461        NVAR det_off = root:SAS:gOffset
1462        NVAR det_width = root:SAS:det_width
1463
1464    theta = atan( ((det_width/2.0)+det_off)/l2s )       //from the instance variables
1465   
1466    return ( 4.0*pi/lambda * sin(theta/2.0) )
1467End
1468
1469// calculate sigma for the resolution function at either limit of q-range
1470Function deltaQ(atQ)
1471        Variable atQ
1472       
1473    Variable k02,lp,l1,l2,sig_02,sigQ2,a1,a2
1474    NVAR l2Diff = root:SAS:L2diff
1475        NVAR lambda = root:SAS:gLambda
1476        NVAR lambda_width = root:SAS:gDeltaLambda
1477        NVAR d_det = root:SAS:d_det
1478        NVAR del_r = root:SAS:del_r
1479       
1480       
1481    l1 = sourceToSampleDist()
1482    l2 = sampleToDetectorDist() + L2diff
1483    a1 = sourceApertureDiam()
1484    a2 = sampleApertureDiam()
1485   
1486    k02 = (6.2832/lambda)*(6.2832/lambda)
1487    lp = 1/(1/l1 + 1/l2)
1488   
1489    sig_02 = (0.25*a1/l1)*(0.25*a1/l1)
1490    sig_02 += (0.25*a2/lp)*(0.25*a2/lp)
1491    sig_02 += (d_det/(2.355*l2))*(d_det/(2.355*l2))
1492    sig_02 += (del_r/l2)*(del_r/l2)/12
1493    sig_02 *= k02
1494   
1495    sigQ2 = sig_02 + (atQ*lambda_width)*(atQ*lambda_width)/6
1496
1497    return(100*sqrt(sigQ2)/atQ)
1498End
1499
1500
1501Function beamIntensity()
1502
1503    Variable alpha,f,t,t4,t5,t6,as,solid_angle,l1,d2_phi
1504    Variable a1,a2,retVal
1505    SetDataFolder root:SAS
1506    NVAR l_gap=l_gap,guide_width =guide_width,ng = gNg
1507    NVAR lambda_t=lambda_t,b=b,c=c
1508    NVAR lambda=gLambda,t1=t1,t2=t2,t3=t3,phi_0=phi_0
1509    NVAR lambda_width=gDeltaLambda
1510   
1511    l1 = sourceToSampleDist()
1512    a1 = sourceApertureDiam()
1513    a2 = sampleApertureDiam()
1514   
1515   
1516    alpha = (a1+a2)/(2*l1)      //angular divergence of beam
1517    f = l_gap*alpha/(2*guide_width)
1518    t4 = (1-f)*(1-f)
1519    t5 = exp(ng*ln(0.96))       // trans losses of guides in pre-sample flight
1520    t6 = 1 - lambda*(b-(ng/8)*(b-c))            //experimental correction factor
1521    t = t1*t2*t3*t4*t5*t6
1522   
1523    as = pi/4*a2*a2             //area of sample in the beam
1524    d2_phi = phi_0/(2*pi)
1525    d2_phi *= exp(4*ln(lambda_t/lambda))
1526    d2_phi *= exp(-1*(lambda_t*lambda_t/lambda/lambda))
1527
1528    solid_angle = pi/4* (a1/l1)*(a1/l1)
1529
1530    retVal = as * d2_phi * lambda_width * solid_angle * t
1531     SetDataFolder root:
1532    return (retVal)
1533end
1534
1535Function figureOfMerit()
1536
1537        Variable bi = beamIntensity()
1538        NVAR lambda = root:SAS:gLambda
1539       
1540    return (lambda*lambda*bi)
1541End
1542
1543//estimate the number of pixels in the beam, and enforce the maximum countrate per pixel (idmax)
1544Function attenuatorTransmission()
1545
1546    Variable num_pixels,i_pix           //i_pix = id in John's notation
1547    Variable bDiam = beamDiameter("horizontal") //!! note that prev calculations used bh (horizontal)
1548    Variable atten,a2
1549    SetDataFolder root:SAS
1550    NVAR a_pixel=a_pixel,idmax=idmax
1551   
1552    a2 = sampleApertureDiam()
1553   
1554    num_pixels = pi/4*(0.5*(a2+bDiam))*(0.5*(a2+bDiam))/a_pixel/a_pixel
1555    i_pix = ( beamIntensity() )/num_pixels
1556   
1557    atten = (i_pix < idmax) ? 1.0 : idmax/i_pix
1558    SetDataFolder root:
1559    return(atten)
1560End
1561
1562Function attenuatorNumber()
1563
1564    Variable atten = attenuatorTransmission()
1565    Variable af,nf,numAtten
1566    SetDataFolder root:SAS
1567    NVAR lambda=gLambda
1568   
1569    af = 0.498 + 0.0792*lambda - 1.66e-3*lambda*lambda
1570    nf = -ln(atten)/af          //floating point
1571   
1572    numAtten = trunc(nf) + 1                    //in c, (int)nf
1573    //correct for larger step thickness at n > 6
1574    if(numAtten > 6)
1575        numAtten = 7 + trunc( (numAtten-6)/2 )          //in c, numAtten = 7 + (int)( (numAtten-6)/2 )
1576    endif
1577   
1578    SetDatafolder root:
1579    return (numAtten)
1580End
Note: See TracBrowser for help on using the repository browser.