*$ 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(2), NAMLAT_Q(6)
      SAVE NAMLAT_B, NAMLAT_Q
      DOUBLE PRECISION B1Angl, B1Len, B0B, QK1,B2Angl,B2Len
      SAVE B1Angl, B1Len, B0B, QK1

C      INTEGER NROT
 

      DIMENSION NROT(10)
      DATA      NROT /10*0 /
      SAVE      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("Q3Vac   ",NAMLAT_Q(3),NROT,IERR)
         CALL GEON2R("Q4Vac   ",NAMLAT_Q(4),NROT,IERR)
         CALL GEON2R("Q5Vac   ",NAMLAT_Q(5),NROT,IERR)
         CALL GEON2R("Q6Vac   ",NAMLAT_Q(6),NROT,IERR)
     
	 CALL GEON2R("B1Vac   ",NAMLAT_B(1),NROT,IERR)
         CALL GEON2R("B2Vac   ",NAMLAT_B(2),NROT,IERR)

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

         PBEAM = 0.2390D0
         B1Angl = 0.35779D0
         B2Angl = 0.69812D0
         B1Len = 0.23D0
         B2Len = 0.40D0
         B0B = B1Angl*PBEAM/299.79D0/B1Len
         B0B1 = B2Angl*PBEAM/299.79D0/B2Len
         WRITE(LUNOUT,*) ' >>> Beam mom=',PBEAM,' Bend B=',B0B

         QK1 = 10.0D0
         QRAP = 5.5
         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
         BTX = ZERZER
         BTY = -ONEONE
         BTZ = ZERZER
         B = 1.243D0
         WRITE(LUNOUT,*) ' >>> Track in B1 ',X,Y,Z
      END IF 
      IF ( NREG.EQ.NAMLAT_B(2) ) THEN
         B = 1.392D+00
         BTX = ZERZER
         BTY = -ONEONE
         BTZ = ZERZER
         WRITE(LUNOUT,*) ' >>> Track in B2 ',X,Y,Z
      END IF 
             Xs = X
             Ys = Y
             Zs = Z
             IF ( MLATTC .GT. 0) THEN
                CALL UNDOTR (1, Xs, Ys, Zs, NROT(MLATTC) )
             END IF

      IF ( NREG.EQ.NAMLAT_Q(1) ) THEN
             B=0.045D0
             BTX=ZERZER
             BTY=ZERZER
             BTZ=ONEONE
C             B0 = QRAP*QK1
             WRITE(LUNOUT,*) ' >>> Track in Q2 ',X,Y,Z
             IF ((Xs*Xs+Ys*Ys).GT.1.2D-12) THEN
                B= B0B*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
             Xs = X
             Ys = Y
             Zs = Z
       IF ( MLATTC .GT. 0) THEN
           CALL UNDOTR (1, Xs, Ys, Zs, NROT(MLATTC) )
       END IF

      IF ( NREG.EQ.NAMLAT_Q(2) ) THEN
             B=0.101D0
             BTX=ZERZER
             BTY=ZERZER
             BTZ=ONEONE
C             B0 = QRAP*QK1
             WRITE(LUNOUT,*) ' >>> Track in Q2 ',X,Y,Z
             IF ((Xs*Xs+Ys*Ys).GT.1.2D-12) THEN
                B= B0B*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        
             Xs = X
             Ys = Y
             Zs = Z
       IF ( MLATTC .GT. 0) THEN
          CALL UNDOTR (1, Xs, Ys, Zs, NROT(MLATTC) )
       END IF

      IF ( NREG.EQ.NAMLAT_Q(3) ) THEN
             B=0.043D0
             BTX=ZERZER
             BTY=ZERZER
             BTZ=ONEONE
C             B0 = QRAP*QK1
             WRITE(LUNOUT,*) ' >>> Track in Q3 ',X,Y,Z
             IF ((Xs*Xs+Ys*Ys).GT.1.2D-12) THEN
                B= B0B*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
             Xs = X
             Ys = Y
             Zs = Z
       IF ( MLATTC .GT. 0) THEN
          CALL UNDOTR (1, Xs, Ys, Zs, NROT(MLATTC) )
       END IF  

      IF ( NREG.EQ.NAMLAT_Q(4) ) THEN
             B=0.0617D0
             BTX=ZERZER
             BTY=ZERZER
             BTZ=ONEONE
C             B0 = QRAP*QK1
             WRITE(LUNOUT,*) ' >>> Track in Q4 ',X,Y,Z
             IF ((Xs*Xs+Ys*Ys).GT.1.2D-12) THEN
                B= B0B*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(5) ) THEN
             B=0.0590D0
             BTX=ZERZER
             BTY=ZERZER
             BTZ=ONEONE
C             B0 = QRAP*QK1
             WRITE(LUNOUT,*) ' >>> Track in Q5 ',X,Y,Z
             IF ((Xs*Xs+Ys*Ys).GT.1.2D-12) THEN
                B= B0B*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
             Xs = X
             Ys = Y
             Zs = Z
       IF ( MLATTC .GT. 0) THEN
          CALL UNDOTR (1, Xs, Ys, Zs, NROT(MLATTC) )
       END IF 
      IF ( NREG.EQ.NAMLAT_Q(6) ) THEN
             B=0.2411D0
             BTX=ZERZER
             BTY=ZERZER
             BTZ=ONEONE
C             B0 = QRAP*QK1
             WRITE(LUNOUT,*) ' >>> Track in Q6 ',X,Y,Z
             IF ((Xs*Xs+Ys*Ys).GT.1.2D-12) THEN
                B= B0B*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 (MLATTC .GT. 0 ) THEN
             CALL DORTNO (1, BTX, BTY, BTZ, NROT(MLATTC) )
             ENDIF 
      RETURN
*=== End of subroutine Magfld =========================================*
      END