*$ 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)'
      DOUBLE PRECISION BTX,BTY,BTZ,B,BR,BZ
      REAL F_MAGN(22,26),R
      INTEGER i1,i2,IR,IZ
      COMMON /MAG_FLD/F_MAGN
*
*  +-------------------------------------------------------------------*
*  |  Earth geomagnetic field:
      IF ( LGMFLD ) THEN
         CALL GEOFLD ( X, Y, Z, BTX, BTY, BTZ, B, NREG, IDISC )
         RETURN
      END IF
*  |
*  +-------------------------------------------------------------------*
      R=SQRT(X**2+Y**2)
      IF((ABS(Z).GE.250.0).OR.(R.GE.100.0))THEN
         BTX=0.0D0
	 BTY=0.0D0
	 BTZ=1.0D0
	 B=0.0001D0
      ELSE
	 IR=R/10.
	 IR=IR*10
	 IZ=ABS(Z)/10.
	 IZ=IZ*10
	 i1=IR/10*2+1
	 i2=IZ/10+1
	 BR=F_MAGN(i1,i2)+(F_MAGN(i1+2,i2)-F_MAGN(i1,i2))*
     &   (R-IR)/10.+(F_MAGN(i1,i2+1)-F_MAGN(i1,i2))*(ABS(Z)-IZ)/10. 
	 BZ=F_MAGN(i1+1,i2)+(F_MAGN(i1+3,i2)-F_MAGN(i1+1,i2))*
     &   (R-IR)/10.+(F_MAGN(i1+1,i2+1)-F_MAGN(i1+1,i2))*(ABS(Z)-IZ)/10.
         B = DSQRT(BR**2+BZ**2)
	 IF(B.LE.0.001)B = 0.001
	 BTZ = BZ/B
	 IF(R.LE.0.0001)THEN
	     BTX = 0.D0
	     BTY = 0.D0
	 ELSE
	     BTY = Y/R*BR
	     BTX = X/R*BR
	 ENDIF
         BTZ = DSQRT(2*BTZ**2+BTY**2+BTX**2-1.D0)
         IF(BZ.LE.0)BTZ=-BTZ
      ENDIF	  
*      print *,X,Y,Z,BTX,BTY,BTZ,B
      RETURN
*=== End of subroutine Magfld =========================================*
      END

