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

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

      LOGICAL LFIRST
      SAVE LFIRST
      DATA LFIRST /.TRUE./
 
      INTEGER INSIDE

      DOUBLE PRECISION BBX, BBY, BBZ, BTT
      
      INTEGER Nx, Ny, Nz

      INTEGER I, J, K, INDX

      DOUBLE PRECISION Xmin, Xmax
      DOUBLE PRECISION Ymin, Ymax
      DOUBLE PRECISION Zmin, Zmax

      DOUBLE PRECISION XX, YY, ZZ

      DOUBLE PRECISION HemiX, HemiY, HemiZ

      DOUBLE PRECISION Xstep, Ystep, Zstep
      
      DOUBLE PRECISION Bx(binsX*binsY*binsZ)
      DOUBLE PRECISION By(binsX*binsY*binsZ)
      DOUBLE PRECISION Bz(binsX*binsY*binsZ)

      DOUBLE PRECISION INTERPOLFIELD

      DOUBLE PRECISION HA(2,2,2)
      DOUBLE PRECISION HB(2,2)
      DOUBLE PRECISION HC(2)  

      COMMON /COORDS/ Xmin, Xmax, Ymin, Ymax, Zmin, Zmax, INSIDE
      COMMON /GRID/ Nx, Ny, Nz
      COMMON /STEPS/ Xstep, Ystep, Zstep
      COMMON /BS/ Bx, By, Bz
      COMMON /HS/ HA, HB, HC, INTERPOLFIELD
      COMMON /HEMI/ HemiX, HemiY, HemiZ

* return message from first call
      IF (LFIRST) THEN
         WRITE(LUNOUT,*) '*************************'
         WRITE(LUNOUT,*) '*************************'
         WRITE(LUNOUT,*) 'Version "Standard" of Magnet Feld'
         WRITE(LUNOUT,*) '*************************'
         WRITE(LUNOUT,*) '*************************'
         WRITE(LUNOUT,*) '*************************'
	 CALL SETFIELD         
*         OPEN(UNIT=10,FILE='field_target.txt',STATUS='NEW')      
*         DO I=1,Nx
*            XX = Xmin + (0.5 + I-1)*Xstep
*            DO J=1,Ny
*               YY = Ymin + (0.5 + J-1)*Ystep
*               DO K=1,Nz
*                  ZZ = Zmin + (0.5 + K-1)*Zstep - 40.
*                  INDX = (I-1)*Ny*Nz + (J-1)*Nz + (K-1) + 1
*                  BBX = Bx(INDX)
*                  BBY = By(INDX)
*                  BBZ = Bz(INDX)
*                  BTT = sqrt(BBX*BBX + BBY*BBY + BBZ*BBZ)
*                  WRITE(LUNOUT,*) ' Z= ', ZZ, ' B= ', BTT
*
*                  IF (ZZ .EQ. -1) THEN
*                     WRITE(10, *) XX, YY, ZZ, BTT
*                  ENDIF
*
*               ENDDO
*            ENDDO
*         ENDDO
*         CLOSE(10)
         LFIRST = .FALSE.
      ENDIF


*      WRITE(LUNOUT,*) 'Xmin = ', Xmin, ', Xmax = ', Xmax
*      WRITE(LUNOUT,*) 'Ymin = ', Ymin, ', Ymax = ', Ymax
*      WRITE(LUNOUT,*) 'Zmin = ', Zmin, ', Zmax = ', Zmax

*      WRITE(LUNOUT,*) 'TARGSHIFT = ', TARGSHIFT

*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1988-2009      by Alberto Fasso` & Alfredo Ferrari *
*     All Rights Reserved.                                             *
*                                                                      *
*                                                                      *
*     Created  in 1988    by     Alberto Fasso`, CERN - TIS            *
*                                                                      *
*     Last change on 15-oct-09     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      *
*                                                                      *
*----------------------------------------------------------------------*
*
*      IDISC = 0
* 2 Tesla uniform field along +z:
*      BTX   = ZERZER
*      BTY   = ZERZER
*      BTZ   = ONEONE
*      B     = TWOTWO

      INSIDE = 0

      CALL GETBX ( X, Y, Z )

      IF (INSIDE .EQ. 1) THEN

      BBX = INTERPOLFIELD * HemiX * HemiY

      CALL GETBY ( X, Y, Z )
      BBY = INTERPOLFIELD

      CALL GETBZ ( X, Y, Z )
      BBZ = INTERPOLFIELD * HemiY * HemiZ

      B = sqrt(BBX*BBX + BBY*BBY + BBZ*BBZ)

      BTX = BBX / B
      BTY = BBY / B
      BTZ = BBZ / B
  
      ELSE

      BTX = 0.
      BTY = 0.
      BTZ = 0.
      B = 0.

      ENDIF

      RETURN
