*$ 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      *
*                                                                      *
*----------------------------------------------------------------------*
*
*     EDIT BY HVS TRYING TO IMPLEMENT DIP AND QUAD MAG FIELDS **********

      CHARACTER*8 REGNAM
      
      INCLUDE '(CMEMFL)'
      INCLUDE '(CSMCRY)'
      INCLUDE '(RTDFCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(BEAMCM)'

      LOGICAL LFIRST
      DATA    LFIRST / .TRUE. /
      SAVE    LFIRST

*     MAGNETIC DATA:
*     DIPOLE FIELD (70 MeV = 0.176 TESLA, 800 MeV = 0.697)
*      DATA DIPFLD / 0.176D+00 /
*      SAVE DIPFLD
      
*     DOUBLET QUAD FIELD (K1 = +/- 0.6377)
*      DATA QD / 0.6377D+00 /
*      SAVE QD
      
*     TRIM QUAD FIELD (K1 = +/-0.103)
*      DATA TQ / 0.103D+00 /
*      SAVE TQ
      
*     SINGLET QUAD FIELD (K1 = -0.7262)
*      DATA QS / -0.7262D+00 /
*      SAVE QS


*     DEFAULT VALUES
      IDISC = 0
      BTX = ZERZER
      BTY = ZERZER
      BTZ = ONEONE
      B   = ZERZER

*     1 GeV/c = 5.344286D-19 kg*m/s
*     ELECTRONIC CHARGE = 1.60217657D-19 C
      GEVC = 5.344286D-19
      QE = 1.60217657D-19
*     BRHO = p/e = (momentum/charge) and conver to Tcm rather than Tm
      BRHO = (PBEAM*GEVC/QE)*100

*     BRHO = 1.0D-09*PBEAM*CLIGHT
      
*     FIRST CALL - INITIALISATION
      IF (LFIRST) THEN
         WRITE (LUNOUT,*)''
         WRITE (LUNOUT,*)'DEDICATED MAGFLD.F ROUTINE'
         WRITE (LUNOUT,*)'INITIALISATION STARTED'
         WRITE (LUNOUT,*)''
         WRITE (LUNOUT,*)'SPEED OF LIGHT: ', CLIGHT, 'cm/s'
         WRITE (LUNOUT,*)'BEAM MOMENTUM: ', PBEAM, 'GeV/c'
         WRITE (LUNOUT,*)'BRHO: ',BRHO,'Tcm'
         WRITE (LUNOUT,*)''
      LFIRST = .FALSE.
      ENDIF
      
*     ASSIGN CORRECT MAGNETIC FIELD TO EACH REGION
      CALL GEOR2N(NREG, REGNAM, IERR)

      DIPFLD = 0.176
*     DIPFLD = 0.697
      QK1 = 0.6377
      TQK1 = 0.103
      QSK1 = -0.7262
           
      IF (REGNAM .EQ. 'DIPVAC') THEN
         BTX = ZERZER
         BTY = DIPFLD
         BTZ = ZERZER
      ELSEIF (REGNAM .EQ. 'QDVAC') THEN       
         BTX = 1.0D-04*QK1*Y*BRHO
         BTY = 1.0D-04*QK1*X*BRHO
         BTZ = ZERZER
      ELSEIF (REGNAM .EQ. 'QFVAC') THEN
         BTX = 1.0D-04*(-1*QK1)*Y*BRHO
         BTY = 1.0D-04*(-1*QK1)*X*BRHO
         BTZ = ZERZER
      ELSEIF (REGNAM .EQ. 'QTDVAC') THEN
         BTX = 1.0D-04*TQK1*Y*BRHO
         BTY = 1.0D-04*TQK1*X*BRHO
         BTZ = ZERZER
      ELSEIF (REGNAM .EQ. 'QTFVAC') THEN
         BTX = 1.0D-04*(-1*TQK1)*Y*BRHO
         BTY = 1.0D-04*(-1*TQK1)*X*BRHO
         BTZ = ZERZER
      ELSEIF (REGNAM .EQ. 'QSVAC') THEN
         BTX = 1.0D-04*QSK1*Y*BRHO
         BTY = 1.0D-04*QSK1*X*BRHO
         BTZ = ZERZER
      ELSE
         CALL FLABRT('MAGFLD', 'NO MAGNETIC FIELD IN THIS REGION !')
      ENDIF

*     GET FIELD STRENGTH:
      B = SQRT(BTX**2 + BTY**2 + BTZ**2)
      IF (B .GT. 1.0D-12) THEN
         BTX = BTX / B
         BTY = BTY / B
         BTZ = BTZ / B
      ELSE
         BTX = ZERZER
         BTY = ZERZER
         BTZ = ONEONE
         B   = ZERZER
      ENDIF

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

