source: sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS/CircSectAve.ipf @ 1212

Last change on this file since 1212 was 1212, checked in by krzywon, 3 years ago

Combine elliptical averaging with circular averaging and wrote the binning routine. Need to verify calculations and fix issues with elliptical drawing routine.

File size: 19.7 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.1
4
5//*********************************************
6//              For AVERAGE and for DRAWING
7//                      DRAWING routines only use a subset of the total list, since saving, naming, etc. don't apply
8//              (10) possible keywords, some numerical, some string values
9//              AVTYPE=string           string from set {Circular,Annular,Rectangular,Sector,2D_ASCII,QxQy_ASCII,PNG_Graphic}
10//              PHI=value                       azimuthal angle (-90,90)
11//              DPHI=value                      +/- angular range around phi for average
12//              WIDTH=value             total width of rectangular section, in pixels
13//              SIDE=string             string from set {left,right,both} **note NOT capitalized
14//              QCENTER=value           q-value (1/A) of center of annulus for annular average
15//              QDELTA=value            total width of annulus centered at QCENTER
16//              PLOT=string             string from set {Yes,No} = truth of generating plot of averaged data
17//              SAVE=string             string from set {Yes,No} = truth of saving averaged data to disk
18//              NAME=string             string from set {Auto,Manual} = Automatic name generation or Manual(dialog)
19//***********************************************
20
21
22// this function also does sector averaging
23//the parameters in the global keyword-string must have already been set somewhere
24//either directly, from the protocol, or from the Average_Panel
25//** the keyword-list has already been "pre-parsed" to send only Circular or Sector
26//averages to this routine. Rectangular or annular averages get done elsewhere
27// TYPE parameter determines which data folder to work from
28//
29//annnulus (step) size is currently fixed at 1 (variable dr, below)
30Function CircularAverageTo1D(type)
31        String type
32       
33        SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr          //this is the list that has it all
34        Variable isCircular = 0, isElliptical = 0
35       
36        if( cmpstr("Circular",StringByKey("AVTYPE",keyListStr,"=",";")) ==0)
37                isCircular = 1          //set a switch for later
38        Endif
39        if( cmpstr("Elliptical",StringByKey("AVTYPE",keyListStr,"=",";")) ==0)
40                isElliptical = 1                //set a switch for later
41        Endif
42       
43        //type is the data type to do the averaging on, and will be set as the current folder
44        //get the current displayed data (so the correct folder is used)
45        String destPath = "root:Packages:NIST:"+type
46       
47        //
48        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist,dr,ddr
49        Variable lambda,trans,raxes
50        WAVE reals = $(destPath + ":RealsRead")
51        WAVE/T textread = $(destPath + ":TextRead")
52//      String fileStr = textread[3]
53       
54        // center of detector, for non-linear corrections
55        NVAR pixelsX = root:myGlobals:gNPixelsX
56        NVAR pixelsY = root:myGlobals:gNPixelsY
57       
58        xcenter = pixelsX/2 + 0.5               // == 64.5 for 128x128 Ordela
59        ycenter = pixelsY/2 + 0.5               // == 64.5 for 128x128 Ordela
60       
61        // beam center, in pixels
62        x0 = reals[16]
63        y0 = reals[17]
64        //detector calibration constants
65        sx = reals[10]          //mm/pixel (x)
66        sx3 = reals[11]         //nonlinear coeff
67        sy = reals[13]          //mm/pixel (y)
68        sy3 = reals[14]         //nonlinear coeff
69       
70        dtsize = 10*reals[20]           //det size in mm
71        dtdist = 1000*reals[18]         // det distance in mm
72       
73        NVAR binWidth=root:Packages:NIST:gBinWidth
74       
75        dr = binWidth           // ***********annulus width set by user, default is one***********
76        ddr = dr*sx             //step size, in mm (this value should be passed to the resolution calculation, not dr 18NOV03)
77               
78        Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,dxbm,dybm,ii
79        Variable phi_rad,dphi_rad,phi_x,phi_y
80        Variable forward,mirror
81       
82        String side = StringByKey("SIDE",keyListStr,"=",";")
83//      Print "side = ",side
84       
85        if(!isCircular)         //must be sector avg (rectangular not sent to this function)
86                //convert from degrees to radians
87                phi_rad = (Pi/180)*NumberByKey("PHI",keyListStr,"=",";")
88                dphi_rad = (Pi/180)*NumberByKey("DPHI",keyListStr,"=",";")
89                //create cartesian values for unit vector in phi direction
90                phi_x = cos(phi_rad)
91                phi_y = sin(phi_rad)
92        Endif
93        if(isElliptical)
94                raxes = NumberByKey("RATIOAXES",keyListStr,"=",";")
95        EndIf
96       
97        /// data wave is data in the current folder which was set at the top of the function
98        WAVE data=$(destPath + ":data")
99        //Check for the existence of the mask, if not, make one (local to this folder) that is null
100       
101        if(WaveExists($"root:Packages:NIST:MSK:data") == 0)
102                Print "There is no mask file loaded (WaveExists)- the data is not masked"
103                Make/O/N=(pixelsX,pixelsY) $(destPath + ":mask")
104                Wave mask = $(destPath + ":mask")
105                mask = 0
106        else
107                Wave mask=$"root:Packages:NIST:MSK:data"
108        Endif
109       
110        //
111        //pixels within rcentr of beam center are broken into 9 parts (units of mm)
112        rcentr = 100            //original
113//      rcentr = 0
114        // values for error if unable to estimate value
115        //large_num = 1e10
116        large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1
117        small_num = 1e-10
118       
119        // output wave are expected to exist (?) initialized to zero, what length?
120        // 200 points on VAX --- use 300 here, or more if SAXS data is used with 1024x1024 detector (1000 pts seems good)
121        Variable defWavePts=500
122        Make/O/N=(defWavePts) $(destPath + ":qval"),$(destPath + ":aveint")
123        Make/O/N=(defWavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave")
124        Make/O/N=(defWavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar")
125
126        WAVE qval = $(destPath + ":qval")
127        WAVE aveint = $(destPath + ":aveint")
128        WAVE ncells = $(destPath + ":ncells")
129        WAVE dsq = $(destPath + ":dsq")
130        WAVE sigave = $(destPath + ":sigave")
131        WAVE qbar = $(destPath + ":QBar")
132        WAVE sigmaq = $(destPath + ":SigmaQ")
133        WAVE fsubs = $(destPath + ":fSubS")
134       
135        qval = 0
136        aveint = 0
137        ncells = 0
138        dsq = 0
139        sigave = 0
140        qbar = 0
141        sigmaq = 0
142        fsubs = 0
143
144        dtdis2 = dtdist^2
145        nq = 1
146        xoffst=0
147        //distance of beam center from detector center
148        dxbm = FX(x0,sx3,xcenter,sx)
149        dybm = FY(y0,sy3,ycenter,sy)
150               
151        //BEGIN AVERAGE **********
152        Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1
153        Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p,rho
154       
155        // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
156        // loop index corresponds to FORTRAN (old code)
157        // and the IGOR array indices must be adjusted (-1) to the correct address
158        ii=1
159        do
160                xi = ii
161                dxi = FX(xi,sx3,xcenter,sx)
162                dx = dxi-dxbm           //dx and dy are in mm
163               
164                jj = 1
165                do
166                        data_pixel = data[ii-1][jj-1]           //assign to local variable
167                        yj = jj
168                        dyj = FY(yj,sy3,ycenter,sy)
169                        dy = dyj - dybm
170                        if(!(mask[ii-1][jj-1]))                 //masked pixels = 1, skip if masked (this way works...)
171                                dr2 = (dx^2 + dy^2)^(0.5)               //distance from beam center NOTE dr2 used here - dr used above
172                                if(dr2>rcentr)          //keep pixel whole
173                                        nd = 1
174                                        fd = 1
175                                else                            //break pixel into 9 equal parts
176                                        nd = 3
177                                        fd = 2
178                                endif
179                                nd2 = nd^2
180                                ll = 1          //"el-el" loop index
181                                do
182                                        dxx = dx + (ll - fd)*sx/3
183                                        kk = 1
184                                        do
185                                                dyy = dy + (kk - fd)*sy/3
186                                                if(isCircular || isElliptical)
187                                                        //use all pixels
188                                                        if (isElliptical)
189                                                                rho = atan(dyy/dxx) - phi_rad
190                                                                nq = IncrementEllipticalPixel(data_pixel,ddr,dxx,dyy,rho,raxes, aveint,dsq,ncells,nq,nd2)
191                                                        else
192                                                                //circular average
193                                                                nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
194                                                        EndIf
195                                                else
196                                                        //a sector average - determine azimuthal angle
197                                                        dphi_p = dphi_pixel(dxx,dyy,phi_x,phi_y)
198                                                        if(dphi_p < dphi_rad)
199                                                                forward = 1                     //within forward sector
200                                                        else
201                                                                forward = 0
202                                                        Endif
203                                                        if((Pi - dphi_p) < dphi_rad)
204                                                                mirror = 1              //within mirror sector
205                                                        else
206                                                                mirror = 0
207                                                        Endif
208                                                        //check if pixel lies within allowed sector(s)
209                                                        if(cmpstr(side,"both")==0)              //both sectors
210                                                                if ( mirror || forward)
211                                                                        //increment
212                                                                        nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
213                                                                Endif
214                                                        else
215                                                                if(cmpstr(side,"right")==0)             //forward sector only
216                                                                        if(forward)
217                                                                                //increment
218                                                                                nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
219                                                                        Endif
220                                                                else                    //mirror sector only
221                                                                        if(mirror)
222                                                                                //increment
223                                                                                nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
224                                                                        Endif
225                                                                Endif
226                                                        Endif           //allowable sectors
227                                                Endif   //circular or sector check
228                                                kk+=1
229                                        while(kk<=nd)
230                                        ll += 1
231                                while(ll<=nd)
232                        Endif           //masked pixel check
233                        jj += 1
234                while (jj<=pixelsY)
235                ii += 1
236        while(ii<=pixelsX)              //end of the averaging
237               
238        //compute q-values and errors
239        Variable ntotal,rr,theta,avesq,aveisq,var
240       
241        lambda = reals[26]
242        ntotal = 0
243        kk = 1
244        do
245                rr = (2*kk-1)*ddr/2
246                theta = 0.5*atan(rr/dtdist)
247                qval[kk-1] = (4*Pi/lambda)*sin(theta)
248                if(ncells[kk-1] == 0)
249                        //no pixels in annuli, data unknown
250                        aveint[kk-1] = 0
251                        sigave[kk-1] = large_num
252                else
253                        if(ncells[kk-1] <= 1)
254                                //need more than one pixel to determine error
255                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
256                                sigave[kk-1] = large_num
257                        else
258                                //assume that the intensity in each pixel in annuli is normally
259                                // distributed about mean...
260                                aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
261                                avesq = aveint[kk-1]^2
262                                aveisq = dsq[kk-1]/ncells[kk-1]
263                                var = aveisq-avesq
264                                if(var<=0)
265//                                      sigave[kk-1] = small_num
266                                        sigave[kk-1] = large_num                //if there are zero counts, make error large_num to be a warning flag
267                                else
268                                        sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1))
269                                endif
270                        endif
271                endif
272                ntotal += ncells[kk-1]
273                kk+=1
274        while(kk<=nq)
275       
276        //Print "NQ = ",nq
277        // data waves were defined as 300 points (=defWavePts), but now have less than that (nq) points
278        // use DeletePoints to remove junk from end of waves
279        //WaveStats would be a more foolproof implementation, to get the # points in the wave
280        Variable startElement,numElements
281        startElement = nq
282        numElements = defWavePts - startElement
283        DeletePoints startElement,numElements, qval,aveint,ncells,dsq,sigave
284       
285        //////////////end of VAX sector_ave()
286               
287        //angle dependent transmission correction
288        Variable uval,arg,cos_th
289        lambda = reals[26]
290        trans = reals[4]
291
292//
293//  The transmission correction is now done at the ADD step, in DetCorr()
294//     
295//      ////this section is the trans_correct() VAX routine
296//      if(trans<0.1)
297//              Print "***transmission is less than 0.1*** and is a significant correction"
298//      endif
299//      if(trans==0)
300//              Print "***transmission is ZERO*** and has been reset to 1.0 for the averaging calculation"
301//              trans = 1
302//      endif
303//      //optical thickness
304//      uval = -ln(trans)               //use natural logarithm
305//      //apply correction to aveint[]
306//      //index from zero here, since only working with IGOR waves
307//      ii=0
308//      do
309//              theta = 2*asin(lambda*qval[ii]/(4*pi))
310//              cos_th = cos(theta)
311//              arg = (1-cos_th)/cos_th
312//              if((uval<0.01) || (cos_th>0.99))                //OR
313//                      //small arg, approx correction
314//                      aveint[ii] /= 1-0.5*uval*arg
315//              else
316//                      //large arg, exact correction
317//                      aveint[ii] /= (1-exp(-uval*arg))/(uval*arg)
318//              endif
319//              ii+=1
320//      while(ii<nq)
321//      //end of transmission/pathlength correction
322
323// ***************************************************************
324//
325// Do the extra 3 columns of resolution calculations starting here.
326//
327// ***************************************************************
328
329        Variable L2 = reals[18]
330        Variable BS = reals[21]
331        Variable S1 = reals[23]
332        Variable S2 = reals[24]
333        Variable L1 = reals[25]
334        lambda = reals[26]
335        Variable lambdaWidth = reals[27]
336        String detStr=textRead[9]
337       
338        Variable usingLenses = reals[28]                //new 2007
339
340        //Two parameters DDET and APOFF are instrument dependent.  Determine
341        //these from the instrument name in the header.
342        //From conversation with JB on 01.06.99 these are the current
343        //good values
344
345        Variable DDet
346        NVAR apOff = root:myGlobals:apOff               //in cm
347       
348//      DDet = DetectorPixelResolution(fileStr,detStr)          //needs detector type and beamline
349        //note that reading the detector pixel size from the header ASSUMES SQUARE PIXELS! - Jan2008
350        DDet = reals[10]/10                     // header value (X) is in mm, want cm here
351       
352       
353        //Width of annulus used for the average is gotten from the
354        //input dialog before.  This also must be passed to the resolution
355        //calculator. Currently the default is dr=1 so just keeping that.
356
357        //Go from 0 to nq doing the calc for all three values at
358        //every Q value
359
360        ii=0
361
362        Variable ret1,ret2,ret3
363        do
364                getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,usingLenses,ret1,ret2,ret3)
365                sigmaq[ii] = ret1       
366                qbar[ii] = ret2
367                fsubs[ii] = ret3       
368                ii+=1
369        while(ii<nq)
370        DeletePoints startElement,numElements, sigmaq, qbar, fsubs
371
372// End of resolution calculations
373// ***************************************************************
374       
375        //Plot the data in the Plot_1d window
376        Avg_1D_Graph(aveint,qval,sigave)
377
378        //get rid of the default mask, if one was created (it is in the current folder)
379        //don't just kill "mask" since it might be pointing to the one in the MSK folder
380        Killwaves/Z $(destPath+":mask")
381       
382        //return to root folder (redundant)
383        SetDataFolder root:
384       
385        Return 0
386End
387
388//returns nq, new number of q-values
389//arrays aveint,dsq,ncells are also changed by this function
390//
391Function IncrementPixel(dataPixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
392        Variable dataPixel,ddr,dxx,dyy
393        Wave aveint,dsq,ncells
394        Variable nq,nd2
395       
396        Variable ir
397       
398        ir = trunc(sqrt(dxx*dxx+dyy*dyy)/ddr)+1
399        if (ir>nq)
400                nq = ir         //resets maximum number of q-values
401        endif
402        aveint[ir-1] += dataPixel/nd2           //ir-1 must be used, since ir is physical
403        dsq[ir-1] += dataPixel*dataPixel/nd2
404        ncells[ir-1] += 1/nd2
405       
406        Return nq
407End
408
409//returns nq, new number of q-values
410//arrays aveint,dsq,ncells are also changed by this function
411//
412Function IncrementEllipticalPixel(dataPixel,ddr,dxx,dyy,rho,raxes,aveint,dsq,ncells,nq,nd2)
413        Variable dataPixel,ddr,dxx,dyy,rho,raxes
414        Wave aveint,dsq,ncells
415        Variable nq,nd2
416       
417        Variable irCircular,ir
418       
419        irCircular = (sqrt(dxx*dxx+dyy*dyy)/ddr)
420        ir = irCircular*sqrt(cos(rho)*cos(rho) + raxes*raxes*sin(rho)*sin(rho))+1
421        if (ir>nq)
422                nq = ir         //resets maximum number of q-values
423        endif
424        aveint[ir-1] += dataPixel/nd2           //ir-1 must be used, since ir is physical
425        dsq[ir-1] += dataPixel*dataPixel/nd2
426        ncells[ir-1] += 1/nd2
427       
428        Return nq
429End
430
431//function determines azimuthal angle dphi that a vector connecting
432//center of detector to pixel makes with respect to vector
433//at chosen azimuthal angle phi -> [cos(phi),sin(phi)] = [phi_x,phi_y]
434//dphi is always positive, varying from 0 to Pi
435//
436Function dphi_pixel(dxx,dyy,phi_x,phi_y)
437        Variable dxx,dyy,phi_x,phi_y
438       
439        Variable val,rr,dot_prod
440       
441        rr = sqrt(dxx^2 + dyy^2)
442        dot_prod = (dxx*phi_x + dyy*phi_y)/rr
443        //? correct for roundoff error? - is this necessary in IGOR, w/ double precision?
444        if(dot_prod > 1)
445                dot_prod =1
446        Endif
447        if(dot_prod < -1)
448                dot_prod = -1
449        Endif
450       
451        val = acos(dot_prod)
452       
453        return val
454
455End
456
457//calculates the x distance from the center of the detector, w/nonlinear corrections
458//
459Function FX(xx,sx3,xcenter,sx)         
460        Variable xx,sx3,xcenter,sx
461       
462        Variable retval
463       
464        retval = sx3*tan((xx-xcenter)*sx/sx3)
465        Return retval
466End
467
468//calculates the y distance from the center of the detector, w/nonlinear corrections
469//
470Function FY(yy,sy3,ycenter,sy)         
471        Variable yy,sy3,ycenter,sy
472       
473        Variable retval
474       
475        retval = sy3*tan((yy-ycenter)*sy/sy3)
476        Return retval
477End
478
479//old function not called anymore, now "ave" button calls routine from AvgGraphics.ipf
480//to get input from panel rather than large prompt for missing parameters
481Function Ave_button(button0) : ButtonControl
482        String button0
483
484        // the button on the graph will average the currently displayed data
485        SVAR type=root:myGlobals:gDataDisplayType
486       
487        //Check for logscale data in "type" folder
488        SetDataFolder "root:Packages:NIST:"+type                //use the full path, so it will always work
489        String dest = "root:Packages:NIST:" + type
490       
491        NVAR isLogScale = $(dest + ":gIsLogScale")
492        if(isLogScale == 1)
493                //data is logscale, convert it back and reset the global
494                Duplicate/O $(dest + ":linear_data") $(dest + ":data")
495//              WAVE vlegend=$(dest + ":vlegend")
496        //  Make the color table linear scale
497//              vlegend = y
498                Variable/G $(dest + ":gIsLogScale") = 0         //copy to keep with the current data folder
499                SetDataFolder root:
500                //rename the button to reflect "isLin" - the displayed name must have been isLog
501                Button bisLog,title="isLin",rename=bisLin
502        Endif
503
504        //set data folder back to root
505        SetDataFolder root:
506       
507        //do the average - ask the user for what type of average
508        //ask the user for averaging paramters
509        Execute "GetAvgInfo()"
510       
511        //dispatch to correct averaging routine
512        //if you want to save the files, see Panel_DoAverageButtonProc(ctrlName) function
513        //for making a fake protocol (needed to write out data)
514        SVAR tempStr = root:myGlobals:Protocols:gAvgInfoStr
515        String choice = StringByKey("AVTYPE",tempStr,"=",";")
516        if(cmpstr("Rectangular",choice)==0)
517                //dispatch to rectangular average
518                RectangularAverageTo1D(type)
519        else
520                if(cmpstr("Annular",choice)==0)
521                        AnnularAverageTo1D(type)
522                else
523                        //circular or sector
524                        CircularAverageTo1D(type)
525                Endif
526        Endif
527       
528        Return 0
529End
530
531
532
533// -- seems to work, now I need to give it a name, add it to the list, and
534// make sure I've thought of all of the cases - then the average can be passed as case "Sector_PlusMinus"
535// and run through the normal average and writing routines.
536//
537//
538// -- depending on what value PHI has - it's [-90,90] "left" and "right" may not be what
539// you expect. so sorting the concatenated values may be necessary (always)
540//
541// -- need documentation of the definition of PHI, left, and right so that it can make better sense
542//              which quadrants of the detector become "negative" depending on the choice of phi. may need a
543//              switch after a little thinking.
544//
545// may want a variation of this where both sides are done, in separate files. but I think that's already
546// called a "sector" average. save them. load them. plot them.
547//
548//
549Function Sector_PlusMinus1D(type)
550        String type
551
552//      do the left side (-)
553// then hold that data in tmp_ waves
554// then do the right (+)
555// then concatenate the data
556
557// the button on the pink panel copies the two strings so they're the same
558        SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr          //this is the list that has it all
559
560        String oldStr = ""
561        String  CurPath="root:myGlobals:Plot_1D:"
562        String destPath = "root:Packages:NIST:"+type+":"
563
564        oldStr = StringByKey("SIDE",keyListStr,"=",";")
565
566// do the left first, and call it negative
567        keyListStr = ReplaceStringByKey("SIDE",keyListStr,"left","=",";")
568
569        CircularAverageTo1D(type)
570       
571        WAVE qval = $(destPath + "qval")
572        WAVE aveint = $(destPath + "aveint")
573        WAVE sigave = $(destPath + "sigave")
574        WAVE qbar = $(destPath + "QBar")
575        WAVE sigmaq = $(destPath + "SigmaQ")
576        WAVE fsubs = $(destPath + "fSubS")
577
578        // copy the average, set the q's negative
579        qval *= -1
580        Duplicate/O qval $(destPath+"tmp_q")
581        Duplicate/O aveint $(destPath+"tmp_i")
582        Duplicate/O sigave $(destPath+"tmp_s")
583        Duplicate/O qbar $(destPath+"tmp_qb")
584        Duplicate/O sigmaq $(destPath+"tmp_sq")
585        Duplicate/O fsubs $(destPath+"tmp_fs")
586       
587       
588// do the right side
589        keyListStr = ReplaceStringByKey("SIDE",keyListStr,"right","=",";")
590
591        CircularAverageTo1D(type)
592       
593        // concatenate
594        WAVE tmp_q = $(destPath + "tmp_q")
595        WAVE tmp_i = $(destPath + "tmp_i")
596        WAVE tmp_s = $(destPath + "tmp_s")
597        WAVE tmp_qb = $(destPath + "tmp_qb")
598        WAVE tmp_sq = $(destPath + "tmp_sq")
599        WAVE tmp_fs = $(destPath + "tmp_fs")
600
601        SetDataFolder destPath          //to get the concatenation in the right folder
602        Concatenate/NP/O {tmp_q,qval},tmp_cat
603        Duplicate/O tmp_cat qval
604        Concatenate/NP/O {tmp_i,aveint},tmp_cat
605        Duplicate/O tmp_cat aveint
606        Concatenate/NP/O {tmp_s,sigave},tmp_cat
607        Duplicate/O tmp_cat sigave
608        Concatenate/NP/O {tmp_qb,qbar},tmp_cat
609        Duplicate/O tmp_cat qbar
610        Concatenate/NP/O {tmp_sq,sigmaq},tmp_cat
611        Duplicate/O tmp_cat sigmaq
612        Concatenate/NP/O {tmp_fs,fsubs},tmp_cat
613        Duplicate/O tmp_cat fsubs
614
615// then sort
616        Sort qval, qval,aveint,sigave,qbar,sigmaq,fsubs
617
618// move these to the Plot_1D folder for plotting
619        Duplicate/O qval $(curPath+"xAxisWave")
620        Duplicate/O aveint $(curPath+"yAxisWave")
621        Duplicate/O sigave $(curPath+"yErrWave")
622       
623        keyListStr = ReplaceStringByKey("SIDE",keyListStr,oldStr,"=",";")
624
625        DoUpdate/W=Plot_1d
626       
627        // clean up
628        KillWaves/Z tmp_q,tmp_i,tmp_s,tmp_qb,tmp_sq,tmp_fs,tmp_cat
629       
630        SetDataFolder root:
631       
632        return(0)
633End
Note: See TracBrowser for help on using the repository browser.