*=== End of subroutine magfld =========================================*
      END


*=== SUBROUTINES & FUNCTIONS =========================================*

********************** SETFIELD

      SUBROUTINE SETFIELD

      INCLUDE 'COMMON'

      INTEGER I, J, K, INDX, INSIDE
      DOUBLE PRECISION b1, b2, b3

      INTEGER Nx, Ny, Nz

      DOUBLE PRECISION Xmin, Xmax
      DOUBLE PRECISION Ymin, Ymax
      DOUBLE PRECISION Zmin, Zmax

      DOUBLE PRECISION Xstep, Ystep, Zstep

      DOUBLE PRECISION Bx(binsX*binsY*binsZ)
      DOUBLE PRECISION By(binsX*binsY*binsZ)
      DOUBLE PRECISION Bz(binsX*binsY*binsZ) 

      COMMON /COORDS/ Xmin, Xmax, Ymin, Ymax, Zmin, Zmax, INSIDE
      COMMON /GRID/ Nx, Ny, Nz
      COMMON /STEPS/ Xstep, Ystep, Zstep
      COMMON /BS/ Bx, By, Bz

      DOUBLE PRECISION factor
      factor = FIELDSCALE

      OPEN(20,FILE='/lustre/cbm/users/anna/field',STATUS='OLD')

      READ(20,*) Xmin, Xmax, Nx 
      READ(20,*) Ymin, Ymax, Ny 
      READ(20,*) Zmin, Zmax, Nz  

*      WRITE(LUNOUT,*) '*************************'
*      WRITE(LUNOUT,*) '*************************'
*      WRITE(LUNOUT,*) '*************************'
*      WRITE(LUNOUT,*) 'Xmin = ', Xmin, ', Xmax = ', Xmax
*      WRITE(LUNOUT,*) 'Ymin = ', Ymin, ', Ymax = ', Ymax
*      WRITE(LUNOUT,*) 'Zmin = ', Zmin, ', Zmax = ', Zmax
*      WRITE(LUNOUT,*) '*************************'
*      WRITE(LUNOUT,*) '*************************'
*      WRITE(LUNOUT,*) '*************************'

      Xstep = ( Xmax - Xmin ) / ( Nx - 1. )
      Ystep = ( Ymax - Ymin ) / ( Ny - 1. )
      Zstep = ( Zmax - Zmin ) / ( Nz - 1. )

      DO I=1,Nx
        DO J=1,Ny
	  DO K=1,Nz
	    READ(20,*) b1, b2, b3
*  (+ 1 - FORTRAN array)
	    INDX = (I-1)*Ny*Nz + (J-1)*Nz + (K-1) + 1
	    Bx(INDX) = b1*factor
	    By(INDX) = b2*factor
	    Bz(INDX) = b3*factor
*	    WRITE(LUNOUT,*) 'BX = ', Bx(INDX), ', BY = ', By(INDX), ', BZ = ', Bz(INDX)
	  ENDDO
        ENDDO
      ENDDO

      WRITE(LUNOUT,*) 'Magnet field is full ', INDX 
      WRITE(LUNOUT,*) '************ FIELDSCALE = ', factor
      WRITE(LUNOUT,*) '*************************'
      WRITE(LUNOUT,*) 'Version "Standard" of Magnet Feld'      
      WRITE(LUNOUT,*) '*************************'
      WRITE(LUNOUT,*) '*************************'
     
      CLOSE(20)

      RETURN
      END

