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

Last change on this file since 564 was 564, checked in by srkline, 13 years ago

A variety of changes...

  • fixed model strings for Sum_Model
  • /Z flag for RTHelp, Transmission help, FIT help, MRED help
  • added help and done buttons to SASCALC
  • added comments in MonteCarlo?
  • restored the 2-line buttons in USANS panels, now supported in Igor 6.1 (COR_Graph and Main_USANS)
  • added simple text notebook of UCALC configuration, and a save of PNG of the whole panel
  • added help and done buttons to UCALC
  • added links in SANS_Reduction_Help to the simluation help file (not added to SVN yet)
File size: 65.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// 27 JAN 2009 SRK - Changed behavior of Lens checkbox. Now, it SETS parameters as needed for proper
20//                                               configuration. 17.2 can be typed in for lens/prism on NG3. Invalid conditions
21//                                               will automatically uncheck the box
22//
23// calculate what q-values you get based on the instruments settings
24// that you would typically input to SASCALC
25// calculation is for (80A radius spheres) * (beam stop shadow factor)
26// or a Debye function (Rg=50A) * (beam stop shadow factor)
27// - NOT true intensity, not counts, just a display
28//
29// To Do:
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)             //will use defaults
62        Endif
63
64// now a checkbox as needed
65//      DoWindow/F MC_SASCALC
66//      if(V_flag==0)
67//              MC_SASCALC()
68//      endif
69End
70
71Proc S_initialize_space()
72        NewDataFolder/O root:Packages
73        NewDataFolder/O root:Packages:NIST
74        NewDataFolder/O root:Packages:NIST:SAS
75       
76        Make/O/D/N=23 root:Packages:NIST:SAS:integersRead
77        Make/O/D/N=52 root:Packages:NIST:SAS:realsRead
78        Make/O/T/N=11 root:Packages:NIST:SAS:textRead
79        // data
80        Make/O/D/N=(128,128) root:Packages:NIST:SAS:data,root:Packages:NIST:SAS:linear_data
81        Make/O/D/N=2 root:Packages:NIST:SAS:aveint,root:Packages:NIST:SAS:qval,root:Packages:NIST:SAS:sigave
82        root:Packages:NIST:SAS:data = 1
83        root:Packages:NIST:SAS:linear_data = 1
84        // fill w/default values
85        S_fillDefaultHeader(root:Packages:NIST:SAS:integersRead,root:Packages:NIST:SAS:realsRead,root:Packages:NIST:SAS:textRead)
86       
87        // other variables
88        // -(hard coded right now - look for NVAR declarations)
89        Variable/G root:Packages:NIST:SAS:gBinWidth=1
90        Variable/G root:Packages:NIST:SAS:gisLogScale=0
91        String/G root:Packages:NIST:SAS:FileList = "SASCALC"
92       
93        // for the panel
94        Variable/G root:Packages:NIST:SAS:gInst=3               //or 7 for NG7
95        Variable/G root:Packages:NIST:SAS:gNg=0
96        Variable/G root:Packages:NIST:SAS:gTable=2              //2=chamber, 1=table
97        Variable/G root:Packages:NIST:SAS:gDetDist=1000         //sample chamber to detector in cm
98        Variable/G root:Packages:NIST:SAS:gSSD=1632             //!!SSD in cm fo 0 guides (derived from Ng)
99        Variable/G root:Packages:NIST:SAS:gOffset=0
100        Variable/G root:Packages:NIST:SAS:gSamAp=1.27           //samAp diameter in cm
101        Variable/G root:Packages:NIST:SAS:gLambda=6
102        Variable/G root:Packages:NIST:SAS:gDeltaLambda=0.125            //default value
103        String/G root:Packages:NIST:SAS:gSourceApString = "1.43 cm;2.54 cm;3.81 cm;"
104        String/G root:Packages:NIST:SAS:gDeltaLambdaStr = "0.109;0.125;0.236;"          //ng3 defaults
105        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;"
106        Variable/G root:Packages:NIST:SAS:gSamApOther = 10              //non-standard aperture diameter, in mm
107        Variable/G root:Packages:NIST:SAS:gUsingLenses = 0              //0=no lenses, 1=lenses(or prisms)
108        Variable/G root:Packages:NIST:SAS:gModelOffsetFactor = 1
109       
110        // for the MC simulation
111        Variable/G root:Packages:NIST:SAS:doSimulation  =0              // == 1 if 1D simulated data, 0 if other from the checkbox
112        Variable/G root:Packages:NIST:SAS:gRanDateTime=datetime
113        Variable/G root:Packages:NIST:SAS:gImon = 10000
114        Variable/G root:Packages:NIST:SAS:gThick = 0.1
115        Variable/G root:Packages:NIST:SAS:gSig_incoh = 0.1
116        String/G root:Packages:NIST:SAS:gFuncStr = ""
117        Variable/G root:Packages:NIST:SAS:gR2 = 2.54/2 
118        Variable/G root:Packages:NIST:SAS:gSamTrans=0.8                 //for 1D, default value
119        Variable/G root:Packages:NIST:SAS:gCntTime = 300
120        Variable/G root:Packages:NIST:SAS:gDoMonteCarlo = 0
121        Variable/G root:Packages:NIST:SAS:gUse_MC_XOP = 1                               //set to zero to use Igor code
122        Variable/G root:Packages:NIST:SAS:gBeamStopIn = 1                       //set to zero for beamstop out (transmission)
123        Variable/G root:Packages:NIST:SAS:gRawCounts = 0
124        Variable/G root:Packages:NIST:SAS:gSaveIndex = 100
125        String/G root:Packages:NIST:SAS:gSavePrefix = "SIMUL"
126        Make/O/D/N=10 root:Packages:NIST:SAS:results = 0
127        Make/O/T/N=10 root:Packages:NIST:SAS:results_desc = {"total X-section (1/cm)","SAS X-section (1/cm)","number that scatter","number that reach detector","avg # times scattered","fraction single coherent","fraction multiple coherent","fraction multiple scattered","fraction transmitted","detector counts w/beamstop"}
128
129        Variable/G root:Packages:NIST:SAS:g_1DTotCts = 0                        //summed counts (simulated)
130        Variable/G root:Packages:NIST:SAS:g_1DEstDetCR = 0              // estimated detector count rate
131        Variable/G root:Packages:NIST:SAS:g_1DFracScatt = 0             // fraction of beam captured on detector
132        Variable/G root:Packages:NIST:SAS:g_1DEstTrans = 0              // estimated transmission of sample
133        Variable/G root:Packages:NIST:SAS:g_1D_DoABS = 1
134        Variable/G root:Packages:NIST:SAS:g_1D_AddNoise = 1
135       
136        //tick labels for SDD slider
137        //userTicks={tvWave,tlblWave }
138        Make/O/D/N=5 root:Packages:NIST:SAS:tickSDDNG3,root:Packages:NIST:SAS:tickSDDNG7
139        Make/O/T/N=5 root:Packages:NIST:SAS:lblSDDNG3,root:Packages:NIST:SAS:lblSDDNG7
140        root:Packages:NIST:SAS:tickSDDNG3 = {133,400,700,1000,1317}
141        root:Packages:NIST:SAS:lblSDDNG3 = {"133","400","700","1000","1317"}
142        root:Packages:NIST:SAS:tickSDDNG7 = {100,450,800,1150,1530}
143        root:Packages:NIST:SAS:lblSDDNG7 = {"100","450","800","1150","1530"}
144       
145        //for the fake dependency
146        Variable/G root:Packages:NIST:SAS:gTouched=0
147        Variable/G root:Packages:NIST:SAS:gCalculate=0
148        //for plotting
149        Variable/G root:Packages:NIST:SAS:gFreezeCount=1                //start the count at 1 to keep Jeff happy
150        Variable/G root:Packages:NIST:SAS:gDoTraceOffset=0              // (1==Yes, offset 2^n), 0==turn off the offset
151       
152End
153
154Function initNG3()
155
156        SetDataFolder root:Packages:NIST:SAS
157       
158        Variable/G instrument = 3
159        Variable/G s12 = 54.8
160        Variable/G d_det = 0.5
161        Variable/G a_pixel = 0.5
162        Variable/G del_r = 0.5
163        Variable/G det_width = 64.0
164        Variable/G phi_0 = 2.95e13
165        Variable/G lambda_t = 5.50
166        Variable/G l2r_lower = 132.3
167        Variable/G l2r_upper =  1317
168        Variable/G lambda_lower = 2.5
169        Variable/G lambda_upper = 20.0
170        Variable/G d_upper = 25.0
171        Variable/G bs_factor = 1.05
172        Variable/G t1 = 0.63
173        Variable/G t2 = 1.0
174        Variable/G t3 = 0.75
175        Variable/G l_gap = 100.0
176        Variable/G guide_width = 6.0
177        Variable/G idmax = 100.0
178        Variable/G b = 0.023
179        Variable/G c = 0.023
180       
181        //fwhm values (new variables) (+3, 0, -3, calibrated 2009)
182        Variable/G fwhm_narrow = 0.109
183        Variable/G fwhm_mid = 0.125
184        Variable/G fwhm_wide = 0.236
185       
186        //source apertures (cm)
187        Variable/G a1_0_0 = 1.43
188        Variable/G a1_0_1 = 2.54
189        Variable/G a1_0_2 = 3.81
190        Variable/G a1_7_0 = 2.5 // after the polarizer         
191        Variable/G a1_7_1 = 5.0
192        Variable/G a1_def = 5.00
193       
194        //default configuration values
195//      ng = 0
196//      a1 = 3.81
197//      pos_table = 2
198//      l2r = 1310.0
199//      a2 = 1.27
200//      det_off = 0.0
201//      lambda = 6.0
202//      lambda_width = 0.15
203        Variable/G      l2diff = 5.0
204//     
205        SetDataFolder root:
206end
207
208Function initNG7()
209
210        SetDataFolder root:Packages:NIST:SAS
211       
212        Variable/G instrument = 7
213        Variable/G s12 = 54.8
214        Variable/G d_det = 0.5
215        Variable/G a_pixel = 0.5
216        Variable/G del_r = 0.5
217        Variable/G det_width = 64.0
218        Variable/G phi_0 = 2.3e13
219        Variable/G lambda_t = 5.50
220        Variable/G l2r_lower = 100
221        Variable/G l2r_upper =  1531
222        Variable/G lambda_lower = 4.0
223        Variable/G lambda_upper = 20.0
224        Variable/G d_upper = 25.0
225        Variable/G bs_factor = 1.05
226        Variable/G t1 = 0.63
227        Variable/G t2 = 0.7
228        Variable/G t3 = 0.75
229        Variable/G l_gap = 188.0
230        Variable/G guide_width = 5.0
231        Variable/G idmax = 100.0
232        Variable/G b = 0.028
233        Variable/G c = 0.028
234       
235        //fwhm values (new variables)
236        Variable/G fwhm_narrow = 0.09
237        Variable/G fwhm_mid = 0.11
238        Variable/G fwhm_wide = 0.22
239       
240        //source apertures (cm)
241        Variable/G a1_0_0 = 1.43
242        Variable/G a1_0_1 = 2.22
243        Variable/G a1_0_2 = 3.81
244        Variable/G a1_7_0 = 5.0         // don't apply to NG7
245        Variable/G a1_7_1 = 5.0
246        Variable/G a1_def = 5.00
247       
248        //default configuration values
249//      ng = 0
250//      a1 = 2.22
251//      pos_table = 2
252//      l2r = 1530.0
253//      a2 = 1.27
254//      det_off = 0.0
255//      lambda = 6.0
256//      lambda_width = 0.11
257        Variable/G      l2diff = 5.0
258//     
259        SetDataFolder root:
260end
261
262Function S_fillDefaultHeader(iW,rW,tW)
263        Wave iW,rW
264        Wave/T tW
265
266        // text wave
267        // don't need anything
268       
269        // integer wave
270        // don't need anything
271       
272        // real wave
273        rw = 0
274        rw[16] = 64             // beamcenter X (pixels)
275        rw[17] = 64             // beamcenter Y
276       
277        rw[10]  = 5.08                  //detector resolution (5mm) and calibration constants (linearity)
278        rw[11] = 10000
279        rw[12] = 0
280        rw[13] = 5.08
281        rw[14] = 10000
282        rw[15] = 0
283       
284        rw[20] = 65             // det size in cm
285        rw[18] = 6              // SDD in meters (=L2)
286       
287        rw[26] = 6              //lambda in Angstroms
288        rw[4] = 1               //transmission
289       
290        rw[21] = 76.2                   //BS diameter in mm
291        rw[23] = 50                     //A1 diameter in mm
292        rw[24] = 12.7                   //A2 diameter in mm
293        rw[25] = 7.02                   //L1 distance in meters (derived from number of guides)
294        rw[27] = 0.11                   //DL/L wavelength spread
295       
296        return(0)
297End
298
299
300Window SASCALC_Panel()
301
302        PauseUpdate; Silent 1           // building window...
303
304// if I make the graph a subwindow in a panel, it breaks the "append 1d" from the wrapper       
305//      NewPanel/W=(5,44,463,570)/K=1 as "SASCALC"
306//      DoWindow/C SASCALC
307//      ModifyPanel cbRGB=(49151,53155,65535)
308//
309//     
310//      String fldrSav0= GetDataFolder(1)
311//      SetDataFolder root:Packages:NIST:SAS:
312//     
313//      Display/HOST=#/W=(5,200,463,570) aveint vs qval
314//      ModifyGraph mode=3
315//      ModifyGraph marker=19
316//      ModifyGraph rgb=(0,0,0)
317//      Modifygraph log=1
318//      Modifygraph grid=1
319//      Modifygraph mirror=2
320//      ModifyGraph msize(aveint)=2
321//      ErrorBars/T=0 aveint Y,wave=(sigave,sigave)
322//      Label bottom, "Q (1/A)"
323//      Label left, "Relative Intensity"
324//      legend
325//     
326//      RenameWindow #,G_aveint
327//      SetActiveSubwindow ##
328
329//// end panel commands
330
331/// draw as a graph
332        String fldrSav0= GetDataFolder(1)
333        SetDataFolder root:Packages:NIST:SAS:
334       
335        Display/W=(5,44,463,570)/K=1  aveint vs qval as "SASCALC"
336        DoWindow/C SASCALC
337        ModifyGraph cbRGB=(49151,53155,65535)
338        ModifyGraph mode=3
339        ModifyGraph marker=19
340        ModifyGraph rgb=(0,0,0)
341        Modifygraph log=1
342        Modifygraph grid=1
343        Modifygraph mirror=2
344        ModifyGraph msize(aveint)=2
345        ErrorBars/T=0 aveint Y,wave=(sigave,sigave)
346        Label bottom, "Q (1/A)"
347        Label left, "Relative Intensity"
348        legend
349       
350        ControlBar/T 200
351
352        SetDataFolder fldrSav0
353       
354        Slider SC_Slider,pos={11,46},size={150,45},proc=GuideSliderProc,live=0
355        Slider SC_Slider,limits={0,8,1},variable= root:Packages:NIST:SAS:gNg,vert= 0//,thumbColor= (1,16019,65535)
356        Slider SC_Slider_1,pos={234,45},size={150,45},proc=DetDistSliderProc,live=0
357        Slider SC_Slider_1,tkLblRot= 90,userTicks={root:Packages:NIST:SAS:tickSDDNG3,root:Packages:NIST:SAS:lblSDDNG3 }
358        Slider SC_Slider_1,limits={133,1317,1},variable= root:Packages:NIST:SAS:gDetDist,vert= 0//,thumbColor= (1,16019,65535)
359        Slider SC_Slider_2,pos={394,21},size={47,65},proc=OffsetSliderProc,live=0,ticks=4
360        Slider SC_Slider_2,limits={0,25,1},variable= root:Packages:NIST:SAS:gOffset//,thumbColor= (1,16019,65535)
361        CheckBox checkNG3,pos={20,19},size={36,14},proc=SelectInstrumentCheckProc,title="NG3"
362        CheckBox checkNG3,value=1,mode=1
363        CheckBox checkNG7,pos={66,19},size={36,14},proc=SelectInstrumentCheckProc,title="NG7"
364        CheckBox checkNG7,value=0,mode=1
365        CheckBox checkChamber,pos={172,48},size={57,14},proc=TableCheckProc,title="Chamber"
366        CheckBox checkChamber,value=1,mode=1
367        CheckBox checkHuber,pos={172,27},size={44,14},proc=TableCheckProc,title="Huber"
368        CheckBox checkHuber,value=0,mode=1
369        PopupMenu popup0,pos={6,94},size={76,20},proc=SourceAperturePopMenuProc
370        PopupMenu popup0,mode=1,popvalue="3.81 cm",value= root:Packages:NIST:SAS:gSourceApString
371        PopupMenu popup0_1,pos={172,72},size={49,20},proc=SampleAperturePopMenuProc
372        PopupMenu popup0_1,mode=8,popvalue="1/2\"",value= root:Packages:NIST:SAS:gApPopStr
373        SetVariable setvar0,pos={301,107},size={130,15},title="Det Dist (cm)",proc=SDDSetVarProc
374        SetVariable setvar0,limits={133,1317,1},value=root:Packages:NIST:SAS:gDetDist
375        SetVariable setvar0_1,pos={321,129},size={110,15},title="Offset (cm)",proc=OffsetSetVarProc
376        SetVariable setvar0_1,limits={0,25,1},value= root:Packages:NIST:SAS:gOffset
377        SetVariable setvar0_2,pos={6,130},size={90,15},title="Lambda",proc=LambdaSetVarProc
378        SetVariable setvar0_2,limits={4,20,0.1},value= root:Packages:NIST:SAS:gLambda
379        PopupMenu popup0_2,pos={108,127},size={55,20},proc=DeltaLambdaPopMenuProc
380        PopupMenu popup0_2,mode=2,value= root:Packages:NIST:SAS:gDeltaLambdaStr
381        Button FreezeButton title="Freeze",size={60,20},pos={180,166}
382        Button FreezeButton proc=FreezeButtonProc
383        Button ClearButton title="Clear",size={60,20},pos={250,166}
384        Button ClearButton proc=S_ClearButtonProc
385        GroupBox group0,pos={6,1},size={108,36},title="Instrument"
386        SetDataFolder fldrSav0
387       
388        SetVariable setvar0_3,pos={140,94},size={110,15},title="Diam (mm)",disable=1
389        SetVariable setvar0_3,limits={0,100,0.1},value= root:Packages:NIST:SAS:gSamApOther,proc=SampleApOtherSetVarProc
390       
391        CheckBox checkLens,pos={6,155},size={44,14},proc=LensCheckProc,title="Lenses?"
392        CheckBox checkLens,value=root:Packages:NIST:SAS:gUsingLenses
393       
394        CheckBox checkSim,pos={6,175},size={44,14},proc=SimCheckProc,title="Simulation?"
395        CheckBox checkSim,value=0
396       
397        // help, done buttons
398        Button SC_helpButton,pos={340,166},size={25,20},proc=showSASCALCHelp,title="?"
399        Button SC_helpButton,help={"Show help file for simulation of SANS Data"}
400        Button SC_DoneButton,pos={380,166},size={50,20},proc=SASCALCDoneButton,title="Done"
401        Button SC_DoneButton,help={"This button will close the panel"}
402
403        // set up a fake dependency to trigger recalculation
404        //root:Packages:NIST:SAS:gCalculate := ReCalculateInten(root:Packages:NIST:SAS:gTouched)
405EndMacro
406
407//clean up
408Proc SASCALCDoneButton(ctrlName): ButtonControl
409        String ctrlName
410        DoWindow/K SASCALC
411        DoWindow/K Trial_Configuration
412        DoWindow/K Saved_Configurations
413        DoWindow/K MC_SASCALC
414        DoWindow/K Sim_1D_Panel
415end
416
417//
418Proc showSASCALCHelp(ctrlName): ButtonControl
419        String ctrlName
420        DisplayHelpTopic/K=1/Z "SASCALC"
421        if(V_flag !=0)
422                DoAlert 0,"The SANS Simulation Help file could not be found"
423        endif
424end
425
426// based on the instrument selected...
427//set the SDD range
428// set the source aperture popup (based on NGx and number of guides)
429// set the wavelength spread popup
430//
431Function UpdateControls()
432        //poll the controls on the panel, and change needed values
433        Variable isNG3,Ng,mode
434        ControlInfo/W=SASCALC checkNG3
435        isNG3=V_value
436        ControlInfo/W=SASCALC SC_slider
437        Ng=V_value
438        SVAR A1str= root:Packages:NIST:SAS:gSourceApString// = "1.43 cm;2.54 cm;3.81 cm;"
439        SVAR dlStr = root:Packages:NIST:SAS:gDeltaLambdaStr
440        if(isNG3)
441                switch(ng)     
442                        case 0:
443                                ControlInfo/W=SASCALC popup0
444                                mode=V_value
445                                A1str="1.43 cm;2.54 cm;3.81 cm;"
446                                break
447                        case 6:
448                                A1str = "! 6 Guides invalid;"
449                                mode=1
450                                break
451                        case 7:                                                 // switched order in 2009 to keep 5 cm as default, 2.5 cm for polarizer
452                                A1Str = "5.00 cm;2.50 cm;"
453                                mode = 1
454                                break
455                        default:
456                                A1str = "5 cm;"
457                                mode=1
458                endswitch
459                //wavelength spread
460                dlStr = "0.109;0.125;0.236;"            //updated calibration 2009
461                //detector limits
462                SetVariable setvar0,limits={133,1317,1}
463                NVAR detDist=root:Packages:NIST:SAS:gDetDist
464                if(detDist < 133 )
465                        detDist = 133
466                elseif (detDist > 1317 )
467                        detDist = 1317
468                endif
469                Slider SC_Slider_1,limits={133,1317,1},userTicks={root:Packages:NIST:SAS:tickSDDNG3,root:Packages:NIST:SAS:lblSDDNG3 }
470                Slider SC_Slider_1,variable=root:Packages:NIST:SAS:gDetDist             //forces update
471        else                    //ng7
472                switch(ng)     
473                        case 0:
474                                ControlInfo/W=SASCALC popup0
475                                mode=V_value
476                                A1str="1.43 cm;2.22 cm;3.81 cm;"
477                                break
478                        default:
479                                A1str = "5.08 cm;"
480                                mode=1
481                endswitch
482               
483                dlStr = "0.09;0.11;0.22;"
484                Slider SC_Slider_1,limits={100,1531,1},userTicks={root:Packages:NIST:SAS:tickSDDNG7,root:Packages:NIST:SAS:lblSDDNG7 }
485                SetVariable setvar0,limits={100,1531,1}
486                Slider SC_Slider_1,variable=root:Packages:NIST:SAS:gDetDist             //forces update
487        endif
488        ControlUpdate popup0
489        PopupMenu popup0,mode=mode              //source Ap
490        ControlInfo/W=SASCALC popup0
491        SourceAperturePopMenuProc("",0,S_Value)                 //send popNum==0 so recalculation won't be done
492       
493        ControlUpdate popup0_2          // delta lambda pop
494        ControlInfo/W=SASCALC popup0_2
495        DeltaLambdaPopMenuProc("",0,S_Value)            //sets the global and the wave
496End
497
498//changing the number of guides changes the SSD
499// the source aperture popup may need to be updated
500//
501Function GuideSliderProc(ctrlName,sliderValue,event)
502        String ctrlName
503        Variable sliderValue
504        Variable event  // bit field: bit 0: value set, 1: mouse down, 2: mouse up, 3: mouse moved
505       
506        Variable recalc=0
507       
508        if(event %& 0x1)        // bit 0, value set
509                if(cmpstr(ctrlName,"") != 0)            //here by direct action, so do LensCheck and recalculate
510                        recalc=1
511                        LensCheckProc("",2)             //make sure lenses are deselected
512                endif
513                sourceToSampleDist()            //updates the SSD global and wave
514                //change the sourceAp popup, SDD range, etc
515                UpdateControls()
516                ReCalculateInten(recalc)
517        endif
518        return 0
519End
520
521// changing the detector position changes the SDD
522//
523Function DetDistSliderProc(ctrlName,sliderValue,event)
524        String ctrlName
525        Variable sliderValue
526        Variable event  // bit field: bit 0: value set, 1: mouse down, 2: mouse up, 3: mouse moved
527
528        Variable recalc=0
529        if(event %& 0x1)        // bit 0, value set
530                if(cmpstr(ctrlName,"") != 0)
531                        recalc=1
532                        LensCheckProc("",2)             //make sure lenses are only selected for valid configurations
533                endif
534                sampleToDetectorDist()          //changes the SDD and wave (DetDist is the global)
535                ReCalculateInten(recalc)
536        endif
537
538        return 0
539End
540
541// change the offset
542// - changes the beamcenter (x,y) position too
543Function OffsetSliderProc(ctrlName,sliderValue,event)
544        String ctrlName
545        Variable sliderValue
546        Variable event  // bit field: bit 0: value set, 1: mouse down, 2: mouse up, 3: mouse moved
547
548        if(event %& 0x1)        // bit 0, value set
549                detectorOffset()
550                ReCalculateInten(1)
551        endif
552
553        return 0
554End
555
556// changing the instrument -
557// re-initialize the global values
558// update the controls to appropriate limits/choices
559//
560Function SelectInstrumentCheckProc(ctrlName,checked) : CheckBoxControl
561        String ctrlName
562        Variable checked
563
564        if(cmpstr(ctrlName,"checkNG3")==0)
565                checkBox checkNG3, value=1
566                checkBox checkNG7, value=0
567                initNG3()
568        else
569                checkBox checkNG3, value=0
570                checkBox checkNG7, value=1
571                initNG7()
572        endif
573        LensCheckProc("",2)             //check if lenses are still valid (they won't be)
574        UpdateControls()
575        ReCalculateInten(1)
576End
577
578//changing the table position
579// changes the SSD and SDD
580Function TableCheckProc(ctrlName,checked) : CheckBoxControl
581        String ctrlName
582        Variable checked
583
584        NVAR table=root:Packages:NIST:SAS:gTable
585        if(cmpstr(ctrlName,"checkHuber")==0)
586                checkBox checkHuber, value=1
587                checkBox checkChamber, value=0
588                table=1         //in Huber position
589        else
590                checkBox checkHuber, value=0
591                checkBox checkChamber, value=1
592                table = 2               //in Sample chamber
593        endif
594        sampleToDetectorDist()
595        sourceToSampleDist()            //update
596        ReCalculateInten(1)
597End
598
599
600//lenses (or prisms) in/out changes resolution
601//
602// passing in a checked == 2 will do a check of the configuration without
603// affecting the box state
604//
605// if an invalid configuration is detected, the box is unchecked
606// and the lens flag is set to zero (=out).
607//
608// When necessary controls are "popped", a ctrlName="" is passed to signal the control
609// to NOT recalculate the intensity, and to NOT (recursively) call LensCheckProc again
610//
611// currently, the 17.2 A for lens/prism @ ng3 must be typed in
612//
613Function LensCheckProc(ctrlName,checked) : CheckBoxControl
614        String ctrlName
615        Variable checked
616
617        // don't let the box get checked if the conditions are wrong
618        // lambda != 8.09,8.4,17.2
619        // Ng != 0
620        NVAR lens = root:Packages:NIST:SAS:gUsingLenses
621        NVAR ng = root:Packages:NIST:SAS:gNg
622        NVAR lam = root:Packages:NIST:SAS:gLambda
623        NVAR dist = root:Packages:NIST:SAS:gDetDist
624        NVAR instrument = root:Packages:NIST:SAS:instrument
625        Wave rw=root:Packages:NIST:SAS:realsRead
626       
627        // directly uncheck the box, just set the flag and get out
628        if(checked == 0)
629                lens = 0
630                CheckBox checkLens,value=0
631                rw[28]=0                //flag for lenses out
632                ReCalculateInten(1)
633                return(0)
634        endif
635       
636        // check the box, enforce the proper conditions
637        if(checked == 1)
638                lens = checked
639                if(instrument == 3)
640                        dist = 1317
641                        DetDistSliderProc("",1317,1)
642                       
643                        lam = 8.4
644                        LambdaSetVarProc("",8.4,"8.4","")
645
646                        ng=0
647                        GuideSliderProc("",0,1)         //this updates the controls to the new # of guides
648                       
649                        PopupMenu popup0,mode=1,popvalue="1.43 cm"              //first item in source aperture menu
650                       
651                        PopupMenu popup0_2,mode=2               //deltaLambda
652                        ControlInfo popup0_2
653                        DeltaLambdaPopMenuProc("",0,S_value)                    //zero as 2nd param skips recalculation
654                else
655                        dist = 1531
656                        DetDistSliderProc("",1531,1)
657                       
658                        lam = 8.09
659                        LambdaSetVarProc("",8.09,"8.09","")
660                       
661                        ng=0
662                        GuideSliderProc("",0,1)
663                        PopupMenu popup0,mode=1,popvalue="1.43 cm"              //first item
664                       
665                        PopupMenu popup0_2,mode=2               //deltaLambda
666                        ControlInfo popup0_2
667                        DeltaLambdaPopMenuProc("",0,S_value)                    //zero as 2nd param skips recalculation
668                endif
669                rw[28]=1                //flag for lenses in (not the true number, but OK)
670                ReCalculateInten(1)
671        endif
672       
673        // this is my internal check to see if conditions are still valid
674        // I'll uncheck as needed
675        if(checked == 2)
676
677                // source aperture must be 1.43 cm diameter
678                // number of guides must be zero
679                Variable a1 = sourceApertureDiam()
680                if(a1 != 1.43  || Ng !=0)
681                        lens = 0
682                        CheckBox checkLens,value=0
683                        rw[28]=0                //flag for lenses out
684                        return(0)
685                endif
686       
687                // instrument specific distance requirements
688                if(instrument == 3 && dist != 1317)
689                        lens = 0
690                        CheckBox checkLens,value=0
691                        rw[28]=0                //flag for lenses out
692                        return(0)
693                endif
694       
695                if(instrument == 7 && dist != 1531)
696                        lens = 0
697                        CheckBox checkLens,value=0
698                        rw[28]=0                //flag for lenses out
699                        return(0)
700                endif
701               
702                // instrument specific wavelength requirements
703                if(instrument == 3 && !(lam == 8.4 || lam == 17.2) )
704                        lens = 0
705                        CheckBox checkLens,value=0
706                        rw[28]=0                //flag for lenses out
707                        return(0)
708                endif
709               
710                if(instrument == 7 && lam != 8.09 )
711                        lens = 0
712                        CheckBox checkLens,value=0
713                        rw[28]=0                //flag for lenses out
714                        return(0)
715                endif
716               
717        endif
718
719        return(1)               //return value not used
720End
721
722//simulation controls as a control bar that toggles on/off to the right
723//
724// depending on the state of the 2D flag, open the 1d or 2d control panel
725Function SimCheckProc(ctrlName,checked) : CheckBoxControl
726        String ctrlName
727        Variable checked
728
729        if(checked)
730                NVAR do2D = root:Packages:NIST:SAS:gDoMonteCarlo
731               
732                if(do2D)
733                        DoWindow/F MC_SASCALC
734                        if(V_flag==0)
735                                Execute "MC_SASCALC()"          //sets the variable
736                                AutoPositionWindow/M=1/R=SASCALC MC_SASCALC
737                        endif
738                else
739                        //draw 1D panel
740                        DoWindow/F Sim_1D_Panel
741                        if(V_flag==0)
742                                Execute "Sim_1D_Panel()"
743                                AutoPositionWindow/M=1/R=SASCALC Sim_1D_Panel
744                        endif
745                endif
746                NVAR doSim = root:Packages:NIST:SAS:doSimulation
747                doSim=1
748        else
749                //get rid of the controls (just try both)
750                DoWindow MC_SASCALC
751                if(V_flag!=0)
752                        NVAR do2D = root:Packages:NIST:SAS:gDoMonteCarlo
753                        do2D = 0
754                        KillWindow MC_SASCALC
755                endif
756                DoWindow Sim_1D_Panel
757                if(V_flag!=0)
758                        KillWindow Sim_1D_Panel
759                endif
760               
761                NVAR doSim = root:Packages:NIST:SAS:doSimulation
762                doSim=0
763        endif
764
765        return(0)
766End
767
768// change the source aperture
769//
770Function SourceAperturePopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
771        String ctrlName
772        Variable popNum
773        String popStr
774
775        Variable a1 = sourceApertureDiam()              //sets the new value in the wave
776       
777        // if LensCheckProc is calling this, don't call it again.
778        // if this is called by a direct pop, then check
779        if(cmpstr(ctrlName,"") != 0)
780                LensCheckProc("",2)             //make sure lenses are only selected for valid configurations
781        endif
782       
783        ReCalculateInten(popnum)                //skip the recalculation if I pass in a zero
784End
785
786// set sample aperture
787//
788Function SampleAperturePopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
789        String ctrlName
790        Variable popNum
791        String popStr
792
793        sampleApertureDiam()
794        ReCalculateInten(1)
795End
796
797// set sample to detector distance
798//
799Function SDDSetVarProc(ctrlName,varNum,varStr,varName) : SetVariableControl
800        String ctrlName
801        Variable varNum
802        String varStr
803        String varName
804       
805        sampleToDetectorDist()
806        if(cmpstr(ctrlName,"") != 0)
807                LensCheckProc("",2)             //make sure lenses are only selected for valid configurations
808        endif
809        ReCalculateInten(1)
810End
811
812// set offset
813// -- also changes the beamcenter (x,y)
814//
815Function OffsetSetVarProc(ctrlName,varNum,varStr,varName) : SetVariableControl
816        String ctrlName
817        Variable varNum
818        String varStr
819        String varName
820
821        detectorOffset()                //sets the offset in the wave and also the new (x,y) beamcenter
822       
823        ReCalculateInten(1)
824        return(0)
825End
826
827// change the wavelength
828Function LambdaSetVarProc(ctrlName,varNum,varStr,varName) : SetVariableControl
829        String ctrlName
830        Variable varNum
831        String varStr
832        String varName
833
834        WAVE rw=root:Packages:NIST:SAS:realsRead
835        Variable recalc=0
836       
837        rw[26] = str2num(varStr)
838        if(cmpstr(ctrlName,"") != 0)
839                recalc=1
840                LensCheckProc("",2)             //make sure lenses are only selected for valid configurations
841        endif
842        ReCalculateInten(recalc)
843        return(0)
844End
845
846
847// change the wavelength spread
848Function DeltaLambdaPopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
849        String ctrlName
850        Variable popNum
851        String popStr
852
853        NVAR dl=root:Packages:NIST:SAS:gDeltaLambda
854        dl=str2num(popStr)
855        WAVE rw=root:Packages:NIST:SAS:realsRead
856        rw[27] = dl
857        ReCalculateInten(popnum)                //skip the calculation if I pass in  zero
858        return(0)
859End
860
861
862// calculate the intensity
863// - either do MC or the straight calculation.
864//
865// *****currently the smeared calculation is turned off, by not prepending "fSmeared" to the FUNCREFs)
866// --- the random deviate can't be calculated for the smeared model since I don't know the resolution
867// function over an infinite q-range, just the detector. Maybe the interpolation is OK, but I
868// don't really have a good way of testing this. Also, the resolution calculation explicitly multiplies
869// by fShad, and this wrecks the random deviate calculation. The 2D looks great, but the probabilities
870// are all wrong. fShad is only appropriate post-simulation.
871//
872//
873// (--NO--) ALWAYS DOES THE RESOLUTION SMEARED CALCULATION
874//  even though the unsmeared model is plotted. this is more realistic to
875//  present to users that are planning base on what they will see in an experiment.
876//
877// some bits of the calculation are in a root:Simulation folder that are needed for the resolution smeared calculation
878// all other bits (as possible) are in the SAS folder (a 2D work folder)
879//
880// passing in one does the calculation, "normal" or MC, depending on the global. Normal calculation is the default
881// passing in zero from a control skips the calculation
882//
883Function ReCalculateInten(doIt)
884        Variable doIt
885       
886        if(doIt==0)                     
887                return(0)
888        endif
889       
890        // update the wave with the beamstop diameter here, since I don't know what
891        // combinations of parameters will change the BS - but anytime the curve is
892        // recalculated, or the text displayed, the right BS must be present
893        beamstopDiam()
894       
895        // generate the resolution waves first, so they are present for a smearing calculation
896        // average the "fake" 2d data now to generate the smearing information
897        S_CircularAverageTo1D("SAS")
898        WAVE aveint=root:Packages:NIST:SAS:aveint
899        WAVE qval=root:Packages:NIST:SAS:qval
900        WAVE sigave=root:Packages:NIST:SAS:sigave
901        WAVE SigmaQ=root:Packages:NIST:SAS:sigmaQ
902        WAVE qbar=root:Packages:NIST:SAS:qbar
903        WAVE fSubS=root:Packages:NIST:SAS:fSubS
904       
905        //generate a "fake" 1d data folder/set named "Simulation"
906        Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,"Simulation")
907       
908        // do the simulation here, or not
909        Variable r1,xCtr,yCtr,sdd,pixSize,wavelength
910        String coefStr,abortStr,str
911
912        // now the cases are: simulation (0|1), 1D (default) or 2D (hidden)
913        NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo                // == 1 if 2D MonteCarlo set by hidden flag
914        NVAR doSimulation = root:Packages:NIST:SAS:doSimulation         // == 1 if 1D simulated data, 0 if other from the checkbox
915        SVAR funcStr = root:Packages:NIST:SAS:gFuncStr          //set by the popup
916
917        if(doSimulation == 1)
918                if(doMonteCarlo == 1)
919                        //2D simulation (in MultiScatter_MonteCarlo_2D.ipf)
920                       
921                        Simulate_2D_MC(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
922                       
923                        //end 2D simulation
924                else
925                        //1D simulation
926                       
927                        if(exists(funcStr) != 0)
928                                FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr)                 //a wrapper for the structure version
929                                FUNCREF SANSModelAAO_MCproto funcUnsmeared=$(funcStr)           //unsmeared
930                                coefStr = MC_getFunctionCoef(funcStr)
931                               
932                                if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
933                                        Abort "Function and coefficients do not match. You must plot the unsmeared function before simulation."
934                                endif
935                               
936                                Wave inten=$"root:Simulation:Simulation_i"              // this will exist and send the smeared calculation to the corect DF
937                               
938                                // the resolution-smeared intensity is calculated, including the incoherent background
939                                func($coefStr,inten,qval)
940
941                                NVAR imon = root:Packages:NIST:SAS:gImon
942                                NVAR ctTime = root:Packages:NIST:SAS:gCntTime
943                                NVAR thick = root:Packages:NIST:SAS:gThick
944                                NVAR trans = root:Packages:NIST:SAS:gSamTrans
945                                NVAR SimDetCts = root:Packages:NIST:SAS:g_1DTotCts                      //summed counts (simulated)
946                                NVAR estDetCR = root:Packages:NIST:SAS:g_1DEstDetCR                     // estimated detector count rate
947                                NVAR fracScat = root:Packages:NIST:SAS:g_1DFracScatt            // fraction of beam captured on detector
948                                NVAR estTrans = root:Packages:NIST:SAS:g_1DEstTrans             // estimated transmission of sample
949                                NVAR SimCountTime = root:Packages:NIST:SAS:gCntTime             //counting time used for simulation
950                               
951                                WAVE rw=root:Packages:NIST:SAS:realsRead
952                                WAVE nCells=root:Packages:NIST:SAS:nCells                               
953                                                               
954                                pixSize = rw[10]/10             // convert pix size in mm to cm
955                                sdd = rw[18]*100                //convert header of [m] to [cm]
956                                wavelength = rw[26]             // in 1/A
957                               
958                                imon = beamIntensity()
959                               
960                                // calculate the scattering cross section simply to be able to estimate the transmission
961                                Variable sig_sas=0
962                               
963                                // remember that the random deviate is the coherent portion ONLY - the incoherent background is
964                                // subtracted before the calculation.
965                                CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
966                               
967//                              if(sig_sas > 100)
968//                                      sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas
969//                              endif
970                                estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
971                                Print "Sig_sas = ",sig_sas
972                               
973                                Duplicate/O qval prob_i,countsInAnnulus
974                               
975                                // not needed - nCells takes care of this when the error is correctly calculated
976//                              Duplicate/O qval circle_fraction,rval,nCells_expected
977//                              rval = sdd*tan(2*asin(qval*wavelength/4/pi))            //radial distance in cm
978//                              nCells_expected = 2*pi*rval/pixSize                                     //does this need to be an integer?
979//                              circle_fraction = nCells / nCells_expected
980                               
981                                                       
982//                              prob_i = trans*thick*nCells*(pixSize/sdd)^2*inten                       //probability of a neutron in q-bin(i) that has nCells
983                                prob_i = trans*thick*(pixSize/sdd)^2*inten                      //probability of a neutron in q-bin(i)
984                               
985                                Variable P_on = sum(prob_i,-inf,inf)
986                                Print "P_on = ",P_on
987                               
988//                              fracScat = P_on
989                                fracScat = 1-estTrans
990                               
991//                              aveint = (Imon*ctTime)*prob_i / circle_fraction / nCells_expected
992                                aveint = (Imon*ctTime)*prob_i
993
994                                countsInAnnulus = aveint*nCells
995                                SimDetCts = sum(countsInAnnulus,-inf,inf)
996                                estDetCR = SimDetCts/SimCountTime
997                               
998                               
999                                NVAR doABS = root:Packages:NIST:SAS:g_1D_DoABS
1000                                NVAR addNoise = root:Packages:NIST:SAS:g_1D_AddNoise
1001                                                       
1002                                sigave = sqrt(aveint/nCells)            // corrected based on John's memo, from 8/9/99
1003                               
1004                                // add in random error in aveint based on the sigave
1005                                if(addNoise)
1006                                        aveint += gnoise(sigave)
1007                                endif
1008
1009                                // convert to absolute scale
1010                                if(doABS)
1011                                        Variable kappa = thick*(pixSize/sdd)^2*trans*iMon*ctTime
1012                                        aveint /= kappa
1013                                        sigave /= kappa
1014                                endif
1015                               
1016                        else
1017                                //no function plotted, no simulation can be done
1018                                DoAlert 0,"No function is selected or plotted, so no simulation is done. The default Debye function is used."
1019               
1020                                aveint = S_Debye(1000,100,0.0,qval)
1021                                aveint *= fSubS         // multiply either estimate by beamstop shadowing
1022                                sigave = 0              //reset for model calculation
1023                        endif
1024                       
1025                endif // end 1D simulation
1026        else
1027                //no simulation
1028               
1029                aveint = S_Debye(1000,100,0.0,qval)
1030                aveint *= fSubS         // multiply either estimate by beamstop shadowing
1031                sigave = 0              //reset for model calculation
1032               
1033                //end no simulation     
1034        endif
1035       
1036
1037        //display the configuration text in a separate notebook
1038        DisplayConfigurationText()
1039       
1040        return(0)
1041End
1042
1043//freezes the current configuration
1044// -1- duplicates the trace on the graph
1045// -2- copies the configuration text to a second notebook window for printing
1046Function FreezeButtonProc(ctrlName) : ButtonControl
1047        String ctrlName
1048
1049        String str=""
1050        NVAR ct=root:Packages:NIST:SAS:gFreezeCount
1051       
1052       
1053        SetDataFolder root:Packages:NIST:SAS
1054       
1055        Duplicate/O aveint,$("aveint_"+num2str(ct))
1056        Duplicate/O qval,$("qval_"+num2str(ct))
1057        Duplicate/O sigave,$("sigave_"+num2str(ct))
1058        Appendtograph $("aveint_"+num2str(ct)) vs $("qval_"+num2str(ct))
1059        ModifyGraph mode($("aveint_"+num2str(ct)))=3
1060        ModifyGraph marker($("aveint_"+num2str(ct)))=19
1061        ModifyGraph msize($("aveint_"+num2str(ct)))=2
1062        ErrorBars/T=0 $("aveint_"+num2str(ct)) Y,wave=($("sigave_"+num2str(ct)),$("sigave_"+num2str(ct)))
1063       
1064        switch(mod(ct,10))      // 10 different colors - black is the unfrozen color
1065                case 0:
1066                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,16385,16385)
1067                        break
1068                case 1:
1069                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(2,39321,1)
1070                        break
1071                case 2:
1072                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(0,0,65535)
1073                        break
1074                case 3:
1075                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(39321,1,31457)
1076                        break
1077                case 4:
1078                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(48059,48059,48059)
1079                        break
1080                case 5:
1081                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,32768,32768)
1082                        break
1083                case 6:
1084                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(0,65535,0)
1085                        break
1086                case 7:
1087                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(16385,65535,65535)
1088                        break
1089                case 8:
1090                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(65535,32768,58981)
1091                        break
1092                case 9:
1093                        ModifyGraph rgb($("aveint_"+num2str(ct)))=(36873,14755,58982)
1094                        break
1095        endswitch
1096       
1097        NVAR doTraceOffset = root:Packages:NIST:SAS:gDoTraceOffset
1098        NVAR offset = root:Packages:NIST:SAS:gModelOffsetFactor
1099        if(doTraceOffset)
1100                offset = 2^ct
1101                //multiply by current offset (>=1)
1102                Wave inten = $("aveint_"+num2str(ct))
1103                inten *= offset
1104                //      Print "new offset = ",offset
1105        endif
1106       
1107        ct +=1
1108        SetDataFolder root:
1109       
1110        // create or append the configuration to a new notebook
1111        if(WinType("Saved_Configurations")==0)
1112                NewNotebook/F=1/N=Saved_Configurations /W=(480,400,880,725)
1113                DoWindow/F SASCALC              //return focus to SASCALC window
1114        endif
1115        //append the text
1116        sprintf str,"\rConfiguration #%d\r",ct-1
1117        Notebook Saved_Configurations showRuler=0,defaultTab=20,selection={endOfFile, endOfFile}
1118        Notebook Saved_Configurations font="Monaco",fSize=10,fstyle=1,text=str          //bold
1119        Notebook Saved_Configurations font="Monaco",fSize=10,fstyle=0,text=SetConfigurationText()
1120        return(0)
1121End
1122
1123//clears the frozen traces on the graph, asks if you want to clear the saved text
1124//
1125Function S_ClearButtonProc(ctrlName) : ButtonControl
1126        String ctrlName
1127
1128        NVAR ct=root:Packages:NIST:SAS:gFreezeCount
1129        Variable ii
1130        Setdatafolder root:Packages:NIST:SAS
1131        for(ii=ct-1;ii>=1;ii-=1)
1132        //remove all traces, replace aveint
1133        // kill all waves _ct
1134                RemoveFromGraph $("aveint_"+num2str(ii))
1135                Killwaves/Z $("aveint_"+num2str(ii)),$("qval_"+num2str(ii))
1136        endfor
1137        ct=1
1138        setdatafolder root:
1139       
1140        DoAlert 1,"Do you also want to clear the \"Saved Configurations\" text?"
1141        if(V_flag == 1)         // yes
1142                DoWindow/K Saved_Configurations
1143        endif
1144       
1145        //reset offset value
1146        NVAR offset = root:Packages:NIST:SAS:gModelOffsetFactor
1147        offset = 1
1148        ReCalculateInten(1)
1149        return(0)
1150End
1151
1152
1153///////////////////////////////////////////////////////////
1154// 19MAR07 uses correction for beamstop diameter projection to get shadow factor correct
1155//
1156Function S_CircularAverageTo1D(type)
1157        String type
1158       
1159        Variable isCircular = 1
1160       
1161        //type is the data type to do the averaging on, and will be set as the current folder
1162        //get the current displayed data (so the correct folder is used)
1163        String destPath = "root:Packages:NIST:"+type
1164       
1165        //
1166        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist,dr,ddr
1167        Variable lambda,trans
1168        WAVE reals = $(destPath + ":RealsRead")
1169//      WAVE/T textread = $(destPath + ":TextRead")
1170//      String fileStr = textread[3]
1171       
1172        // center of detector, for non-linear corrections
1173        Variable pixelsX=128,pixelsY=128
1174       
1175        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
1176        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
1177       
1178        // beam center, in pixels
1179        x0 = reals[16]
1180        y0 = reals[17]
1181        //detector calibration constants
1182        sx = reals[10]          //mm/pixel (x)
1183        sx3 = reals[11]         //nonlinear coeff
1184        sy = reals[13]          //mm/pixel (y)
1185        sy3 = reals[14]         //nonlinear coeff
1186       
1187        dtsize = 10*reals[20]           //det size in mm
1188        dtdist = 1000*reals[18]         // det distance in mm
1189       
1190        NVAR binWidth=root:Packages:NIST:SAS:gBinWidth
1191//      Variable binWidth = 1
1192       
1193        dr = binWidth           // ***********annulus width set by user, default is one***********
1194        ddr = dr*sx             //step size, in mm (this value should be passed to the resolution calculation, not dr 18NOV03)
1195               
1196        Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,dxbm,dybm,ii
1197        Variable phi_rad,dphi_rad,phi_x,phi_y
1198        Variable forward,mirror
1199       
1200        String side = "both"
1201//      side = StringByKey("SIDE",keyListStr,"=",";")
1202//      Print "side = ",side
1203       
1204        /// data wave is data in the current folder which was set at the top of the function
1205        WAVE data=$(destPath + ":data")
1206
1207// fake mask that uses all of the detector
1208        Make/D/O/N=(pixelsX,pixelsY) $(destPath + ":mask")
1209        Wave mask = $(destPath + ":mask")
1210        mask = 0
1211        //two pixels all around
1212        mask[0,1][] = 1
1213        mask[126,127][] = 1
1214        mask[][0,1] = 1
1215        mask[][126,127] = 1
1216        //
1217        //pixels within rcentr of beam center are broken into 9 parts (units of mm)
1218        rcentr = 100            //original
1219//      rcentr = 0
1220        // values for error if unable to estimate value
1221        //large_num = 1e10
1222        large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1
1223        small_num = 1e-10
1224       
1225        // output wave are expected to exist (?) initialized to zero, what length?
1226        // 200 points on VAX --- use 300 here, or more if SAXS data is used with 1024x1024 detector (1000 pts seems good)
1227        Variable defWavePts=500
1228        Make/O/D/N=(defWavePts) $(destPath + ":qval"),$(destPath + ":aveint")
1229        Make/O/D/N=(defWavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave")
1230        Make/O/D/N=(defWavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar")
1231
1232        WAVE qval = $(destPath + ":qval")
1233        WAVE aveint = $(destPath + ":aveint")
1234        WAVE ncells = $(destPath + ":ncells")
1235        WAVE dsq = $(destPath + ":dsq")
1236        WAVE sigave = $(destPath + ":sigave")
1237        WAVE qbar = $(destPath + ":QBar")
1238        WAVE sigmaq = $(destPath + ":SigmaQ")
1239        WAVE fsubs = $(destPath + ":fSubS")
1240       
1241        qval = 0
1242        aveint = 0
1243        ncells = 0
1244        dsq = 0
1245        sigave = 0
1246        qbar = 0
1247        sigmaq = 0
1248        fsubs = 0
1249
1250        dtdis2 = dtdist^2
1251        nq = 1
1252        xoffst=0
1253        //distance of beam center from detector center
1254        dxbm = S_FX(x0,sx3,xcenter,sx)
1255        dybm = S_FY(y0,sy3,ycenter,sy)
1256               
1257        //BEGIN AVERAGE **********
1258        Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1
1259        Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p
1260       
1261        // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
1262        // loop index corresponds to FORTRAN (old code)
1263        // and the IGOR array indices must be adjusted (-1) to the correct address
1264        ii=1
1265        do
1266                xi = ii
1267                dxi = S_FX(xi,sx3,xcenter,sx)
1268                dx = dxi-dxbm           //dx and dy are in mm
1269               
1270                jj = 1
1271                do
1272                        data_pixel = data[ii-1][jj-1]           //assign to local variable
1273                        yj = jj
1274                        dyj = S_FY(yj,sy3,ycenter,sy)
1275                        dy = dyj - dybm
1276                        if(!(mask[ii-1][jj-1]))                 //masked pixels = 1, skip if masked (this way works...)
1277                                dr2 = (dx^2 + dy^2)^(0.5)               //distance from beam center NOTE dr2 used here - dr used above
1278                                if(dr2>rcentr)          //keep pixel whole
1279                                        nd = 1
1280                                        fd = 1
1281                                else                            //break pixel into 9 equal parts
1282                                        nd = 3
1283                                        fd = 2
1284                                endif
1285                                nd2 = nd^2
1286                                ll = 1          //"el-el" loop index
1287                                do
1288                                        dxx = dx + (ll - fd)*sx/3
1289                                        kk = 1
1290                                        do
1291                                                dyy = dy + (kk - fd)*sy/3
1292                                                if(isCircular)
1293                                                        //circular average, use all pixels
1294                                                        //(increment)
1295                                                        nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1296                                                else
1297                                                        //a sector average - determine azimuthal angle
1298                                                        dphi_p = S_dphi_pixel(dxx,dyy,phi_x,phi_y)
1299                                                        if(dphi_p < dphi_rad)
1300                                                                forward = 1                     //within forward sector
1301                                                        else
1302                                                                forward = 0
1303                                                        Endif
1304                                                        if((Pi - dphi_p) < dphi_rad)
1305                                                                mirror = 1              //within mirror sector
1306                                                        else
1307                                                                mirror = 0
1308                                                        Endif
1309                                                        //check if pixel lies within allowed sector(s)
1310                                                        if(cmpstr(side,"both")==0)              //both sectors
1311                                                                if ( mirror || forward)
1312                                                                        //increment
1313                                                                        nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1314                                                                Endif
1315                                                        else
1316                                                                if(cmpstr(side,"right")==0)             //forward sector only
1317                                                                        if(forward)
1318                                                                                //increment
1319                                                                                nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1320                                                                        Endif
1321                                                                else                    //mirror sector only
1322                                                                        if(mirror)
1323                                                                                //increment
1324                                                                                nq = S_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1325                                                                        Endif
1326                                                                Endif
1327                                                        Endif           //allowable sectors
1328                                                Endif   //circular or sector check
1329                                                kk+=1
1330                                        while(kk<=nd)
1331                                        ll += 1
1332                                while(ll<=nd)
1333                        Endif           //masked pixel check
1334                        jj += 1
1335                while (jj<=pixelsY)
1336                ii += 1
1337        while(ii<=pixelsX)              //end of the averaging
1338               
1339        //compute q-values and errors
1340        Variable ntotal,rr,theta,avesq,aveisq,var
1341       
1342        lambda = reals[26]
1343        ntotal = 0
1344        kk = 1
1345        do
1346                rr = (2*kk-1)*ddr/2
1347                theta = 0.5*atan(rr/dtdist)
1348                qval[kk-1] = (4*Pi/lambda)*sin(theta)
1349                if(ncells[kk-1] == 0)
1350                        //no pixels in annuli, data unknown
1351                        aveint[kk-1] = 0
1352                        sigave[kk-1] = large_num
1353                else
1354                        if(ncells[kk-1] <= 1)
1355                                //need more than one pixel to determine error
1356                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
1357                                sigave[kk-1] = large_num
1358                        else
1359                                //assume that the intensity in each pixel in annuli is normally
1360                                // distributed about mean...
1361                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
1362                                avesq = aveint[kk-1]^2
1363                                aveisq = dsq[kk-1]/ncells[kk-1]
1364                                var = aveisq-avesq
1365                                if(var<=0)
1366                                        sigave[kk-1] = small_num
1367                                else
1368                                        sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1))
1369                                endif
1370                        endif
1371                endif
1372                ntotal += ncells[kk-1]
1373                kk+=1
1374        while(kk<=nq)
1375       
1376        //Print "NQ = ",nq
1377        // data waves were defined as 300 points (=defWavePts), but now have less than that (nq) points
1378        // use DeletePoints to remove junk from end of waves
1379        //WaveStats would be a more foolproof implementation, to get the # points in the wave
1380        Variable startElement,numElements
1381        startElement = nq
1382        numElements = defWavePts - startElement
1383        DeletePoints startElement,numElements, qval,aveint,ncells,dsq,sigave
1384       
1385        //////////////end of VAX sector_ave()
1386               
1387
1388
1389// ***************************************************************
1390//
1391// Do the extra 3 columns of resolution calculations starting here.
1392//
1393// ***************************************************************
1394
1395        Variable L2 = reals[18]
1396        Variable BS = reals[21]         //this the diameter is stored in mm
1397//  SRK - why was I using the projected diameter of the beam stop?? I added a step at the beginning of every recalculation
1398// of the intensity to get the right beamstop diameter into RealsRead...
1399//      Variable BS = beamstopDiamProjection(1) * 10            //calculated projection in cm *10 = mm
1400        Variable S1 = reals[23]
1401        Variable S2 = reals[24]
1402        Variable L1 = reals[25]
1403        lambda = reals[26]
1404        Variable lambdaWidth = reals[27]
1405
1406        Variable DDet, apOff
1407        //typical value for NG3 and NG7 - distance between sample aperture and sample in (cm)
1408        apOff=5.0
1409        // hard wire value for Ordela detectors
1410        DDet = 0.5              // resolution in cm
1411        //      String detStr=textRead[9]
1412        //      DDet = DetectorPixelResolution(fileStr,detStr)          //needs detector type and beamline
1413
1414        //Go from 0 to nq doing the calc for all three values at
1415        //every Q value
1416
1417        ii=0
1418        Variable ret1,ret2,ret3
1419        do
1420                S_getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,ret1,ret2,ret3)
1421                sigmaq[ii] = ret1       //res_wave[0]
1422                qbar[ii] = ret2         //res_wave[1]
1423                fsubs[ii] = ret3                //res_wave[2]
1424                ii+=1
1425        while(ii<nq)
1426        DeletePoints startElement,numElements, sigmaq, qbar, fsubs
1427
1428        fsubs += 1e-8           //keep the values from being too small
1429
1430// End of resolution calculations
1431// ***************************************************************
1432       
1433        //get rid of the default mask, if one was created (it is in the current folder)
1434        //don't just kill "mask" since it might be pointing to the one in the MSK folder
1435        Killwaves/Z $(destPath+":mask")
1436       
1437        //return to root folder (redundant)
1438        SetDataFolder root:
1439       
1440        Return 0
1441End
1442
1443//returns nq, new number of q-values
1444//arrays aveint,dsq,ncells are also changed by this function
1445//
1446Function S_IncrementPixel(dataPixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1447        Variable dataPixel,ddr,dxx,dyy
1448        Wave aveint,dsq,ncells
1449        Variable nq,nd2
1450       
1451        Variable ir
1452       
1453        ir = trunc(sqrt(dxx*dxx+dyy*dyy)/ddr)+1
1454        if (ir>nq)
1455                nq = ir         //resets maximum number of q-values
1456        endif
1457        aveint[ir-1] += dataPixel/nd2           //ir-1 must be used, since ir is physical
1458        dsq[ir-1] += dataPixel*dataPixel/nd2
1459        ncells[ir-1] += 1/nd2
1460       
1461        Return nq
1462End
1463
1464//function determines azimuthal angle dphi that a vector connecting
1465//center of detector to pixel makes with respect to vector
1466//at chosen azimuthal angle phi -> [cos(phi),sin(phi)] = [phi_x,phi_y]
1467//dphi is always positive, varying from 0 to Pi
1468//
1469Function S_dphi_pixel(dxx,dyy,phi_x,phi_y)
1470        Variable dxx,dyy,phi_x,phi_y
1471       
1472        Variable val,rr,dot_prod
1473       
1474        rr = sqrt(dxx^2 + dyy^2)
1475        dot_prod = (dxx*phi_x + dyy*phi_y)/rr
1476        //? correct for roundoff error? - is this necessary in IGOR, w/ double precision?
1477        if(dot_prod > 1)
1478                dot_prod =1
1479        Endif
1480        if(dot_prod < -1)
1481                dot_prod = -1
1482        Endif
1483       
1484        val = acos(dot_prod)
1485       
1486        return val
1487
1488End
1489
1490//calculates the x distance from the center of the detector, w/nonlinear corrections
1491//
1492Function S_FX(xx,sx3,xcenter,sx)               
1493        Variable xx,sx3,xcenter,sx
1494       
1495        Variable retval
1496       
1497        retval = sx3*tan((xx-xcenter)*sx/sx3)
1498        Return retval
1499End
1500
1501//calculates the y distance from the center of the detector, w/nonlinear corrections
1502//
1503Function S_FY(yy,sy3,ycenter,sy)               
1504        Variable yy,sy3,ycenter,sy
1505       
1506        Variable retval
1507       
1508        retval = sy3*tan((yy-ycenter)*sy/sy3)
1509        Return retval
1510End
1511
1512//**********************
1513// Resolution calculation - used by the averaging routines
1514// to calculate the resolution function at each q-value
1515// - the return value is not used
1516//
1517// equivalent to John's routine on the VAX Q_SIGMA_AVE.FOR
1518// Incorporates eqn. 3-15 from J. Appl. Cryst. (1995) v. 28 p105-114
1519//
1520Function/S S_getResolution(inQ,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,SigmaQ,QBar,fSubS)
1521        Variable inQ, lambda, lambdaWidth, DDet, apOff, S1, S2, L1, L2, BS, del_r
1522        Variable &fSubS, &QBar, &SigmaQ         //these are the output quantities at the input Q value
1523       
1524        //lots of calculation variables
1525        Variable a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g
1526        Variable r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r
1527
1528        //Constants
1529        //Variable del_r = .1
1530        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
1531        Variable g = 981.0                              //gravity acceleration [cm/s^2]
1532
1533        String results
1534        results ="Failure"
1535       
1536        NVAR usingLenses = root:Packages:NIST:SAS:gUsingLenses
1537
1538        //rename for working variables,  these must be gotten from global
1539        //variables
1540
1541//      Variable wLam, wLW, wL1, wL2, wS1, wS2
1542//      Variable wBS, wDDet, wApOff
1543//      wLam = lambda
1544//      wLW = lambdaWidth
1545//      wDDet = DDet
1546//      wApOff = apOff
1547        S1 *= 0.5*0.1                   //convert to radius and [cm]
1548        S2 *= 0.5*0.1
1549
1550        L1 *= 100.0                     // [cm]
1551        L1 -= apOff                             //correct the distance
1552
1553        L2 *= 100.0
1554        L2 += apOff
1555
1556        BS *= 0.5*0.1                   //convert to radius and [cm]
1557        del_r *= 0.1                            //width of annulus, convert mm to [cm]
1558       
1559        //Start resolution calculation
1560        a2 = S1*L2/L1 + S2*(L1+L2)/L1
1561        q_small = 2.0*Pi*(BS-a2)*(1.0-lambdaWidth)/(lambda*L2)
1562        lp = 1.0/( 1.0/L1 + 1.0/L2)
1563
1564        v_lambda = lambdaWidth^2/6.0
1565       
1566        if(usingLenses==1)                      //SRK 2007
1567                v_b = 0.25*(S1*L2/L1)^2 +0.25*(2/3)*(lambdaWidth/lambda)^2*(S2*L2/lp)^2         //correction to 2nd term
1568        else
1569                v_b = 0.25*(S1*L2/L1)^2 +0.25*(S2*L2/lp)^2              //original form
1570        endif
1571       
1572        v_d = (DDet/2.3548)^2 + del_r^2/12.0
1573        vz = vz_1 / lambda
1574        yg = 0.5*g*L2*(L1+L2)/vz^2
1575        v_g = 2.0*(2.0*yg^2*v_lambda)                                   //factor of 2 correction, B. Hammouda, 2007
1576
1577        r0 = L2*tan(2.0*asin(lambda*inQ/(4.0*Pi) ))
1578        delta = 0.5*(BS - r0)^2/v_d
1579
1580        if (r0 < BS)
1581                inc_gamma=exp(gammln(1.5))*(1-gammp(1.5,delta))
1582        else
1583                inc_gamma=exp(gammln(1.5))*(1+gammp(1.5,delta))
1584        endif
1585
1586        fSubS = 0.5*(1.0+erf( (r0-BS)/sqrt(2.0*v_d) ) )
1587        if (fSubS <= 0.0)
1588                fSubS = 1.e-10
1589        endif
1590        fr = 1.0 + sqrt(v_d)*exp(-1.0*delta) /(r0*fSubS*sqrt(2.0*Pi))
1591        fv = inc_gamma/(fSubS*sqrt(Pi)) - r0^2*(fr-1.0)^2/v_d
1592
1593        rmd = fr*r0
1594        v_r1 = v_b + fv*v_d +v_g
1595
1596        rm = rmd + 0.5*v_r1/rmd
1597        v_r = v_r1 - 0.5*(v_r1/rmd)^2
1598        if (v_r < 0.0)
1599                v_r = 0.0
1600        endif
1601        QBar = (4.0*Pi/lambda)*sin(0.5*atan(rm/L2))
1602        SigmaQ = QBar*sqrt(v_r/rmd^2 +v_lambda)
1603       
1604        results = "success"
1605        Return results
1606End
1607
1608Function S_Debye(scale,rg,bkg,x)
1609        Variable scale,rg,bkg
1610        Variable x
1611       
1612        // variables are:
1613        //[0] scale factor
1614        //[1] radius of gyration [A]
1615        //[2] background        [cm-1]
1616       
1617        // calculates (scale*debye)+bkg
1618        Variable Pq,qr2
1619       
1620        qr2=(x*rg)^2
1621        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
1622       
1623        //scale
1624        Pq *= scale
1625        // then add in the background
1626        return (Pq+bkg)
1627End
1628
1629
1630Function S_SphereForm(scale,radius,delrho,bkg,x)                               
1631        Variable scale,radius,delrho,bkg
1632        Variable x
1633       
1634        // variables are:                                                       
1635        //[0] scale
1636        //[1] radius (A)
1637        //[2] delrho (A-2)
1638        //[3] background (cm-1)
1639       
1640//      Variable scale,radius,delrho,bkg                               
1641//      scale = w[0]
1642//      radius = w[1]
1643//      delrho = w[2]
1644//      bkg = w[3]
1645       
1646       
1647        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
1648        // and is rescaled to give [=] cm^-1
1649       
1650        Variable bes,f,vol,f2
1651        //
1652        //handle q==0 separately
1653        If(x==0)
1654                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
1655                return(f)
1656        Endif
1657       
1658        bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
1659        vol = 4*pi/3*radius^3
1660        f = vol*bes*delrho              // [=] A
1661        // normalize to single particle volume, convert to 1/cm
1662        f2 = f * f / vol * 1.0e8                // [=] 1/cm
1663       
1664        return (scale*f2+bkg)   // Scale, then add in the background
1665       
1666End
1667
1668Function/S SetConfigurationText()
1669
1670        String str="",temp
1671
1672        SetDataFolder root:Packages:NIST:SAS
1673       
1674        NVAR numberOfGuides=gNg
1675        NVAR gTable=gTable              //2=chamber, 1=table
1676        NVAR wavelength=gLambda
1677        NVAR lambdaWidth=gDeltaLambda
1678        NVAR instrument = instrument
1679        NVAR L2diff = L2diff
1680        NVAR lens = root:Packages:NIST:SAS:gUsingLenses
1681        SVAR aStr = root:myGlobals:gAngstStr
1682       
1683        sprintf temp,"Source Aperture Diameter =\t\t%6.2f cm\r",sourceApertureDiam()
1684        str += temp
1685        sprintf temp,"Source to Sample =\t\t\t\t%6.0f cm\r",sourceToSampleDist()
1686        str += temp
1687        sprintf temp,"Sample Aperture to Detector =\t%6.0f cm\r",sampleToDetectorDist() + L2diff
1688        str += temp
1689        sprintf temp,"Beam diameter =\t\t\t\t\t%6.2f cm\r",beamDiameter("maximum")
1690        str += temp
1691        sprintf temp,"Beamstop diameter =\t\t\t\t%6.2f inches\r",beamstopDiam()/2.54
1692        str += temp
1693        sprintf temp,"Minimum Q-value =\t\t\t\t%6.4f 1/%s (sigQ/Q = %3.1f %%)\r",qMin(),aStr,deltaQ(qMin())
1694        str += temp
1695        sprintf temp,"Maximum Horizontal Q-value =\t%6.4f 1/%s\r",qMaxHoriz(),aStr
1696        str += temp
1697        sprintf temp,"Maximum Vertical Q-value =\t\t%6.4f 1/%s\r",qMaxVert(),aStr
1698        str += temp
1699        sprintf temp,"Maximum Q-value =\t\t\t\t%6.4f 1/%s (sigQ/Q = %3.1f %%)\r",qMaxCorner(),aStr,deltaQ(qMaxCorner())
1700        str += temp
1701        sprintf temp,"Beam Intensity =\t\t\t\t%.0f counts/s\r",beamIntensity()
1702        str += temp
1703        sprintf temp,"Figure of Merit =\t\t\t\t%3.3g %s^2/s\r",figureOfMerit(),aStr
1704        str += temp
1705        sprintf temp,"Attenuator transmission =\t\t%3.3g = Atten # %d\r",attenuatorTransmission(),attenuatorNumber()
1706        str += temp
1707//     
1708//      // add text of the user-edited values
1709//      //
1710        sprintf temp,"***************** NG %d *****************\r",instrument
1711        str += temp
1712        sprintf temp,"Sample Aperture Diameter =\t\t\t\t%.2f cm\r",sampleApertureDiam()
1713        str += temp
1714        sprintf temp,"Number of Guides =\t\t\t\t\t\t%d \r", numberOfGuides
1715        str += temp
1716        sprintf temp,"Sample Chamber to Detector =\t\t\t%.1f cm\r", chamberToDetectorDist()
1717        str += temp
1718        if(gTable==1)
1719                sprintf temp,"Sample Position is \t\t\t\t\t\tHuber\r"
1720        else
1721                sprintf temp,"Sample Position is \t\t\t\t\t\tChamber\r"
1722        endif
1723        str += temp
1724        sprintf temp,"Detector Offset =\t\t\t\t\t\t%.1f cm\r", detectorOffset()
1725        str += temp
1726        sprintf temp,"Neutron Wavelength =\t\t\t\t\t%.2f %s\r", wavelength,aStr
1727        str += temp
1728        sprintf temp,"Wavelength Spread, FWHM =\t\t\t\t%.3f\r", lambdaWidth
1729        str += temp
1730        sprintf temp,"Sample Aperture to Sample Position =\t%.2f cm\r", L2Diff
1731        str += temp
1732        if(lens==1)
1733                sprintf temp,"Lenses are IN\r"
1734        else
1735                sprintf temp,"Lenses are OUT\r"
1736        endif
1737        str += temp
1738       
1739   setDataFolder root:
1740   return str                   
1741End
1742
1743Function DisplayConfigurationText()
1744
1745        if(WinType("Trial_Configuration")==0)
1746                NewNotebook/F=0/K=1/N=Trial_Configuration /W=(480,44,880,369)
1747        endif
1748        //replace the text
1749        Notebook Trial_Configuration selection={startOfFile, endOfFile}
1750        Notebook Trial_Configuration font="Monaco",fSize=10,text=SetConfigurationText()
1751        return(0)
1752end
1753
1754//parses the control for A1 diam
1755// updates the wave
1756Function sourceApertureDiam()
1757        ControlInfo/W=SASCALC popup0
1758        Variable diam
1759        sscanf S_Value, "%g cm", diam
1760        WAVE rw=root:Packages:NIST:SAS:realsRead
1761        rw[23] = diam*10                        //sample aperture diameter in mm
1762        return(diam)
1763end
1764
1765// change the sample aperture to a non-standard value
1766Function SampleApOtherSetVarProc(ctrlName,varNum,varStr,varName) : SetVariableControl
1767        String ctrlName
1768        Variable varNum
1769        String varStr
1770        String varName
1771               
1772        WAVE rw=root:Packages:NIST:SAS:realsRead
1773        rw[24] = varNum                 //sample aperture diameter in mm
1774        ReCalculateInten(1)
1775        return(0)
1776End
1777
1778//parses the control for A2 diam
1779// updates the wave and global
1780// returns a2 in cm
1781Function sampleApertureDiam()
1782        ControlInfo/W=SASCALC popup0_1
1783       
1784//      Print "In sampleApertureDiam()"
1785        //set the global
1786        NVAR a2=root:Packages:NIST:SAS:gSamAp
1787       
1788        if(cmpstr(S_Value,"other") == 0)                // "other" selected
1789                //enable the setvar, diameter in mm!
1790                SetVariable setvar0_3 disable=0
1791                // read its value (a global)
1792                NVAR a2other = root:Packages:NIST:SAS:gSamApOther
1793                a2=a2other/10                           //a2 in cm
1794        else
1795                SetVariable setvar0_3 disable=1
1796                //1st item is 1/16", popup steps by 1/16"
1797                a2 = 2.54/16.0 * (V_Value)                      //convert to cm         
1798        endif
1799        WAVE rw=root:Packages:NIST:SAS:realsRead
1800        rw[24] = a2*10                  //sample aperture diameter in mm
1801       
1802        return(a2)
1803end
1804
1805//1=huber 2=chamber
1806Function tableposition()
1807        ControlInfo/W=SASCALC checkHuber
1808        if(V_Value)
1809                return(1)               //huber
1810        else
1811                return(2)               //chamber
1812        endif
1813End
1814
1815//compute SSD and update both the global and the wave
1816//
1817Function sourceToSampleDist()
1818
1819        NVAR NG=root:Packages:NIST:SAS:gNg
1820        NVAR S12 = root:Packages:NIST:SAS:S12
1821        NVAR L2Diff = root:Packages:NIST:SAS:L2Diff
1822        NVAR SSD = root:Packages:NIST:SAS:gSSD
1823       
1824        SSD = 1632 - 155*NG - s12*(2-tableposition()) - L2Diff
1825       
1826        WAVE rw=root:Packages:NIST:SAS:realsRead
1827        rw[25] = SSD/100                // in meters
1828        return(SSD)
1829End
1830
1831
1832// not part of SASCALC, but can be used to convert the SSD to number of guides
1833//
1834// SSD in meters
1835Function numGuides(SSD)
1836        variable SSD
1837       
1838        Variable Ng
1839        Ng = SSD*100 + 5 - 1632
1840        Ng /= -155
1841       
1842        Ng = round(Ng)
1843        return(Ng)
1844End
1845
1846
1847//returns the offset value
1848// slider and setVar are linked to the same global
1849// updates the wave and changes the beamcenter (x,y) in the wave
1850Function detectorOffset()
1851       
1852        WAVE rw=root:Packages:NIST:SAS:RealsRead
1853        NVAR val = root:Packages:NIST:SAS:gOffset
1854        rw[19] = val            // already in cm
1855        //move the beamcenter, make it an integer value for the MC simulation
1856        rw[16] = 64 + round(2*rw[19])           //approximate beam X is 64 w/no offset, 114 w/25 cm offset
1857        rw[17] = 64             //typical value
1858       
1859        return(val)
1860end
1861
1862//returns the detector distance (slider and setVar are linked by the global)
1863//
1864Function chamberToDetectorDist()
1865
1866        NVAR val = root:Packages:NIST:SAS:gDetDist
1867        return(val)
1868End
1869
1870//sets the SDD (slider and setVar are linked by the global and is the detector position
1871//  relative to the chamber)
1872// updates the wave
1873Function sampleToDetectorDist()
1874
1875        NVAR detDist = root:Packages:NIST:SAS:gDetDist
1876        NVAR S12 = root:Packages:NIST:SAS:S12
1877        WAVE rw=root:Packages:NIST:SAS:RealsRead       
1878        Variable SDD   
1879       
1880        SDD = detDist + s12*(2-tableposition())
1881        rw[18] = SDD/100                // convert to meters for header
1882        return(SDD)
1883End
1884
1885//direction = one of "vertical;horizontal;maximum;"
1886// all of this is bypassed if the lenses are in
1887//
1888Function beamDiameter(direction)
1889        String direction
1890
1891        NVAR lens = root:Packages:NIST:SAS:gUsingLenses
1892        if(lens)
1893                return sourceApertureDiam()
1894        endif
1895       
1896    Variable l1 = sourceToSampleDist()
1897    Variable l2 //= sampleAperToDetDist()
1898    Variable d1,d2,bh,bv,bm,umbra,a1,a2
1899   
1900    NVAR L2diff = root:Packages:NIST:SAS:L2diff
1901    NVAR lambda = root:Packages:NIST:SAS:gLambda
1902    NVAR lambda_width = root:Packages:NIST:SAS:gDeltaLambda
1903    NVAR bs_factor = root:Packages:NIST:SAS:bs_factor
1904   
1905    l2 = sampleToDetectorDist() + L2diff
1906    a1 = sourceApertureDiam()
1907    a2 = sampleApertureDiam()
1908   
1909    d1 = a1*l2/l1
1910    d2 = a2*(l1+l2)/l1
1911    bh = d1+d2          //beam size in horizontal direction
1912    umbra = abs(d1-d2)
1913    //vertical spreading due to gravity
1914    bv = bh + 1.25e-8*(l1+l2)*l2*lambda*lambda*lambda_width
1915    bm = (bs_factor*bh > bv) ? bs_factor*bh : bv //use the larger of horiz*safety or vertical
1916   
1917    strswitch(direction)        // string switch
1918        case "vertical":                // execute if case matches expression
1919                return(bv)
1920                break                                           // exit from switch
1921        case "horizontal":              // execute if case matches expression
1922                return(bh)
1923                break
1924        case "maximum":         // execute if case matches expression
1925                return(bm)
1926                break
1927        default:                                                        // optional default expression executed
1928                return(bm)                                              // when no case matches
1929    endswitch
1930End
1931
1932//on NG3 and NG7, allowable sizes are 1,2,3,4" diameter
1933//will return values larger than 4.0*2.54 if a larger beam is needed
1934//
1935// - in an approximate way, account for lenses
1936Function beamstopDiam()
1937
1938        NVAR yesLens = root:Packages:NIST:SAS:gUsingLenses
1939        Variable bm=0
1940        Variable bs=0.0
1941   
1942        if(yesLens)
1943                //bm = sourceApertureDiam()             //ideal result, not needed
1944                bs = 1                                                          //force the diameter to 1"
1945        else
1946                bm = beamDiameter("maximum")
1947                do
1948                bs += 1
1949           while ( (bs*2.54 < bm) || (bs > 30.0))                       //30 = ridiculous limit to avoid inf loop
1950        endif
1951
1952        //update the wave
1953        WAVE rw=root:Packages:NIST:SAS:realsRead
1954        rw[21] = bs*25.4                //store the BS diameter in mm
1955       
1956    return (bs*2.54)            //return diameter in cm, not inches for txt
1957End
1958
1959//returns the projected diameter of the beamstop at the anode plane.
1960// most noticeable at short SDD
1961//if flag == 0 use conservative estimate = largest diameter (for SASCALC, default)
1962//if flag != 0 use point aperture = average diameter (for resolution calculation)
1963Function beamstopDiamProjection(flag)
1964        Variable flag
1965       
1966        NVAR L2diff = root:Packages:NIST:SAS:L2diff
1967        Variable a2 = sampleApertureDiam()
1968        Variable bs = beamstopDiam()
1969        Variable l2, LB, BS_P
1970   
1971        l2 = sampleToDetectorDist() + L2diff
1972        LB = 20.1 + 1.61*BS                     //distance in cm from beamstop to anode plane (empirical)
1973        if(flag==0)
1974                BS_P = bs + (bs+a2)*lb/(l2-lb)          //diameter of shadow from parallax
1975        else
1976                BS_P = bs + bs*lb/(l2-lb)               //diameter of shadow, point A2
1977        endif
1978        return (bs_p)           //return projected diameter in cm
1979End
1980
1981// 19MAR07 - using correction from John for an estimate of the shadow of the beamstop
1982// at the detection plane. This is a noticeable effect at short SDD, where the projected
1983// diameter of the beamstop is much larger than the physical diameter.
1984Function qMin()
1985
1986    Variable l2s = sampleToDetectorDist()       //distance from sample to detector in cm
1987//    Variable bs = beamstopDiam()              //beamstop diameter in cm
1988    Variable bs_p = beamstopDiamProjection(0)           //projected beamstop diameter in cm
1989    NVAR lambda = root:Packages:NIST:SAS:gLambda
1990    NVAR d_det = root:Packages:NIST:SAS:d_det           //cm
1991    NVAR a_pixel = root:Packages:NIST:SAS:a_pixel       //cm
1992
1993    return( (pi/lambda)*(bs_p + d_det + a_pixel)/l2s )          //use bs_p rather than bs
1994   // return( (pi/lambda)*(bs + d_det + a_pixel)/l2s )          //use bs (incorrect)
1995End
1996
1997Function qMaxVert()
1998
1999    Variable theta
2000    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
2001    NVAR lambda = root:Packages:NIST:SAS:gLambda
2002        NVAR det_width = root:Packages:NIST:SAS:det_width
2003       
2004    theta = atan( (det_width/2.0)/l2s )
2005   
2006    return ( 4.0*pi/lambda * sin(theta/2.0) )
2007end
2008
2009Function qMaxCorner()
2010
2011    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
2012    Variable radial
2013    NVAR lambda = root:Packages:NIST:SAS:gLambda
2014        NVAR det_off = root:Packages:NIST:SAS:gOffset
2015        NVAR det_width = root:Packages:NIST:SAS:det_width
2016
2017    radial=sqrt( (0.5*det_width)*(0.5*det_width)+(0.5*det_width+det_off)*(0.5*det_width+det_off) )
2018   
2019    return ( 4*pi/lambda*sin(0.5*atan(radial/l2s)) )
2020End
2021
2022Function qMaxHoriz()
2023
2024    Variable theta
2025    Variable l2s = sampleToDetectorDist()       //distance from sample to detector
2026    NVAR lambda = root:Packages:NIST:SAS:gLambda
2027        NVAR det_off = root:Packages:NIST:SAS:gOffset
2028        NVAR det_width = root:Packages:NIST:SAS:det_width
2029
2030    theta = atan( ((det_width/2.0)+det_off)/l2s )       //from the instance variables
2031   
2032    return ( 4.0*pi/lambda * sin(theta/2.0) )
2033End
2034
2035// calculate sigma for the resolution function at either limit of q-range
2036Function deltaQ(atQ)
2037        Variable atQ
2038       
2039    Variable k02,lp,l1,l2,sig_02,sigQ2,a1,a2
2040    NVAR l2Diff = root:Packages:NIST:SAS:L2diff
2041        NVAR lambda = root:Packages:NIST:SAS:gLambda
2042        NVAR lambda_width = root:Packages:NIST:SAS:gDeltaLambda
2043        NVAR d_det = root:Packages:NIST:SAS:d_det
2044        NVAR del_r = root:Packages:NIST:SAS:del_r
2045       
2046       
2047    l1 = sourceToSampleDist()
2048    l2 = sampleToDetectorDist() + L2diff
2049    a1 = sourceApertureDiam()
2050    a2 = sampleApertureDiam()
2051   
2052    k02 = (6.2832/lambda)*(6.2832/lambda)
2053    lp = 1/(1/l1 + 1/l2)
2054   
2055    sig_02 = (0.25*a1/l1)*(0.25*a1/l1)
2056    sig_02 += (0.25*a2/lp)*(0.25*a2/lp)
2057    sig_02 += (d_det/(2.355*l2))*(d_det/(2.355*l2))
2058    sig_02 += (del_r/l2)*(del_r/l2)/12
2059    sig_02 *= k02
2060   
2061    sigQ2 = sig_02 + (atQ*lambda_width)*(atQ*lambda_width)/6
2062
2063    return(100*sqrt(sigQ2)/atQ)
2064End
2065
2066
2067Function beamIntensity()
2068
2069    Variable alpha,f,t,t4,t5,t6,as,solid_angle,l1,d2_phi
2070    Variable a1,a2,retVal
2071    SetDataFolder root:Packages:NIST:SAS
2072    NVAR l_gap=l_gap,guide_width =guide_width,ng = gNg
2073    NVAR lambda_t=lambda_t,b=b,c=c
2074    NVAR lambda=gLambda,t1=t1,t2=t2,t3=t3,phi_0=phi_0
2075    NVAR lambda_width=gDeltaLambda
2076   
2077    l1 = sourceToSampleDist()
2078    a1 = sourceApertureDiam()
2079    a2 = sampleApertureDiam()
2080   
2081   
2082    alpha = (a1+a2)/(2*l1)      //angular divergence of beam
2083    f = l_gap*alpha/(2*guide_width)
2084    t4 = (1-f)*(1-f)
2085    t5 = exp(ng*ln(0.96))       // trans losses of guides in pre-sample flight
2086    t6 = 1 - lambda*(b-(ng/8)*(b-c))            //experimental correction factor
2087    t = t1*t2*t3*t4*t5*t6
2088   
2089    as = pi/4*a2*a2             //area of sample in the beam
2090    d2_phi = phi_0/(2*pi)
2091    d2_phi *= exp(4*ln(lambda_t/lambda))
2092    d2_phi *= exp(-1*(lambda_t*lambda_t/lambda/lambda))
2093
2094    solid_angle = pi/4* (a1/l1)*(a1/l1)
2095
2096    retVal = as * d2_phi * lambda_width * solid_angle * t
2097     SetDataFolder root:
2098    return (retVal)
2099end
2100
2101Function figureOfMerit()
2102
2103        Variable bi = beamIntensity()
2104        NVAR lambda = root:Packages:NIST:SAS:gLambda
2105       
2106    return (lambda*lambda*bi)
2107End
2108
2109//estimate the number of pixels in the beam, and enforce the maximum countrate per pixel (idmax)
2110Function attenuatorTransmission()
2111
2112    Variable num_pixels,i_pix           //i_pix = id in John's notation
2113    Variable bDiam = beamDiameter("horizontal") //!! note that prev calculations used bh (horizontal)
2114    Variable atten,a2
2115    SetDataFolder root:Packages:NIST:SAS
2116    NVAR a_pixel=a_pixel,idmax=idmax
2117   
2118    a2 = sampleApertureDiam()
2119   
2120    num_pixels = pi/4*(0.5*(a2+bDiam))*(0.5*(a2+bDiam))/a_pixel/a_pixel
2121    i_pix = ( beamIntensity() )/num_pixels
2122   
2123    atten = (i_pix < idmax) ? 1.0 : idmax/i_pix
2124    SetDataFolder root:
2125    return(atten)
2126End
2127
2128Function attenuatorNumber()
2129
2130    Variable atten = attenuatorTransmission()
2131    Variable af,nf,numAtten
2132    SetDataFolder root:Packages:NIST:SAS
2133    NVAR lambda=gLambda
2134   
2135    af = 0.498 + 0.0792*lambda - 1.66e-3*lambda*lambda
2136    nf = -ln(atten)/af          //floating point
2137   
2138    numAtten = trunc(nf) + 1                    //in c, (int)nf
2139    //correct for larger step thickness at n > 6
2140    if(numAtten > 6)
2141        numAtten = 7 + trunc( (numAtten-6)/2 )          //in c, numAtten = 7 + (int)( (numAtten-6)/2 )
2142    endif
2143   
2144    SetDatafolder root:
2145    return (numAtten)
2146End
Note: See TracBrowser for help on using the repository browser.