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

Last change on this file since 84 was 84, checked in by srkline, 15 years ago

Added default sample aperture to sample offset of 5.0 cm to match the VAX implementation.

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