source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_DetectorBinning.ipf @ 935

Last change on this file since 935 was 935, checked in by srkline, 9 years ago

Adding a directory with the beginning procedures for "VCALC" -- the VSANS version of SASCALC.

Including the VSANS_Includes file is all that is needed - no other SANS procedures are used.
File size: 53.4 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3
4///// Procedures for:
5//      Generating the detector panels
6//      Filling the panels with Qtot, QxQyQz
7// Filling the "data" with a model function
8// Averaging the panels (independently) into I(Q)
9//
10//
11//There are some things in the current circular averaging that don't make any sense
12//and don't seem to really do anything at all, so trim them out.
13//1) subdividing pixels near the beam stop into 9 sub-pixels
14//2) non-linear correction (maybe keep this as a dummy). I may need to put some sort
15//of correction here for the offset, or it may already be done as I receive it.
16//
17//
18//
19//Do I separate out the circular, sector, rectangular, annular averaging into
20//separate routines?
21//
22//
23
24
25Proc PlotFrontPanels()
26        fPlotFrontPanels()
27End
28
29// to plot I(q) for the 4 front panels
30//
31// *** Call this function when front panels are adjusted, or wavelength, etc. changed
32//
33Function fPlotFrontPanels()
34
35        // space is allocated for all of the detectors and Q's on initialization
36        // calculate Qtot, qxqyqz arrays from geometry
37        V_CalculateQFrontPanels()
38       
39        // fill the panels with fake sphere scattering data
40        // TODO: am I in the right data folder??
41        SetDataFolder root:Packages:NIST:VSANS:VCALC
42
43        WAVE det_FL = det_FL
44        WAVE det_FR = det_FR
45        WAVE det_FT = det_FT
46        WAVE det_FB = det_FB
47       
48        WAVE qTot_FL = qTot_FL
49        WAVE qTot_FR = qTot_FR
50        WAVE qTot_FT = qTot_FT
51        WAVE qTot_FB = qTot_FB
52
53        FillPanel_wModelData(det_FL,qTot_FL,"FL")
54        FillPanel_wModelData(det_FR,qTot_FR,"FR")
55        FillPanel_wModelData(det_FT,qTot_FT,"FT")
56        FillPanel_wModelData(det_FB,qTot_FB,"FB")
57//      det_FL = V_SphereForm(1,60,1e-6,1,qTot_FL[p][q])                               
58//      det_FR = V_SphereForm(1,60,1e-6,1,qTot_FR[p][q])                               
59//      det_FT = V_SphereForm(1,60,1e-6,1,qTot_FT[p][q])                               
60//      det_FB = V_SphereForm(1,60,1e-6,1,qTot_FB[p][q])                       
61
62        SetDataFolder root:
63               
64        // set any "shadowed" area of the T/B detectors to NaN to get a realitic
65        // view of how much of the detectors are actually collecting data
66        // -- I can get the separation L/R from the panel - only this "open" width is visible.
67        //TODO - make this a proper shadow - TB extent of the LR panels matters too, not just the LR separation
68        V_SetShadow_TopBottom("","FT")          // TODO: -- be sure the data folder is properly set (within the function...)
69        V_SetShadow_TopBottom("","FB")
70       
71        // do the q-binning for each of the panels to get I(Q)
72        Execute "BinAllFrontPanels()"
73
74        // plot the results
75        Execute "Front_IQ_Graph()"
76        Execute "FrontPanels_AsQ()"
77End
78
79// TODO: hard wired for a sphere - change this to allow minimal selections and altering of coefficients
80// TODO: add the "fake" 2D simulation to fill the panels which are then later averaged as I(Q)
81Function FillPanel_wModelData(det,qTot,type)
82        Wave det,qTot
83        String type
84
85        SetDataFolder root:Packages:NIST:VSANS:VCALC
86
87        // q-values and detector arrays already allocated and calculated
88        Duplicate/O det tmpInten,tmpSig,prob_i
89               
90        Variable imon,trans,thick,sdd,pixSizeX,pixSizeY,sdd_offset
91
92        //imon = V_BeamIntensity()*CountTime
93        imon = VCALC_getImon()          //TODO: currently from the panel, not calculated
94        trans = 0.8
95        thick = 0.1
96       
97        // need SDD
98        // need pixel dimensions
99        // nominal sdd in meters, offset in mm, want result in cm !
100        sdd = VCALC_getSDD(type)*100    +  VSANS_getTopBottomSDDOffset(type) / 10               // result is sdd in [cm]
101
102        pixSizeX = VCALC_getPixSizeX(type)              // cm
103        pixSizeY = VCALC_getPixSizeY(type)
104       
105       
106        //?? pick the function from a popup on the panel? (bypass the analysis panel, or maybe it's better to
107        //  keep the panel to keep people used to using it.)
108        String funcStr = VCALC_getModelFunctionStr()
109        strswitch(funcStr)
110                case "Big Debye":
111                        tmpInten = V_Debye(10,3000,0.0001,qTot[p][q])
112                        break
113                case "Big Sphere":
114                        tmpInten = V_SphereForm(1,900,1e-6,0.01,qTot[p][q])     
115                        break
116                case "Debye":
117                        tmpInten = V_Debye(10,300,0.0001,qTot[p][q])
118                        break
119                case "Sphere":
120                        tmpInten = V_SphereForm(1,60,1e-6,0.001,qTot[p][q])     
121                        break
122                default:
123                        tmpInten = V_Debye(10,300,0.1,qTot[p][q])
124        endswitch
125//      tmpInten = V_SphereForm(1,100,1e-6,1,qTot[p][q])       
126       
127//      tmpInten = V_Debye(scale,rg,bkg,x)                     
128//      tmpInten = V_Debye(10,300,0.1,qTot[p][q])
129
130
131///////////////
132//      // calculate the scattering cross section simply to be able to estimate the transmission
133//      Variable sig_sas=0
134//     
135//      // remember that the random deviate is the coherent portion ONLY - the incoherent background is
136//      // subtracted before the calculation.
137//      CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
138//
139//      if(sig_sas > 100)
140//              DoAlert 0,"SAS cross section > 100. Estimates of multiple scattering are unreliable. Choosing a model with a well-defined Rg may help"
141//      endif           
142//
143//      // calculate the multiple scattering fraction for display (10/2009)
144//      Variable ii,nMax=10,tau
145//      mScat=0
146//      tau = thick*sig_sas
147//      // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
148//      // neutrons out of those that were scattered
149//      for(ii=2;ii<nMax;ii+=1)
150//              mScat += tau^(ii)/factorial(ii)
151////            print tau^(ii)/factorial(ii)
152//      endfor
153//      estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
154//      mscat *= (estTrans)/(1-estTrans)
155//
156////    if(mScat > 0.1)         //  Display warning
157//
158//      Print "Sig_sas = ",sig_sas
159////////////////////
160       
161        prob_i = trans*thick*pixSizeX*pixSizeY/(sdd)^2*tmpInten                 //probability of a neutron in q-bin(i)
162               
163        tmpInten = (imon)*prob_i                //tmpInten is not the model calculation anymore!!
164
165
166/// **** can I safely assume a Gaussian error in the count rate??
167        tmpSig = sqrt(tmpInten)         // corrected based on John's memo, from 8/9/99
168
169        tmpInten += gnoise(tmpSig)
170        tmpInten = (tmpInten[p][q] < 0) ? 0 : tmpInten[p][q]                    // MAR 2013 -- is this the right thing to do
171        tmpInten = trunc(tmpInten)
172               
173       
174        det = tmpInten
175
176// if I want "absolute" scale -- then I lose the integer nature of the detector (but keep the random)
177        det /= trans*thick*pixSizeX*pixSizeY/(sdd)^2*imon
178       
179        KillWaves/Z tmpInten,tmpSig,prob_i     
180        SetDataFolder root:
181
182        return(0)
183End
184
185
186
187// works for Left, works for Right... works for T/B too.
188//
189// - TODO: be sure that the Q's are calculated correctly even when the beam is off of the
190//     detector, and on different sides of the detector (or T/B) - since it will be in a different
191//     relative postion to 0,0 on the detector. If the postions are symmetric, then the Q's should be identical.
192//     --- test this...
193// TODO -- be sure I'm in the right data folder. nothing is set correctly right now
194//
195// TODO: make all detector parameters global, not hard-wired
196//
197//
198// --- Panels are all allocated in the initialization. Here, only the q-values are calculated
199//     when anything changes
200//
201Function V_CalculateQFrontPanels()
202
203        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
204        Variable F_LR_sep,F_TB_sep,F_offset,F_sdd_offset
205
206// get the values from the panel + constants   
207        F_LR_sep = VCALC_getPanelSeparation("FLR")
208        F_TB_sep = VCALC_getPanelSeparation("FTB")
209        F_offset = VCALC_getLateralOffset("FL")
210       
211        SDD = VCALC_getSDD("FL")                //nominal SDD - need offset for TB
212        lam = VCALC_getWavelength()
213
214//separations are in mm -- need to watch the units, convert to cm
215        F_LR_sep /= 10
216        F_TB_sep /= 10
217// TODO - I'm treating the separation as the TOTAL width - so the difference
218//      from the "center" to the edge is 1/2 of the separation
219
220// TODO (make the N along the tube length a variable, since this can be reset @ acquisition)
221
222        F_sdd_offset = VSANS_getTopBottomSDDOffset("FT")        //T/B are 300 mm farther back  //TODO: make all detector parameters global, not hard-wired
223
224        SetDataFolder root:Packages:NIST:VSANS:VCALC   
225        Wave det_FL,det_FR                      // these are (48,256)
226        Wave det_FT,det_FB                      // these are (128,48)
227
228//FRONT/LEFT   
229        WAVE qTot_FL,qx_FL,qy_FL,qz_FL
230        qTot_FL = 0
231        qx_FL = 0
232        qy_FL = 0
233        qz_FL = 0       
234       
235// TODO - these are to be set from globals, not hard-wired. N and pixelSixze will be known (or pre-measured)
236// pixel sizes are in cm
237        pixSizeX = VCALC_getPixSizeX("FL")
238        pixSizeY = VCALC_getPixSizeY("FL")
239//      pixSizeX = 0.8                  // 0.8 cm/pixel along width
240//      pixSizeY = 0.4                  // approx 0.4 cm/pixel along length
241       
242        xCtr = 48+(F_LR_sep/2/pixSizeX)         // TODO  -- check -- starting from 47 rather than 48 (but I'm in pixel units for centers)??
243        yCtr = 127     
244        V_Detector_2Q(det_FL,qTot_FL,qx_FL,qy_FL,qz_FL,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
245//      Print "xy for FL = ",xCtr,yCtr
246       
247        //set the wave scaling for the detector image so that it can be plotted in q-space
248        // TODO: this is only approximate - since the left "edge" is not the same from top to bottom, so I crudely
249        // take the middle value. At very small angles, OK, at 1m, this is a crummy approximation.
250        // since qTot is magnitude only, I need to put in the (-ve)
251        SetScale/I x WaveMin(qx_FL),WaveMax(qx_FL),"", det_FL           //this sets the left and right ends of the data scaling
252        SetScale/I y WaveMin(qy_FL),WaveMax(qy_FL),"", det_FL
253//////////////////
254
255//FRONT/RIGHT
256        SetDataFolder root:Packages:NIST:VSANS:VCALC
257        WAVE qTot_FR,qx_FR,qy_FR,qz_FR
258        qTot_FR = 0
259        qx_FR = 0
260        qy_FR = 0
261        qz_FR = 0
262
263// TODO - these are to be set from globals, not hard-wired
264// pixel sizes are in cm
265        pixSizeX = VCALC_getPixSizeX("FR")
266        pixSizeY = VCALC_getPixSizeY("FR")
267       
268        xCtr = -(F_LR_sep/2/pixSizeX)-1         
269        yCtr = 127
270        V_Detector_2Q(det_FR,qTot_FR,qx_FR,qy_FR,qz_FR,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
271//      Print "xy for FR = ",xCtr,yCtr
272        SetScale/I x WaveMin(qx_FR),WaveMax(qx_FR),"", det_FR           //this sets the left and right ends of the data scaling
273        SetScale/I y WaveMin(qy_FR),WaveMax(qy_FR),"", det_FR
274/////////////////
275
276//FRONT/TOP
277        SetDataFolder root:Packages:NIST:VSANS:VCALC
278        WAVE qTot_FT,qx_FT,qy_FT,qz_FT
279        qTot_FT = 0
280        qx_FT = 0
281        qy_FT = 0
282        qz_FT = 0
283
284// TODO - these are to be set from globals, not hard-wired
285// pixel sizes are in cm
286        pixSizeX = VCALC_getPixSizeX("FT")
287        pixSizeY = VCALC_getPixSizeY("FT")
288
289        xCtr = 64
290        yCtr = -(F_TB_sep/2/pixSizeY)-1   
291        // global sdd_offset is in (mm), convert to meters here for the Q-calculation
292        V_Detector_2Q(det_FT,qTot_FT,qx_FT,qy_FT,qz_FT,xCtr,yCtr,sdd+F_sdd_offset/1000,lam,pixSizeX,pixSizeY)
293//      Print "xy for FT = ",xCtr,yCtr
294        SetScale/I x WaveMin(qx_FT),WaveMax(qx_FT),"", det_FT           //this sets the left and right ends of the data scaling
295        SetScale/I y WaveMin(qy_FT),WaveMax(qy_FT),"", det_FT
296//////////////////
297
298//FRONT/BOTTOM
299        SetDataFolder root:Packages:NIST:VSANS:VCALC
300        WAVE qTot_FB,qx_FB,qy_FB,qz_FB
301        qTot_FB = 0
302        qx_FB = 0
303        qy_FB = 0
304        qz_FB = 0
305
306// TODO - these are to be set from globals, not hard-wired
307// pixel sizes are in cm
308        pixSizeX = VCALC_getPixSizeX("FB")
309        pixSizeY = VCALC_getPixSizeY("FB")
310       
311        xCtr = 64
312        yCtr = 48+(F_TB_sep/2/pixSizeY)                 // TODO  -- check -- starting from 47 rather than 48 (but I'm in pixel units for centers)??
313        // global sdd_offset is in (mm), convert to meters here for the Q-calculation
314        V_Detector_2Q(det_FB,qTot_FB,qx_FB,qy_FB,qz_FB,xCtr,yCtr,sdd+F_sdd_offset/1000,lam,pixSizeX,pixSizeY)
315//      Print "xy for FB = ",xCtr,yCtr
316        SetScale/I x WaveMin(qx_FB),WaveMax(qx_FB),"", det_FB           //this sets the left and right ends of the data scaling
317        SetScale/I y WaveMin(qy_FB),WaveMax(qy_FB),"", det_FB
318/////////////////
319
320        SetDataFolder root:
321               
322        return(0)
323End
324
325
326// TODO" - this doesn't quite mask things out as they should be (too much is masked L/R of center)
327// and the outer edges of the detector are masked even if the TB panels extend past the TB of the LR panels.
328// ? skip the masking? but then I bin the detector data directly to get I(q), skipping the masked NaN values...
329//
330Function V_SetShadow_TopBottom(folderStr,type)
331        String folderStr,type
332       
333        Variable LR_sep,nPix,xCtr,ii,jj,numCol,pixSizeX,pixSizeY
334
335/// !! type passed in will be FT, FB, MT, MB, so I can't ask for the panel separation -- or I'll get the TB separation...
336        if(cmpstr(type[0],"F")==0)
337                //front
338                ControlInfo/W=VCALC VCALCCtrl_2a
339                LR_sep = V_Value       
340        else
341                //middle
342                ControlInfo/W=VCALC VCALCCtrl_3a
343                LR_sep = V_Value       
344        endif           
345//separations on panel are in mm -- need to watch the units, convert to cm
346        LR_sep /= 10
347
348//detector data
349        Wave det = $("root:Packages:NIST:VSANS:VCALC:"+"det_"+type)
350
351// TODO - these are to be set from globals, not hard-wired
352// pixel sizes are in cm for T/B detector
353// TODO: the "FT" check is hard wired for FRONT -- get rid of this...
354
355        pixSizeX = VCALC_getPixSizeX(type)
356        pixSizeY = VCALC_getPixSizeY(type)
357
358        //TODO -- get this from a global
359        xCtr = 64
360        nPix = (LR_sep/2/pixSizeX)              // # of pixels Left/right of center that are not obscured by L/R panels
361       
362        numCol = DimSize(det,0)         // x dim (columns)
363        for(ii=0;ii<(xCtr-nPix-2);ii+=1)
364                det[ii][] = NaN
365        endfor
366        for(ii=(xCtr+nPix+2);ii<numCol;ii+=1)
367                det[ii][] = NaN
368        endfor
369       
370        return(0)
371end
372
373
374Window FrontPanels_AsQ() : Graph
375//      DoWindow/F FrontPanels_AsQ
376//      if(V_flag == 0)
377//      PauseUpdate; Silent 1           // building window...
378//      Display /W=(1477,44,1978,517)
379
380        SetDataFolder root:Packages:NIST:VSANS:VCALC
381
382        CheckDisplayed /W=VCALC#Panels_Q det_FB
383        if(V_flag == 0)
384                AppendImage/W=VCALC#Panels_Q det_FB
385                ModifyImage/W=VCALC#Panels_Q det_FB ctab= {*,*,ColdWarm,0}
386                AppendImage/W=VCALC#Panels_Q det_FT
387                ModifyImage/W=VCALC#Panels_Q det_FT ctab= {*,*,ColdWarm,0}
388                AppendImage/W=VCALC#Panels_Q det_FL
389                ModifyImage/W=VCALC#Panels_Q det_FL ctab= {*,*,ColdWarm,0}
390                AppendImage/W=VCALC#Panels_Q det_FR
391                ModifyImage/W=VCALC#Panels_Q det_FR ctab= {*,*,ColdWarm,0}
392        endif
393
394        Variable dval
395        ControlInfo/W=VCALC setVar_b
396        dval = V_Value
397
398        SetAxis/W=VCALC#Panels_Q left -dval,dval
399        SetAxis/W=VCALC#Panels_Q bottom -dval,dval     
400
401        ControlInfo/W=VCALC check_0a
402// V_Value == 1 if checked
403        ModifyImage/W=VCALC#Panels_Q det_FB log=V_Value
404        ModifyImage/W=VCALC#Panels_Q det_FT log=V_Value
405        ModifyImage/W=VCALC#Panels_Q det_FL log=V_Value
406        ModifyImage/W=VCALC#Panels_Q det_FR log=V_Value
407
408
409        SetDataFolder root:
410       
411//      ModifyGraph width={Aspect,1},height={Aspect,1},gbRGB=(56797,56797,56797)
412//      ModifyGraph grid=2
413//      ModifyGraph mirror=2
414//      SetAxis left -0.2,0.2
415//      SetAxis bottom -0.2,0.2
416//      endif
417EndMacro
418
419
420// For a given detector panel, calculate the q-values
421// -work with everything as arrays
422// Input needed:
423// detector data
424// detector type (LRTB?)
425// beam center (may be off the detector)
426// SDD
427// lambda
428//
429// pixel dimensions for detector type (global constants)
430// - data dimensions read directly from array
431//
432// --What is calculated:
433// array of Q
434// array of qx,qy,qz
435// array of error already exists
436//
437//
438// -- sdd in meters
439// -- lambda in Angstroms
440Function V_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
441        Wave data,qTot,qx,qy,qz
442        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
443               
444        // loop over the array and calculate the values - this is done as a wave assignment
445// TODO -- be sure that it's p,q -- or maybe p+1,q+1 as used in WriteQIS.ipf   
446        qTot = V_CalcQval(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
447        qx = V_CalcQX(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
448        qy = V_CalcQY(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
449        qz = V_CalcQZ(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
450       
451        return(0)
452End
453
454
455//////////////////////
456// NOTE: The Q calculations are different than what is in GaussUtils in that they take into
457// accout the different x/y pixel sizes and the beam center not being on the detector -
458// off a different edge for each LRTB type
459/////////////////////
460
461//function to calculate the overall q-value, given all of the necesary trig inputs
462//NOTE: detector locations passed in are pixels = 0.5cm real space on the detector
463//and are in detector coordinates (1,128) rather than axis values
464//the pixel locations need not be integers, reals are ok inputs
465//sdd is in meters
466//wavelength is in Angstroms
467//
468//returned magnitude of Q is in 1/Angstroms
469//
470Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
471        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
472       
473        Variable dx,dy,qval,two_theta,dist
474               
475        sdd *=100               //convert to cm
476        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
477        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
478        dist = sqrt(dx^2 + dy^2)
479       
480        two_theta = atan(dist/sdd)
481
482        qval = 4*Pi/lam*sin(two_theta/2)
483       
484        return qval
485End
486
487//calculates just the q-value in the x-direction on the detector
488//input/output is the same as CalcQval()
489//ALL inputs are in detector coordinates
490//
491//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector
492//sdd is in meters
493//wavelength is in Angstroms
494//
495// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
496// now properly accounts for qz
497//
498Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
499        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
500
501        Variable qx,qval,phi,dx,dy,dist,two_theta
502       
503        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
504       
505        sdd *=100               //convert to cm
506        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
507        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
508        phi = V_FindPhi(dx,dy)
509       
510        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
511        dist = sqrt(dx^2 + dy^2)
512        two_theta = atan(dist/sdd)
513
514        qx = qval*cos(two_theta/2)*cos(phi)
515       
516        return qx
517End
518
519//calculates just the q-value in the y-direction on the detector
520//input/output is the same as CalcQval()
521//ALL inputs are in detector coordinates
522//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector
523//sdd is in meters
524//wavelength is in Angstroms
525//
526// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
527// now properly accounts for qz
528//
529Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
530        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
531       
532        Variable dy,qval,dx,phi,qy,dist,two_theta
533       
534        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
535       
536        sdd *=100               //convert to cm
537        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
538        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
539        phi = V_FindPhi(dx,dy)
540       
541        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
542        dist = sqrt(dx^2 + dy^2)
543        two_theta = atan(dist/sdd)
544       
545        qy = qval*cos(two_theta/2)*sin(phi)
546       
547        return qy
548End
549
550//calculates just the z-component of the q-vector, not measured on the detector
551//input/output is the same as CalcQval()
552//ALL inputs are in detector coordinates
553//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector
554//sdd is in meters
555//wavelength is in Angstroms
556//
557// not actually used, but here for completeness if anyone asks
558//
559Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
560        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
561       
562        Variable dy,qval,dx,phi,qz,dist,two_theta
563       
564        qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
565       
566        sdd *=100               //convert to cm
567       
568        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
569        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
570        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
571        dist = sqrt(dx^2 + dy^2)
572        two_theta = atan(dist/sdd)
573       
574        qz = qval*sin(two_theta/2)
575       
576        return qz
577End
578
579//phi is defined from +x axis, proceeding CCW around [0,2Pi]
580Threadsafe Function V_FindPhi(vx,vy)
581        variable vx,vy
582       
583        variable phi
584       
585        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
586       
587        // special cases
588        if(vx==0 && vy > 0)
589                return(pi/2)
590        endif
591        if(vx==0 && vy < 0)
592                return(3*pi/2)
593        endif
594        if(vx >= 0 && vy == 0)
595                return(0)
596        endif
597        if(vx < 0 && vy == 0)
598                return(pi)
599        endif
600       
601       
602        if(vx > 0 && vy > 0)
603                return(phi)
604        endif
605        if(vx < 0 && vy > 0)
606                return(phi + pi)
607        endif
608        if(vx < 0 && vy < 0)
609                return(phi + pi)
610        endif
611        if( vx > 0 && vy < 0)
612                return(phi + 2*pi)
613        endif
614       
615        return(phi)
616end
617
618Function V_SphereForm(scale,radius,delrho,bkg,x)                               
619        Variable scale,radius,delrho,bkg
620        Variable x
621       
622        // variables are:                                                       
623        //[0] scale
624        //[1] radius (A)
625        //[2] delrho (A-2)
626        //[3] background (cm-1)
627       
628//      Variable scale,radius,delrho,bkg                               
629//      scale = w[0]
630//      radius = w[1]
631//      delrho = w[2]
632//      bkg = w[3]
633       
634       
635        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
636        // and is rescaled to give [=] cm^-1
637       
638        Variable bes,f,vol,f2
639        //
640        //handle q==0 separately
641        If(x==0)
642                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
643                return(f)
644        Endif
645       
646//      bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
647       
648        bes = 3*sqrt(pi/(2*x*radius))*BesselJ(1.5,x*radius)/(x*radius)
649       
650        vol = 4*pi/3*radius^3
651        f = vol*bes*delrho              // [=] A
652        // normalize to single particle volume, convert to 1/cm
653        f2 = f * f / vol * 1.0e8                // [=] 1/cm
654       
655        return (scale*f2+bkg)   // Scale, then add in the background
656       
657End
658
659Function V_Debye(scale,rg,bkg,x)
660        Variable scale,rg,bkg
661        Variable x
662       
663        // variables are:
664        //[0] scale factor
665        //[1] radius of gyration [A]
666        //[2] background        [cm-1]
667       
668        // calculates (scale*debye)+bkg
669        Variable Pq,qr2
670       
671        qr2=(x*rg)^2
672        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
673       
674        //scale
675        Pq *= scale
676        // then add in the background
677        return (Pq+bkg)
678End
679
680
681
682
683
684
685
686
687
688//
689// these routines bin the 2D q data to 1D I(q). Currently the Qtot is magnitude only, no sign (since
690// it's being binned to I(Q), having a sign makes no sense. If you want the sign, work from qxqyqz
691//
692// first - the DeltaQ step is set as the smaller detector resolution (along tube)
693//       which is different for LR / TB geometry. This is not set in stone.
694//
695// second - each detector is binned separately
696//
697// -- like the routines in CircSectAve, start with 500 points, and trim after binning is done.
698//      you'l end up with < 200 points.
699//
700// the results are in iBin_qxqy, qBin_qxqy, and eBin_qxqy, in the folder passed
701//
702Proc BinAllFrontPanels()
703
704        SetDeltaQ("","FL")
705        SetDeltaQ("","FT")
706
707        Variable binType       
708        ControlInfo/W=VCALC popup_b
709        binType = V_Value               // V_value counts menu items from 1, so 1=1, 2=2, 3=4
710
711        if(binType == 1)
712                V_BinQxQy_to_1D("","FL")
713                V_BinQxQy_to_1D("","FR")
714                V_BinQxQy_to_1D("","FT")
715                V_BinQxQy_to_1D("","FB")
716        endif
717       
718        if(binType == 2)       
719                V_BinQxQy_to_1D("","FLR")
720                V_BinQxQy_to_1D("","FTB")
721        endif
722
723        if(binType == 3)
724                V_BinQxQy_to_1D("","FLRTB")
725        endif
726
727// TODO -- this is only a temporary fix for slit mode   
728        if(binType == 4)
729                /// this is for a tall, narrow slit mode       
730                V_fBinDetector_byRows("FL")
731                V_fBinDetector_byRows("FR")
732                V_fBinDetector_byRows("FT")
733                V_fBinDetector_byRows("FB")
734        endif
735               
736End
737
738
739//TODO -- folderStr is ignored in this function
740Function SetDeltaQ(folderStr,type)
741        String folderStr,type
742       
743        WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "det_"+type)         // 2D detector data
744       
745        Variable xDim,yDim,delQ
746       
747        xDim=DimSize(inten,0)
748        yDim=DimSize(inten,1)
749       
750        if(xDim<yDim)
751                WAVE qx = $("root:Packages:NIST:VSANS:VCALC:" + "qx_"+type)
752                delQ = abs(qx[0][0] - qx[1][0])/2
753        else
754                WAVE qy = $("root:Packages:NIST:VSANS:VCALC:" + "qy_"+type)
755                delQ = abs(qy[0][1] - qy[0][0])/2
756        endif
757       
758        // set the global
759        Variable/G $("root:Packages:NIST:VSANS:VCALC:" + "delQ_"+type) = delQ
760//      Print "SET delQ = ",delQ," for ",type
761       
762        return(0)
763end
764
765
766//TODO -- need a switch here to dispatch to the averaging type
767Proc V_BinQxQy_to_1D(folderStr,type)
768        String folderStr
769        String type
770//      Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
771//      Prompt type,"detector identifier"
772
773
774        V_fDoBinning_QxQy2D("", type)
775
776
777/// this is for a tall, narrow slit mode       
778//      V_fBinDetector_byRows(type)
779       
780End
781
782Proc V_Graph_1D_detType(folderStr,type)
783        String folderStr,type
784       
785        SetDataFolder root:Packages:NIST:VSANS:VCALC
786       
787        Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
788        ModifyGraph mirror=2,grid=1,log=1
789        ModifyGraph mode=4,marker=19,msize=2
790//      ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy)                // for simulations, I don't have 2D uncertainty
791        ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
792        legend
793       
794        SetDataFolder root:
795
796End
797
798
799// see the equivalent function in PlotUtils2D_v40.ipf
800//
801//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
802//
803// this has been modeified to accept different detector panels and to take arrays
804// -- type = FL or FR or...other panel identifiers
805//
806// TODO "iErr" is all messed up since it doesn't really apply here for data that is not 2D simulation
807//
808//
809Function V_fDoBinning_QxQy2D(folderStr,type)
810        String folderStr,type
811
812        // TODO: folderStr is ignored here
813        folderStr = ""
814       
815        Variable nSets = 0
816        Variable xDim,yDim
817        Variable ii,jj
818        Variable qVal,nq,var,avesq,aveisq
819        Variable binIndex,val
820       
821       
822        SetDataFolder root:Packages:NIST:VSANS:VCALC
823       
824// now switch on the type to determine which waves to declare and create
825// since there may be more than one panel to step through. There may be two, there may be four
826//
827
828        strswitch(type) // string switch
829                case "FL":              // execute if case matches expression
830                case "FR":
831                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "delQ_FL")
832                        WAVE inten = $("det_"+type)             // 2D detector data
833                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
834                        Wave qTotal = $("qTot_"+type)                   // 2D q-values
835                        nSets = 1
836                        break   
837                                                               
838                case "FT":             
839                case "FB":
840                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "delQ_FT")
841                        WAVE inten = $("det_"+type)             // 2D detector data
842                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
843                        Wave qTotal = $("qTot_"+type)                   // 2D q-values
844                        nSets = 1
845                        break
846                       
847                case "ML":             
848                case "MR":             
849                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "delQ_ML")
850                        WAVE inten = $("det_"+type)             // 2D detector data
851                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
852                        Wave qTotal = $("qTot_"+type)                   // 2D q-values
853                        nSets = 1
854                        break   
855                                       
856                case "MT":             
857                case "MB":             
858                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "delQ_MT")
859                        WAVE inten = $("det_"+type)             // 2D detector data
860                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
861                        Wave qTotal = $("qTot_"+type)                   // 2D q-values
862                        nSets = 1
863                        break   
864                                       
865                case "B":               
866                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "delQ_B")
867                        WAVE inten = $("det_"+type)             // 2D detector data
868                        WAVE/Z iErr = $("iErr_"+type)                   // 2D errors -- may not exist, especially for simulation
869                        Wave qTotal = $("qTot_"+type)                   // 2D q-values
870                        nSets = 1
871                        break   
872                       
873                case "FLR":
874                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "delQ_FL")
875                        WAVE inten = $("det_"+"FL")             // 2D detector data
876                        WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation
877                        Wave qTotal = $("qTot_"+"FL")                   // 2D q-values
878                        WAVE inten2 = $("det_"+"FR")            // 2D detector data
879                        WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation
880                        Wave qTotal2 = $("qTot_"+"FR")                  // 2D q-values
881                        nSets = 2
882                        break                   
883               
884                case "FTB":
885                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "delQ_FT")
886                        WAVE inten = $("det_"+"FT")             // 2D detector data
887                        WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation
888                        Wave qTotal = $("qTot_"+"FT")                   // 2D q-values
889                        WAVE inten2 = $("det_"+"FB")            // 2D detector data
890                        WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation
891                        Wave qTotal2 = $("qTot_"+"FB")                  // 2D q-values
892                        nSets = 2
893                        break           
894               
895                case "FLRTB":
896                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "delQ_FL")
897                        WAVE inten = $("det_"+"FL")             // 2D detector data
898                        WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation
899                        Wave qTotal = $("qTot_"+"FL")                   // 2D q-values
900                        WAVE inten2 = $("det_"+"FR")            // 2D detector data
901                        WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation
902                        Wave qTotal2 = $("qTot_"+"FR")                  // 2D q-values
903                        WAVE inten3 = $("det_"+"FT")            // 2D detector data
904                        WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation
905                        Wave qTotal3 = $("qTot_"+"FT")                  // 2D q-values
906                        WAVE inten4 = $("det_"+"FB")            // 2D detector data
907                        WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation
908                        Wave qTotal4 = $("qTot_"+"FB")                  // 2D q-values
909                        nSets = 4
910                        break           
911                       
912
913                case "MLR":
914                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "delQ_ML")
915                        WAVE inten = $("det_"+"ML")             // 2D detector data
916                        WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation
917                        Wave qTotal = $("qTot_"+"ML")                   // 2D q-values
918                        WAVE inten2 = $("det_"+"MR")            // 2D detector data
919                        WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation
920                        Wave qTotal2 = $("qTot_"+"MR")                  // 2D q-values
921                        nSets = 2
922                        break                   
923               
924                case "MTB":
925                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "delQ_MT")
926                        WAVE inten = $("det_"+"MT")             // 2D detector data
927                        WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation
928                        Wave qTotal = $("qTot_"+"MT")                   // 2D q-values
929                        WAVE inten2 = $("det_"+"MB")            // 2D detector data
930                        WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation
931                        Wave qTotal2 = $("qTot_"+"MB")                  // 2D q-values
932                        nSets = 2
933                        break                           
934               
935                case "MLRTB":
936                        NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "delQ_ML")
937                        WAVE inten = $("det_"+"ML")             // 2D detector data
938                        WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation
939                        Wave qTotal = $("qTot_"+"ML")                   // 2D q-values
940                        WAVE inten2 = $("det_"+"MR")            // 2D detector data
941                        WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation
942                        Wave qTotal2 = $("qTot_"+"MR")                  // 2D q-values
943                        WAVE inten3 = $("det_"+"MT")            // 2D detector data
944                        WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation
945                        Wave qTotal3 = $("qTot_"+"MT")                  // 2D q-values
946                        WAVE inten4 = $("det_"+"MB")            // 2D detector data
947                        WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation
948                        Wave qTotal4 = $("qTot_"+"MB")                  // 2D q-values
949                        nSets = 4
950                        break                                                                   
951                                       
952                default:
953                        nSets = 0                                                       // optional default expression executed
954                        Print "ERROR   ---- type is not recognized "
955        endswitch
956
957//      Print "delQ = ",delQ," for ",type
958
959        if(nSets == 0)
960                return(0)
961        endif
962
963
964//TODO: properly define the errors here - I'll have this if I do the simulation
965        if(WaveExists(iErr)==0)
966                Duplicate/O inten,iErr
967                Wave iErr=iErr
968//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
969                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
970        endif
971
972        nq = 600
973
974        // note that the back panel of 320x320 (1mm res) results in 447 data points!
975        // - so I upped nq to 600
976       
977//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
978        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin_qxqy"+"_"+type)
979        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"qBin_qxqy"+"_"+type)
980        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"nBin_qxqy"+"_"+type)
981        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin2_qxqy"+"_"+type)
982        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin_qxqy"+"_"+type)
983        Make/O/D/N=(nq)  $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin2D_qxqy"+"_"+type)
984       
985        Wave iBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin_qxqy_"+type)
986        Wave qBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"qBin_qxqy"+"_"+type)
987        Wave nBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"nBin_qxqy"+"_"+type)
988        Wave iBin2_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin2_qxqy"+"_"+type)
989        Wave eBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin_qxqy"+"_"+type)
990        Wave eBin2D_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin2D_qxqy"+"_"+type)
991       
992       
993//      delQ = abs(sqrt(qx[2]^2+qy[2]^2+qz[2]^2) - sqrt(qx[1]^2+qy[1]^2+qz[1]^2))               //use bins of 1 pixel width
994// TODO: not sure if I want to so dQ in x or y direction...
995        // the short dimension is the 8mm tubes, use this direction as dQ?
996        // but don't use the corner of the detector, since dQ will be very different on T/B or L/R due to the location of [0,0]
997        // WRT the beam center. use qx or qy directly. Still not happy with this way...
998
999
1000        qBin_qxqy[] =  p*delQ   
1001        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1002
1003        iBin_qxqy = 0
1004        iBin2_qxqy = 0
1005        eBin_qxqy = 0
1006        eBin2D_qxqy = 0
1007        nBin_qxqy = 0   //number of intensities added to each bin
1008
1009// now there are situations of:
1010// 1 panel
1011// 2 panels
1012// 4 panels
1013//
1014// this needs to be a double loop now...
1015
1016// use set 1 (no number) only
1017        if(nSets >= 1)
1018                xDim=DimSize(inten,0)
1019                yDim=DimSize(inten,1)
1020       
1021                for(ii=0;ii<xDim;ii+=1)
1022                        for(jj=0;jj<yDim;jj+=1)
1023                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1024                                qVal = qTotal[ii][jj]
1025                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1026                                val = inten[ii][jj]
1027                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
1028                                        iBin_qxqy[binIndex] += val
1029                                        iBin2_qxqy[binIndex] += val*val
1030                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1031                                        nBin_qxqy[binIndex] += 1
1032                                endif
1033                        endfor
1034                endfor
1035               
1036        endif
1037
1038// add in set 2 (set 1 already done)
1039        if(nSets >= 2)
1040                xDim=DimSize(inten2,0)
1041                yDim=DimSize(inten2,1)
1042       
1043                for(ii=0;ii<xDim;ii+=1)
1044                        for(jj=0;jj<yDim;jj+=1)
1045                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1046                                qVal = qTotal2[ii][jj]
1047                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1048                                val = inten2[ii][jj]
1049                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
1050                                        iBin_qxqy[binIndex] += val
1051                                        iBin2_qxqy[binIndex] += val*val
1052                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1053                                        nBin_qxqy[binIndex] += 1
1054                                endif
1055                        endfor
1056                endfor
1057               
1058        endif
1059
1060// add in set 3 and 4 (set 1 and 2already done)
1061        if(nSets == 4)
1062                xDim=DimSize(inten3,0)
1063                yDim=DimSize(inten3,1)
1064       
1065                for(ii=0;ii<xDim;ii+=1)
1066                        for(jj=0;jj<yDim;jj+=1)
1067                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1068                                qVal = qTotal3[ii][jj]
1069                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1070                                val = inten3[ii][jj]
1071                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
1072                                        iBin_qxqy[binIndex] += val
1073                                        iBin2_qxqy[binIndex] += val*val
1074                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1075                                        nBin_qxqy[binIndex] += 1
1076                                endif
1077                        endfor
1078                endfor
1079               
1080               
1081                xDim=DimSize(inten4,0)
1082                yDim=DimSize(inten4,1)
1083       
1084                for(ii=0;ii<xDim;ii+=1)
1085                        for(jj=0;jj<yDim;jj+=1)
1086                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1087                                qVal = qTotal4[ii][jj]
1088                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1089                                val = inten4[ii][jj]
1090                                if (numType(val)==0)            //count only the good points, ignore Nan or Inf
1091                                        iBin_qxqy[binIndex] += val
1092                                        iBin2_qxqy[binIndex] += val*val
1093                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1094                                        nBin_qxqy[binIndex] += 1
1095                                endif
1096                        endfor
1097                endfor
1098               
1099        endif
1100
1101
1102// after looping through all of the data on the panels, calculate errors on I(q),
1103// just like in CircSectAve.ipf
1104        for(ii=0;ii<nq;ii+=1)
1105                if(nBin_qxqy[ii] == 0)
1106                        //no pixels in annuli, data unknown
1107                        iBin_qxqy[ii] = 0
1108                        eBin_qxqy[ii] = 1
1109                        eBin2D_qxqy[ii] = NaN
1110                else
1111                        if(nBin_qxqy[ii] <= 1)
1112                                //need more than one pixel to determine error
1113                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1114                                eBin_qxqy[ii] = 1
1115                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1116                        else
1117                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
1118                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1119                                avesq = iBin_qxqy[ii]^2
1120                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1121                                var = aveisq-avesq
1122                                if(var<=0)
1123                                        eBin_qxqy[ii] = 1e-6
1124                                else
1125                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1126                                endif
1127                                // and calculate as it is propagated pixel-by-pixel
1128                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1129                        endif
1130                endif
1131        endfor
1132       
1133        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
1134       
1135        // find the last non-zero point, working backwards
1136        val=nq
1137        do
1138                val -= 1
1139        while((nBin_qxqy[val] == 0) && val > 0)
1140       
1141//      print val, nBin_qxqy[val]
1142        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1143
1144        if(val == 0)
1145                // all the points were deleted
1146                return(0)
1147        endif
1148       
1149       
1150        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1151        // find the first non-zero point, working forwards
1152        val = -1
1153        do
1154                val += 1
1155        while(nBin_qxqy[val] == 0)     
1156        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1157
1158        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1159        val = numpnts(nBin_qxqy)-1
1160        do
1161                if(nBin_qxqy[val] == 0)
1162                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1163                endif
1164                val -= 1
1165        while(val>0)
1166       
1167        SetDataFolder root:
1168       
1169        return(0)
1170End
1171
1172////////////to plot the (4) 2D panels and to plot the I(Q) data on the same plot
1173//
1174// ** but now I need to check and see if these waves exist before trying to append them
1175// since the panels may bave been combined when binned - rather than all separate.
1176//
1177// TODO
1178// -- so maybe I want to clear the traces from the graph?
1179// -- set a flag on the panel to know how the binning is applied?
1180//
1181Window Front_IQ_Graph() : Graph
1182
1183        Variable binType
1184       
1185        ControlInfo/W=VCALC popup_b
1186        binType = V_Value               // V_value counts menu items from 1, so 1=1, 2=2, 3=4
1187
1188        SetDataFolder root:Packages:NIST:VSANS:VCALC
1189
1190        if(binType==1)
1191                ClearIQIfDisplayed("FLRTB")
1192                ClearIQIfDisplayed("FLR")
1193                ClearIQIfDisplayed("FTB")
1194               
1195                SetDataFolder root:Packages:NIST:VSANS:VCALC
1196                CheckDisplayed/W=VCALC#Panels_IQ iBin_qxqy_FL
1197               
1198                if(V_flag==0)
1199                        AppendtoGraph/W=VCALC#Panels_IQ iBin_qxqy_FL vs qBin_qxqy_FL
1200                        AppendToGraph/W=VCALC#Panels_IQ iBin_qxqy_FR vs qBin_qxqy_FR
1201                        AppendToGraph/W=VCALC#Panels_IQ iBin_qxqy_FT vs qBin_qxqy_FT
1202                        AppendToGraph/W=VCALC#Panels_IQ iBin_qxqy_FB vs qBin_qxqy_FB
1203                        ModifyGraph/W=VCALC#Panels_IQ mode=4
1204                        ModifyGraph/W=VCALC#Panels_IQ marker=19
1205                        ModifyGraph/W=VCALC#Panels_IQ rgb(iBin_qxqy_FL)=(39321,26208,1),rgb(iBin_qxqy_FB)=(2,39321,1),rgb(iBin_qxqy_FR)=(39321,26208,1),rgb(iBin_qxqy_FT)=(2,39321,1)
1206                        ModifyGraph/W=VCALC#Panels_IQ msize=2
1207                        ModifyGraph/W=VCALC#Panels_IQ muloffset(iBin_qxqy_FL)={0,4},muloffset(iBin_qxqy_FB)={0,2},muloffset(iBin_qxqy_FR)={0,8}
1208                        ModifyGraph/W=VCALC#Panels_IQ grid=1
1209                        ModifyGraph/W=VCALC#Panels_IQ log=1
1210                        ModifyGraph/W=VCALC#Panels_IQ mirror=2
1211                        Label/W=VCALC#Panels_IQ left "Intensity (1/cm)"
1212                        Label/W=VCALC#Panels_IQ bottom "Q (1/A)"
1213                endif   
1214                               
1215        endif
1216
1217        if(binType==2)
1218                ClearIQIfDisplayed("FLRTB")
1219                ClearIQIfDisplayed("FT")       
1220                ClearIQIfDisplayed("FL")       
1221                ClearIQIfDisplayed("FR")       
1222                ClearIQIfDisplayed("FB")
1223       
1224
1225                SetDataFolder root:Packages:NIST:VSANS:VCALC
1226                CheckDisplayed/W=VCALC#Panels_IQ iBin_qxqy_FLR
1227               
1228                if(V_flag==0)
1229                        AppendtoGraph/W=VCALC#Panels_IQ iBin_qxqy_FLR vs qBin_qxqy_FLR
1230                        AppendToGraph/W=VCALC#Panels_IQ iBin_qxqy_FTB vs qBin_qxqy_FTB
1231                        ModifyGraph/W=VCALC#Panels_IQ mode=4
1232                        ModifyGraph/W=VCALC#Panels_IQ marker=19
1233                        ModifyGraph/W=VCALC#Panels_IQ rgb(iBin_qxqy_FLR)=(39321,26208,1),rgb(iBin_qxqy_FTB)=(2,39321,1)
1234                        ModifyGraph/W=VCALC#Panels_IQ msize=2
1235                        ModifyGraph/W=VCALC#Panels_IQ muloffset(iBin_qxqy_FLR)={0,2}
1236                        ModifyGraph/W=VCALC#Panels_IQ grid=1
1237                        ModifyGraph/W=VCALC#Panels_IQ log=1
1238                        ModifyGraph/W=VCALC#Panels_IQ mirror=2
1239                        Label/W=VCALC#Panels_IQ left "Intensity (1/cm)"
1240                        Label/W=VCALC#Panels_IQ bottom "Q (1/A)"
1241                endif   
1242                       
1243        endif
1244       
1245        if(binType==3)
1246                ClearIQIfDisplayed("FLR")
1247                ClearIQIfDisplayed("FTB")       
1248                ClearIQIfDisplayed("FT")       
1249                ClearIQIfDisplayed("FL")       
1250                ClearIQIfDisplayed("FR")       
1251                ClearIQIfDisplayed("FB")       
1252       
1253                SetDataFolder root:Packages:NIST:VSANS:VCALC
1254                CheckDisplayed/W=VCALC#Panels_IQ iBin_qxqy_FLRTB
1255               
1256                if(V_flag==0)
1257                        AppendtoGraph/W=VCALC#Panels_IQ iBin_qxqy_FLRTB vs qBin_qxqy_FLRTB
1258                        ModifyGraph/W=VCALC#Panels_IQ mode=4
1259                        ModifyGraph/W=VCALC#Panels_IQ marker=19
1260                        ModifyGraph/W=VCALC#Panels_IQ rgb(iBin_qxqy_FLRTB)=(39321,26208,1)
1261                        ModifyGraph/W=VCALC#Panels_IQ msize=2
1262                        ModifyGraph/W=VCALC#Panels_IQ grid=1
1263                        ModifyGraph/W=VCALC#Panels_IQ log=1
1264                        ModifyGraph/W=VCALC#Panels_IQ mirror=2
1265                        Label/W=VCALC#Panels_IQ left "Intensity (1/cm)"
1266                        Label/W=VCALC#Panels_IQ bottom "Q (1/A)"
1267                endif   
1268                       
1269        endif
1270       
1271        SetDataFolder root:
1272       
1273EndMacro
1274
1275Function        ClearIQIfDisplayed(type)
1276        String type
1277       
1278        SetDataFolder root:Packages:NIST:VSANS:VCALC
1279        CheckDisplayed/W=VCALC#Panels_IQ $("iBin_qxqy_"+type)
1280        if(V_flag==1)
1281                RemoveFromGraph/W=VCALC#Panels_IQ $("iBin_qxqy_"+type)
1282        endif
1283        SetDataFolder root:
1284       
1285        return(0)
1286end
1287
1288Window Table_of_QBins() : Table
1289        PauseUpdate; Silent 1           // building window...
1290        String fldrSav0= GetDataFolder(1)
1291        SetDataFolder root:Packages:NIST:VSANS:VCALC:
1292        Edit/W=(5,44,771,898) qBin_qxqy_FL,qBin_qxqy_FR,qBin_qxqy_FT,qBin_qxqy_FB,qBin_qxqy_FLR
1293        AppendToTable qBin_qxqy_FTB,qBin_qxqy_FLRTB
1294        ModifyTable format(Point)=1,width(qBin_qxqy_FLR)=136,width(qBin_qxqy_FLRTB)=120
1295        SetDataFolder fldrSav0
1296EndMacro
1297
1298
1299
1300
1301
1302
1303
1304//////NOTE///
1305//// some chunks of the code here have been trimmed out for fun
1306//Function V_CircularAverageTo1D(type)
1307//      String type
1308//             
1309//////// get information about the detector (in type folder) that is needed for reduction
1310//// pixel dimensions
1311//// beam center
1312//// distances
1313//// wavelength
1314//// wavelength spread
1315////
1316//
1317//      String destPath = "root:"
1318//
1319////    NVAR pixelsX = root:myGlobals:gNPixelsX
1320////    NVAR pixelsY = root:myGlobals:gNPixelsY
1321//     
1322////    pixelsX = 48
1323////    pixelsY = 250
1324//     
1325//      // this is for non-linear corrections not applicable?
1326////    xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
1327////    ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
1328//     
1329//      // beam center, in pixels
1330//      x0 = 60         //reals[16]
1331//      y0 = 125                //reals[17]
1332////    //detector calibration constants
1333////    sx = reals[10]          //mm/pixel (x)
1334////    sx3 = reals[11]         //nonlinear coeff
1335////    sy = reals[13]          //mm/pixel (y)
1336////    sy3 = reals[14]         //nonlinear coeff
1337//     
1338////    dtsize = 10*reals[20]           //det size in mm
1339//      dtdist = 1000           //1000*reals[18]                // det distance in mm
1340//
1341//
1342/////// decide how the binning width is to be determined
1343//     
1344////    NVAR binWidth=root:Packages:NIST:gBinWidth
1345//     
1346//      dr = 1          //binWidth              // ***********annulus width set by user, default is one***********
1347//      ddr = 4         //dr*sx         //step size, in mm (this value should be passed to the resolution calculation, not dr 18NOV03)
1348//             
1349//      Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,dxbm,dybm,ii
1350//      Variable phi_rad,dphi_rad,phi_x,phi_y
1351//      Variable forward,mirror
1352//     
1353////// do I need to pick sides?? probably, for consistency, but confusing nomenclature now
1354////    String side = StringByKey("SIDE",keyListStr,"=",";")
1355//
1356///////// keep the sector calculations, I'll want to do this...
1357////// 
1358////    if(!isCircular)         //must be sector avg (rectangular not sent to this function)
1359////            //convert from degrees to radians
1360////            phi_rad = (Pi/180)*NumberByKey("PHI",keyListStr,"=",";")
1361////            dphi_rad = (Pi/180)*NumberByKey("DPHI",keyListStr,"=",";")
1362////            //create cartesian values for unit vector in phi direction
1363////            phi_x = cos(phi_rad)
1364////            phi_y = sin(phi_rad)
1365////    Endif
1366//     
1367//      /// data wave is data in the current folder which was set at the top of the function
1368////    WAVE data=$(destPath + ":data")
1369//      Make/O/N=(pixelsX,pixelsY)  $(destPath + ":data")
1370//      Wave data = $(destPath + ":data")
1371//      data = 1
1372//
1373//      //Check for the existence of the mask, if not, make one (local to this folder) that is null
1374//     
1375//      if(WaveExists($"root:Packages:NIST:MSK:data") == 0)
1376//              Print "There is no mask file loaded (WaveExists)- the data is not masked"
1377//              Make/O/N=(pixelsX,pixelsY) $(destPath + ":mask")
1378//              Wave mask = $(destPath + ":mask")
1379//              mask = 0
1380//      else
1381//              Wave mask=$"root:Packages:NIST:MSK:data"
1382//      Endif
1383//     
1384//      //
1385//      //pixels within rcentr of beam center are broken into 9 parts (units of mm)
1386//      rcentr = 100            //original
1387////    rcentr = 0
1388//      // values for error if unable to estimate value
1389//      //large_num = 1e10
1390//      large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1
1391//      small_num = 1e-10
1392//     
1393//      // output wave are expected to exist (?) initialized to zero, what length?
1394//      // 200 points on VAX --- use 300 here, or more if SAXS data is used with 1024x1024 detector (1000 pts seems good)
1395//      Variable defWavePts=500
1396//      Make/O/N=(defWavePts) $(destPath + ":qval"),$(destPath + ":aveint")
1397//      Make/O/N=(defWavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave")
1398//      Make/O/N=(defWavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar")
1399//
1400//      WAVE qval = $(destPath + ":qval")
1401//      WAVE aveint = $(destPath + ":aveint")
1402//      WAVE ncells = $(destPath + ":ncells")
1403//      WAVE dsq = $(destPath + ":dsq")
1404//      WAVE sigave = $(destPath + ":sigave")
1405//      WAVE qbar = $(destPath + ":QBar")
1406//      WAVE sigmaq = $(destPath + ":SigmaQ")
1407//      WAVE fsubs = $(destPath + ":fSubS")
1408//     
1409//      qval = 0
1410//      aveint = 0
1411//      ncells = 0
1412//      dsq = 0
1413//      sigave = 0
1414//      qbar = 0
1415//      sigmaq = 0
1416//      fsubs = 0
1417//
1418//      dtdis2 = dtdist^2
1419//      nq = 1
1420//      xoffst=0
1421//      //distance of beam center from detector center
1422//////  // the linearity corrections for the 2D Ordela detectors are applied from the center of the detector,
1423//      // so figure where that is, relative to the beam center.
1424////    dxbm = V_FX(x0,sx3,xcenter,sx)
1425////    dybm = V_FY(y0,sy3,ycenter,sy)
1426//             
1427//      //BEGIN AVERAGE **********
1428//      Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1
1429//      Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p
1430//     
1431//      // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
1432//      // loop index corresponds to FORTRAN (old code)
1433//      // and the IGOR array indices must be adjusted (-1) to the correct address
1434//     
1435//      /////
1436//      //// need to add in here a step that calculates the q-values, and bins based on q-value
1437//      //// since the 4 panels are not in the same plane (but relatively close)
1438//      //
1439//      //    if we're always using pairs in the same plane, then binning in r is OK
1440//      ////
1441//     
1442//      ii=1
1443//      do
1444//              xi = ii
1445//              dxi = V_FX(xi,sx3,xcenter,sx)
1446//              dx = dxi-dxbm           //dx and dy are in mm
1447//             
1448//              jj = 1
1449//              do
1450//                      data_pixel = data[ii-1][jj-1]           //assign to local variable
1451//                      yj = jj
1452//                      dyj = V_FY(yj,sy3,ycenter,sy)
1453//                      dy = dyj - dybm
1454//                      if(!(mask[ii-1][jj-1]))                 //masked pixels = 1, skip if masked (this way works...)
1455//                              dr2 = (dx^2 + dy^2)^(0.5)               //distance from beam center NOTE dr2 used here - dr used above
1456//                              if(dr2>rcentr)          //keep pixel whole
1457//                                      nd = 1
1458//                                      fd = 1
1459//                              else                            //break pixel into 9 equal parts
1460//                                      nd = 3
1461//                                      fd = 2
1462//                              endif
1463//                              nd2 = nd^2
1464//                              ll = 1          //"el-el" loop index
1465//                              do
1466//                                      dxx = dx + (ll - fd)*sx/3
1467//                                      kk = 1
1468//                                      do
1469//                                              dyy = dy + (kk - fd)*sy/3
1470//                                              if(isCircular)
1471//                                                      //circular average, use all pixels
1472//                                                      //(increment)
1473//                                                      nq = V_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1474////                                            else
1475////                                                    //a sector average - determine azimuthal angle
1476////                                                    dphi_p = V_dphi_pixel(dxx,dyy,phi_x,phi_y)
1477////                                                    if(dphi_p < dphi_rad)
1478////                                                            forward = 1                     //within forward sector
1479////                                                    else
1480////                                                            forward = 0
1481////                                                    Endif
1482////                                                    if((Pi - dphi_p) < dphi_rad)
1483////                                                            mirror = 1              //within mirror sector
1484////                                                    else
1485////                                                            mirror = 0
1486////                                                    Endif
1487////                                                    //check if pixel lies within allowed sector(s)
1488////                                                    if(cmpstr(side,"both")==0)              //both sectors
1489////                                                            if ( mirror || forward)
1490////                                                                    //increment
1491////                                                                    nq = V_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1492////                                                            Endif
1493////                                                    else
1494////                                                            if(cmpstr(side,"right")==0)             //forward sector only
1495////                                                                    if(forward)
1496////                                                                            //increment
1497////                                                                            nq = V_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1498////                                                                    Endif
1499////                                                            else                    //mirror sector only
1500////                                                                    if(mirror)
1501////                                                                            //increment
1502////                                                                            nq = V_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1503////                                                                    Endif
1504////                                                            Endif
1505////                                                    Endif           //allowable sectors
1506//                                              Endif   //circular or sector check
1507//                                              kk+=1
1508//                                      while(kk<=nd)
1509//                                      ll += 1
1510//                              while(ll<=nd)
1511//                      Endif           //masked pixel check
1512//                      jj += 1
1513//              while (jj<=pixelsY)
1514//              ii += 1
1515//      while(ii<=pixelsX)              //end of the averaging
1516//             
1517//      //compute q-values and errors
1518//      Variable ntotal,rr,theta,avesq,aveisq,var
1519//     
1520//      lambda = reals[26]
1521//      ntotal = 0
1522//      kk = 1
1523//      do
1524//              rr = (2*kk-1)*ddr/2
1525//              theta = 0.5*atan(rr/dtdist)
1526//              qval[kk-1] = (4*Pi/lambda)*sin(theta)
1527//              if(ncells[kk-1] == 0)
1528//                      //no pixels in annuli, data unknown
1529//                      aveint[kk-1] = 0
1530//                      sigave[kk-1] = large_num
1531//              else
1532//                      if(ncells[kk-1] <= 1)
1533//                              //need more than one pixel to determine error
1534//                              aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
1535//                              sigave[kk-1] = large_num
1536//                      else
1537//                              //assume that the intensity in each pixel in annuli is normally
1538//                              // distributed about mean...
1539//                              aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
1540//                              avesq = aveint[kk-1]^2
1541//                              aveisq = dsq[kk-1]/ncells[kk-1]
1542//                              var = aveisq-avesq
1543//                              if(var<=0)
1544//                                      sigave[kk-1] = small_num
1545//                              else
1546//                                      sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1))
1547//                              endif
1548//                      endif
1549//              endif
1550//              ntotal += ncells[kk-1]
1551//              kk+=1
1552//      while(kk<=nq)
1553//     
1554//      //Print "NQ = ",nq
1555//      // data waves were defined as 300 points (=defWavePts), but now have less than that (nq) points
1556//      // use DeletePoints to remove junk from end of waves
1557//      //WaveStats would be a more foolproof implementation, to get the # points in the wave
1558//      Variable startElement,numElements
1559//      startElement = nq
1560//      numElements = defWavePts - startElement
1561//      DeletePoints startElement,numElements, qval,aveint,ncells,dsq,sigave
1562//     
1563//      //////////////end of VAX sector_ave()
1564//             
1565//      //angle dependent transmission correction
1566//      Variable uval,arg,cos_th
1567//      lambda = reals[26]
1568//      trans = reals[4]
1569//
1570////
1571////  The transmission correction is now done at the ADD step, in DetCorr()
1572////   
1573////    ////this section is the trans_correct() VAX routine
1574////    if(trans<0.1)
1575////            Print "***transmission is less than 0.1*** and is a significant correction"
1576////    endif
1577////    if(trans==0)
1578////            Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
1579////            trans = 1
1580////    endif
1581////    //optical thickness
1582////    uval = -ln(trans)               //use natural logarithm
1583////    //apply correction to aveint[]
1584////    //index from zero here, since only working with IGOR waves
1585////    ii=0
1586////    do
1587////            theta = 2*asin(lambda*qval[ii]/(4*pi))
1588////            cos_th = cos(theta)
1589////            arg = (1-cos_th)/cos_th
1590////            if((uval<0.01) || (cos_th>0.99))                //OR
1591////                    //small arg, approx correction
1592////                    aveint[ii] /= 1-0.5*uval*arg
1593////            else
1594////                    //large arg, exact correction
1595////                    aveint[ii] /= (1-exp(-uval*arg))/(uval*arg)
1596////            endif
1597////            ii+=1
1598////    while(ii<nq)
1599////    //end of transmission/pathlength correction
1600//
1601//// ***************************************************************
1602////
1603//// Do the extra 3 columns of resolution calculations starting here.
1604////
1605//// ***************************************************************
1606//
1607//      Variable L2 = reals[18]
1608//      Variable BS = reals[21]
1609//      Variable S1 = reals[23]
1610//      Variable S2 = reals[24]
1611//      Variable L1 = reals[25]
1612//      lambda = reals[26]
1613//      Variable lambdaWidth = reals[27]
1614//      String detStr=textRead[9]
1615//     
1616//      Variable usingLenses = reals[28]                //new 2007
1617//
1618//      //Two parameters DDET and APOFF are instrument dependent.  Determine
1619//      //these from the instrument name in the header.
1620//      //From conversation with JB on 01.06.99 these are the current
1621//      //good values
1622//
1623//      Variable DDet
1624//      NVAR apOff = root:myGlobals:apOff               //in cm
1625//     
1626////    DDet = DetectorPixelResolution(fileStr,detStr)          //needs detector type and beamline
1627//      //note that reading the detector pixel size from the header ASSUMES SQUARE PIXELS! - Jan2008
1628//      DDet = reals[10]/10                     // header value (X) is in mm, want cm here
1629//     
1630//     
1631//      //Width of annulus used for the average is gotten from the
1632//      //input dialog before.  This also must be passed to the resolution
1633//      //calculator. Currently the default is dr=1 so just keeping that.
1634//
1635//      //Go from 0 to nq doing the calc for all three values at
1636//      //every Q value
1637//
1638//      ii=0
1639//
1640//      Variable ret1,ret2,ret3
1641//      do
1642//      // commented out for compiler
1643////            getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,usingLenses,ret1,ret2,ret3)
1644//              sigmaq[ii] = ret1       
1645//              qbar[ii] = ret2
1646//              fsubs[ii] = ret3       
1647//              ii+=1
1648//      while(ii<nq)
1649//      DeletePoints startElement,numElements, sigmaq, qbar, fsubs
1650//
1651//// End of resolution calculations
1652//// ***************************************************************
1653//     
1654//      //Plot the data in the Plot_1d window
1655////    Avg_1D_Graph(aveint,qval,sigave)                //commented out for compiler
1656//
1657//      //get rid of the default mask, if one was created (it is in the current folder)
1658//      //don't just kill "mask" since it might be pointing to the one in the MSK folder
1659//      Killwaves/Z $(destPath+":mask")
1660//     
1661//      //return to root folder (redundant)
1662//      SetDataFolder root:
1663//     
1664//      Return 0
1665//End
1666
1667////returns nq, new number of q-values
1668////arrays aveint,dsq,ncells are also changed by this function
1669////
1670//Function V_IncrementPixel(dataPixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1671//      Variable dataPixel,ddr,dxx,dyy
1672//      Wave aveint,dsq,ncells
1673//      Variable nq,nd2
1674//     
1675//      Variable ir
1676//     
1677//      ir = trunc(sqrt(dxx*dxx+dyy*dyy)/ddr)+1
1678//      if (ir>nq)
1679//              nq = ir         //resets maximum number of q-values
1680//      endif
1681//      aveint[ir-1] += dataPixel/nd2           //ir-1 must be used, since ir is physical
1682//      dsq[ir-1] += dataPixel*dataPixel/nd2
1683//      ncells[ir-1] += 1/nd2
1684//     
1685//      Return nq
1686//End
1687//
1688////function determines azimuthal angle dphi that a vector connecting
1689////center of detector to pixel makes with respect to vector
1690////at chosen azimuthal angle phi -> [cos(phi),sin(phi)] = [phi_x,phi_y]
1691////dphi is always positive, varying from 0 to Pi
1692////
1693//Function V_dphi_pixel(dxx,dyy,phi_x,phi_y)
1694//      Variable dxx,dyy,phi_x,phi_y
1695//     
1696//      Variable val,rr,dot_prod
1697//     
1698//      rr = sqrt(dxx^2 + dyy^2)
1699//      dot_prod = (dxx*phi_x + dyy*phi_y)/rr
1700//      //? correct for roundoff error? - is this necessary in IGOR, w/ double precision?
1701//      if(dot_prod > 1)
1702//              dot_prod =1
1703//      Endif
1704//      if(dot_prod < -1)
1705//              dot_prod = -1
1706//      Endif
1707//     
1708//      val = acos(dot_prod)
1709//     
1710//      return val
1711//
1712//End
1713//
1714////calculates the x distance from the center of the detector, w/nonlinear corrections
1715////
1716//Function V_FX(xx,sx3,xcenter,sx)             
1717//      Variable xx,sx3,xcenter,sx
1718//     
1719//      Variable retval
1720//     
1721//      retval = sx3*tan((xx-xcenter)*sx/sx3)
1722//      Return retval
1723//End
1724//
1725////calculates the y distance from the center of the detector, w/nonlinear corrections
1726////
1727//Function V_FY(yy,sy3,ycenter,sy)             
1728//      Variable yy,sy3,ycenter,sy
1729//     
1730//      Variable retval
1731//     
1732//      retval = sy3*tan((yy-ycenter)*sy/sy3)
1733//      Return retval
1734//End
Note: See TracBrowser for help on using the repository browser.