source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SASCALC.ipf @ 409

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

Removed Angstrom symbol (Mac code) and replaced it with simply a capital "A" so strange symbols won't appear anymore on Win.
File size: 50.2 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
55Proc 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/W=SASCALC checkNG3
336        isNG3=V_value
337        ControlInfo/W=SASCALC 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/W=SASCALC 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/W=SASCALC 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/W=SASCALC 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/W=SASCALC 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(nCts)
600        variable nCts = 1e6
601       
602        test2DSphere(nCts)
603       
604        root:Packages:NIST:SAS:linear_data = root:sim_lin
605       
606        S_CircularAverageTo1D("Packages:NIST:SAS")
607       
608        //multiply by sphere FF (S_SphereForm(scale,radius,delrho,bkg,x))
609        // or Debye Function S_Debye(scale,rg,bkg,x)
610       
611        //root:Packages:NIST:SAS:aveint = S_SphereForm(1,80,1e-6,0,root:Packages:NIST:SAS:qval)
612        //root:Packages:NIST:SAS:aveint = S_Debye(1000,100,0,qval)
613       
614        // multiply by beamstop shadowing
615        //root:Packages:NIST:SAS:aveint *= root:Packages:NIST:SAS:fSubS
616       
617        //multiply by current offset (>=1)
618        //root:Packages:NIST:SAS:aveint *= root:Packages:NIST:SAS:gModelOffsetFactor
619       
620End
621
622// passing in zero from a control skips the calculation
623Function ReCalculateInten(doIt)
624        Variable doIt
625       
626        if(doIt==0)                     
627                return(0)
628        endif
629       
630        // update the wave with the beamstop diameter here, since I don't know what
631        // combinations of parameters will change the BS - but anytime the curve is
632        // recalculated, or the text displayed, the right BS must be present
633        beamstopDiam()
634       
635        S_CircularAverageTo1D("Packages:NIST:SAS")
636        WAVE aveint=root:Packages:NIST:SAS:aveint
637        WAVE qval=root:Packages:NIST:SAS:qval
638        WAVE fSubS=root:Packages:NIST:SAS:fSubS
639       
640        //multiply by sphere FF (S_SphereForm(scale,radius,delrho,bkg,x))
641        // or Debye Function S_Debye(scale,rg,bkg,x)
642       
643        //aveint = S_SphereForm(1,80,1e-6,0,qval)
644        aveint = S_Debye(1000,100,0.0,qval)
645
646        // multiply by beamstop shadowing
647        aveint *= fSubS
648       
649        //display the configuration text in a separate notebook
650        DisplayConfigurationText()
651       
652        return(0)
653End
654
655//freezes the current configuration
656// -1- duplicates the trace on the graph
657// -2- copies the configuration text to a second notebook window for printing
658Function FreezeButtonProc(ctrlName) : ButtonControl
659        String ctrlName
660
661        String str=""
662        NVAR ct=root:Packages:NIST:SAS:gFreezeCount
663       
664       
665        SetDataFolder root:Packages:NIST:SAS
666       
667        Duplicate/O aveint,$("aveint_"+num2str(ct))
668        Duplicate/O qval,$("qval_"+num2str(ct))
669        Appendtograph $("aveint_"+num2str(ct)) vs $("qval_"+num2str(ct))
670        ModifyGraph mode=3
671        ModifyGraph marker=19
672        switch(mod(ct,10))      // 10 different colors - black is the unfrozen color
673                case 0:
674                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,16385,16385)
675                        break
676                case 1:
677                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(2,39321,1)
678                        break
679                case 2:
680                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(0,0,65535)
681                        break
682                case 3:
683                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(39321,1,31457)
684                        break
685                case 4:
686                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(48059,48059,48059)
687                        break
688                case 5:
689                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,32768,32768)
690                        break
691                case 6:
692                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(0,65535,0)
693                        break
694                case 7:
695                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(16385,65535,65535)
696                        break
697                case 8:
698                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,32768,58981)
699                        break
700                case 9:
701                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(36873,14755,58982)
702                        break
703        endswitch
704       
705        NVAR offset = root:Packages:NIST:SAS:gModelOffsetFactor
706        offset = 2^ct
707        //multiply by current offset (>=1)
708        Wave inten = $("aveint_"+num2str(ct))
709        inten *= offset
710//      Print "new offset = ",offset
711       
712        ct +=1
713        SetDataFolder root:
714       
715        // create or append the configuration to a new notebook
716        if(WinType("Saved_Configurations")==0)
717                NewNotebook/F=1/N=Saved_Configurations /W=(480,400,880,725)
718                DoWindow/F SASCALC              //return focus to SASCALC window
719        endif
720        //append the text
721        sprintf str,"\rConfiguration #%d\r",ct-1
722        Notebook Saved_Configurations showRuler=0,defaultTab=20,selection={endOfFile, endOfFile}
723        Notebook Saved_Configurations font="Monaco",fSize=10,fstyle=1,text=str          //bold
724        Notebook Saved_Configurations font="Monaco",fSize=10,fstyle=0,text=SetConfigurationText()
725        return(0)
726End
727
728//clears the frozen traces on the graph, asks if you want to clear the saved text
729//
730Function ClearButtonProc(ctrlName) : ButtonControl
731        String ctrlName
732
733        NVAR ct=root:Packages:NIST:SAS:gFreezeCount
734        Variable ii
735        Setdatafolder root:Packages:NIST:SAS
736        for(ii=ct-1;ii>=1;ii-=1)
737        //remove all traces, replace aveint
738        // kill all waves _ct
739                RemoveFromGraph $("aveint_"+num2str(ii))
740                Killwaves/Z $("aveint_"+num2str(ii)),$("qval_"+num2str(ii))
741        endfor
742        ct=1
743        setdatafolder root:
744       
745        DoAlert 1,"Do you also want to clear the \"Saved Configurations\" text?"
746        if(V_flag == 1)         // yes
747                DoWindow/K Saved_Configurations
748        endif
749       
750        //reset offset value
751        NVAR offset = root:Packages:NIST:SAS:gModelOffsetFactor
752        offset = 1
753        ReCalculateInten(1)
754        return(0)
755End
756
757
758///////////////////////////////////////////////////////////
759// 19MAR07 uses correction for beamstop diameter projection to get shadow factor correct
760//
761Function S_CircularAverageTo1D(type)
762        String type
763       
764        Variable isCircular = 1
765       
766        //type is the data type to do the averaging on, and will be set as the current folder
767        //get the current displayed data (so the correct folder is used)
768        String destPath = "root:"+type
769       
770        //
771        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist,dr,ddr
772        Variable lambda,trans
773        WAVE reals = $(destPath + ":RealsRead")
774//      WAVE/T textread = $(destPath + ":TextRead")
775//      String fileStr = textread[3]
776       
777        // center of detector, for non-linear corrections
778        Variable pixelsX=128,pixelsY=128
779       
780        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
781        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
782       
783        // beam center, in pixels
784        x0 = reals[16]
785        y0 = reals[17]
786        //detector calibration constants
787        sx = reals[10]          //mm/pixel (x)
788        sx3 = reals[11]         //nonlinear coeff
789        sy = reals[13]          //mm/pixel (y)
790        sy3 = reals[14]         //nonlinear coeff
791       
792        dtsize = 10*reals[20]           //det size in mm
793        dtdist = 1000*reals[18]         // det distance in mm
794       
795        NVAR binWidth=root:Packages:NIST:SAS:gBinWidth
796//      Variable binWidth = 1
797       
798        dr = binWidth           // ***********annulus width set by user, default is one***********
799        ddr = dr*sx             //step size, in mm (this value should be passed to the resolution calculation, not dr 18NOV03)
800               
801        Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,dxbm,dybm,ii
802        Variable phi_rad,dphi_rad,phi_x,phi_y
803        Variable forward,mirror
804       
805        String side = "both"
806//      side = StringByKey("SIDE",keyListStr,"=",";")
807//      Print "side = ",side
808       
809        /// data wave is data in the current folder which was set at the top of the function
810        WAVE data=$(destPath + ":data")
811
812// fake mask that uses all of the detector
813        Make/O/N=(pixelsX,pixelsY) $(destPath + ":mask")
814        Wave mask = $(destPath + ":mask")
815        mask = 0
816        //two pixels all around
817        mask[0,1][] = 1
818        mask[126,127][] = 1
819        mask[][0,1] = 1
820        mask[][126,127] = 1
821        //
822        //pixels within rcentr of beam center are broken into 9 parts (units of mm)
823        rcentr = 100            //original
824//      rcentr = 0
825        // values for error if unable to estimate value
826        //large_num = 1e10
827        large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1
828        small_num = 1e-10
829       
830        // output wave are expected to exist (?) initialized to zero, what length?
831        // 200 points on VAX --- use 300 here, or more if SAXS data is used with 1024x1024 detector (1000 pts seems good)
832        Variable defWavePts=500
833        Make/O/N=(defWavePts) $(destPath + ":qval"),$(destPath + ":aveint")
834        Make/O/N=(defWavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave")
835        Make/O/N=(defWavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar")
836
837        WAVE qval = $(destPath + ":qval")
838        WAVE aveint = $(destPath + ":aveint")
839        WAVE ncells = $(destPath + ":ncells")
840        WAVE dsq = $(destPath + ":dsq")
841        WAVE sigave = $(destPath + ":sigave")
842        WAVE qbar = $(destPath + ":QBar")
843        WAVE sigmaq = $(destPath + ":SigmaQ")
844        WAVE fsubs = $(destPath + ":fSubS")
845       
846        qval = 0
847        aveint = 0
848        ncells = 0
849        dsq = 0
850        sigave = 0
851        qbar = 0
852        sigmaq = 0
853        fsubs = 0
854
855        dtdis2 = dtdist^2
856        nq = 1
857        xoffst=0
858        //distance of beam center from detector center
859        dxbm = S_FX(x0,sx3,xcenter,sx)
860        dybm = S_FY(y0,sy3,ycenter,sy)
861               
862        //BEGIN AVERAGE **********
863        Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1
864        Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p
865       
866        // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
867        // loop index corresponds to FORTRAN (old code)
868        // and the IGOR array indices must be adjusted (-1) to the correct address
869        ii=1
870        do
871                xi = ii
872                dxi = S_FX(xi,sx3,xcenter,sx)
873                dx = dxi-dxbm           //dx and dy are in mm
874               
875                jj = 1
876                do
877                        data_pixel = data[ii-1][jj-1]           //assign to local variable
878                        yj = jj
879                        dyj = S_FY(yj,sy3,ycenter,sy)
880                        dy = dyj - dybm
881                        if(!(mask[ii-1][jj-1]))                 //masked pixels = 1, skip if masked (this way works...)
882                                dr2 = (dx^2 + dy^2)^(0.5)               //distance from beam center NOTE dr2 used here - dr used above
883                                if(dr2>rcentr)          //keep pixel whole
884                                        nd = 1
885                                        fd = 1
886                                else                            //break pixel into 9 equal parts
887                                        nd = 3
888                                        fd = 2
889                                endif
890                                nd2 = nd^2
891                                ll = 1          //"el-el" loop index
892                                do
893                                        dxx = dx + (ll - fd)*sx/3
894                                        kk = 1
895                                        do
896                                                dyy = dy + (kk - fd)*sy/3
897                                                if(isCircular)
898                                                        //circular average, use all pixels
899                                                        //(increment)
900                                                        nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
901                                                else
902                                                        //a sector average - determine azimuthal angle
903                                                        dphi_p = S_dphi_pixel(dxx,dyy,phi_x,phi_y)
904                                                        if(dphi_p < dphi_rad)
905                                                                forward = 1                     //within forward sector
906                                                        else
907                                                                forward = 0
908                                                        Endif
909                                                        if((Pi - dphi_p) < dphi_rad)
910                                                                mirror = 1              //within mirror sector
911                                                        else
912                                                                mirror = 0
913                                                        Endif
914                                                        //check if pixel lies within allowed sector(s)
915                                                        if(cmpstr(side,"both")==0)              //both sectors
916                                                                if ( mirror || forward)
917                                                                        //increment
918                                                                        nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
919                                                                Endif
920                                                        else
921                                                                if(cmpstr(side,"right")==0)             //forward sector only
922                                                                        if(forward)
923                                                                                //increment
924                                                                                nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
925                                                                        Endif
926                                                                else                    //mirror sector only
927                                                                        if(mirror)
928                                                                                //increment
929                                                                                nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
930                                                                        Endif
931                                                                Endif
932                                                        Endif           //allowable sectors
933                                                Endif   //circular or sector check
934                                                kk+=1
935                                        while(kk<=nd)
936                                        ll += 1
937                                while(ll<=nd)
938                        Endif           //masked pixel check
939                        jj += 1
940                while (jj<=pixelsY)
941                ii += 1
942        while(ii<=pixelsX)              //end of the averaging
943               
944        //compute q-values and errors
945        Variable ntotal,rr,theta,avesq,aveisq,var
946       
947        lambda = reals[26]
948        ntotal = 0
949        kk = 1
950        do
951                rr = (2*kk-1)*ddr/2
952                theta = 0.5*atan(rr/dtdist)
953                qval[kk-1] = (4*Pi/lambda)*sin(theta)
954                if(ncells[kk-1] == 0)
955                        //no pixels in annuli, data unknown
956                        aveint[kk-1] = 0
957                        sigave[kk-1] = large_num
958                else
959                        if(ncells[kk-1] <= 1)
960                                //need more than one pixel to determine error
961                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
962                                sigave[kk-1] = large_num
963                        else
964                                //assume that the intensity in each pixel in annuli is normally
965                                // distributed about mean...
966                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
967                                avesq = aveint[kk-1]^2
968                                aveisq = dsq[kk-1]/ncells[kk-1]
969                                var = aveisq-avesq
970                                if(var<=0)
971                                        sigave[kk-1] = small_num
972                                else
973                                        sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1))
974                                endif
975                        endif
976                endif
977                ntotal += ncells[kk-1]
978                kk+=1
979        while(kk<=nq)
980       
981        //Print "NQ = ",nq
982        // data waves were defined as 300 points (=defWavePts), but now have less than that (nq) points
983        // use DeletePoints to remove junk from end of waves
984        //WaveStats would be a more foolproof implementation, to get the # points in the wave
985        Variable startElement,numElements
986        startElement = nq
987        numElements = defWavePts - startElement
988        DeletePoints startElement,numElements, qval,aveint,ncells,dsq,sigave
989       
990        //////////////end of VAX sector_ave()
991               
992
993
994// ***************************************************************
995//
996// Do the extra 3 columns of resolution calculations starting here.
997//
998// ***************************************************************
999
1000        Variable L2 = reals[18]
1001//      Variable BS = reals[21]         //this the diameter is stored in mm
1002        Variable BS = beamstopDiamProjection(1) * 10            //calculated projection in cm *10 = mm
1003        Variable S1 = reals[23]
1004        Variable S2 = reals[24]
1005        Variable L1 = reals[25]
1006        lambda = reals[26]
1007        Variable lambdaWidth = reals[27]
1008
1009        Variable DDet, apOff
1010        //typical value for NG3 and NG7 - distance between sample aperture and sample in (cm)
1011        apOff=5.0
1012        // hard wire value for Ordela detectors
1013        DDet = 0.5              // resolution in cm
1014        //      String detStr=textRead[9]
1015        //      DDet = DetectorPixelResolution(fileStr,detStr)          //needs detector type and beamline
1016
1017        //Go from 0 to nq doing the calc for all three values at
1018        //every Q value
1019
1020        ii=0
1021        Variable ret1,ret2,ret3
1022        do
1023                S_getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,ret1,ret2,ret3)
1024                sigmaq[ii] = ret1       //res_wave[0]
1025                qbar[ii] = ret2         //res_wave[1]
1026                fsubs[ii] = ret3                //res_wave[2]
1027                ii+=1
1028        while(ii<nq)
1029        DeletePoints startElement,numElements, sigmaq, qbar, fsubs
1030
1031        fsubs += 1e-8           //keep the values from being too small
1032
1033// End of resolution calculations
1034// ***************************************************************
1035       
1036        //get rid of the default mask, if one was created (it is in the current folder)
1037        //don't just kill "mask" since it might be pointing to the one in the MSK folder
1038        Killwaves/Z $(destPath+":mask")
1039       
1040        //return to root folder (redundant)
1041        SetDataFolder root:
1042       
1043        Return 0
1044End
1045
1046//returns nq, new number of q-values
1047//arrays aveint,dsq,ncells are also changed by this function
1048//
1049Function S_IncrementPixel(dataPixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1050        Variable dataPixel,ddr,dxx,dyy
1051        Wave aveint,dsq,ncells
1052        Variable nq,nd2
1053       
1054        Variable ir
1055       
1056        ir = trunc(sqrt(dxx*dxx+dyy*dyy)/ddr)+1
1057        if (ir>nq)
1058                nq = ir         //resets maximum number of q-values
1059        endif
1060        aveint[ir-1] += dataPixel/nd2           //ir-1 must be used, since ir is physical
1061        dsq[ir-1] += dataPixel*dataPixel/nd2
1062        ncells[ir-1] += 1/nd2
1063       
1064        Return nq
1065End
1066
1067//function determines azimuthal angle dphi that a vector connecting
1068//center of detector to pixel makes with respect to vector
1069//at chosen azimuthal angle phi -> [cos(phi),sin(phi)] = [phi_x,phi_y]
1070//dphi is always positive, varying from 0 to Pi
1071//
1072Function S_dphi_pixel(dxx,dyy,phi_x,phi_y)
1073        Variable dxx,dyy,phi_x,phi_y
1074       
1075        Variable val,rr,dot_prod
1076       
1077        rr = sqrt(dxx^2 + dyy^2)
1078        dot_prod = (dxx*phi_x + dyy*phi_y)/rr
1079        //? correct for roundoff error? - is this necessary in IGOR, w/ double precision?
1080        if(dot_prod > 1)
1081                dot_prod =1
1082        Endif
1083        if(dot_prod < -1)
1084                dot_prod = -1
1085        Endif
1086       
1087        val = acos(dot_prod)
1088       
1089        return val
1090
1091End
1092
1093//calculates the x distance from the center of the detector, w/nonlinear corrections
1094//
1095Function S_FX(xx,sx3,xcenter,sx)               
1096        Variable xx,sx3,xcenter,sx
1097       
1098        Variable retval
1099       
1100        retval = sx3*tan((xx-xcenter)*sx/sx3)
1101        Return retval
1102End
1103
1104//calculates the y distance from the center of the detector, w/nonlinear corrections
1105//
1106Function S_FY(yy,sy3,ycenter,sy)               
1107        Variable yy,sy3,ycenter,sy
1108       
1109        Variable retval
1110       
1111        retval = sy3*tan((yy-ycenter)*sy/sy3)
1112        Return retval
1113End
1114
1115//**********************
1116// Resolution calculation - used by the averaging routines
1117// to calculate the resolution function at each q-value
1118// - the return value is not used
1119//
1120// equivalent to John's routine on the VAX Q_SIGMA_AVE.FOR
1121// Incorporates eqn. 3-15 from J. Appl. Cryst. (1995) v. 28 p105-114
1122//
1123Function/S S_getResolution(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,SigmaQ,QBar,fSubS)
1124        Variable inQ, lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r
1125        Variable &fSubS, &QBar, &SigmaQ         //these are the output quantities at the input Q value
1126       
1127        //lots of calculation variables
1128        Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g
1129        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
1130
1131        //Constants
1132        //Variable del_r = .1
1133        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
1134        Variable g = 981.0                              //gravity acceleration [cm/s^2]
1135
1136        String results
1137        results ="Failure"
1138       
1139        NVAR usingLenses = root:Packages:NIST:SAS:gUsingLenses
1140
1141        //rename for working variables,  these must be gotten from global
1142        //variables
1143
1144//      Variable wLam, wLW, wL1, wL2, wS1, wS2
1145//      Variable wBS, wDDet, wApOff
1146//      wLam = lambda
1147//      wLW = lambdaWidth
1148//      wDDet = DDet
1149//      wApOff = apOff
1150        S1 *= 0.5*0.1                   //convert to radius and [cm]
1151        S2 *= 0.5*0.1
1152
1153        L1 *= 100.0                     // [cm]
1154        L1 -= apOff                             //correct the distance
1155
1156        L2 *= 100.0
1157        L2 += apOff
1158
1159        BS *= 0.5*0.1                   //convert to radius and [cm]
1160        del_r *= 0.1                            //width of annulus, convert mm to [cm]
1161       
1162        //Start resolution calculation
1163        a2 = S1*L2/L1 + S2*(L1+L2)/L1
1164        q_small = 2.0*Pi*(BS-a2)*(1.0-lambdaWidth)/(lambda*L2)
1165        lp = 1.0/( 1.0/L1 + 1.0/L2)
1166
1167        v_lambda = lambdaWidth^2/6.0
1168       
1169        if(usingLenses==1)                      //SRK 2007
1170                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term
1171        else
1172                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form
1173        endif
1174       
1175        v_d = (DDet/2.3548)^2 + del_r^2/12.0
1176        vz = vz_1 / lambda
1177        yg = 0.5*g*L2*(L1+L2)/vz^2
1178        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007
1179
1180        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
1181        delta = 0.5*(BS - r0)^2/v_d
1182
1183        if (r0 < BS)
1184                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
1185        else
1186                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
1187        endif
1188
1189        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
1190        if (fSubS <= 0.0)
1191                fSubS = 1.e-10
1192        endif
1193        fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
1194        fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
1195
1196        rmd = fr*r0
1197        v_r1 = v_b + fv*v_d +v_g
1198
1199        rm = rmd + 0.5*v_r1/rmd
1200        v_r = v_r1 - 0.5*(v_r1/rmd)^2
1201        if (v_r < 0.0)
1202                v_r = 0.0
1203        endif
1204        QBar = (4.0*Pi/lambda)*sin(0.5*atan(rm/L2))
1205        SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda)
1206
1207        results = "success"
1208        Return results
1209End
1210
1211Function S_Debye(scale,rg,bkg,x)
1212        Variable scale,rg,bkg
1213        Variable x
1214       
1215        // variables are:
1216        //[0] scale factor
1217        //[1] radius of gyration [A]
1218        //[2] background        [cm-1]
1219       
1220        // calculates (scale*debye)+bkg
1221        Variable Pq,qr2
1222       
1223        qr2=(x*rg)^2
1224        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
1225       
1226        //scale
1227        Pq *= scale
1228        // then add in the background
1229        return (Pq+bkg)
1230End
1231
1232
1233Function S_SphereForm(scale,radius,delrho,bkg,x)                               
1234        Variable scale,radius,delrho,bkg
1235        Variable x
1236       
1237        // variables are:                                                       
1238        //[0] scale
1239        //[1] radius (A)
1240        //[2] delrho (A-2)
1241        //[3] background (cm-1)
1242       
1243//      Variable scale,radius,delrho,bkg                               
1244//      scale = w[0]
1245//      radius = w[1]
1246//      delrho = w[2]
1247//      bkg = w[3]
1248       
1249       
1250        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
1251        // and is rescaled to give [=] cm^-1
1252       
1253        Variable bes,f,vol,f2
1254        //
1255        //handle q==0 separately
1256        If(x==0)
1257                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
1258                return(f)
1259        Endif
1260       
1261        bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
1262        vol = 4*pi/3*radius^3
1263        f = vol*bes*delrho              // [=] A
1264        // normalize to single particle volume, convert to 1/cm
1265        f2 = f * f / vol * 1.0e8                // [=] 1/cm
1266       
1267        return (scale*f2+bkg)   // Scale, then add in the background
1268       
1269End
1270
1271Function/S SetConfigurationText()
1272
1273        String str="",temp
1274
1275        SetDataFolder root:Packages:NIST:SAS
1276       
1277        NVAR numberOfGuides=gNg
1278        NVAR gTable=gTable              //2=chamber, 1=table
1279        NVAR wavelength=gLambda
1280        NVAR lambdaWidth=gDeltaLambda
1281        NVAR instrument = instrument
1282        NVAR L2diff = L2diff
1283        NVAR lens = root:Packages:NIST:SAS:gUsingLenses
1284        SVAR aStr = root:myGlobals:gAngstStr
1285       
1286        sprintf temp,"Source Aperture Diameter =\t\t%6.2f cm\r",sourceApertureDiam()
1287        str += temp
1288        sprintf temp,"Source to Sample =\t\t\t\t%6.0f cm\r",sourceToSampleDist()
1289        str += temp
1290        sprintf temp,"Sample Aperture to Detector =\t%6.0f cm\r",sampleToDetectorDist() + L2diff
1291        str += temp
1292        sprintf temp,"Beam diameter =\t\t\t\t\t%6.2f cm\r",beamDiameter("maximum")
1293        str += temp
1294        sprintf temp,"Beamstop diameter =\t\t\t\t%6.2f inches\r",beamstopDiam()/2.54
1295        str += temp
1296        sprintf temp,"Minimum Q-value =\t\t\t\t%6.4f 1/%s (sigQ/Q = %3.1f %%)\r",qMin(),aStr,deltaQ(qMin())
1297        str += temp
1298        sprintf temp,"Maximum Horizontal Q-value =\t%6.4f 1/%s\r",qMaxHoriz(),aStr
1299        str += temp
1300        sprintf temp,"Maximum Vertical Q-value =\t\t%6.4f 1/%s\r",qMaxVert(),aStr
1301        str += temp
1302        sprintf temp,"Maximum Q-value =\t\t\t\t%6.4f 1/%s (sigQ/Q = %3.1f %%)\r",qMaxCorner(),aStr,deltaQ(qMaxCorner())
1303        str += temp
1304        sprintf temp,"Beam Intensity =\t\t\t\t%.0f counts/s\r",beamIntensity()
1305        str += temp
1306        sprintf temp,"Figure of Merit =\t\t\t\t%3.3g %s^2/s\r",figureOfMerit(),aStr
1307        str += temp
1308        sprintf temp,"Attenuator transmission =\t\t%3.3g = Atten # %d\r",attenuatorTransmission(),attenuatorNumber()
1309        str += temp
1310//     
1311//      // add text of the user-edited values
1312//      //
1313        sprintf temp,"***************** NG %d *****************\r",instrument
1314        str += temp
1315        sprintf temp,"Sample Aperture Diameter =\t\t\t\t%.2f cm\r",sampleApertureDiam()
1316        str += temp
1317        sprintf temp,"Number of Guides =\t\t\t\t\t\t%d \r", numberOfGuides
1318        str += temp
1319        sprintf temp,"Sample Chamber to Detector =\t\t\t%.1f cm\r", chamberToDetectorDist()
1320        str += temp
1321        if(gTable==1)
1322                sprintf temp,"Sample Position is \t\t\t\t\t\tHuber\r"
1323        else
1324                sprintf temp,"Sample Position is \t\t\t\t\t\tChamber\r"
1325        endif
1326        str += temp
1327        sprintf temp,"Detector Offset =\t\t\t\t\t\t%.1f cm\r", detectorOffset()
1328        str += temp
1329        sprintf temp,"Neutron Wavelength =\t\t\t\t\t%.2f %s\r", wavelength,aStr
1330        str += temp
1331        sprintf temp,"Wavelength Spread, FWHM =\t\t\t\t%.3f\r", lambdaWidth
1332        str += temp
1333        sprintf temp,"Sample Aperture to Sample Position =\t%.2f cm\r", L2Diff
1334        str += temp
1335        if(lens==1)
1336                sprintf temp,"Lenses are IN\r"
1337        else
1338                sprintf temp,"Lenses are OUT\r"
1339        endif
1340        str += temp
1341       
1342   setDataFolder root:
1343   return str                   
1344End
1345
1346Function DisplayConfigurationText()
1347
1348        if(WinType("Trial_Configuration")==0)
1349                NewNotebook/F=0/K=1/N=Trial_Configuration /W=(480,44,880,369)
1350        endif
1351        //replace the text
1352        Notebook Trial_Configuration selection={startOfFile, endOfFile}
1353        Notebook Trial_Configuration font="Monaco",fSize=10,text=SetConfigurationText()
1354        return(0)
1355end
1356
1357//parses the control for A1 diam
1358// updates the wave
1359Function sourceApertureDiam()
1360        ControlInfo/W=SASCALC popup0
1361        Variable diam
1362        sscanf S_Value, "%g cm", diam
1363        WAVE rw=root:Packages:NIST:SAS:realsRead
1364        rw[23] = diam*10                        //sample aperture diameter in mm
1365        return(diam)
1366end
1367
1368// change the sample aperture to a non-standard value
1369Function SourceApOtherSetVarProc(ctrlName,varNum,varStr,varName) : SetVariableControl
1370        String ctrlName
1371        Variable varNum
1372        String varStr
1373        String varName
1374               
1375        WAVE rw=root:Packages:NIST:SAS:realsRead
1376        rw[24] = varNum                 //sample aperture diameter in mm
1377        ReCalculateInten(1)
1378        return(0)
1379End
1380
1381//parses the control for A2 diam
1382// updates the wave and global
1383// returns a2 in cm
1384Function sampleApertureDiam()
1385        ControlInfo/W=SASCALC popup0_1
1386       
1387//      Print "In sampleApertureDiam()"
1388        //set the global
1389        NVAR a2=root:Packages:NIST:SAS:gSamAp
1390       
1391        if(cmpstr(S_Value,"other") == 0)                // "other" selected
1392                //enable the setvar, diameter in mm!
1393                SetVariable setvar0_3 disable=0
1394                // read its value (a global)
1395                NVAR a2other = root:Packages:NIST:SAS:gSamApOther
1396                a2=a2other/10                           //a2 in cm
1397        else
1398                SetVariable setvar0_3 disable=1
1399                //1st item is 1/16", popup steps by 1/16"
1400                a2 = 2.54/16.0 * (V_Value)                      //convert to cm         
1401        endif
1402        WAVE rw=root:Packages:NIST:SAS:realsRead
1403        rw[24] = a2*10                  //sample aperture diameter in mm
1404       
1405        return(a2)
1406end
1407
1408//1=huber 2=chamber
1409Function tableposition()
1410        ControlInfo/W=SASCALC checkHuber
1411        if(V_Value)
1412                return(1)               //huber
1413        else
1414                return(2)               //chamber
1415        endif
1416End
1417
1418//compute SSD and update both the global and the wave
1419//
1420Function sourceToSampleDist()
1421
1422        NVAR NG=root:Packages:NIST:SAS:gNg
1423        NVAR S12 = root:Packages:NIST:SAS:S12
1424        NVAR L2Diff = root:Packages:NIST:SAS:L2Diff
1425        NVAR SSD = root:Packages:NIST:SAS:gSSD
1426       
1427        SSD = 1632 - 155*NG - s12*(2-tableposition()) - L2Diff
1428       
1429        WAVE rw=root:Packages:NIST:SAS:realsRead
1430        rw[25] = SSD/100                // in meters
1431        return(SSD)
1432End
1433
1434//returns the offset value
1435// slider and setVar are linked to the same global
1436// updates the wave and changes the beamcenter (x,y) in the wave
1437Function detectorOffset()
1438       
1439        WAVE rw=root:Packages:NIST:SAS:RealsRead
1440        NVAR val = root:Packages:NIST:SAS:gOffset
1441        rw[19] = val            // already in cm
1442        //move the beamcenter
1443        rw[16] = 64 + 2*rw[19]          //approximate beam X is 64 w/no offset, 114 w/25 cm offset
1444        rw[17] = 64             //typical value
1445       
1446        return(val)
1447end
1448
1449//returns the detector distance (slider and setVar are linked by the global)
1450//
1451Function chamberToDetectorDist()
1452
1453        NVAR val = root:Packages:NIST:SAS:gDetDist
1454        return(val)
1455End
1456
1457//sets the SDD (slider and setVar are linked by the global and is the detector position
1458//  relative to the chamber)
1459// updates the wave
1460Function sampleToDetectorDist()
1461
1462        NVAR detDist = root:Packages:NIST:SAS:gDetDist
1463        NVAR S12 = root:Packages:NIST:SAS:S12
1464        WAVE rw=root:Packages:NIST:SAS:RealsRead       
1465        Variable SDD   
1466       
1467        SDD = detDist + s12*(2-tableposition())
1468        rw[18] = SDD/100                // convert to meters for header
1469        return(SDD)
1470End
1471
1472//direction = one of "vertical;horizontal;maximum;"
1473// all of this is bypassed if the lenses are in
1474//
1475Function beamDiameter(direction)
1476        String direction
1477
1478        NVAR lens = root:Packages:NIST:SAS:gUsingLenses
1479        if(lens)
1480                return sourceApertureDiam()
1481        endif
1482       
1483    Variable l1 = sourceToSampleDist()
1484    Variable l2 //= sampleAperToDetDist()
1485    Variable d1,d2,bh,bv,bm,umbra,a1,a2
1486   
1487    NVAR L2diff = root:Packages:NIST:SAS:L2diff
1488    NVAR lambda = root:Packages:NIST:SAS:gLambda
1489    NVAR lambda_width = root:Packages:NIST:SAS:gDeltaLambda
1490    NVAR bs_factor = root:Packages:NIST:SAS:bs_factor
1491   
1492    l2 = sampleToDetectorDist() + L2diff
1493    a1 = sourceApertureDiam()
1494    a2 = sampleApertureDiam()
1495   
1496    d1 = a1*l2/l1
1497    d2 = a2*(l1+l2)/l1
1498    bh = d1+d2          //beam size in horizontal direction
1499    umbra = abs(d1-d2)
1500    //vertical spreading due to gravity
1501    bv = bh + 1.25e-8*(l1+l2)*l2*lambda*lambda*lambda_width
1502    bm = (bs_factor*bh > bv) ? bs_factor*bh : bv //use the larger of horiz*safety or vertical
1503   
1504    strswitch(direction)        // string switch
1505        case "vertical":                // execute if case matches expression
1506                return(bv)
1507                break                                           // exit from switch
1508        case "horizontal":              // execute if case matches expression
1509                return(bh)
1510                break
1511        case "maximum":         // execute if case matches expression
1512                return(bm)
1513                break
1514        default:                                                        // optional default expression executed
1515                return(bm)                                              // when no case matches
1516    endswitch
1517End
1518
1519//on NG3 and NG7, allowable sizes are 1,2,3,4" diameter
1520//will return values larger than 4.0*2.54 if a larger beam is needed
1521//
1522// - in an approximate way, account for lenses
1523Function beamstopDiam()
1524
1525        NVAR yesLens = root:Packages:NIST:SAS:gUsingLenses
1526        Variable bm=0
1527        Variable bs=0.0
1528   
1529        if(yesLens)
1530                //bm = sourceApertureDiam()             //ideal result, not needed
1531                bs = 1                                                          //force the diameter to 1"
1532        else
1533                bm = beamDiameter("maximum")
1534                do
1535                bs += 1
1536           while ( (bs*2.54 < bm) || (bs > 30.0))                       //30 = ridiculous limit to avoid inf loop
1537        endif
1538
1539        //update the wave
1540        WAVE rw=root:Packages:NIST:SAS:realsRead
1541        rw[21] = bs*25.4                //store the BS diameter in mm
1542       
1543    return (bs*2.54)            //return diameter in cm, not inches for txt
1544End
1545
1546//returns the projected diameter of the beamstop at the anode plane.
1547// most noticeable at short SDD
1548//if flag == 0 use conservative estimate = largest diameter (for SASCALC, default)
1549//if flag != 0 use point aperture = average diameter (for resolution calculation)
1550Function beamstopDiamProjection(flag)
1551        Variable flag
1552       
1553        NVAR L2diff = root:Packages:NIST:SAS:L2diff
1554        Variable a2 = sampleApertureDiam()
1555        Variable bs = beamstopDiam()
1556        Variable l2, LB, BS_P
1557   
1558        l2 = sampleToDetectorDist() + L2diff
1559        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical)
1560        if(flag==0)
1561                BS_P = bs + (bs+a2)*lb/(l2-lb)          //diameter of shadow from parallax
1562        else
1563                BS_P = bs + bs*lb/(l2-lb)               //diameter of shadow, point A2
1564        endif
1565        return (bs_p)           //return projected diameter in cm
1566End
1567
1568// 19MAR07 - using correction from John for an estimate of the shadow of the beamstop
1569// at the detection plane. This is a noticeable effect at short SDD, where the projected
1570// diameter of the beamstop is much larger than the physical diameter.
1571Function qMin()
1572
1573    Variable l2s = sampleToDetectorDist()       //distance from sample to detector in cm
1574//    Variable bs = beamstopDiam()              //beamstop diameter in cm
1575    Variable bs_p = beamstopDiamProjection(0)           //projected beamstop diameter in cm
1576    NVAR lambda = root:Packages:NIST:SAS:gLambda
1577    NVAR d_det = root:Packages:NIST:SAS:d_det           //cm
1578    NVAR a_pixel = root:Packages:NIST:SAS:a_pixel       //cm
1579
1580    return( (pi/lambda)*(bs_p + d_det + a_pixel)/l2s )          //use bs_p rather than bs
1581   // return( (pi/lambda)*(bs + d_det + a_pixel)/l2s )          //use bs (incorrect)
1582End
1583
1584Function qMaxVert()
1585
1586    Variable theta
1587    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
1588    NVAR lambda = root:Packages:NIST:SAS:gLambda
1589        NVAR det_width = root:Packages:NIST:SAS:det_width
1590       
1591    theta = atan( (det_width/2.0)/l2s )
1592   
1593    return ( 4.0*pi/lambda * sin(theta/2.0) )
1594end
1595
1596Function qMaxCorner()
1597
1598    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
1599    Variable radial
1600    NVAR lambda = root:Packages:NIST:SAS:gLambda
1601        NVAR det_off = root:Packages:NIST:SAS:gOffset
1602        NVAR det_width = root:Packages:NIST:SAS:det_width
1603
1604    radial=sqrt( (0.5*det_width)*(0.5*det_width)+(0.5*det_width+det_off)*(0.5*det_width+det_off) )
1605   
1606    return ( 4*pi/lambda*sin(0.5*atan(radial/l2s)) )
1607End
1608
1609Function qMaxHoriz()
1610
1611    Variable theta
1612    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
1613    NVAR lambda = root:Packages:NIST:SAS:gLambda
1614        NVAR det_off = root:Packages:NIST:SAS:gOffset
1615        NVAR det_width = root:Packages:NIST:SAS:det_width
1616
1617    theta = atan( ((det_width/2.0)+det_off)/l2s )       //from the instance variables
1618   
1619    return ( 4.0*pi/lambda * sin(theta/2.0) )
1620End
1621
1622// calculate sigma for the resolution function at either limit of q-range
1623Function deltaQ(atQ)
1624        Variable atQ
1625       
1626    Variable k02,lp,l1,l2,sig_02,sigQ2,a1,a2
1627    NVAR l2Diff = root:Packages:NIST:SAS:L2diff
1628        NVAR lambda = root:Packages:NIST:SAS:gLambda
1629        NVAR lambda_width = root:Packages:NIST:SAS:gDeltaLambda
1630        NVAR d_det = root:Packages:NIST:SAS:d_det
1631        NVAR del_r = root:Packages:NIST:SAS:del_r
1632       
1633       
1634    l1 = sourceToSampleDist()
1635    l2 = sampleToDetectorDist() + L2diff
1636    a1 = sourceApertureDiam()
1637    a2 = sampleApertureDiam()
1638   
1639    k02 = (6.2832/lambda)*(6.2832/lambda)
1640    lp = 1/(1/l1 + 1/l2)
1641   
1642    sig_02 = (0.25*a1/l1)*(0.25*a1/l1)
1643    sig_02 += (0.25*a2/lp)*(0.25*a2/lp)
1644    sig_02 += (d_det/(2.355*l2))*(d_det/(2.355*l2))
1645    sig_02 += (del_r/l2)*(del_r/l2)/12
1646    sig_02 *= k02
1647   
1648    sigQ2 = sig_02 + (atQ*lambda_width)*(atQ*lambda_width)/6
1649
1650    return(100*sqrt(sigQ2)/atQ)
1651End
1652
1653
1654Function beamIntensity()
1655
1656    Variable alpha,f,t,t4,t5,t6,as,solid_angle,l1,d2_phi
1657    Variable a1,a2,retVal
1658    SetDataFolder root:Packages:NIST:SAS
1659    NVAR l_gap=l_gap,guide_width =guide_width,ng = gNg
1660    NVAR lambda_t=lambda_t,b=b,c=c
1661    NVAR lambda=gLambda,t1=t1,t2=t2,t3=t3,phi_0=phi_0
1662    NVAR lambda_width=gDeltaLambda
1663   
1664    l1 = sourceToSampleDist()
1665    a1 = sourceApertureDiam()
1666    a2 = sampleApertureDiam()
1667   
1668   
1669    alpha = (a1+a2)/(2*l1)      //angular divergence of beam
1670    f = l_gap*alpha/(2*guide_width)
1671    t4 = (1-f)*(1-f)
1672    t5 = exp(ng*ln(0.96))       // trans losses of guides in pre-sample flight
1673    t6 = 1 - lambda*(b-(ng/8)*(b-c))            //experimental correction factor
1674    t = t1*t2*t3*t4*t5*t6
1675   
1676    as = pi/4*a2*a2             //area of sample in the beam
1677    d2_phi = phi_0/(2*pi)
1678    d2_phi *= exp(4*ln(lambda_t/lambda))
1679    d2_phi *= exp(-1*(lambda_t*lambda_t/lambda/lambda))
1680
1681    solid_angle = pi/4* (a1/l1)*(a1/l1)
1682
1683    retVal = as * d2_phi * lambda_width * solid_angle * t
1684     SetDataFolder root:
1685    return (retVal)
1686end
1687
1688Function figureOfMerit()
1689
1690        Variable bi = beamIntensity()
1691        NVAR lambda = root:Packages:NIST:SAS:gLambda
1692       
1693    return (lambda*lambda*bi)
1694End
1695
1696//estimate the number of pixels in the beam, and enforce the maximum countrate per pixel (idmax)
1697Function attenuatorTransmission()
1698
1699    Variable num_pixels,i_pix           //i_pix = id in John's notation
1700    Variable bDiam = beamDiameter("horizontal") //!! note that prev calculations used bh (horizontal)
1701    Variable atten,a2
1702    SetDataFolder root:Packages:NIST:SAS
1703    NVAR a_pixel=a_pixel,idmax=idmax
1704   
1705    a2 = sampleApertureDiam()
1706   
1707    num_pixels = pi/4*(0.5*(a2+bDiam))*(0.5*(a2+bDiam))/a_pixel/a_pixel
1708    i_pix = ( beamIntensity() )/num_pixels
1709   
1710    atten = (i_pix < idmax) ? 1.0 : idmax/i_pix
1711    SetDataFolder root:
1712    return(atten)
1713End
1714
1715Function attenuatorNumber()
1716
1717    Variable atten = attenuatorTransmission()
1718    Variable af,nf,numAtten
1719    SetDataFolder root:Packages:NIST:SAS
1720    NVAR lambda=gLambda
1721   
1722    af = 0.498 + 0.0792*lambda - 1.66e-3*lambda*lambda
1723    nf = -ln(atten)/af          //floating point
1724   
1725    numAtten = trunc(nf) + 1                    //in c, (int)nf
1726    //correct for larger step thickness at n > 6
1727    if(numAtten > 6)
1728        numAtten = 7 + trunc( (numAtten-6)/2 )          //in c, numAtten = 7 + (int)( (numAtten-6)/2 )
1729    endif
1730   
1731    SetDatafolder root:
1732    return (numAtten)
1733End
Note: See TracBrowser for help on using the repository browser.