PROGRAM newx 
************************************************************************
* Lumps and interpolates an x-ordered x-y-z file.                        *
* From standard input it reads:                                        *
* input file name                                                      *
* filename containing array of x-values to used on output              *
* The resulting x y z file flushed to standard output.                   *
* It interpolates linearly over data-points around the output x-values.*
* Points that fall outside the available data-set are omitted.         *
* If the new grid is actually smaller than the original one,           *
* or if some datapoints are missing, a linear  interpolation is made   *
* between the two neighboring original data-points.                    *
************************************************************************
      parameter(nphys=100000)
      real eps,sy,sz,x,y,z
     *    ,xx(nphys),yy(nphys),zz(nphys)
     *    ,xmx,xmn,xlast,xnext,ylast,ynext,zlast,znext,wlow,whigh
     *    ,a,b,xnw(nphys),ynw(nphys),znw(nphys),tmx,tmn
      CHARACTER*40 flin,flix
      INTEGER I,J,n,mm,nn,ilast                                          
      
      read(*,'(a40)') flin
      read(*,'(a40)') flix
      open(23,FILE=flin)                                      
      mm=10000 
      xmx=-1e20
      xmn=1e20
      do 10 i=1,mm                                                  
       READ(23,*,END=11) xx(i),yy(i),zz(i)
       if (xx(i).lt.xmn) xmn=xx(i)
       if (xx(i).gt.xmx) xmx=xx(i)
10    continue   
11    mm=i-1
      close(23)

      open(24,FILE=flix)                                      
      nn=10000 
      do 20 i=1,nn                                                  
       READ(24,*,END=21) xnw(i)
20    continue   
21    nn=i-1
      close(24)

******lump to tsp intervals 
      ilast=1
      eps=1e-8
      do 31 i=1,nn
       if (i.le.(nn-1)) then
        tmx=xnw(i+1)
       else
        tmx=xnw(i)
       endif
       if (i.ge.2) then
        tmn=xnw(i-1)
       else
        tmn=xnw(i)
       endif
       sy=0
       sz=0
       n=0
       do 200 j=ilast,mm
	if (xx(j).lt.tmn) goto 200
        if (xx(j).ge.tmx) goto 210
        sy=sy+yy(j)
	sz=sz+zz(j)
        n=n+1
200     continue
210    ilast=j
       if (n.eq.0) then
*****  first omit all points that are outside interval       
        if ((xnw(i).lt.xmn).or.(xnw(i).gt.xmx)) then
         y=1
*****   otherwise take appropriate weighted average between
*****   last point and first subsequent point       
        else 
         xlast=xx(ilast-1) 
         xnext=xx(ilast)
         ylast=yy(ilast-1)
         ynext=yy(ilast)
	 zlast=zz(ilast-1)
	 znext=zz(ilast)
         whigh=(xnw(i)-xlast)/(xnext-xlast)
         wlow=1.-whigh
         y=ylast*wlow+ynext*whigh
	 z=zlast*wlow+znext*whigh
        endif
       endif
       if (n.eq.1) then
	y=sy
	z=sz
       endif
       if (n.ge.2) then
        y=sy/n
	z=sz/n
       endif
       write(*,*) xnw(i),y,z
31    continue                   

      end