*$ 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      *
*                                                                      *
*----------------------------------------------------------------------*
*
*     	lfirst, to flag first call of routine
		LOGICAL LFIRST
		DATA    LFIRST / .TRUE. /
		SAVE    LFIRST
*     magnetic regions:
         SAVE tar
	    IF ( LFIRST ) THEN		 
          CALL GEON2R( "TARGET", tar, IERR )
         WRITE (LUNOUT,*) ' "TARGET" region is number: ',tar
*        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
*  +-------------------------------------------------------------------*
*  |  Earth geomagnetic field:
*      IF ( LGMFLD ) THEN
*         CALL GEOFLD ( X, Y, Z, BTX, BTY, BTZ, B, NREG, IDISC )
*         RETURN
*      END IF
*  |
*  +-------------------------------------------------------------------*
      IDISC = 0
* 2 Tesla uniform field along +z:
		IF ( NREG .EQ. tar) THEN
		R0=0
		RE=1
		B0=10
	   R =SQRT(X**2 + Y**2) 
       B =B0 * EXP(-ABS(R - R0)/RE) 	
       BTX =ZERZER 
       BTY = ONEONE 
       BTZ = ZERZER 	
       END IF
	   RETURN
*=== End of subroutine Magfld =========================================*
      END