********************** GETBX

      SUBROUTINE GETBX(X, Y, Z) 

      INCLUDE 'COMMON'

      INTEGER IX, IY, IZ
      DOUBLE PRECISION DX, DY, DZ
      DOUBLE PRECISION X, Y, Z

      INTEGER Nx, Ny, Nz

      DOUBLE PRECISION HemiX, HemiY, HemiZ

      DOUBLE PRECISION INTERPOLFIELD

      DOUBLE PRECISION Bx(binsX*binsY*binsZ)
      DOUBLE PRECISION By(binsX*binsY*binsZ)
      DOUBLE PRECISION Bz(binsX*binsY*binsZ)

      DOUBLE PRECISION HA(2,2,2)
      DOUBLE PRECISION HB(2,2)
      DOUBLE PRECISION HC(2)  

      COMMON /GRID/ Nx, Ny, Nz
      COMMON /BS/ Bx, By, Bz
      COMMON /HS/ HA, HB, HC, INTERPOLFIELD
      COMMON /HEMI/ HemiX, HemiY, HemiZ

      IX = 0
      IY = 0
      IZ = 0
      DX = 0.
      DY = 0.
      DZ = 0.

      IF ( ISINSIDE ( X, Y, Z, IX, IY, IZ, DX, DY, DZ ).EQ.1 ) THEN

* Get Bx field values at grid cell corners (+ 1 - FORTRAN array)
	HA(1, 1, 1) = Bx(IX    *Ny*Nz + IY    *Nz + IZ  + 1 )
	HA(2, 1, 1) = Bx((IX+1)*Ny*Nz + IY    *Nz + IZ  + 1 )
	HA(1, 2, 1) = Bx(IX    *Ny*Nz + (IY+1)*Nz + IZ  + 1 )
	HA(2, 2, 1) = Bx((IX+1)*Ny*Nz + (IY+1)*Nz + IZ  + 1 )
	HA(1, 1, 2) = Bx(IX    *Ny*Nz + IY    *Nz + (IZ+1)  + 1 )
	HA(2, 1, 2) = Bx((IX+1)*Ny*Nz + IY    *Nz + (IZ+1)  + 1 )
	HA(1, 2, 2) = Bx(IX    *Ny*Nz + (IY+1)*Nz + (IZ+1)  + 1 )
	HA(2, 2, 2) = Bx((IX+1)*Ny*Nz + (IY+1)*Nz + (IZ+1)  + 1 )

	CALL INTERPOLATE ( DX, DY, DZ )
* Return interpolated field value
	RETURN 
      ENDIF

      RETURN 
      END 

********************** GETBY

      SUBROUTINE GETBY(X, Y, Z) 

      INCLUDE 'COMMON'

      INTEGER IX, IY, IZ
      DOUBLE PRECISION DX, DY, DZ
      DOUBLE PRECISION X, Y, Z

      INTEGER Nx, Ny, Nz

      DOUBLE PRECISION INTERPOLFIELD

      DOUBLE PRECISION Bx(binsX*binsY*binsZ)
      DOUBLE PRECISION By(binsX*binsY*binsZ)
      DOUBLE PRECISION Bz(binsX*binsY*binsZ)

      DOUBLE PRECISION HA(2,2,2)
      DOUBLE PRECISION HB(2,2)
      DOUBLE PRECISION HC(2)  

      COMMON /GRID/ Nx, Ny, Nz
      COMMON /BS/ Bx, By, Bz
      COMMON /HS/ HA, HB, HC, INTERPOLFIELD

      IX = 0
      IY = 0
      IZ = 0
      DX = 0.
      DY = 0.
      DZ = 0.

      IF ( ISINSIDE ( X, Y, Z, IX, IY, IZ, DX, DY, DZ ).EQ.1 ) THEN

* Get By field values at grid cell corners (+ 1 - FORTRAN array)
	HA(1, 1, 1) = By(IX    *Ny*Nz + IY    *Nz + IZ  + 1 )
	HA(2, 1, 1) = By((IX+1)*Ny*Nz + IY    *Nz + IZ  + 1 )
	HA(1, 2, 1) = By(IX    *Ny*Nz + (IY+1)*Nz + IZ  + 1 )
	HA(2, 2, 1) = By((IX+1)*Ny*Nz + (IY+1)*Nz + IZ  + 1 )
	HA(1, 1, 2) = By(IX    *Ny*Nz + IY    *Nz + (IZ+1)  + 1 )
	HA(2, 1, 2) = By((IX+1)*Ny*Nz + IY    *Nz + (IZ+1)  + 1 )
	HA(1, 2, 2) = By(IX    *Ny*Nz + (IY+1)*Nz + (IZ+1)  + 1 )
	HA(2, 2, 2) = By((IX+1)*Ny*Nz + (IY+1)*Nz + (IZ+1)  + 1 )
