*$ 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 Yang Yao trying to implent DIP and QUAD field to HFRS
      INCLUDE '(CMEMFL)'
      INCLUDE '(CSMCRY)'
      INCLUDE '(RTDFCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(BEAMCM)'
      INCLUDE '(TRACKR)'
*     Define the string length and allocate space  
      CHARACTER*8 REGNAM 
*     LFIRST , to flag first call of routine
      LOGICAL LFIRST 
      DATA LFIRST /.TRUE./
      SAVE LFIRST
*     Magnetic filed data
*     (1)Dipole field(T)
      DATA DIPFLD /1.5923567D+00/
      SAVE DIPFLD
*     (2)Quadrupole field(T/m)
*        Q1
      DATA QUAFLD1 /4.33355D+00/
      SAVE QUAFLD1
*        Q2
      DATA QUAFLD2 /4.665797D+00/
      SAVE QUAFLD2
*        Q3
      DATA QUAFLD3 /2.51449D+00/
      SAVE QUAFLD3
*     Magnetic Regions:F:First;S:Second;T:Third
      SAVE QUADF, QUADS, QUADT, DIPF
* Charge Correction for the magnetic field
*     CHFACT: Quadrupoles charge factor
*     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 for U beam in HFRS'
         WRITE (LUNOUT,*) 'Initialisation started'
         WRITE (LUNOUT,*) ''
         WRITE (LUNOUT,*) ''
         WRITE (LUNOUT,*) 'Effective charge for U: ', Zfrttk,'+'
         WRITE (LUNOUT,*) 'Dipo field: ',DIPFLD,'T;'
         WRITE (LUNOUT,*) 'Quad field: ',QUAFLD1,'T/m;'
         WRITE (LUNOUT,*) 'Quad field: ',QUAFLD2,'T/m;'
         WRITE (LUNOUT,*) 'Quad field: ',QUAFLD3,'T/m;'
         WRITE (LUNOUT,*) ''
*     Retrieve the indices of the magnetic regions starting from their
*     names and store them
      CALL GEON2R("Q1Vac",QUADF,IERR)
      CALL GEON2R("Q2Vac",QUADS,IERR)
      CALL GEON2R("Q3Vac",QUADT,IERR)
      CALL GEON2R("Dip1Vac",DIPF,IERR)
*     Flush the .out file, in order to free the buffer:
      CALL FFLUSH
      CALL FLUSH(LUNOUT)
*     Avoid to initialise again the routine
      LFIRST =.FALSE.
      ENDIF
*  ASSIGN correct magnetic field to each REGION
         CALL GEOR2N (NREG,REGNAM,IERR)
         CHFACT =  34/Zfrttk
         IF(REGNAM .EQ. 'Q1Vac')THEN
           BTX = QUAFLD1*CHFACT*Y/1D+02
           BTY = QUAFLD1*CHFACT*X/1D+02
           BTZ = ZERZER
         ELSEIF(REGNAM .EQ. 'Q2Vac')THEN
           BTX = -QUAFLD2*CHFACT*Y/1D+02
           BTY = -QUAFLD2*CHFACT*X/1D+02
           BTZ = ZERZER
         ELSEIF(REGNAM .EQ. 'Q3Vac')THEN
           BTX = QUAFLD3*CHFACT*Y/1D+02
           BTY = QUAFLD3*CHFACT*X/1D+02
           BTZ = ZERZER
         ELSEIF(REGNAM .EQ. 'Dip1Vac')THEN
     	     BTX = ZERZER
     	     BTY = DIPFLD*CHFACT
           BTY = ZERZER
        ELSE
           CALL FLABRT('MAGFLD','NO MAGFLD 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
*  Normalised field components
                BTX = BTX/B
                BTY = BTY/B
                BTZ = BTZ/B
        ELSE
*  Restore default values
                BTX = ZERZER
                BTY = ZERZER
                BTZ = ONEONE
                B   = ZERZER
        ENDIF

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

