*$ 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 '(RTDFCM)'
      INCLUDE '(LTCLCM)'
*
      DIMENSION IROTIN(10)
      LOGICAL LFIRST
      SAVE LFIRST
      SAVE  DIPFLD1, NPIPE2, NPIPE4, IROTIN
*
      DATA LFIRST / .TRUE. /, IROTIN / 10*0 /
*
      IDISC = 0

      IF (LFIRST) THEN 
        CALL GEON2R( 'pipenei2', NPIPE2, IERR)
        IF (IERR .EQ. 1) THEN
          CALL FLABRT ( 'GEON2R', 'CANT FIND pipenei2 REGION NUMBER' )
        END IF
        CALL GEON2R( 'pipenei4', NPIPE4, IERR)
        IF (IERR .EQ. 1) THEN
          CALL FLABRT ( 'GEON2R', 'CANT FIND pipenei4 REGION NUMBER' )
        END IF

        WRITE (LUNOUT,*) ' pipenei2 region is number: ', NPIPE2
        WRITE (LUNOUT,*) ' pipenei4 region is number: ', NPIPE4

        LFIRST  = .FALSE.
      END IF
      IF ( NREG .EQ. NPIPE2 ) THEN
        R = SQRT(X**2 + Y**2)
        BTX = -Y/R
        BTY = X/R
        BTZ = ZERZER
        B = 1.716d-2 * R
      ELSE IF ( NREG .EQ. NPIPE4 ) THEN
        R = SQRT(X**2 + Y**2)
        BTX = Y/R
        BTY = -X/R
        BTZ = ZERZER
        B = 1.716d-2 * R
      ELSE
        CALL FLABRT ( 'MAGFLD', 'NO MAGNETIC FIELD IN THIS REGION!' )
      END IF
      RETURN
*=== End of subroutine magfld =========================================*
      END

