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

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

Changed upper limit of SDD @ NG3 to be 1317 cm, not 1319 cm (too keep Cedric happy)

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