*$ CREATE MAGFLD.FOR
*COPY MAGFLD
*
*===magfld=============================================================*
*
      SUBROUTINE MAGFLD ( X, Y, Z, BTX, BTY, BTZ, B, NREG, IDISC )

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
* *
* Copyright (C) 1988-2005 by Alberto Fasso` & Alfredo Ferrari *
* All Rights Reserved. *
* *
* *
* Created in 1988 by Alberto Fasso`, CERN - TIS *
* *
* Last change on 11-dec-92 by Alfredo Ferrari *
* *
* Input variables: *
* x,y,z = current position *
* nreg = current region *
* Output variables: *
* btx,bty,btz = cosines of the magn. field vector *
* B = magnetic field intensity (Tesla) *
* idisc = set to 1 if the particle has to be discarded *
* *
*----------------------------------------------------------------------*
*
      PARAMETER (NX=20059) !TT: There are 20059 meshing nodes hence data points
      PARAMETER (NY=20059)
      PARAMETER (NZ=20059)
      DIMENSION BX(NX,NY,NZ)
      DIMENSION BY(NX,NY,NZ)
      DIMENSION BZ(NX,NY,NZ)
      DIMENSION XB(NX), YB(NY), ZB(NZ)

      LOGICAL LFIRST
      DATA LFIRST /.TRUE./
      SAVE LFIRST

      IDISC = 0

      IF (LFIRST) THEN
         
* Here, at the first time in which MAGFLD is called
* open and read a file containing a map of the 3 magnetic field components
* BX(I,J,K), BY(I,J,K), BZ(I,J,K)
* in a lattice of points in 3 dimensions XB(I), YB(J), ZB(K). There are
* NX points in X with running index I
* NY points in Y with running index J
* NZ points in Z with running index K
         OPEN (1, file = '3DMagFieldMap_Final.txt')
         DO i=1, NX !TT: loop over the entire data points
	 READ (1,*) BTX(i), BTY(i), BTZ(i)
	 END DO !TT
         CLOSE(1) !TT
         LFIRST = .FALSE.
      ENDIF

* Optional check on the Region at run time, to avoid calculation in
* regions where magnetic field is expected to be null..
**      IF(NREG .NE. ...) THEN
**         WRITE (LUNOUT,*)
**     & ' Magfld called in zone ',NREG,' where there should'
**         WRITE (LUNOUT,*)
**     & ' be no magnetic field. Something is wrong'
**         IDISC = 1
**         RETURN
**      ELSE
* Each time the routine is called, activate an
* algorithm of your choice that, given the particle coordinates X, Y
* and Z, finds the closest XB(), YB() and ZB() coordinates and
* interpolates the local value of BX, BY and BZ...
* Then:
         B = SQRT(BX**2 + BY**2 + BZ**2)
         BTX = X / B
         BTY = Y / B
         BTZ = Z / B
         RETURN
      ENDIF
*=== End of subroutine magfld =========================================*
      END 