*        WRITE(LUNOUT,*) '*************************'
*        WRITE(LUNOUT,*) 'HA1 = ', HA(1, 1, 1), ' HA2 = ', HA(2, 1, 1)
*        WRITE(LUNOUT,*) 'HA3 = ', HA(1, 2, 1), ' HA4 = ', HA(2, 2, 1)
*        WRITE(LUNOUT,*) 'HA5 = ', HA(1, 1, 2), ' HA6 = ', HA(2, 1, 2)
*        WRITE(LUNOUT,*) 'HA7 = ', HA(1, 2, 2), ' HA8 = ', HA(2, 2, 2)
*        WRITE(LUNOUT,*) '*************************'

	CALL INTERPOLATE ( DX, DY, DZ )
* Return interpolated field value
*        WRITE(LUNOUT,*) '*************************'*
*	 WRITE(LUNOUT,*) 'BY = ', INTERPOLFIELD
*        WRITE(LUNOUT,*) '*************************'
	RETURN 
      ENDIF

      RETURN 
      END 

********************** GETBZ

      SUBROUTINE GETBZ(X, Y, Z) 

      INCLUDE 'COMMON'

      INTEGER IX, IY, IZ
      DOUBLE PRECISION DX, DY, DZ
      DOUBLE PRECISION X, Y, Z

      INTEGER Nx, Ny, Nz

      DOUBLE PRECISION HemiX, HemiY, HemiZ

      DOUBLE PRECISION INTERPOLFIELD

      DOUBLE PRECISION Bx(binsX*binsY*binsZ)
      DOUBLE PRECISION By(binsX*binsY*binsZ)
      DOUBLE PRECISION Bz(binsX*binsY*binsZ)

      DOUBLE PRECISION HA(2,2,2)
      DOUBLE PRECISION HB(2,2)
      DOUBLE PRECISION HC(2)  

      COMMON /GRID/ Nx, Ny, Nz
      COMMON /BS/ Bx, By, Bz
      COMMON /HS/ HA, HB, HC, INTERPOLFIELD
      COMMON /HEMI/ HemiX, HemiY, HemiZ

      IX = 0
      IY = 0
      IZ = 0
      DX = 0.
      DY = 0.
      DZ = 0.

      IF ( ISINSIDE ( X, Y, Z, IX, IY, IZ, DX, DY, DZ ).EQ.1 ) THEN

* Get Bz field values at grid cell corners (+ 1 - FORTRAN array)
	HA(1, 1, 1) = Bz(IX    *Ny*Nz + IY    *Nz + IZ  + 1 )
	HA(2, 1, 1) = Bz((IX+1)*Ny*Nz + IY    *Nz + IZ  + 1 )
	HA(1, 2, 1) = Bz(IX    *Ny*Nz + (IY+1)*Nz + IZ  + 1 )
	HA(2, 2, 1) = Bz((IX+1)*Ny*Nz + (IY+1)*Nz + IZ  + 1 )
	HA(1, 1, 2) = Bz(IX    *Ny*Nz + IY    *Nz + (IZ+1)  + 1 )
	HA(2, 1, 2) = Bz((IX+1)*Ny*Nz + IY    *Nz + (IZ+1)  + 1 )
	HA(1, 2, 2) = Bz(IX    *Ny*Nz + (IY+1)*Nz + (IZ+1)  + 1 )
	HA(2, 2, 2) = Bz((IX+1)*Ny*Nz + (IY+1)*Nz + (IZ+1)  + 1 )

	CALL INTERPOLATE ( DX, DY, DZ )
* Return interpolated field value
	RETURN
      ENDIF

      RETURN 
      END 

********************** ISINSIDE

      FUNCTION ISINSIDE ( X, Y, Z, IX, IY, IZ, DX, DY, DZ )

      INCLUDE 'COMMON'

      INTEGER IX, IY, IZ, INSIDE
      DOUBLE PRECISION DX, DY, DZ
      DOUBLE PRECISION X, Y, Z

      DOUBLE PRECISION Xmin, Xmax
      DOUBLE PRECISION Ymin, Ymax
      DOUBLE PRECISION Zmin, Zmax

      DOUBLE PRECISION HemiX, HemiY, HemiZ

      DOUBLE PRECISION Xstep, Ystep, Zstep  

      COMMON /COORDS/ Xmin, Xmax, Ymin, Ymax, Zmin, Zmax, INSIDE
      COMMON /STEPS/ Xstep, Ystep, Zstep
      COMMON /HEMI/ HemiX, HemiY, HemiZ
 
* Local coordinate system
      DOUBLE PRECISION XL, YL, ZL

