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

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

initial import of SASCALC version 1.00 files

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