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

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

incorporated beamstop diameter projection to the qMin and shadow
factor calculations. very important at short SDD, negligible at long SDD.

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