*      WRITE(LUNOUT,*) '*************************'
*      WRITE(LUNOUT,*) 'X = ', X, ', Y = ', Y, ', Z = ', Z
*      WRITE(LUNOUT,*) '*************************'

      INSIDE = 0

* Transform into local coordinate system
      XL = X - PosX
      YL = Y - PosY
      ZL = Z - PosZ

*  Reflect w.r.t. symmetry axes
      HemiX = 1.
      HemiY = 1.
      HemiZ = 1.

      IF ( XL .LT. 0. ) THEN
	HemiX = -1.
	XL    = -1. * XL
      ENDIF

      IF ( YL .LT. 0. ) THEN
	HemiY = -1.
	YL    = -1. * YL
      ENDIF

      IF ( sym .EQ. 3 ) THEN      
	IF ( ZL .LT. 0. ) THEN
	  HemiZ = -1.
	  ZL    = -1. * ZL
	ENDIF
      ENDIF

* Check for being outside the map range
      IF ((XL.LT.Xmin .OR. XL.GT.Xmax) 
     & .OR. (YL.LT.Ymin .OR. YL.GT.Ymax) 
     & .OR. (ZL.LT.Zmin .OR. ZL.GT.Zmax)) THEN
	IX = 0
	IY = 0
	IZ = 0
	DX = 0.
	DY = 0.
	DZ = 0.
	ISINSIDE = 0
      INSIDE = 0
      ELSE
* Determine grid cell
	IX = INT( (XL-Xmin) / Xstep )
	IY = INT( (YL-Ymin) / Ystep )
	IZ = INT( (ZL-Zmin) / Zstep )	
* Relative distance from grid point (in units of cell size)
	DX = (XL-Xmin) / Xstep - FLOAT(IX)
	DY = (YL-Ymin) / Ystep - FLOAT(IY)
	DZ = (ZL-Zmin) / Zstep - FLOAT(IZ)
	ISINSIDE = 1
      INSIDE = 1
      ENDIF

      RETURN 
      END

********************** INTERPOLATE

      SUBROUTINE INTERPOLATE (DX, DY, DZ) 

      INCLUDE 'COMMON'

      DOUBLE PRECISION DX, DY, DZ

      DOUBLE PRECISION INTERPOLFIELD

      DOUBLE PRECISION HA(2,2,2)
      DOUBLE PRECISION HB(2,2)
      DOUBLE PRECISION HC(2)  

      COMMON /HS/ HA, HB, HC, INTERPOLFIELD


*      WRITE(LUNOUT,*) '*************************'
*      WRITE(LUNOUT,*) 'DX = ', DX, ' DY = ', DY, ' DZ = ', DZ
*      WRITE(LUNOUT,*) '*************************'

* Interpolate in x coordinate
      HB(1, 1) = HA(1, 1, 1) + ( HA(2, 1, 1)-HA(1, 1, 1) ) * DX
      HB(2, 1) = HA(1, 2, 1) + ( HA(2, 2, 1)-HA(1, 2, 1) ) * DX
      HB(1, 2) = HA(1, 1, 2) + ( HA(2, 1, 2)-HA(1, 1, 2) ) * DX
      HB(2, 2) = HA(1, 2, 2) + ( HA(2, 2, 2)-HA(1, 2, 2) ) * DX

*      WRITE(LUNOUT,*) '*************************'
*      WRITE(LUNOUT,*) 'HB1 = ', HB(1, 1), ' HB2 = ', HB(2, 1)
*      WRITE(LUNOUT,*) 'HB3 = ', HB(1, 2), ' HB4 = ', HB(2, 2)
*      WRITE(LUNOUT,*) '*************************'

* Interpolate in y coordinate
      HC(1) = HB(1, 1) + ( HB(2, 1) - HB(1, 1) ) * DY
      HC(2) = HB(1, 2) + ( HB(2, 2) - HB(1, 2) ) * DY
*      WRITE(LUNOUT,*) '*************************'
*      WRITE(LUNOUT,*) 'HC1 = ', HC(1), ' HC2 = ', HC(2)
*      WRITE(LUNOUT,*) '*************************'

* Interpolate in z coordinate
      INTERPOLFIELD = HC(1) + ( HC(2) - HC(1) ) * DZ
*      WRITE(LUNOUT,*) '*************************'
*      WRITE(LUNOUT,*) 'INTERPOLATE = ', INTERPOLFIELD
*      WRITE(LUNOUT,*) '*************************'
      RETURN 
      END

*=== END =========================================*
