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

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

now uses canSANS approved data folders for internal storage

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