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

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

changed to count frozen traces from 1, rather than zero, to keep Jeff happy

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