*$ 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-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      *
*                                                                      *
*----------------------------------------------------------------------*
*
*     Last change on 5-oct-10     by    Advanced FLUKA course teacher  *
*

*
      LOGICAL LFIRST
      SAVE LFIRST
      SAVE GRADONE, GRADTWO, NQUADRA
*
      DATA LFIRST / .TRUE. /
*
      IDISC = 0
      IF (LFIRST) THEN
*     gradient in tesla per meter for two quadropoles
         GRADONE = 5.0694D+00
         GRADTWO = -9.2982D+00
         CALL GEON2R("QuadrA ",NQUADRA,IERR)
         CALL GEON2R("QuadrAA ",NQUADRA,IERR)
         LFIRST  = .FALSE.
      END IF
*
      IF (NREG.EQ.NQUADRA) THEN 
         XNEW = X
         YNEW = Y
         ZNEW = Z
         
c     gradient in tesla per meter
         IF (NQUADRA.EQ.QuadrA) THEN
            BTX = GRADONE * YNEW/1D+02
            BTY = GRADONE * XNEW/1D+02
         ELSE 
            BTX = GRADTWO * YNEW/1D+02
            BTY = GRADTWO * XNEW/1D+02
         END IF 
         
         BTZ=ZERZER
         B=SQRT(BTX**2+BTY**2)
         IF (B.GT.1D-12) THEN
            BTX = BTX / B
            BTY = BTY / B
         END IF
            
*=== End of subroutine magfld =========================================*
      END
