subroutine mag_by(x0,y0,z0,byave)

      implicit none

      integer i, j, k, n, repno
      parameter(n=1500)
      double precision x0, y0, z0
      double precision x(n), y(n), z(n), by(n), tmp
      double precision dr(8),rat(8),byp(8), rattot, bymrat, byave
      integer xrep, yrep, zrep
      parameter(xrep=3)
      parameter(yrep=1)
      parameter(zrep=48)
c
      open(11,file='y.dat',status='old')
      i=1
 20   read(11,*,end=30) x(i), y(i), z(i), by(i)
      i=i+1
      goto 20
 30   tmp = 1
c
      byave=0
      repno=0
      do k=1,i-1
         if(  (x0.ge.x(k).and.x0.le.x(k+xrep) .and.
     &        abs(x(k)-x(k+xrep)).lt.1. )     .and.
     &        (y0.ge.y(k).and.y0.le.y(k+yrep) .and.
     &        abs(y(k)-y(k+yrep)).lt.3. )     .and.
     &        (z0.ge.z(k).and.z0.le.z(k+zrep) .and.
     &        abs(z(k)-z(k+zrep)).lt.2. )   ) then
            repno=k
         endif
      enddo
c
      do k=1,i-1
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno)) byp(1)=by(k)
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno)) byp(2)=by(k)
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno+zrep)) byp(3)=by(k)
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno+zrep)) byp(4)=by(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno)) byp(5)=by(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno)) byp(6)=by(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno+zrep)) byp(7)=by(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno+zrep)) byp(8)=by(k)
      enddo
c
      dr(1)=sqrt((x0-x(repno))**2      + (y0-y(repno))**2      +
     &     (z0-z(repno))**2)
      dr(2)=sqrt((x0-x(repno))**2      + (y0-y(repno+yrep))**2 +
     &     (z0-z(repno))**2)
      dr(3)=sqrt((x0-x(repno))**2      + (y0-y(repno))**2      +
     &     (z0-z(repno+zrep))**2)
      dr(4)=sqrt((x0-x(repno))**2      + (y0-y(repno+yrep))**2 +
     &     (z0-z(repno+zrep))**2)
      dr(5)=sqrt((x0-x(repno+xrep))**2 + (y0-y(repno))**2      +
     &     (z0-z(repno))**2)
      dr(6)=sqrt((x0-x(repno+xrep))**2 + (y0-y(repno+yrep))**2 +
     &     (z0-z(repno))**2)
      dr(7)=sqrt((x0-x(repno+xrep))**2 + (y0-y(repno))**2      +
     &     (z0-z(repno+zrep))**2)
      dr(8)=sqrt((x0-x(repno+xrep))**2 + (y0-y(repno+yrep))**2 +
     &     (z0-z(repno+zrep))**2)
c
      do j=1,8
         if(dr(j).eq.0) dr(j)=1e-10
         rat(j)=1./dr(j)
      enddo
c
      rattot=0.
      bymrat=0.
      do j=1,8
         rattot=rattot+rat(j)
         bymrat=bymrat+rat(j)*byp(j)
      enddo
      byave=bymrat/rattot

      close(11)

 100  return
      end