      SUBROUTINE MAGFLD (X, Y, Z, BTX, BTY, BTZ, B, NREG, IDISC)

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'

      LOGICAL LFIRST
      CHARACTER BUFF
      INTEGER XI, YI, ZI

c X and Y are indexed from -30 to 30, with 61 points, while
c Z is indexed from 0 to 600 in steps of 5, with 121 points,
c so each point taken from Z must be divided by 5, then have
c 1 added to it before used as an index.  X, Y must have 31 
c added to them.  From there,
c all three axes must use the nint() function to round to the
c nearest integer, then passed into the three B-field arrays.

      PARAMETER(NX=61)
      PARAMETER(NY=61)
      PARAMETER(NZ=121)
      DIMENSION BX(NX,NY,NZ)
      DIMENSION BY(NX,NY,NZ)
      DIMENSION BZ(NX,NY,NZ)

      DATA LFIRST /.TRUE./
      SAVE LFIRST

      IDISC = 0

      IF (LFIRST) THEN

c make sure you remove the first 6 lines of the file first!
c these just say the units and quantities in the columns, so
c should be removed so the file starts immediately with the
c B-field data

         OPEN(U,FILE='magfld.txt',STATUS='OLD')

c X,Y,Z do not need to be read, simply passed over by 
c "read" commands to a buffer.  The array indeces will 
c take care of the rest.
c B-field values must be:
c a) put into double precision using dble() - done below
c b) converted from gauss to teslas - done here

         DO 30 I=1,61
            DO 20 J=1,61
               DO 10 K=1,121
                  READ(U,*) BUFF, BUFF, BUFF, BX(I,J,K),
     &             BY(I,J,K), BZ(I,J,K)
   10          CONTINUE
   20       CONTINUE
   30    CONTINUE

         CLOSE(U)

         BX = BX*1.0E-4
         BY = BY*1.0E-4
         BZ = BZ*1.0E-4

         LFIRST = .FALSE.
      ENDIF

c Here is where the routine is actually called, and Y/X/Z
c are nint()ed to find the nearest point, and given a
c B-field magnitude.

      XI = NINT(X)+31
      YI = NINT(Y)+31
      ZI = NINT(Z/5.0)+1

      B = DBLE(SQRT(BX(XI,YI,ZI)**2 + BY(XI,YI,ZI)**2 
     & + BZ(XI,YI,ZI)**2))
      BTX = BX/B
      BTY = BY/B
      BTZ = BZ/B

      RETURN
      END
