c      subroutine mag_b(x0,y0,z0,bave,zrep,flag)
      program mag_b

      implicit none

      real HRNDM1, vect(6)

      include 'common.inc'

c
      zrep=48
c
      call rluxgo(3,4756874,0,0)  ! RANLUX luxury level 3

      do jj=1,1e3
      call ranlux(vect,4)
      x0= 0.
*9.0*vect(1)-4.5
      y0 = 0.
*4.0*vect(2)-2.0
      z0 = 18.0*vect(3)-9.0

      open(11,file='../z.dat',status='old')
      i=1
 20   read(11,*,end=30) x(i), y(i), z(i), b(i)
      i=i+1
      goto 20
 30   tmp = 1
c
c      bave=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)) bp(1)=b(k)
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno)) bp(2)=b(k)
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno+zrep)) bp(3)=b(k)
         if(x(k).eq.x(repno)     .and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno+zrep)) bp(4)=b(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno)) bp(5)=b(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno)) bp(6)=b(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno)     .and.
     &        z(k).eq.z(repno+zrep)) bp(7)=b(k)
         if(x(k).eq.x(repno+xrep).and.y(k).eq.y(repno+yrep).and.
     &        z(k).eq.z(repno+zrep)) bp(8)=b(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)**2)+(0.1/dr(j))+(0.1/dr(j)**3)
      enddo
c
      rattot=0.
      bmrat=0.
      do j=1,8
         rattot=rattot+rat(j)
         bmrat=bmrat+rat(j)*bp(j)
      enddo
      bave=bmrat/rattot

c      if(abs(x0).le.0.1.and.abs(y0).le.0.1)
c     &
      write(15,*) x0, y0, z0, bave

      close(11)

      enddo

 100  return
      end