*$ 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-2010      by Alberto Fasso` & Alfredo Ferrari *
*     All Rights Reserved.                                             *
*                                                                      *
*                                                                      *
*     Created  in     1988         by    Alberto Fasso`                *
*                                                                      *
*                                                                      *
*     Last change on 06-Nov-10     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      *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(CMEMFL)'
      INCLUDE '(CSMCRY)'
*
*     Variables storing region number of the magnetic regions:
      SAVE NTURN1, NTURN2, NTURN3, NTURN4, NTURN5
*
*     lfirst, to flag first call of routine
      LOGICAL LFIRST
      DATA    LFIRST / .TRUE. /
      SAVE    LFIRST
*
*  +-------------------------------------------------------------------*
*	Initialize Core Routine                                          *
*  +-------------------------------------------------------------------*
      IDISC = 0
*
* Default vaules - No magnetic field
      BTX = ZERZER
      BTY = ONEONE
      BTZ = ZERZER
      B   = ZERZER

*----------------------------------------------------------------------*
*     first call: initialisation
*----------------------------------------------------------------------*
*
      IF ( LFIRST ) THEN
*
*        retrieve the indices of the magnetic regions starting from their
*        names, and store them
         CALL GEON2R( "TurnVac1", NTURN1, IERR )
         CALL GEON2R( "TurnVac2", NTURN2, IERR )
         CALL GEON2R( "BMvac8  ", NTURN3, IERR )
         CALL GEON2R( "BMvac1  ", NTURN4, IERR )
         CALL GEON2R( "BMvac2  ", NTURN5, IERR )
         CALL GEON2R( "BMvac3  ", NTURN6, IERR )
         CALL GEON2R( "BMvac4  ", NTURN7, IERR )
         CALL GEON2R( "BMvac5  ", NTURN8, IERR )
         CALL GEON2R( "BMvac6  ", NTURN9, IERR )
         CALL GEON2R( "BMvac7  ", NTURN10, IERR )
         CALL GEON2R( "InjB3Vac", NTURN11, IERR )
         CALL GEON2R( "InjB4Vac", NTURN12, IERR )
         CALL GEON2R( "ThinVac ", NTURN13, IERR )

*
*         avoid to initialise again the routine:
          LFIRST  = .FALSE.
      END IF

* 0.793388 T fields in 2 regions TurnVac1 BTY = 1.0, TurnVac2 BTY = -1.0  *
      IF ( NREG .EQ. NTURN1 ) THEN
*        First linac turn, 0.793388 T in positive y-axis
         BTX = ZERZER
         BTY = ONEONE
         BTZ = ZERZER
         B   = 0.793388
      ELSEIF ( NREG .EQ. NTURN2 ) THEN
*        Second linac turn, 0.793388 T in negative y-axis
         BTX = ZERZER
         BTY = -ONEONE
         BTZ = ZERZER
         B   = 0.793388
      ELSEIF ( NREG .EQ. NTURN3 ) THEN
*        Bending Magnet 8 turn, 0.20455 T in negative x-axis
         BTX = -ONEONE
         BTY = ZERZER
         BTZ = ZERZER
         B   = 0.20455
      ELSEIF ( NREG .EQ. NTURN4 ) THEN
*        Bending Magnet 1 turn, 0.2063 T in negative x-axis
         BTX = -ONEONE
         BTY = ZERZER
         BTZ = ZERZER
         B   = 0.2063
      ELSEIF ( NREG .EQ. NTURN5 ) THEN
*        Bending Magnet 2 turn, 0.197 T in negative x-axis
         BTX = -ONEONE
         BTY = ZERZER
         BTZ = ZERZER
         B   = 0.197	
      ELSEIF ( NREG .EQ. NTURN6 ) THEN
*        Bending Magnet 3 turn, 0.213 T in negative x-axis
         BTX = -ONEONE
         BTY = ZERZER
         BTZ = ZERZER
         B   = 0.213
      ELSEIF ( NREG .EQ. NTURN7 ) THEN
*        Bending Magnet 4 turn, 0.204 T in negative x-axis
         BTX = -ONEONE
         BTY = ZERZER
         BTZ = ZERZER
         B   = 0.204
      ELSEIF ( NREG .EQ. NTURN8 ) THEN
*        Bending Magnet 5 turn, 0.2043 T in negative x-axis
         BTX = -ONEONE
         BTY = ZERZER
         BTZ = ZERZER
         B   = 0.2043
      ELSEIF ( NREG .EQ. NTURN9 ) THEN
*        Bending Magnet 6 turn, 0.2108 T in negative x-axis
         BTX = -ONEONE
         BTY = ZERZER
         BTZ = ZERZER
         B   = 0.2108
      ELSEIF ( NREG .EQ. NTURN10 ) THEN
*        Bending Magnet 7 turn, 0.19889 T in negative x-axis
         BTX = -ONEONE
         BTY = ZERZER
         BTZ = ZERZER
         B   = 0.19889
      ELSEIF ( NREG .EQ. NTURN11 ) THEN
*        Injection Magnet B3, 1.128 T in negative x-axis
         BTX = -ONEONE
         BTY = ZERZER
         BTZ = ZERZER
         B   = 1.128
      ELSEIF ( NREG .EQ. NTURN12 ) THEN
*        Injection Magnet B4, 0.193 T in negative x-axis
         BTX = -ONEONE
         BTY = ZERZER
         BTZ = ZERZER
         B   = 0.193
      ELSEIF ( NREG .EQ. NTURN13 ) THEN
*        Thin Septum, 0.19889 T in negative x-axis
         BTX = -ONEONE
         BTY = ZERZER
         BTZ = ZERZER
         B   = 17.8
      ENDIF

      RETURN
*=== End of subroutine Magfld =========================================*
      END
