*$ 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 '(LTCLCM)'
      INCLUDE '(TRACKR)'
      INCLUDE '(BEAMCM)'
      INCLUDE '(CMEMFL)'
      INCLUDE '(CSMCRY)'

      LOGICAL LFIRST
      SAVE LFIRST
      DATA LFIRST/.TRUE./

      DOUBLE PRECISION Xs,Ys,Zs
      INTEGER NAMLAT_B(1), NAMLAT_Q(2)
      SAVE NAMLAT_B, NAMLAT_Q
      DOUBLE PRECISION B1Angl, B1Len, B0B, QK1
      SAVE B1Angl, B1Len, B0B, QK1

      INTEGER NROT

      IF (LFIRST) THEN
         WRITE(LUNOUT,*) ' >>> Version 1.0 of routine MAGFLD.F called '
         LFIRST = .FALSE.

         CALL GEON2R("Q1Vac   ",NAMLAT_Q(1),NROT,IERR)
         CALL GEON2R("Q2Vac   ",NAMLAT_Q(2),NROT,IERR)
        
	 CALL GEON2R("B1Vac   ",NAMLAT_B(1),NROT,IERR)
      

         WRITE(LUNOUT,*) ' >>> Magnet IDs: BENDS=',NAMLAT_B,
     +        ' QUADS=',NAMLAT_Q

C        PBEAM = 0.030D0
C        B1Angl = 0.35779D0
C        B1Len = 0.23D0
C        B0B = B1Angl*PBEAM/299.79D0/B1Len

         WRITE(LUNOUT,*) ' >>> Beam mom=',PBEAM,' Bend B=',B0B

         QK1 = 10.0D0
         QRAP = 11.0D0
         WRITE(LUNOUT,*) ' >>> Quad app=',QRAP

      ENDIF
      
C --- Identify the current lattice and reference region each time

      WRITE(LUNOUT,*) ' >>> Lattice MLATTC=',MLATTC,' NREG=',NREG

      IF ( NREG.EQ.NAMLAT_B(1) ) THEN
          B = -1.243D0
          BTX = ZERZER
          BTY = ONEONE
          BTZ = ZERZER

         WRITE(LUNOUT,*) ' >>> Track in B1 ',X,Y,Z
      END IF  

      IF ( NREG.EQ.NAMLAT_Q(1) ) THEN
             Xs = X
             Ys = Y
             Zs = Z
             B=0.045D0
             BTX=ZERZER
             BTY=ZERZER
             BTZ=ONEONE
             B0 = QRAP*QK1
             WRITE(LUNOUT,*) ' >>> Track in Q1 ',X,Y,Z
             IF ((Xs*Xs+Ys*Ys).GT.0.0D0) THEN
                B= B0*SQRT(Xs*Xs+Ys*Ys)/QRAP
                BTX=Ys/SQRT(Xs*Xs+Ys*Ys)
                BTY=-Xs/SQRT(Xs*Xs+Ys*Ys)
                BTZ=0.0D0
             END IF
      END IF
      IF ( NREG.EQ.NAMLAT_Q(2) ) THEN
             Xs = X
             Ys = Y
             Zs = Z
             B=0.045D0
             BTX=ZERZER
             BTY=ZERZER
             BTZ=ONEONE
             B0 = QRAP*QK1
             WRITE(LUNOUT,*) ' >>> Track in Q2 ',X,Y,Z
             IF ((Xs*Xs+Ys*Ys).GT.0.0D0) THEN
                B= B0*SQRT(Xs*Xs+Ys*Ys)/QRAP
                BTX=Ys/SQRT(Xs*Xs+Ys*Ys)
                BTY=-Xs/SQRT(Xs*Xs+Ys*Ys)
                BTZ=0.0D0
             END IF
      END IF
      RETURN
*=== End of subroutine Magfld =========================================*
      END
