*$ 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 12-sep-12    by    Advanced FLUKA course teacher  *
*
      INCLUDE '(RTDFCM)'
      INCLUDE '(LTCLCM)'

*     rotation indices
      DIMENSION IROTIN(10)
      DATA      IROTIN / 10*0 /
      SAVE      IROTIN

*     lfirst, to flag first call of routine
      LOGICAL LFIRST
      DATA    LFIRST / .TRUE. /
      SAVE    LFIRST

*     magnetic data:
*     - dipole field (T):
      DATA DIPFLD / 1.66782D+00 /
      SAVE DIPFLD
*     - quadrupole field (T/m):
      DATA QUAFLD / 4.9D+00 /
      SAVE QUAFLD

*     magnetic regions:
      SAVE NUMQUA, NUMDIP

*     default values
      IDISC = 0
      BTX = ZERZER
      BTY = ZERZER
      BTZ = ONEONE
      B   = ZERZER

*----------------------------------------------------------------------*
*     first call: initialisation
*----------------------------------------------------------------------*

      IF ( LFIRST ) THEN

*        some useful printout
         WRITE (LUNOUT,*) ''
         WRITE (LUNOUT,*) ' dedicated magfld.f routine for FLUKA course'
         WRITE (LUNOUT,*) ' initialisation started'
         WRITE (LUNOUT,*) ''
         WRITE (LUNOUT,*) ''
         WRITE (LUNOUT,*) ' dipole field: ',DIPFLD,' T;'
         WRITE (LUNOUT,*) ' quad field:   ',QUAFLD,' T/m;'
         WRITE (LUNOUT,*) ''

*        retrieve the indices of the magnetic regions starting from their
*            names, and store them
         CALL GEON2R( "DIP_BPVC", NUMDIP, IERR )
         CALL GEON2R( "QUADBPVC", NUMQUA, IERR )
         WRITE (LUNOUT,*) ''
         WRITE (LUNOUT,*) ' "DIP_BPVC" region is number: ',NUMDIP
         WRITE (LUNOUT,*) ' "QUADBPVC" region is number: ',NUMQUA
         WRITE (LUNOUT,*) ''

*        retrieve the index of the lattices and the associated rotations,
*        and store them in an array indexed with the lattice index

         CALL GEON2L("QUAD___1", NLATT, IRTLAT, IERR)
         IF ( IERR .GT. 0 ) THEN
            CALL FLABRT('MAGFLD','no "QUAD___1" lattice')
         ENDIF
         IROTIN( NLATT ) = IRTLAT
         WRITE (LUNOUT,*) ' "QUAD___1" is lattice number: ',NLATT,
     &'; associated transformation index: ',IROTIN( NLATT )

         CALL GEON2L("QUAD___2", NLATT, IRTLAT, IERR)
         IF ( IERR .GT. 0 ) THEN
            CALL FLABRT('MAGFLD','no "QUAD___2" lattice')
         ENDIF
         IROTIN( NLATT ) = IRTLAT
         WRITE (LUNOUT,*) ' "QUAD___2" is lattice number: ',NLATT,
     &'; associated transformation index: ',IROTIN( NLATT )

         CALL GEON2L("QUAD___3", NLATT, IRTLAT, IERR)
         IF ( IERR .GT. 0 ) THEN
            CALL FLABRT('MAGFLD','no "QUAD___3" lattice')
         ENDIF
         IROTIN( NLATT ) = IRTLAT
         WRITE (LUNOUT,*) ' "QUAD___3" is lattice number: ',NLATT,
     &'; associated transformation index: ',IROTIN( NLATT )

         CALL GEON2L("QUAD___4", NLATT, IRTLAT, IERR)
         IF ( IERR .GT. 0 ) THEN
            CALL FLABRT('MAGFLD','no "QUAD___4" lattice')
         ENDIF
         IROTIN( NLATT ) = IRTLAT
         WRITE (LUNOUT,*) ' "QUAD___4" is lattice number: ',NLATT,
     &'; associated transformation index: ',IROTIN( NLATT )
         WRITE (LUNOUT,*) ''

*        flush the .out file, in order to free the buffer:
         CALL FFLUSH
         CALL FLUSH(LUNOUT)

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

      END IF

*----------------------------------------------------------------------*
*     core of the routine
*----------------------------------------------------------------------*

*     if in a lattice, go to model:
      XNEW = X
      YNEW = Y
      ZNEW = Z
      IF ( MLATTC .GT. 0 ) THEN
*        update XNEW, YNEW, ZNEW with the correct transformation
         CALL DOTRSF ( 1, XNEW, YNEW, ZNEW, IROTIN(MLATTC) )
      ENDIF

*     disentangle the magnetic field type on the basis of the current 
*        region, and compute field components
      IF ( NREG .EQ. NUMDIP ) THEN
*        dipole field: conventionally along the y-axis, upward oriented,
*            no matter if the particle is in lattice or not, if there is
*            an associated transformation or not
         BTX = ZERZER
         BTY = DIPFLD
         BTZ = ZERZER
      ELSEIF ( NREG .EQ. NUMQUA ) THEN
*        quadrupole field:
*        the analytical expression is a function of the transverse position
*            of the particle with respect to the magnetic centre
         BTX = QUAFLD * YNEW/1D+02
         BTY = QUAFLD * XNEW/1D+02
         BTZ = ZERZER
      ELSE
*        no other region is magnetic for the moment
         CALL FLABRT('MAGFLD','NO MAGNETIC FIELD IN THIS REGION !')
      ENDIF

*     get field module:
      B   = SQRT( BTX**2 + BTY**2 + BTZ**2 )

*     check against numerical precision
      IF ( B .GT. 1.0D-12 ) THEN
*        normalise field components
         BTX = BTX / B
         BTY = BTY / B
         BTZ = BTZ / B
      ELSE
*        restore default values
         BTX = ZERZER
         BTY = ZERZER
         BTZ = ONEONE
         B   = ZERZER
      END IF

*     if in a lattice, transform back the field components in the global
*        FLUKA reference system
      IF ( MLATTC .GT. 0 ) THEN
         CALL UNDRTO ( 1, BTX, BTY, BTZ, IROTIN(MLATTC) )
      ENDIF

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