source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/VC_DetectorBinning_Utils.ipf @ 1098

Last change on this file since 1098 was 1098, checked in by srkline, 5 years ago

renamed white beam smearing models so they would be grouped together on the file list.

Re-worked the logic and flow of the averaging/plotting/saving steps of the reduction protocol so that it would flow cleanly and leave room for changes for the multitude of different collimation conditions. The averaging routines are now aware of the collimation conditions so that the appropriate resolution can be calculated. The collimation string is also written out to the averaged data file as element[9] of the protocol. The hope is that one could key on this collimation string to decide how to proceed with the analysis.

File size: 50.8 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2
3/////////////////////////
4//
5// Utility functions to:
6//              calculate Q, Qx, Qy, Qz
7//              fill the detector panels with simulated data (the model functions are here)
8//              bin the 2D detector to 1D I(Q) based on Q and deltaQ (bin width)
9//
10/////////////////////////
11
12
13
14
15// x- hard wired for a sphere - change this to allow minimal selections and altering of coefficients
16// x- add the "fake" 2D simulation to fill the panels which are then later averaged as I(Q)
17//
18// NOTE - this is a VCALC only routine, so it's not been made completely generic
19//
20Function FillPanel_wModelData(det,qTot,type)
21        Wave det,qTot
22        String type
23
24//      SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
25
26        // q-values and detector arrays already allocated and calculated
27        Duplicate/O det tmpInten,tmpSig,prob_i
28               
29        Variable imon,trans,thick,sdd,pixSizeX,pixSizeY,sdd_offset
30
31        //imon = VC_BeamIntensity()*CountTime
32        imon = VCALC_getImon()          //TODO: currently from the panel, not calculated
33        trans = 0.8
34        thick = 0.1
35       
36        // need SDD
37        // need pixel dimensions
38        // nominal sdd in cm, offset in cm, want result in cm !
39        sdd = VCALC_getSDD(type)        +  VCALC_getTopBottomSDDSetback(type)           // result is sdd in [cm]
40
41        pixSizeX = VCALC_getPixSizeX(type)              // cm
42        pixSizeY = VCALC_getPixSizeY(type)
43       
44       
45        //?? pick the function from a popup on the panel? (bypass the analysis panel, or maybe it's better to
46        //  keep the panel to keep people used to using it.)
47        // peak @ 0.1 ~ AgBeh
48        //      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}               
49        //
50        // peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
51        //      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}     
52        //
53        //
54        Variable addEmpBgd=0
55       
56               
57        String funcStr = VCALC_getModelFunctionStr()
58        strswitch(funcStr)
59                case "Big Debye":
60                        tmpInten = VC_Debye(100,3000,0.0001,qTot[p][q])
61                        break
62                case "Big Sphere":
63                        tmpInten = VC_SphereForm(1,900,1e-6,0.01,qTot[p][q])   
64                        break
65                case "Debye":
66                        tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
67                        break
68                case "Sphere":
69                        tmpInten = VC_SphereForm(1,60,1e-6,0.001,qTot[p][q])   
70                        break
71                case "AgBeh":
72                        tmpInten = VC_BroadPeak(1e-11,3,20,100.0,0.1,3,0.1,qTot[p][q])
73                        break
74                case "Vycor":
75                        tmpInten = VC_BroadPeak(1e-9,3,20,500.0,0.015,3,0.1,qTot[p][q])
76                        break   
77                case "Empty Cell":
78                        tmpInten = VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
79                        break
80                case "Blocked Beam":
81                        tmpInten = VC_BlockedBeam(0.01,qTot[p][q])
82                        break
83                case "Debye +":
84                        tmpInten = VC_Debye(10,300,0.0001,qTot[p][q])
85                        addEmpBgd = 1
86                        break
87                case "AgBeh +":
88                        tmpInten = VC_BroadPeak(1e-11,3,20,100.0,0.1,3,0.1,qTot[p][q])
89                        addEmpBgd = 1
90                        break
91                case "Empty Cell +":
92                        tmpInten = VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
93                        tmpInten += VC_BlockedBeam(0.01,qTot[p][q])
94                        break
95                default:
96                        tmpInten = VC_Debye(10,300,0.1,qTot[p][q])
97        endswitch
98
99
100        if(addEmpBgd == 1)
101                tmpInten += VC_EC_Empirical(2.2e-12,3.346,0.0065,9.0,0.016,qTot[p][q])
102                tmpInten += VC_BlockedBeam(0.01,qTot[p][q])
103        endif
104
105       
106// x- this is faked to get around the singularity at the center of the back detector
107//
108//
109        if(cmpstr(type,"B") == 0)
110                Variable nx,ny,px,py
111                nx = VCALC_get_nPix_X(type)
112                ny = VCALC_get_nPix_Y(type)
113                px = trunc(nx/2)
114                py = trunc(ny/2)
115               
116                tmpInten[px][py] = (tmpInten[px][py+1] + tmpInten[px][py-1])/2
117        endif
118
119
120
121///////////////
122//      // calculate the scattering cross section simply to be able to estimate the transmission
123//      Variable sig_sas=0
124//     
125//      // remember that the random deviate is the coherent portion ONLY - the incoherent background is
126//      // subtracted before the calculation.
127//      CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
128//
129//      if(sig_sas > 100)
130//              DoAlert 0,"SAS cross section > 100. Estimates of multiple scattering are unreliable. Choosing a model with a well-defined Rg may help"
131//      endif           
132//
133//      // calculate the multiple scattering fraction for display (10/2009)
134//      Variable ii,nMax=10,tau
135//      mScat=0
136//      tau = thick*sig_sas
137//      // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
138//      // neutrons out of those that were scattered
139//      for(ii=2;ii<nMax;ii+=1)
140//              mScat += tau^(ii)/factorial(ii)
141////            print tau^(ii)/factorial(ii)
142//      endfor
143//      estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
144//      mscat *= (estTrans)/(1-estTrans)
145//
146////    if(mScat > 0.1)         //  Display warning
147//
148//      Print "Sig_sas = ",sig_sas
149////////////////////
150       
151        prob_i = trans*thick*pixSizeX*pixSizeY/(sdd)^2*tmpInten                 //probability of a neutron in q-bin(i)
152               
153        tmpInten = (imon)*prob_i                //tmpInten is not the model calculation anymore!!
154
155
156/// **** can I safely assume a Gaussian error in the count rate??
157        tmpSig = sqrt(tmpInten)         // corrected based on John's memo, from 8/9/99
158
159        tmpInten += gnoise(tmpSig)
160        tmpInten = (tmpInten[p][q] < 0) ? 0 : tmpInten[p][q]                    // MAR 2013 -- is this the right thing to do
161        tmpInten = trunc(tmpInten)
162               
163       
164        det = tmpInten
165
166// if I want "absolute" scale -- then I lose the integer nature of the detector (but keep the random)
167//      det /= trans*thick*pixSizeX*pixSizeY/(sdd)^2*imon
168
169       
170        KillWaves/Z tmpInten,tmpSig,prob_i     
171        SetDataFolder root:
172
173        return(0)
174End
175
176
177// For a given detector panel, calculate the q-values
178// -work with everything as arrays
179// Input needed:
180// detector data
181// detector type (LRTB?)
182// beam center (may be off the detector)
183// SDD
184// lambda
185//
186// pixel dimensions for detector type (global constants)
187// - data dimensions read directly from array
188//
189// --What is calculated:
190// array of Q
191// array of qx,qy,qz
192// array of error already exists
193//
194//
195// -- sdd in cm
196// -- lambda in Angstroms
197Function VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
198        Wave data,qTot,qx,qy,qz
199        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
200               
201        // loop over the array and calculate the values - this is done as a wave assignment
202// TODO -- be sure that it's p,q -- or maybe p+1,q+1 as used in WriteQIS.ipf   
203        qTot = VC_CalcQval(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
204        qx = VC_CalcQX(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
205        qy = VC_CalcQY(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
206        qz = VC_CalcQZ(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
207       
208        return(0)
209End
210
211
212// for testing, a version that will calculate the q-arrays for VCALC based on whatever nonlinear coefficients
213// exist in the RAW data folder
214//
215// reverts to the "regular" linear detector if waves not found or a flag is set
216//
217// need to call the VSANS V_CalcQval routines (these use the real-space distance, not pixel dims)
218//
219// ***** everything passed in is [cm], except for wavelength [A]
220//
221// ****  TODO :: calibration constants are still in [mm]
222//
223//
224// TODO:
225// -- tube width is hard-wired in
226//
227//
228Function VC_Detector_2Q_NonLin(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY,detStr)
229        Wave data,qTot,qx,qy,qz
230        Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
231        String detStr
232       
233        String destPath = "root:Packages:NIST:VSANS:VCALC"
234       
235        // be sure that the real distance waves exist
236        // TODO -- this may not be the best location?
237
238// calibration waves do not exist yet, so make some fake ones   '
239        // do I count on the orientation as an input, or do I just figure it out on my own?
240        String orientation
241        Variable dimX,dimY
242        dimX = DimSize(data,0)
243        dimY = DimSize(data,1)
244        if(dimX > dimY)
245                orientation = "horizontal"
246        else
247                orientation = "vertical"
248        endif
249       
250        if(cmpstr(orientation,"vertical")==0)
251                Make/O/D/N=(3,48) tmpCalib
252                // for the "tall" L/R banks
253                tmpCalib[0][] = -512
254                tmpCalib[1][] = 8
255                tmpCalib[2][] = 0
256        else
257                Make/O/D/N=(3,48) tmpCalib
258                // for the "short" T/B banks
259                tmpCalib[0][] = -256
260                tmpCalib[1][] = 4
261                tmpCalib[2][] = 0
262        endif
263        // override if back panel
264        if(cmpstr(detStr,"B") == 0)
265                // and for the back detector "B"
266                Make/O/D/N=3 tmpCalibX,tmpCalibY
267                tmpCalibX[0] = VCALC_getPixSizeX(detStr)*10                     // pixel size in mm  VCALC_getPixSizeX(detStr) is [cm]
268                tmpCalibX[1] = 1
269                tmpcalibX[2] = 10000
270                tmpCalibY[0] = VCALC_getPixSizeY(detStr)*10                     // pixel size in mm  VCALC_getPixSizeX(detStr) is [cm]
271                tmpCalibY[1] = 1
272                tmpcalibY[2] = 10000
273        endif
274       
275//      Wave w_calib = V_getDetTube_spatialCalib("VCALC",detStr)
276        Variable tube_width = 8.4                       // TODO: UNITS!!! Hard-wired value in [mm]
277        if(cmpstr(detStr,"B") == 0)
278                V_NonLinearCorrection_B("VCALC",data,tmpCalibX,tmpCalibY,detStr,destPath)
279        else
280                V_NonLinearCorrection("VCALC",data,tmpCalib,tube_width,detStr,destPath)
281        endif
282                               
283        Wave/Z data_realDistX = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistX")
284        Wave/Z data_realDistY = $(destPath + ":entry:instrument:detector_"+detStr+":data_realDistY")
285        NVAR gUseNonLinearDet = root:Packages:NIST:VSANS:VCALC:gUseNonLinearDet
286
287        if(kBCTR_CM)
288                if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
289                        // beam ctr is in cm already
290
291                        // calculate all of the q-values
292                        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
293                        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
294                        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
295                        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
296               
297        //              Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
298        //              Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
299                else
300                        // do the q-calculation using linear detector
301                        //VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
302                        qTot = V_CalcQval(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
303                        qx = V_CalcQX(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
304                        qy = V_CalcQY(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
305                        qz = V_CalcQZ(p,q,xCtr,yCtr,sdd,lam,data_realDistX,data_realDistY)
306                endif   
307       
308       
309        else
310        // using the old calculation with beam center in pixels
311                if(gUseNonLinearDet && WaveExists(data_realDistX) && WaveExists(data_realDistY))
312                        // convert the beam centers to mm
313//                      String orientation
314                        Variable newX,newY
315                        dimX = DimSize(data_realDistX,0)
316                        dimY = DimSize(data_realDistX,1)
317                        if(dimX > dimY)
318                                orientation = "horizontal"
319                        else
320                                orientation = "vertical"
321                        endif
322                       
323               
324                //
325                        if(cmpstr(orientation,"vertical")==0)
326                                //      this is data dimensioned as (Ntubes,Npix)
327                                newX = tube_width*xCtr
328                                newY = data_realDistY[0][yCtr]
329                        else
330                                //      this is data (horizontal) dimensioned as (Npix,Ntubes)
331                                newX = data_realDistX[xCtr][0]
332                                newY = tube_width*yCtr
333                        endif   
334       
335                        //if detector "B", different calculation for the centers (not tubes)
336                        if(cmpstr(detStr,"B")==0)
337                                newX = data_realDistX[xCtr][0]
338                                newY = data_realDistY[0][yCtr]
339                                //newX = xCtr
340                                //newY = yCtr
341                        endif           
342                                       
343                        // calculate all of the q-values
344                        qTot = V_CalcQval(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
345                        qx = V_CalcQX(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
346                        qy = V_CalcQY(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
347                        qz = V_CalcQZ(p,q,newX,newY,sdd,lam,data_realDistX,data_realDistY)
348               
349        //              Print "det, x_mm, y_mm ",detStr,num2str(newX),num2str(newY)
350        //              Print "det, x_pix, y_pix ",detStr,num2str(xCtr),num2str(yCtr)
351                else
352                        // do the q-calculation using linear detector
353                        VC_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
354                endif   
355       
356        endif
357       
358        KillWaves/Z tmpCalib,tmpCalibX,tmpCalibY
359       
360        return(0)
361End
362
363
364//////////////////////
365// NOTE: The Q calculations are different than what is in GaussUtils in that they take into
366// accout the different x/y pixel sizes and the beam center not being on the detector -
367// off a different edge for each LRTB type
368/////////////////////
369
370//function to calculate the overall q-value, given all of the necesary trig inputs
371//and are in detector coordinates (1,128) rather than axis values
372//the pixel locations need not be integers, reals are ok inputs
373//sdd is in [cm]
374//wavelength is in Angstroms
375//
376//returned magnitude of Q is in 1/Angstroms
377//
378//
379Function VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
380        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
381       
382        Variable dx,dy,qval,two_theta,dist
383               
384
385        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
386        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
387        dist = sqrt(dx^2 + dy^2)
388       
389        two_theta = atan(dist/sdd)
390
391        qval = 4*Pi/lam*sin(two_theta/2)
392       
393        return qval
394End
395
396//calculates just the q-value in the x-direction on the detector
397//input/output is the same as CalcQval()
398//ALL inputs are in detector coordinates
399//
400//sdd is in [cm]
401//wavelength is in Angstroms
402//
403// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
404// now properly accounts for qz
405//
406Function VC_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
407        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
408
409        Variable qx,qval,phi,dx,dy,dist,two_theta
410       
411        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
412       
413//      sdd *=100               //convert to cm
414        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
415        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
416        phi = V_FindPhi(dx,dy)
417       
418        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
419        dist = sqrt(dx^2 + dy^2)
420        two_theta = atan(dist/sdd)
421
422        qx = qval*cos(two_theta/2)*cos(phi)
423       
424        return qx
425End
426
427//calculates just the q-value in the y-direction on the detector
428//input/output is the same as CalcQval()
429//ALL inputs are in detector coordinates
430//sdd is in [cm]
431//wavelength is in Angstroms
432//
433// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
434// now properly accounts for qz
435//
436Function VC_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
437        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
438       
439        Variable dy,qval,dx,phi,qy,dist,two_theta
440       
441        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
442       
443//      sdd *=100               //convert to cm
444        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
445        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
446        phi = V_FindPhi(dx,dy)
447       
448        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
449        dist = sqrt(dx^2 + dy^2)
450        two_theta = atan(dist/sdd)
451       
452        qy = qval*cos(two_theta/2)*sin(phi)
453       
454        return qy
455End
456
457//calculates just the z-component of the q-vector, not measured on the detector
458//input/output is the same as CalcQval()
459//ALL inputs are in detector coordinates
460//sdd is in [cm]
461//wavelength is in Angstroms
462//
463// not actually used, but here for completeness if anyone asks
464//
465Function VC_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
466        Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
467       
468        Variable dy,qval,dx,phi,qz,dist,two_theta
469       
470        qval = VC_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
471       
472//      sdd *=100               //convert to cm
473       
474        //get scattering angle to project onto flat detector => Qr = qval*cos(theta)
475        dx = (xaxval - xctr)*pixSizeX           //delta x in cm
476        dy = (yaxval - yctr)*pixSizeY           //delta y in cm
477        dist = sqrt(dx^2 + dy^2)
478        two_theta = atan(dist/sdd)
479       
480        qz = qval*sin(two_theta/2)
481       
482        return qz
483End
484
485//phi is defined from +x axis, proceeding CCW around [0,2Pi]
486Threadsafe Function V_FindPhi(vx,vy)
487        variable vx,vy
488       
489        variable phi
490       
491        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
492       
493        // special cases
494        if(vx==0 && vy > 0)
495                return(pi/2)
496        endif
497        if(vx==0 && vy < 0)
498                return(3*pi/2)
499        endif
500        if(vx >= 0 && vy == 0)
501                return(0)
502        endif
503        if(vx < 0 && vy == 0)
504                return(pi)
505        endif
506       
507       
508        if(vx > 0 && vy > 0)
509                return(phi)
510        endif
511        if(vx < 0 && vy > 0)
512                return(phi + pi)
513        endif
514        if(vx < 0 && vy < 0)
515                return(phi + pi)
516        endif
517        if( vx > 0 && vy < 0)
518                return(phi + 2*pi)
519        endif
520       
521        return(phi)
522end
523
524Function VC_SphereForm(scale,radius,delrho,bkg,x)                               
525        Variable scale,radius,delrho,bkg
526        Variable x
527       
528        // variables are:                                                       
529        //[0] scale
530        //[1] radius (A)
531        //[2] delrho (A-2)
532        //[3] background (cm-1)
533       
534//      Variable scale,radius,delrho,bkg                               
535//      scale = w[0]
536//      radius = w[1]
537//      delrho = w[2]
538//      bkg = w[3]
539       
540       
541        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
542        // and is rescaled to give [=] cm^-1
543       
544        Variable bes,f,vol,f2
545        //
546        //handle q==0 separately
547        If(x==0)
548                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
549                return(f)
550        Endif
551       
552//      bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
553       
554        bes = 3*sqrt(pi/(2*x*radius))*BesselJ(1.5,x*radius)/(x*radius)
555       
556        vol = 4*pi/3*radius^3
557        f = vol*bes*delrho              // [=] A
558        // normalize to single particle volume, convert to 1/cm
559        f2 = f * f / vol * 1.0e8                // [=] 1/cm
560       
561        return (scale*f2+bkg)   // Scale, then add in the background
562       
563End
564
565Function VC_Debye(scale,rg,bkg,x)
566        Variable scale,rg,bkg
567        Variable x
568       
569        // variables are:
570        //[0] scale factor
571        //[1] radius of gyration [A]
572        //[2] background        [cm-1]
573       
574        // calculates (scale*debye)+bkg
575        Variable Pq,qr2
576       
577        qr2=(x*rg)^2
578        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
579       
580        //scale
581        Pq *= scale
582        // then add in the background
583        return (Pq+bkg)
584End
585
586// a sum of a power law and debye to approximate the scattering from a real empty cell
587//
588//      make/O/D coef_ECEmp = {2.2e-8,3.346,0.0065,9.0,0.016}
589//
590Function VC_EC_Empirical(aa,mm,scale,rg,bkg,x)
591        Variable aa,mm,scale,rg,bkg
592        Variable x
593       
594        // variables are:
595        //[0] = A
596        //[1] = power m
597        //[2] scale factor
598        //[3] radius of gyration [A]
599        //[4] background        [cm-1]
600       
601        Variable Iq
602       
603        // calculates (scale*debye)+bkg
604        Variable Pq,qr2
605       
606//      if(x*Rg < 1e-3)         //added Oct 2008 to avoid numerical errors at low arg values
607//              return(scale+bkg)
608//      endif
609       
610        Iq = aa*x^-mm
611       
612        qr2=(x*rg)^2
613        Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
614       
615        //scale
616        Pq *= scale
617        // then add the terms up
618        return (Iq + Pq + bkg)
619End
620
621// blocked beam
622//
623Function VC_BlockedBeam(bkg,x)
624        Variable bkg
625        Variable x
626       
627        return (bkg)
628End
629
630
631//
632// a broad peak to simulate silver behenate or vycor
633//
634// peak @ 0.1 ~ AgBeh
635//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 100.0, 0.1,3,0.1}               
636//
637//
638// peak @ 0.015 in middle of middle detector, maybe not "real" vycor, but that is to be resolved
639//      Make/O/D coef_BroadPeak = {1e-9, 3, 20, 500.0, 0.015,3,0.1}             
640//
641//
642Function VC_BroadPeak(aa,nn,cc,LL,Qzero,mm,bgd,x)
643        Variable aa,nn,cc,LL,Qzero,mm,bgd
644        Variable x
645       
646        // variables are:                                                       
647        //[0] Porod term scaling
648        //[1] Porod exponent
649        //[2] Lorentzian term scaling
650        //[3] Lorentzian screening length [A]
651        //[4] peak location [1/A]
652        //[5] Lorentzian exponent
653        //[6] background
654       
655//      local variables
656        Variable inten, qval
657//      x is the q-value for the calculation
658        qval = x
659//      do the calculation and return the function value
660       
661        inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd
662
663        Return (inten)
664       
665End
666
667//
668// updated to new folder structure Feb 2016
669// folderStr = RAW,SAM, VCALC or other
670// detStr is the panel identifer "ML", etc.
671//
672Function V_SetDeltaQ(folderStr,detStr)
673        String folderStr,detStr
674
675        Variable isVCALC
676        if(cmpstr(folderStr,"VCALC") == 0)
677                isVCALC = 1
678        endif
679       
680        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
681        String instPath = ":entry:instrument:detector_"
682               
683        if(isVCALC)
684                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)               // 2D detector data
685        else
686                Wave inten = V_getDetectorDataW(folderStr,detStr)
687        endif
688
689        Wave qx = $(folderPath+instPath+detStr+":qx_"+detStr)
690        Wave qy = $(folderPath+instPath+detStr+":qy_"+detStr)
691       
692        Variable xDim,yDim,delQ
693       
694        xDim=DimSize(inten,0)
695        yDim=DimSize(inten,1)
696       
697        if(xDim<yDim)
698                delQ = abs(qx[0][0] - qx[1][0])/2
699        else
700                delQ = abs(qy[0][1] - qy[0][0])/2
701        endif
702       
703        // set the global
704        Variable/G $(folderPath+instPath+detStr+":gDelQ_"+detStr) = delQ
705//      Print "SET delQ = ",delQ," for ",type
706       
707        return(delQ)
708end
709
710
711//TODO -- need a switch here to dispatch to the averaging type
712Proc VC_BinQxQy_to_1D(folderStr,type)
713        String folderStr
714        String type
715//      Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
716//      Prompt type,"detector identifier"
717
718
719        VC_fDoBinning_QxQy2D(folderStr, type)
720
721
722/// this is for a tall, narrow slit mode       
723//      VC_fBinDetector_byRows(folderStr,type)
724       
725End
726
727
728// folderStr is RAW, VCALC, SAM, etc.
729// type is "B", "FL" for single binning, "FLR", or "MLRTB" or similar if multiple panels are combined
730//
731Proc VC_Graph_1D_detType(folderStr,type)
732        String folderStr,type
733       
734        SetDataFolder $("root:Packages:NIST:VSANS:"+folderStr)
735       
736        Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
737        ModifyGraph mirror=2,grid=1,log=1
738        ModifyGraph mode=4,marker=19,msize=2
739//      ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy)                // for simulations, I don't have 2D uncertainty
740        ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
741        legend
742       
743        SetDataFolder root:
744
745End
746
747
748
749//////////
750//
751//              Function that bins a 2D detctor panel into I(q) based on the q-value of the pixel
752//              - each pixel QxQyQz has been calculated beforehand
753//              - if multiple panels are selected to be combined, it is done here during the binning
754//              - the setting of deltaQ step is still a little suspect (TODO)
755//
756//
757// see the equivalent function in PlotUtils2D_v40.ipf
758//
759//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
760//
761// this has been modified to accept different detector panels and to take arrays
762// -- type = FL or FR or...other panel identifiers
763//
764// TODO "iErr" is not always defined correctly since it doesn't really apply here for data that is not 2D simulation
765//
766//
767// updated Feb2016 to take new folder structure
768// TODO
769// -- VERIFY
770// -- figure out what the best location is to put the averaged data? currently @ top level of WORK folder
771//    but this is a lousy choice.
772// x- binning is now Mask-aware. If mask is not present, all data is used. If data is from VCALC, all data is used
773// x- Where do I put the solid angle correction? In here as a weight for each point, or later on as
774//    a blanket correction (matrix multiply) for an entire panel? (Solid Angle correction is done in the
775//    step where data is added to a WORK file (see Raw_to_Work())
776//
777//
778// TODO:
779// -- some of the input parameters for the resolution calcuation are either assumed (apOff) or are currently
780//    hard-wired. these need to be corrected before even the pinhole resolution is correct
781// x- resolution calculation is in the correct place. The calculation is done per-panel (specified by TYPE),
782//    and then the unwanted points can be discarded (all 6 columns) as the data is trimmed and concatenated
783//    is separate functions that are resolution-aware.
784//
785//
786// folderStr = WORK folder, type = the binning type (may include multiple detectors)
787Function VC_fDoBinning_QxQy2D(folderStr,type,collimationStr)
788        String folderStr,type,collimationStr
789       
790        Variable nSets = 0
791        Variable xDim,yDim
792        Variable ii,jj
793        Variable qVal,nq,var,avesq,aveisq
794        Variable binIndex,val,isVCALC=0,maskMissing
795
796        String folderPath = "root:Packages:NIST:VSANS:"+folderStr
797        String instPath = ":entry:instrument:detector_"
798        String detStr
799               
800        if(cmpstr(folderStr,"VCALC") == 0)
801                isVCALC = 1
802        endif
803       
804        detStr = type
805
806// now switch on the type to determine which waves to declare and create
807// since there may be more than one panel to step through. There may be two, there may be four
808//
809
810// TODO:
811// -- Solid_Angle -- waves will be present for WORK data other than RAW, but not for RAW
812//
813// assume that the mask files are missing unless we can find them. If VCALC data,
814//  then the Mask is missing by definition
815        maskMissing = 1
816
817        strswitch(type) // string switch
818//              case "FL":              // execute if case matches expression
819//              case "FR":
820//                      detStr = type
821//                      if(isVCALC)
822//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
823//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation
824//                      else
825//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
826//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
827//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
828//                              if(WaveExists(mask) == 1)
829//                                      maskMissing = 0
830//                              endif
831//                             
832//                      endif
833//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
834//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
835//                      nSets = 1
836//                      break   
837                                                               
838//              case "FT":             
839//              case "FB":
840//                      detStr = type
841//                      if(isVCALC)
842//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
843//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
844//                      else
845//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
846//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
847//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
848//                              if(WaveExists(mask) == 1)
849//                                      maskMissing = 0
850//                              endif
851//                      endif
852//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
853//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
854//                      nSets = 1
855//                      break
856                       
857//              case "ML":             
858//              case "MR":
859//                      detStr = type
860//                      if(isVCALC)
861//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
862//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
863//                      else
864//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
865//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
866//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
867//                              if(WaveExists(mask) == 1)
868//                                      maskMissing = 0
869//                              endif
870//                      endif   
871//                      //TODO:
872//                      // -- decide on the proper deltaQ for binning. either nominal value for LR, or one
873//                      //    determined specifically for that panel (currently using one tube width as deltaQ)
874//                      // -- this is repeated multiple times in this switch
875//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
876//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
877//                      nSets = 1
878//                      break   
879                                       
880//              case "MT":             
881//              case "MB":
882//                      detStr = type
883//                      if(isVCALC)
884//                              WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
885//                              WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
886//                      else
887//                              Wave inten = V_getDetectorDataW(folderStr,detStr)
888//                              Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
889//                              Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
890//                              if(WaveExists(mask) == 1)
891//                                      maskMissing = 0
892//                              endif
893//                      endif   
894//                      NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
895//                      Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values
896//                      nSets = 1
897//                      break   
898
899// only one panel, simply pick that panel and move on out of the switch
900                case "FL":
901                case "FR":
902                case "FT":
903                case "FB":
904                case "ML":
905                case "MR":
906                case "MT":
907                case "MB":                     
908                case "B":       
909                        if(isVCALC)
910                                WAVE inten = $(folderPath+instPath+detStr+":det_"+detStr)
911                                WAVE/Z iErr = $("iErr_"+detStr)                 // 2D errors -- may not exist, especially for simulation               
912                        else
913                                Wave inten = V_getDetectorDataW(folderStr,detStr)
914                                Wave iErr = V_getDetectorDataErrW(folderStr,detStr)
915                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+detStr+":data")
916                                if(WaveExists(mask) == 1)
917                                        maskMissing = 0
918                                endif
919                        endif   
920                        NVAR delQ = $(folderPath+instPath+detStr+":gDelQ_"+detStr)
921                        Wave qTotal = $(folderPath+instPath+detStr+":qTot_"+detStr)                     // 2D q-values 
922                        nSets = 1
923                        break   
924                       
925                case "FLR":
926                // detStr has multiple values now, so unfortuntely, I'm hard-wiring things...
927                // TODO
928                // -- see if I can un-hard-wire some of this below when more than one panel is combined
929                        if(isVCALC)
930                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
931                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
932                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
933                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
934                        else
935                                Wave inten = V_getDetectorDataW(folderStr,"FL")
936                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
937                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
938                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
939                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
940                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
941                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
942                                        maskMissing = 0
943                                endif
944                        endif   
945                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
946                       
947                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
948                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
949               
950                        nSets = 2
951                        break                   
952               
953                case "FTB":
954                        if(isVCALC)
955                                WAVE inten = $(folderPath+instPath+"FT"+":det_"+"FT")
956                                WAVE/Z iErr = $("iErr_"+"FT")                   // 2D errors -- may not exist, especially for simulation               
957                                WAVE inten2 = $(folderPath+instPath+"FB"+":det_"+"FB")
958                                WAVE/Z iErr2 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
959                        else
960                                Wave inten = V_getDetectorDataW(folderStr,"FT")
961                                Wave iErr = V_getDetectorDataErrW(folderStr,"FT")
962                                Wave inten2 = V_getDetectorDataW(folderStr,"FB")
963                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FB")
964                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
965                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
966                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
967                                        maskMissing = 0
968                                endif
969                        endif   
970                        NVAR delQ = $(folderPath+instPath+"FT"+":gDelQ_FT")
971                       
972                        Wave qTotal = $(folderPath+instPath+"FT"+":qTot_"+"FT")                 // 2D q-values 
973                        Wave qTotal2 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
974       
975                        nSets = 2
976                        break           
977               
978                case "FLRTB":
979                        if(isVCALC)
980                                WAVE inten = $(folderPath+instPath+"FL"+":det_"+"FL")
981                                WAVE/Z iErr = $("iErr_"+"FL")                   // 2D errors -- may not exist, especially for simulation               
982                                WAVE inten2 = $(folderPath+instPath+"FR"+":det_"+"FR")
983                                WAVE/Z iErr2 = $("iErr_"+"FR")                  // 2D errors -- may not exist, especially for simulation       
984                                WAVE inten3 = $(folderPath+instPath+"FT"+":det_"+"FT")
985                                WAVE/Z iErr3 = $("iErr_"+"FT")                  // 2D errors -- may not exist, especially for simulation               
986                                WAVE inten4 = $(folderPath+instPath+"FB"+":det_"+"FB")
987                                WAVE/Z iErr4 = $("iErr_"+"FB")                  // 2D errors -- may not exist, especially for simulation       
988                        else
989                                Wave inten = V_getDetectorDataW(folderStr,"FL")
990                                Wave iErr = V_getDetectorDataErrW(folderStr,"FL")
991                                Wave inten2 = V_getDetectorDataW(folderStr,"FR")
992                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"FR")
993                                Wave inten3 = V_getDetectorDataW(folderStr,"FT")
994                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"FT")
995                                Wave inten4 = V_getDetectorDataW(folderStr,"FB")
996                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"FB")
997                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FL"+":data")
998                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FR"+":data")
999                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FT"+":data")
1000                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"FB"+":data")
1001                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1002                                        maskMissing = 0
1003                                endif
1004                        endif   
1005                        NVAR delQ = $(folderPath+instPath+"FL"+":gDelQ_FL")
1006                       
1007                        Wave qTotal = $(folderPath+instPath+"FL"+":qTot_"+"FL")                 // 2D q-values 
1008                        Wave qTotal2 = $(folderPath+instPath+"FR"+":qTot_"+"FR")                        // 2D q-values 
1009                        Wave qTotal3 = $(folderPath+instPath+"FT"+":qTot_"+"FT")                        // 2D q-values 
1010                        Wave qTotal4 = $(folderPath+instPath+"FB"+":qTot_"+"FB")                        // 2D q-values 
1011               
1012                        nSets = 4
1013                        break           
1014                       
1015                case "MLR":
1016                        if(isVCALC)
1017                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1018                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
1019                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1020                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
1021                        else
1022                                Wave inten = V_getDetectorDataW(folderStr,"ML")
1023                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1024                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1025                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1026                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1027                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1028                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1029                                        maskMissing = 0
1030                                endif
1031                        endif   
1032                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1033                       
1034                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
1035                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
1036               
1037                        nSets = 2
1038                        break                   
1039               
1040                case "MTB":
1041                        if(isVCALC)
1042                                WAVE inten = $(folderPath+instPath+"MT"+":det_"+"MT")
1043                                WAVE/Z iErr = $("iErr_"+"MT")                   // 2D errors -- may not exist, especially for simulation               
1044                                WAVE inten2 = $(folderPath+instPath+"MB"+":det_"+"MB")
1045                                WAVE/Z iErr2 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
1046                        else
1047                                Wave inten = V_getDetectorDataW(folderStr,"MT")
1048                                Wave iErr = V_getDetectorDataErrW(folderStr,"MT")
1049                                Wave inten2 = V_getDetectorDataW(folderStr,"MB")
1050                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MB")
1051                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1052                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1053                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1)
1054                                        maskMissing = 0
1055                                endif
1056                        endif   
1057                        NVAR delQ = $(folderPath+instPath+"MT"+":gDelQ_MT")
1058                       
1059                        Wave qTotal = $(folderPath+instPath+"MT"+":qTot_"+"MT")                 // 2D q-values 
1060                        Wave qTotal2 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
1061               
1062                        nSets = 2
1063                        break                           
1064               
1065                case "MLRTB":
1066                        if(isVCALC)
1067                                WAVE inten = $(folderPath+instPath+"ML"+":det_"+"ML")
1068                                WAVE/Z iErr = $("iErr_"+"ML")                   // 2D errors -- may not exist, especially for simulation               
1069                                WAVE inten2 = $(folderPath+instPath+"MR"+":det_"+"MR")
1070                                WAVE/Z iErr2 = $("iErr_"+"MR")                  // 2D errors -- may not exist, especially for simulation       
1071                                WAVE inten3 = $(folderPath+instPath+"MT"+":det_"+"MT")
1072                                WAVE/Z iErr3 = $("iErr_"+"MT")                  // 2D errors -- may not exist, especially for simulation               
1073                                WAVE inten4 = $(folderPath+instPath+"MB"+":det_"+"MB")
1074                                WAVE/Z iErr4 = $("iErr_"+"MB")                  // 2D errors -- may not exist, especially for simulation       
1075                        else
1076                                Wave inten = V_getDetectorDataW(folderStr,"ML")
1077                                Wave iErr = V_getDetectorDataErrW(folderStr,"ML")
1078                                Wave inten2 = V_getDetectorDataW(folderStr,"MR")
1079                                Wave iErr2 = V_getDetectorDataErrW(folderStr,"MR")
1080                                Wave inten3 = V_getDetectorDataW(folderStr,"MT")
1081                                Wave iErr3 = V_getDetectorDataErrW(folderStr,"MT")
1082                                Wave inten4 = V_getDetectorDataW(folderStr,"MB")
1083                                Wave iErr4 = V_getDetectorDataErrW(folderStr,"MB")
1084                                Wave/Z mask = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"ML"+":data")
1085                                Wave/Z mask2 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MR"+":data")
1086                                Wave/Z mask3 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MT"+":data")
1087                                Wave/Z mask4 = $("root:Packages:NIST:VSANS:MSK:entry:instrument:detector_"+"MB"+":data")
1088                                if(WaveExists(mask) == 1 && WaveExists(mask2) == 1 && WaveExists(mask3) == 1 && WaveExists(mask4) == 1)
1089                                        maskMissing = 0
1090                                endif
1091                        endif   
1092                        NVAR delQ = $(folderPath+instPath+"ML"+":gDelQ_ML")
1093                       
1094                        Wave qTotal = $(folderPath+instPath+"ML"+":qTot_"+"ML")                 // 2D q-values 
1095                        Wave qTotal2 = $(folderPath+instPath+"MR"+":qTot_"+"MR")                        // 2D q-values 
1096                        Wave qTotal3 = $(folderPath+instPath+"MT"+":qTot_"+"MT")                        // 2D q-values 
1097                        Wave qTotal4 = $(folderPath+instPath+"MB"+":qTot_"+"MB")                        // 2D q-values 
1098               
1099                        nSets = 4
1100                        break                                                                   
1101                                       
1102                default:
1103                        nSets = 0                                                       
1104                        Print "ERROR   ---- type is not recognized "
1105        endswitch
1106
1107//      Print "delQ = ",delQ," for ",type
1108
1109        if(nSets == 0)
1110                SetDataFolder root:
1111                return(0)
1112        endif
1113
1114
1115// RAW data is currently read in and the 2D error wave is correctly generated
1116// 2D error is propagated through all reduction steps, but I have not
1117// verified that it is an exact duplication of the 1D error
1118//
1119//
1120//
1121// IF ther is no 2D error wave present for some reason, make a fake one
1122        if(WaveExists(iErr)==0  && WaveExists(inten) != 0)
1123                Duplicate/O inten,iErr
1124                Wave iErr=iErr
1125//              iErr = 1+sqrt(inten+0.75)                       // can't use this -- it applies to counts, not intensity (already a count rate...)
1126                iErr = sqrt(inten+0.75)                 // TODO -- here I'm just using some fictional value
1127        endif
1128        if(WaveExists(iErr2)==0 && WaveExists(inten2) != 0)
1129                Duplicate/O inten2,iErr2
1130                Wave iErr2=iErr2
1131//              iErr2 = 1+sqrt(inten2+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1132                iErr2 = sqrt(inten2+0.75)                       // TODO -- here I'm just using some fictional value
1133        endif
1134        if(WaveExists(iErr3)==0  && WaveExists(inten3) != 0)
1135                Duplicate/O inten3,iErr3
1136                Wave iErr3=iErr3
1137//              iErr3 = 1+sqrt(inten3+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1138                iErr3 = sqrt(inten3+0.75)                       // TODO -- here I'm just using some fictional value
1139        endif
1140        if(WaveExists(iErr4)==0  && WaveExists(inten4) != 0)
1141                Duplicate/O inten4,iErr4
1142                Wave iErr4=iErr4
1143//              iErr4 = 1+sqrt(inten4+0.75)                     // can't use this -- it applies to counts, not intensity (already a count rate...)
1144                iErr4 = sqrt(inten4+0.75)                       // TODO -- here I'm just using some fictional value
1145        endif
1146
1147        // TODO -- nq will need to be larger, once the back detector is installed
1148        //
1149        // note that the back panel of 320x320 (1mm res) results in 447 data points!
1150        // - so I upped nq to 600
1151
1152        if(cmpstr(type,"B") == 0)
1153                nq = 8000
1154        else
1155                nq=600
1156        endif
1157
1158//******TODO****** -- where to put the averaged data -- right now, folderStr is forced to ""   
1159//      SetDataFolder $("root:"+folderStr)              //should already be here, but make sure...     
1160        Make/O/D/N=(nq)  $(folderPath+":"+"iBin_qxqy"+"_"+type)
1161        Make/O/D/N=(nq)  $(folderPath+":"+"qBin_qxqy"+"_"+type)
1162        Make/O/D/N=(nq)  $(folderPath+":"+"nBin_qxqy"+"_"+type)
1163        Make/O/D/N=(nq)  $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1164        Make/O/D/N=(nq)  $(folderPath+":"+"eBin_qxqy"+"_"+type)
1165        Make/O/D/N=(nq)  $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1166       
1167        Wave iBin_qxqy = $(folderPath+":"+"iBin_qxqy_"+type)
1168        Wave qBin_qxqy = $(folderPath+":"+"qBin_qxqy"+"_"+type)
1169        Wave nBin_qxqy = $(folderPath+":"+"nBin_qxqy"+"_"+type)
1170        Wave iBin2_qxqy = $(folderPath+":"+"iBin2_qxqy"+"_"+type)
1171        Wave eBin_qxqy = $(folderPath+":"+"eBin_qxqy"+"_"+type)
1172        Wave eBin2D_qxqy = $(folderPath+":"+"eBin2D_qxqy"+"_"+type)
1173       
1174       
1175//      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
1176// TODO: not sure if I want to set dQ in x or y direction...
1177        // the short dimension is the 8mm tubes, use this direction as dQ?
1178        // 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]
1179        // WRT the beam center. use qx or qy directly. Still not happy with this way...
1180
1181
1182        qBin_qxqy[] =  p*delQ   
1183        SetScale/P x,0,delQ,"",qBin_qxqy                //allows easy binning
1184
1185        iBin_qxqy = 0
1186        iBin2_qxqy = 0
1187        eBin_qxqy = 0
1188        eBin2D_qxqy = 0
1189        nBin_qxqy = 0   //number of intensities added to each bin
1190
1191// now there are situations of:
1192// 1 panel
1193// 2 panels
1194// 4 panels
1195//
1196// this needs to be a double loop now...
1197// TODO:
1198// -- the iErr (=2D) wave and accumulation of error is NOT CALCULATED CORRECTLY YET
1199// -- verify the 2D error propagation by reducing it to 1D error
1200//
1201//
1202// The 1D error does not use iErr, and IS CALCULATED CORRECTLY
1203//
1204// x- the solid angle per pixel will be present for WORK data other than RAW, but not for RAW
1205
1206//
1207// if any of the masks don't exist, display the error, and proceed with the averaging, using all data
1208        if(maskMissing == 1)
1209                Print "Mask file not found for at least one detector - so all data is used"
1210        endif
1211       
1212        NVAR gIgnoreDetB = root:Packages:NIST:VSANS:Globals:gIgnoreDetB
1213        if(gIgnoreDetB && cmpstr(type,"B") == 0)
1214                maskMissing = 1
1215                Print "Mask skipped for B due to possible mismatch (Panel B ignored in preferences)"
1216        endif
1217
1218        Variable mask_val
1219// use set 1 (no number) only
1220        if(nSets >= 1)
1221                xDim=DimSize(inten,0)
1222                yDim=DimSize(inten,1)
1223       
1224                for(ii=0;ii<xDim;ii+=1)
1225                        for(jj=0;jj<yDim;jj+=1)
1226                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1227                                qVal = qTotal[ii][jj]
1228                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1229                                val = inten[ii][jj]
1230                               
1231                                if(isVCALC || maskMissing)              // mask_val == 0 == keep, mask_val == 1 = YES, mask out the point
1232                                        mask_val = 0
1233                                else
1234                                        mask_val = mask[ii][jj]
1235                                endif
1236                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1237                                        iBin_qxqy[binIndex] += val
1238                                        iBin2_qxqy[binIndex] += val*val
1239                                        eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1240                                        nBin_qxqy[binIndex] += 1
1241                                endif
1242                        endfor
1243                endfor
1244               
1245        endif
1246
1247// add in set 2 (set 1 already done)
1248        if(nSets >= 2)
1249                xDim=DimSize(inten2,0)
1250                yDim=DimSize(inten2,1)
1251       
1252                for(ii=0;ii<xDim;ii+=1)
1253                        for(jj=0;jj<yDim;jj+=1)
1254                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1255                                qVal = qTotal2[ii][jj]
1256                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1257                                val = inten2[ii][jj]
1258                               
1259                                if(isVCALC || maskMissing)
1260                                        mask_val = 0
1261                                else
1262                                        mask_val = mask2[ii][jj]
1263                                endif
1264                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1265                                        iBin_qxqy[binIndex] += val
1266                                        iBin2_qxqy[binIndex] += val*val
1267                                        eBin2D_qxqy[binIndex] += iErr2[ii][jj]*iErr2[ii][jj]
1268                                        nBin_qxqy[binIndex] += 1
1269                                endif
1270                        endfor
1271                endfor
1272               
1273        endif
1274
1275// add in set 3 and 4 (set 1 and 2 already done)
1276        if(nSets == 4)
1277                xDim=DimSize(inten3,0)
1278                yDim=DimSize(inten3,1)
1279       
1280                for(ii=0;ii<xDim;ii+=1)
1281                        for(jj=0;jj<yDim;jj+=1)
1282                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1283                                qVal = qTotal3[ii][jj]
1284                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1285                                val = inten3[ii][jj]
1286                               
1287                                if(isVCALC || maskMissing)
1288                                        mask_val = 0
1289                                else
1290                                        mask_val = mask3[ii][jj]
1291                                endif
1292                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1293                                        iBin_qxqy[binIndex] += val
1294                                        iBin2_qxqy[binIndex] += val*val
1295                                        eBin2D_qxqy[binIndex] += iErr3[ii][jj]*iErr3[ii][jj]
1296                                        nBin_qxqy[binIndex] += 1
1297                                endif
1298                        endfor
1299                endfor
1300               
1301               
1302                xDim=DimSize(inten4,0)
1303                yDim=DimSize(inten4,1)
1304       
1305                for(ii=0;ii<xDim;ii+=1)
1306                        for(jj=0;jj<yDim;jj+=1)
1307                                //qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1308                                qVal = qTotal4[ii][jj]
1309                                binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1310                                val = inten4[ii][jj]
1311                               
1312                                if(isVCALC || maskMissing)
1313                                        mask_val = 0
1314                                else
1315                                        mask_val = mask4[ii][jj]
1316                                endif
1317                                if (numType(val)==0 && mask_val == 0)           //count only the good points, ignore Nan or Inf
1318                                        iBin_qxqy[binIndex] += val
1319                                        iBin2_qxqy[binIndex] += val*val
1320                                        eBin2D_qxqy[binIndex] += iErr4[ii][jj]*iErr4[ii][jj]
1321                                        nBin_qxqy[binIndex] += 1
1322                                endif
1323                        endfor
1324                endfor
1325               
1326        endif
1327
1328
1329// after looping through all of the data on the panels, calculate errors on I(q),
1330// just like in CircSectAve.ipf
1331// TODO:
1332// -- 2D Errors were (maybe) properly acculumated through reduction, so this loop of calculations is NOT VERIFIED (yet)
1333// x- the error on the 1D intensity, is correctly calculated as the standard error of the mean.
1334        for(ii=0;ii<nq;ii+=1)
1335                if(nBin_qxqy[ii] == 0)
1336                        //no pixels in annuli, data unknown
1337                        iBin_qxqy[ii] = 0
1338                        eBin_qxqy[ii] = 1
1339                        eBin2D_qxqy[ii] = NaN
1340                else
1341                        if(nBin_qxqy[ii] <= 1)
1342                                //need more than one pixel to determine error
1343                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1344                                eBin_qxqy[ii] = 1
1345                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1346                        else
1347                                //assume that the intensity in each pixel in annuli is normally distributed about mean...
1348                                //  -- this is correctly calculating the error as the standard error of the mean, as
1349                                //    was always done for SANS as well.
1350                                iBin_qxqy[ii] /= nBin_qxqy[ii]
1351                                avesq = iBin_qxqy[ii]^2
1352                                aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1353                                var = aveisq-avesq
1354                                if(var<=0)
1355                                        eBin_qxqy[ii] = 1e-6
1356                                else
1357                                        eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1358                                endif
1359                                // and calculate as it is propagated pixel-by-pixel
1360                                eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1361                        endif
1362                endif
1363        endfor
1364       
1365        eBin2D_qxqy = sqrt(eBin2D_qxqy)         // as equation (3) of John's memo
1366       
1367        // find the last non-zero point, working backwards
1368        val=nq
1369        do
1370                val -= 1
1371        while((nBin_qxqy[val] == 0) && val > 0)
1372       
1373//      print val, nBin_qxqy[val]
1374        DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1375
1376        if(val == 0)
1377                // all the points were deleted
1378                return(0)
1379        endif
1380       
1381       
1382        // since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1383        // find the first non-zero point, working forwards
1384        val = -1
1385        do
1386                val += 1
1387        while(nBin_qxqy[val] == 0)     
1388        DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1389
1390        // ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1391        val = numpnts(nBin_qxqy)-1
1392        do
1393                if(nBin_qxqy[val] == 0)
1394                        DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1395                endif
1396                val -= 1
1397        while(val>0)
1398       
1399        // TODO:
1400        // -- This is where I calculate the resolution in SANS (see CircSectAve)
1401        // -- use the isVCALC flag to exclude VCALC from the resolution calculation if necessary
1402        // -- from the top of the function, folderStr = work folder, type = "FLRTB" or other type of averaging
1403        //
1404        nq = numpnts(qBin_qxqy)
1405        Make/O/D/N=(nq)  $(folderPath+":"+"sigmaQ_qxqy"+"_"+type)
1406        Make/O/D/N=(nq)  $(folderPath+":"+"qBar_qxqy"+"_"+type)
1407        Make/O/D/N=(nq)  $(folderPath+":"+"fSubS_qxqy"+"_"+type)
1408        Wave sigmaq = $(folderPath+":"+"sigmaQ_qxqy_"+type)
1409        Wave qbar = $(folderPath+":"+"qBar_qxqy_"+type)
1410        Wave fsubs = $(folderPath+":"+"fSubS_qxqy_"+type)
1411                               
1412
1413        Variable ret1,ret2,ret3
1414        Variable lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses
1415
1416// TODO: check the units of all of the inputs
1417
1418// lambda = wavelength [A]
1419        lambda = V_getWavelength(folderStr)
1420       
1421// lambdaWidth = [dimensionless]
1422        lambdaWidth = V_getWavelength_spread(folderStr)
1423       
1424// DDet = detector pixel resolution [cm]        **assumes square pixel
1425        // V_getDet_pixel_fwhm_x(folderStr,detStr)
1426        // V_getDet_pixel_fwhm_y(folderStr,detStr)
1427//      DDet = 0.8              // TODO -- this is hard-wired
1428
1429        if(strlen(type) == 1)
1430                // it's "B"
1431                DDet = V_getDet_pixel_fwhm_x(folderStr,type)            // value is already in cm
1432        else
1433                DDet = V_getDet_pixel_fwhm_x(folderStr,type[0,1])               // value is already in cm
1434        endif
1435               
1436// apOff = sample aperture to sample distance [cm]
1437        apOff = 10              // TODO -- this is hard-wired
1438       
1439// S1 = source aperture diameter [mm]
1440// may be either circle or rectangle
1441        String s1_shape="",bs_shape=""
1442        Variable width,height,equiv_S1,equiv_bs
1443       
1444       
1445        s1_shape = V_getSourceAp_shape(folderStr)
1446        if(cmpstr(s1_shape,"CIRCLE") == 0)
1447                S1 = str2num(V_getSourceAp_size(folderStr))
1448        else
1449                S1 = V_getSourceAp_height(folderStr)            // TODO: need the width or at least an equivalent diameter
1450        endif
1451       
1452       
1453// S2 = sample aperture diameter [cm]
1454// as of 3/2018, the "internal" sample aperture is not in use, only the external
1455// TODO : verify the units on the Ap2 (external)
1456// sample aperture 1(internal) is set to report "12.7 mm" as a STRING
1457// sample aperture 2(external) reports the number typed in...
1458//
1459// so I'm trusting [cm] is in the file
1460        S2 = V_getSampleAp2_size(folderStr)*10          // sample ap 1 or 2? 2 = the "external", convert to [mm]
1461       
1462// L1 = source to sample distance [m]
1463        L1 = V_getSourceAp_distance(folderStr)/100
1464
1465// L2 = sample to detector distance [m]
1466// take the first two characters of the "type" to get the correct distance.
1467// if the type is say, MLRTB, then the implicit assumption in combining all four panels is that the resolution
1468// is not an issue for the slightly different distances.
1469        if(strlen(type) == 1)
1470                // it's "B"
1471                L2 = V_getDet_ActualDistance(folderStr,type)/100                //convert cm to m
1472        else
1473                L2 = V_getDet_ActualDistance(folderStr,type[0,1])/100           //convert cm to m
1474        endif
1475       
1476// BS = beam stop diameter [mm]
1477//TODO:? which BS is in? carr2, carr3, none?
1478// -- need to check the detector, num_beamstops field, then description, then shape/size or shape/height and shape/width
1479//
1480// TODO: the values in the file are incorrect!!! BS = 1000 mm diameter!!!
1481//      BS = V_getBeamStopC2_size(folderStr)            // Units are [mm]
1482        BS = 25.4                       //TODO hard-wired value
1483       
1484//      bs_shape = V_getBeamStopC2_shape(folderStr)
1485//      if(cmpstr(s1_shape,"CIRCLE") == 0)
1486//              bs = V_getBeamStopC2_size(folderStr)
1487//      else
1488//              bs = V_getBeamStopC2_height(folderStr) 
1489//      endif
1490
1491
1492       
1493// del_r = step size [mm] = binWidth*(mm/pixel)
1494        del_r = 1*8                     // TODO: 8mm/pixel hard-wired
1495
1496// usingLenses = flag for lenses = 0 if no lenses, non-zero if lenses are in-beam
1497        usingLenses = 0
1498
1499if(cmpstr(detStr,"FL")==0)
1500        Print "(FL) Resolution lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses"
1501        Print lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses
1502endif
1503
1504
1505// TODO:
1506// this is the point where I need to switch on the different collimation types (white beam, slit, Xtal, etc)
1507// to calculate the correct resolution, or fill the waves with the correct "flags"
1508//
1509
1510// For white beam data, the wavelength distribution can't be represented as a gaussian, but all of the other
1511//  geometric corrections still apply. Passing zero for the lambdaWidth will return the geometry contribution,
1512//  as long as the wavelength can be handled separately. It appears to be correct to do as a double integral,
1513//  with the inner(lambda) calculated first, then the outer(geometry).
1514//
1515
1516// possible values are:
1517//
1518// pinhole
1519// pinhole_whiteBeam
1520// convergingPinholes
1521//
1522// *slit data should be reduced using the slit routine, not here, proceed but warn
1523// narrowSlit
1524// narrowSlit_whiteBeam
1525
1526        if(cmpstr(collimationStr,"pinhole") == 0)
1527
1528                ii=0
1529                do
1530                        V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3)
1531                        sigmaq[ii] = ret1       
1532                        qbar[ii] = ret2
1533                        fsubs[ii] = ret3       
1534                        ii+=1
1535                while(ii<nq)
1536       
1537        endif
1538       
1539
1540        if(cmpstr(collimationStr,"pinhole_whiteBeam") == 0)
1541
1542//              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution.
1543// the white beam distribution will need to be flagged some other way
1544//
1545                lambdaWidth = 0
1546               
1547                ii=0
1548                do
1549                        V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3)
1550                        sigmaq[ii] = ret1       
1551                        qbar[ii] = ret2
1552                        fsubs[ii] = ret3       
1553                        ii+=1
1554                while(ii<nq)
1555       
1556        endif
1557
1558        if(cmpstr(collimationStr,"convergingPinholes") == 0)
1559
1560//              set usingLenses == 1 so that the Gaussian resolution calculation will be for a focus condition
1561//
1562                usingLenses = 1
1563               
1564                ii=0
1565                do
1566                        V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3)
1567                        sigmaq[ii] = ret1       
1568                        qbar[ii] = ret2
1569                        fsubs[ii] = ret3       
1570                        ii+=1
1571                while(ii<nq)
1572       
1573        endif
1574
1575
1576// should not end up here, except for odd testing cases
1577        if(cmpstr(collimationStr,"narrowSlit") == 0)
1578
1579                Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???"
1580                ii=0
1581                do
1582                        V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3)
1583                        sigmaq[ii] = ret1       
1584                        qbar[ii] = ret2
1585                        fsubs[ii] = ret3       
1586                        ii+=1
1587                while(ii<nq)
1588       
1589        endif
1590       
1591// should not end up here, except for odd testing cases
1592        if(cmpstr(collimationStr,"narrowSlit_whiteBeam") == 0)
1593
1594//              set lambdaWidth == 0 so that the gaussian resolution calculates only the geometry contribution.
1595// the white beam distribution will need to be flagged some other way
1596//
1597                Print "??? Slit data is being averaged as pinhole - reset the AVERAGE parameters in the protocol ???"
1598
1599                lambdaWidth = 0
1600               
1601                ii=0
1602                do
1603                        V_getResolution(qBin_qxqy[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,ret1,ret2,ret3)
1604                        sigmaq[ii] = ret1       
1605                        qbar[ii] = ret2
1606                        fsubs[ii] = ret3       
1607                        ii+=1
1608                while(ii<nq)
1609       
1610        endif
1611
1612
1613
1614               
1615        SetDataFolder root:
1616       
1617        return(0)
1618End
1619
1620
Note: See TracBrowser for help on using the repository browser.